Changeset 2696 for palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
- Timestamp:
- Dec 14, 2017 5:12:51 PM (6 years ago)
- Location:
- palm/trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk
-
palm/trunk/SOURCE
-
palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
r2547 r2696 1 1 !> @file surface_layer_fluxes_mod.f90 2 2 !------------------------------------------------------------------------------! 3 ! This file is part of PALM.3 ! This file is part of the PALM model system. 4 4 ! 5 5 ! PALM is free software: you can redistribute it and/or modify it under the … … 18 18 ! 19 19 !------------------------------------------------------------------------------! 20 ! 20 21 ! Current revisions: 21 22 ! ------------------ … … 25 26 ! ----------------- 26 27 ! $Id$ 28 ! - Implementation of chemistry module (FK) 29 ! - Added calculation of pt1 and qv1 for all surface types. Added calculation of 30 ! pt_surface for default-type surfaces (BM) 31 ! - Add flag to disable computation of qsws in case of urban surface (MS) 32 ! 33 ! 2547 2017-10-16 12:41:56Z schwenkel 27 34 ! extended by cloud_droplets option 28 35 ! … … 188 195 !> @todo (re)move large_scale_forcing actions 189 196 !> @todo check/optimize OpenMP directives 197 !> @todo simplify if conditions (which flux need to be computed in which case) 190 198 !------------------------------------------------------------------------------! 191 199 MODULE surface_layer_fluxes_mod … … 195 203 drho_air_zw, rho_air_zw 196 204 205 #if defined( __chem ) 206 USE chem_modules, & 207 ONLY: constant_csflux, nvar 208 #endif 209 197 210 USE cloud_parameters, & 198 211 ONLY: l_d_cp, pt_d_t … … 204 217 205 218 USE control_parameters, & 206 ONLY: cloud_droplets, cloud_physics, constant_heatflux, & 207 constant_scalarflux, constant_waterflux, coupling_mode, g, & 208 humidity, ibc_e_b, ibc_pt_b, initializing_actions, kappa, & 219 ONLY: air_chemistry, cloud_droplets, cloud_physics, & 220 constant_heatflux, constant_scalarflux, & 221 constant_waterflux, coupling_mode, g, humidity, ibc_e_b, & 222 ibc_pt_b, initializing_actions, kappa, & 209 223 intermediate_timestep_count, intermediate_timestep_count_max, & 210 224 land_surface, large_scale_forcing, lsf_surf, & … … 291 305 downward = .FALSE. 292 306 ! 293 !-- In case cloud physics is used, it is required to derive potential 294 !-- temperature and specific humidity at first grid level from the fields pt 295 !-- and q 296 IF ( cloud_physics .OR. cloud_droplets ) THEN 297 ! 298 !-- First call for horizontal default-type surfaces (l=0 - upward facing, 299 !-- l=1 - downward facing) 300 DO l = 0, 1 301 IF ( surf_def_h(l)%ns >= 1 ) THEN 302 surf => surf_def_h(l) 303 CALL calc_pt_q 304 ENDIF 305 ENDDO 306 ! 307 !-- Call for natural-type surfaces 308 IF ( surf_lsm_h%ns >= 1 ) THEN 309 surf => surf_lsm_h 307 !-- Derive potential temperature and specific humidity at first grid level 308 !-- from the fields pt and q 309 ! 310 !-- First call for horizontal default-type surfaces (l=0 - upward facing, 311 !-- l=1 - downward facing) 312 DO l = 0, 1 313 IF ( surf_def_h(l)%ns >= 1 ) THEN 314 surf => surf_def_h(l) 310 315 CALL calc_pt_q 311 ENDIF 312 ! 313 !-- Call for urban-type surfaces 314 IF ( surf_usm_h%ns >= 1 ) THEN 315 surf => surf_usm_h 316 IF ( .NOT. neutral ) CALL calc_pt_surface 317 ENDIF 318 ENDDO 319 ! 320 !-- Call for natural-type horizontal surfaces 321 IF ( surf_lsm_h%ns >= 1 ) THEN 322 surf => surf_lsm_h 323 CALL calc_pt_q 324 ENDIF 325 326 ! 327 !-- Call for urban-type horizontal surfaces 328 IF ( surf_usm_h%ns >= 1 ) THEN 329 surf => surf_usm_h 330 CALL calc_pt_q 331 ENDIF 332 333 ! 334 !-- Call for natural-type vertical surfaces 335 DO l = 0, 3 336 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 337 surf => surf_lsm_v(l) 316 338 CALL calc_pt_q 317 339 ENDIF 318 ENDIF 340 341 ! 342 !-- Call for urban-type vertical surfaces 343 IF ( surf_usm_v(l)%ns >= 1 ) THEN 344 surf => surf_usm_v(l) 345 CALL calc_pt_q 346 ENDIF 347 ENDDO 319 348 320 349 ! … … 923 952 !-- However, this do not apply for downward-facing walls. To mask this, 924 953 !-- use ibit, which checks for upward/downward-facing surfaces. 925 926 954 surf%uvw_abs(m) = SQRT( & 927 955 ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) & … … 1538 1566 surf%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i) 1539 1567 surf%qv1(m) = q(k,j,i) 1568 ELSE 1569 surf%pt1(m) = pt(k,j,i) 1570 IF ( humidity ) THEN 1571 surf%qv1(m) = q(k,j,i) 1572 ELSE 1573 surf%qv1(m) = 0.0_wp 1574 ENDIF 1540 1575 ENDIF 1541 1576 … … 1543 1578 1544 1579 END SUBROUTINE calc_pt_q 1580 1581 1582 ! 1583 !-- Calculate potential temperature and specific humidity at first grid level 1584 !-- ( only for upward-facing surfs ) 1585 SUBROUTINE calc_pt_surface 1586 1587 IMPLICIT NONE 1588 1589 INTEGER(iwp) :: koff !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls) 1590 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1591 1592 koff = surf%koff 1593 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1594 DO m = 1, surf%ns 1595 1596 i = surf%i(m) 1597 j = surf%j(m) 1598 k = surf%k(m) 1599 1600 surf%pt_surface(m) = pt(k+koff,j,i) 1601 1602 ENDDO 1603 1604 END SUBROUTINE calc_pt_surface 1545 1605 1546 1606 ! … … 1552 1612 1553 1613 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1554 1614 INTEGER(iwp) :: lsp !< running index for chemical species 1555 1615 ! 1556 1616 !-- Compute theta* at horizontal surfaces … … 1592 1652 ENDIF 1593 1653 1594 IF ( cloud_physics .OR. cloud_droplets ) THEN 1595 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1596 DO m = 1, surf%ns 1597 1598 i = surf%i(m) 1599 j = surf%j(m) 1600 k = surf%k(m) 1601 1602 z_mo = surf%z_mo(m) 1603 1604 surf%ts(m) = kappa * ( surf%pt1(m) - pt(k-1,j,i) ) & 1605 / ( LOG( z_mo / surf%z0h(m) ) & 1606 - psi_h( z_mo / surf%ol(m) ) & 1607 + psi_h( surf%z0h(m) / surf%ol(m) ) ) 1608 1609 ENDDO 1610 ELSE 1611 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1612 DO m = 1, surf%ns 1613 1614 i = surf%i(m) 1615 j = surf%j(m) 1616 k = surf%k(m) 1617 1618 z_mo = surf%z_mo(m) 1619 1620 surf%ts(m) = kappa * ( pt(k,j,i) - pt(k-1,j,i) ) & 1621 / ( LOG( z_mo / surf%z0h(m) ) & 1622 - psi_h( z_mo / surf%ol(m) ) & 1623 + psi_h( surf%z0h(m) / surf%ol(m) ) ) 1624 ENDDO 1625 ENDIF 1654 !$OMP PARALLEL DO PRIVATE( z_mo ) 1655 DO m = 1, surf%ns 1656 1657 z_mo = surf%z_mo(m) 1658 1659 surf%ts(m) = kappa * ( surf%pt1(m) - surf%pt_surface(m) ) & 1660 / ( LOG( z_mo / surf%z0h(m) ) & 1661 - psi_h( z_mo / surf%ol(m) ) & 1662 + psi_h( surf%z0h(m) / surf%ol(m) ) ) 1663 1664 ENDDO 1665 1626 1666 ENDIF 1627 1667 ! … … 1769 1809 1770 1810 ! 1811 !-- If required compute cs* (chemical species) 1812 #if defined( __chem ) 1813 IF ( air_chemistry ) THEN 1814 ! 1815 !-- At horizontal surfaces 1816 DO lsp = 1, nvar 1817 IF ( constant_csflux(lsp) .AND. .NOT. surf_vertical ) THEN 1818 !-- For a given chemical species' flux in the surface layer 1819 !$OMP PARALLEL DO PRIVATE( i, j ) 1820 DO m = 1, surf%ns 1821 i = surf%i(m) 1822 j = surf%j(m) 1823 surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp ) 1824 ENDDO 1825 ENDIF 1826 ENDDO 1827 ! 1828 !-- At vertical surfaces 1829 IF ( surf_vertical ) THEN 1830 DO lsp = 1, nvar 1831 !$OMP PARALLEL DO PRIVATE( i, j ) 1832 DO m = 1, surf%ns 1833 i = surf%i(m) 1834 j = surf%j(m) 1835 surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp ) 1836 ENDDO 1837 ENDDO 1838 ENDIF 1839 ENDIF 1840 #endif 1841 1842 ! 1771 1843 !-- If required compute qc* and nc* 1772 1844 IF ( cloud_physics .AND. microphysics_morrison .AND. & … … 1830 1902 IMPLICIT NONE 1831 1903 1832 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1904 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1905 INTEGER(iwp) :: lsp !< running index for chemical species 1833 1906 1834 1907 REAL(wp) :: dum !< dummy to precalculate logarithm … … 1883 1956 1884 1957 surf%usws(m) = kappa * u(k,j,i) / LOG( z_mo / surf%z0(m) ) 1885 !1886 !-- To Do: Is the sign correct???1887 1958 surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k) 1888 1959 … … 1929 2000 1930 2001 surf%vsws(m) = kappa * v(k,j,i) / LOG( z_mo / surf%z0(m) ) 1931 1932 2002 surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k) 1933 2003 ENDDO … … 1937 2007 IF ( .NOT. constant_heatflux .AND. ( ( time_since_reference_point& 1938 2008 <= skip_time_do_lsm .AND. simulated_time > 0.0_wp ) .OR. & 1939 .NOT. land_surface ) .AND. .NOT. urban_surface .AND.&2009 .NOT. land_surface ) .AND. .NOT. urban_surface .AND. & 1940 2010 .NOT. downward ) THEN 1941 2011 !$OMP PARALLEL DO PRIVATE( i, j, k ) … … 1951 2021 IF ( .NOT. constant_waterflux .AND. humidity .AND. & 1952 2022 ( ( time_since_reference_point <= skip_time_do_lsm .AND. & 1953 simulated_time > 0.0_wp ) .OR. .NOT. land_surface ) .AND.&1954 .NOT. downward ) THEN2023 simulated_time > 0.0_wp ) .OR. .NOT. land_surface ) .AND. & 2024 .NOT. urban_surface .AND. .NOT. downward ) THEN 1955 2025 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1956 2026 DO m = 1, surf%ns … … 1974 2044 ENDDO 1975 2045 ENDIF 2046 ! 2047 !-- Compute the vertical chemical species' flux 2048 #if defined( __chem ) 2049 DO lsp = 1, nvar 2050 IF ( .NOT. constant_csflux(lsp) .AND. air_chemistry .AND. & 2051 .NOT. downward ) THEN 2052 !$OMP PARALLEL DO PRIVATE( i, j ) 2053 DO m = 1, surf%ns 2054 i = surf%i(m) 2055 j = surf%j(m) 2056 surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m) 2057 ENDDO 2058 ENDIF 2059 ENDDO 2060 #endif 2061 1976 2062 ! 1977 2063 !-- Compute (turbulent) fluxes of cloud water content and cloud drop conc.
Note: See TracChangeset
for help on using the changeset viewer.