Changeset 4593 for palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
- Timestamp:
- Jul 9, 2020 12:48:18 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/surface_layer_fluxes_mod.f90 (modified) (38 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
r4562 r4593 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 23 ! 22 ! 23 ! 24 24 ! Former revisions: 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - Pre-calculate ln(z/z0) at each timestep in order to reduce the number of log-calculations 28 ! - Bugfix - add missing density to fluxes of passive-scalars, chemistry and cloud-phyiscal 29 ! quantities at upward-facing surfaces 30 ! - Move if-statement out of inner loop 31 ! - Remove unnecessary index referencing 32 ! 33 ! 4562 2020-06-12 08:38:47Z raasch 27 34 ! File re-formatted to follow the PALM coding standard 28 35 ! … … 268 275 downward = .FALSE. !< flag indicating downward-facing surface elements 269 276 ! 277 !-- First, precalculate ln(z/z0) for all surfaces. This is done each timestep, in order to account 278 !-- for time-dependent roughness or user-modifications. 279 DO l = 0, 1 280 IF ( surf_def_h(l)%ns >= 1 ) THEN 281 surf => surf_def_h(l) 282 CALL calc_ln 283 ENDIF 284 ENDDO 285 IF ( surf_lsm_h%ns >= 1 ) THEN 286 surf => surf_lsm_h 287 CALL calc_ln 288 ENDIF 289 IF ( surf_usm_h%ns >= 1 ) THEN 290 surf => surf_usm_h 291 CALL calc_ln 292 ENDIF 293 294 DO l = 0, 3 295 IF ( surf_def_v(l)%ns >= 1 ) THEN 296 surf => surf_def_v(l) 297 CALL calc_ln 298 ENDIF 299 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 300 surf => surf_lsm_v(l) 301 CALL calc_ln 302 ENDIF 303 IF ( surf_usm_v(l)%ns >= 1 ) THEN 304 surf => surf_usm_v(l) 305 CALL calc_ln 306 ENDIF 307 ENDDO 308 ! 270 309 !-- Derive potential temperature and specific humidity at first grid level from the fields pt and q 271 !272 310 !-- First call for horizontal default-type surfaces (l=0 - upward facing, l=1 - downward facing) 273 311 DO l = 0, 1 … … 314 352 315 353 ! 316 !-- First, calculate the new Obukhov length, then new friction velocity, followed by the new scaling 354 !-- First, calculate the new Obukhov length from precalculated values of log(z/z0) and wind speeds. 355 !-- As a second step, then calculate new friction velocity, followed by the new scaling 317 356 !-- parameters (th*, q*, etc.), and the new surface fluxes, if required. Note, each routine is called 318 357 !-- for different surface types. First call for default-type horizontal surfaces, for natural- and 319 358 !-- urban-type surfaces. Note, here only upward-facing horizontal surfaces are treated. 320 321 ! 322 !-- Default-type upward-facing horizontal surfaces359 !-- Note, calculation of log(z/z0) is redone each timestep, in order to account for time-dependent 360 !-- values. 361 !-- Start with default-type upward-facing horizontal surfaces 323 362 IF ( surf_def_h(0)%ns >= 1 ) THEN 324 363 surf => surf_def_h(0) … … 645 684 ! Description: 646 685 ! ------------ 686 !> Compute ln(z/z0). 687 !--------------------------------------------------------------------------------------------------! 688 SUBROUTINE calc_ln 689 690 INTEGER(iwp) :: m !< running index surface elements 691 692 ! 693 !-- Note, ln(z/z0h) and ln(z/z0q) is also calculated even if neural simulations are applied. 694 !-- This is because the scalar coefficients are also used for other scalars such as passive scalars, 695 !-- chemistry and aerosols. 696 !$OMP PARALLEL DO PRIVATE( z_mo ) 697 !$ACC PARALLEL LOOP PRIVATE(z_mo) & 698 !$ACC PRESENT(surf) 699 DO m = 1, surf%ns 700 z_mo = surf%z_mo(m) 701 surf%ln_z_z0(m) = LOG( z_mo / surf%z0(m) ) 702 surf%ln_z_z0h(m) = LOG( z_mo / surf%z0h(m) ) 703 surf%ln_z_z0q(m) = LOG( z_mo / surf%z0q(m) ) 704 ENDDO 705 706 END SUBROUTINE calc_ln 707 708 !--------------------------------------------------------------------------------------------------! 709 ! Description: 710 ! ------------ 647 711 !> Compute the absolute value of the horizontal velocity (relative to the surface) for horizontal 648 712 !> surface elements. This is required by all methods. … … 663 727 ibit = MERGE( 1, 0, .NOT. downward ) 664 728 665 !$OMP PARALLEL DO PRIVATE(i, j, k, w_lfc)666 !$ACC PARALLEL LOOP PRIVATE(i, j, k, w_lfc) &667 !$ACC PRESENT(surf, u, v)668 DO m = 1, surf%ns669 i = surf%i(m)670 j = surf%j(m)671 k = surf%k(m)672 673 ! 674 ! -- Calculate free convection velocity scale w_lfc is use_free_convection_scaling = .T.. This675 !-- will maintain a horizontal velocity even for very weak wind convective conditions. SIGNis676 !-- used to set w_lfc to zero under stable conditions.677 IF ( use_free_convection_scaling ) THEN 729 IF ( use_free_convection_scaling ) THEN 730 !$OMP PARALLEL DO PRIVATE(i, j, k, w_lfc) 731 !$ACC PARALLEL LOOP PRIVATE(i, j, k, w_lfc) & 732 !$ACC PRESENT(surf, u, v) 733 DO m = 1, surf%ns 734 i = surf%i(m) 735 j = surf%j(m) 736 k = surf%k(m) 737 738 ! 739 !-- Calculate free convection velocity scale w_lfc is use_free_convection_scaling = .T.. This 740 !-- will maintain a horizontal velocity even for very weak wind convective conditions. SIGN is 741 !-- used to set w_lfc to zero under stable conditions. 678 742 w_lfc = ABS(g / surf%pt1(m) * surf%z_mo(m) * surf%shf(m)) 679 743 w_lfc = ( 0.5_wp * ( w_lfc + SIGN(w_lfc,surf%shf(m)) ) )**(0.33333_wp) 680 ELSE 681 w_lfc = 0.0_wp 682 ENDIF 683 684 ! 685 !-- Compute the absolute value of the horizontal velocity. (relative to the surface in case the 686 !-- lower surface is the ocean). Please note, in new surface modelling concept the index values 687 !-- changed, i.e. the reference grid point is not the surface-grid point itself but the first 688 !-- grid point outside of the topography. Note, in case of coupled ocean-atmosphere simulations 689 !-- relative velocity with respect to the ocean surface is used, hence, (k-1,j,i) values are used 690 !-- to calculate the absolute velocity. However, this does not apply for downward-facing walls. 691 !-- To mask this, use ibit, which checks for upward/downward-facing surfaces. 692 surf%uvw_abs(m) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) - ( u(k-1,j,i) + u(k-1,j,i+1) ) & 693 * ibit ) )**2 & 694 + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) - ( v(k-1,j,i) + v(k-1,j+1,i) & 695 ) * ibit ) & 744 ! 745 !-- Compute the absolute value of the horizontal velocity. (relative to the surface in case the 746 !-- lower surface is the ocean). Please note, in new surface modelling concept the index values 747 !-- changed, i.e. the reference grid point is not the surface-grid point itself but the first 748 !-- grid point outside of the topography. Note, in case of coupled ocean-atmosphere simulations 749 !-- relative velocity with respect to the ocean surface is used, hence, (k-1,j,i) values are used 750 !-- to calculate the absolute velocity. However, this does not apply for downward-facing walls. 751 !-- To mask this, use ibit, which checks for upward/downward-facing surfaces. 752 surf%uvw_abs(m) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) - ( u(k-1,j,i) + u(k-1,j,i+1) & 753 ) * ibit ) & 754 )**2 & 755 + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) - ( v(k-1,j,i) + v(k-1,j+1,i) & 756 ) * ibit ) & 696 757 )**2 + w_lfc**2 ) 697 ENDDO 758 ENDDO 759 ELSE 760 !$OMP PARALLEL DO PRIVATE(i, j, k) 761 !$ACC PARALLEL LOOP PRIVATE(i, j, k) & 762 !$ACC PRESENT(surf, u, v) 763 DO m = 1, surf%ns 764 i = surf%i(m) 765 j = surf%j(m) 766 k = surf%k(m) 767 ! 768 !-- Compute the absolute value of the horizontal velocity. (relative to the surface in case the 769 !-- lower surface is the ocean). Please note, in new surface modelling concept the index values 770 !-- changed, i.e. the reference grid point is not the surface-grid point itself but the first 771 !-- grid point outside of the topography. Note, in case of coupled ocean-atmosphere simulations 772 !-- relative velocity with respect to the ocean surface is used, hence, (k-1,j,i) values are used 773 !-- to calculate the absolute velocity. However, this does not apply for downward-facing walls. 774 !-- To mask this, use ibit, which checks for upward/downward-facing surfaces. 775 surf%uvw_abs(m) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) - ( u(k-1,j,i) + u(k-1,j,i+1) & 776 ) * ibit ) & 777 )**2 & 778 + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) - ( v(k-1,j,i) + v(k-1,j+1,i) & 779 ) * ibit ) & 780 )**2 ) 781 ENDDO 782 ENDIF 698 783 699 784 END SUBROUTINE calc_uvw_abs … … 980 1065 ! 981 1066 !-- Calculate f = Ri - [...]/[...]^2 = 0 982 f = surf%rib(m) - ( z_mo / ol_m ) * ( LOG( z_mo / surf%z0h(m) )&1067 f = surf%rib(m) - ( z_mo / ol_m ) * ( surf%ln_z_z0h(m) & 983 1068 - psi_h( z_mo / ol_m ) & 984 1069 + psi_h( surf%z0h(m) / ol_m ) ) / & 985 ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / ol_m )&986 + psi_m( surf%z0(m) / ol_m ) )**21070 ( surf%ln_z_z0(m) - psi_m( z_mo / ol_m ) & 1071 + psi_m( surf%z0(m) / ol_m ) )**2 987 1072 988 1073 ! 989 1074 !-- Calculate df/dL 990 f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / surf%z0h(m) )&1075 f_d_ol = ( - ( z_mo / ol_u ) * ( surf%ln_z_z0h(m) & 991 1076 - psi_h( z_mo / ol_u ) & 992 1077 + psi_h( surf%z0h(m) / ol_u ) ) / & 993 ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / ol_u )&994 + psi_m( surf%z0(m) / ol_u ) )**2&995 + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) ) - psi_h( z_mo / ol_l )&996 + psi_h( surf%z0h(m) / ol_l ) ) /&997 ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / ol_l )&998 + psi_m( surf%z0(m) / ol_l ) )**2 ) / ( ol_u - ol_l )1078 ( surf%ln_z_z0(m) - psi_m( z_mo / ol_u ) & 1079 + psi_m( surf%z0(m) / ol_u ) )**2 & 1080 + ( z_mo / ol_l ) * ( surf%ln_z_z0h(m) - psi_h( z_mo / ol_l ) & 1081 + psi_h( surf%z0h(m) / ol_l ) ) /& 1082 ( surf%ln_z_z0(m) - psi_m( z_mo / ol_l ) & 1083 + psi_m( surf%z0(m) / ol_l ) )**2 ) / ( ol_u - ol_l ) 999 1084 ELSE 1000 1085 ! 1001 1086 !-- Calculate f = Ri - 1 /[...]^3 = 0 1002 1087 f = surf%rib(m) - ( z_mo / ol_m ) / & 1003 ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / ol_m ) + psi_m( surf%z0(m) / ol_m ) & 1004 )**3 1088 ( surf%ln_z_z0(m) - psi_m( z_mo / ol_m ) + psi_m( surf%z0(m) / ol_m ) )**3 1005 1089 1006 1090 ! 1007 1091 !-- Calculate df/dL 1008 f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / surf%z0(m) )&1092 f_d_ol = ( - ( z_mo / ol_u ) / ( surf%ln_z_z0(m) & 1009 1093 - psi_m( z_mo / ol_u ) & 1010 1094 + psi_m( surf%z0(m) / ol_u ) & 1011 1095 )**3 & 1012 + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )&1096 + ( z_mo / ol_l ) / ( surf%ln_z_z0(m) & 1013 1097 - psi_m( z_mo / ol_l ) & 1014 1098 + psi_m( surf%z0(m) / ol_l ) & … … 1095 1179 ! 1096 1180 !-- Calculate f = Ri - [...]/[...]^2 = 0 1097 f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) * ( LOG( z_mo_vec(m) / surf%z0h(m) )&1181 f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) * ( surf%ln_z_z0(m) & 1098 1182 - psi_h( z_mo_vec(m) / ol_m ) & 1099 1183 + psi_h( surf%z0h(m) / ol_m ) & 1100 1184 ) / & 1101 ( LOG( z_mo_vec(m) / surf%z0(m) )&1185 ( surf%ln_z_z0(m) & 1102 1186 - psi_m( z_mo_vec(m) / ol_m ) & 1103 1187 + psi_m( surf%z0(m) / ol_m ) & … … 1106 1190 ! 1107 1191 !-- Calculate df/dL 1108 f_d_ol = ( - ( z_mo_vec(m) / ol_u ) * ( LOG( z_mo_vec(m) / surf%z0h(m) )&1192 f_d_ol = ( - ( z_mo_vec(m) / ol_u ) * ( surf%ln_z_z0h(m) & 1109 1193 - psi_h( z_mo_vec(m) / ol_u ) & 1110 1194 + psi_h( surf%z0h(m) / ol_u ) & 1111 1195 ) / & 1112 ( LOG( z_mo_vec(m) / surf%z0(m) )&1196 ( surf%ln_z_z0(m) & 1113 1197 - psi_m( z_mo_vec(m) / ol_u ) & 1114 1198 + psi_m( surf%z0(m) / ol_u ) & 1115 1199 )**2 & 1116 + ( z_mo_vec(m) / ol_l ) * ( LOG( z_mo_vec(m) / surf%z0h(m) )&1200 + ( z_mo_vec(m) / ol_l ) * ( surf%ln_z_z0h(m) & 1117 1201 - psi_h( z_mo_vec(m) / ol_l ) & 1118 1202 + psi_h( surf%z0h(m) / ol_l ) & 1119 1203 ) / & 1120 ( LOG( z_mo_vec(m) / surf%z0(m) )&1204 ( surf%ln_z_z0(m) & 1121 1205 - psi_m( z_mo_vec(m) / ol_l ) & 1122 1206 + psi_m( surf%z0(m) / ol_l ) & … … 1126 1210 ! 1127 1211 !-- Calculate f = Ri - 1 /[...]^3 = 0 1128 f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) / ( LOG( z_mo_vec(m) / surf%z0(m) )&1212 f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) / ( surf%ln_z_z0(m) & 1129 1213 - psi_m( z_mo_vec(m) / ol_m ) & 1130 1214 + psi_m( surf%z0(m) / ol_m ) & … … 1133 1217 ! 1134 1218 !-- Calculate df/dL 1135 f_d_ol = ( - ( z_mo_vec(m) / ol_u ) / ( LOG( z_mo_vec(m) / surf%z0(m) )&1219 f_d_ol = ( - ( z_mo_vec(m) / ol_u ) / ( surf%ln_z_z0(m) & 1136 1220 - psi_m( z_mo_vec(m) / ol_u ) & 1137 1221 + psi_m( surf%z0(m) / ol_u ) & 1138 1222 )**3 & 1139 + ( z_mo_vec(m) / ol_l ) / ( LOG( z_mo_vec(m) / surf%z0(m) )&1223 + ( z_mo_vec(m) / ol_l ) / ( surf%ln_z_z0(m) & 1140 1224 - psi_m( z_mo_vec(m) / ol_l ) & 1141 1225 + psi_m( surf%z0(m) / ol_l ) & … … 1210 1294 ! 1211 1295 !-- Compute u* at the scalars' grid points 1212 surf%us(m) = kappa * surf%uvw_abs(m) / ( LOG( z_mo / surf%z0(m) )&1296 surf%us(m) = kappa * surf%uvw_abs(m) / ( surf%ln_z_z0(m) & 1213 1297 - psi_m( z_mo / surf%ol(m) ) & 1214 1298 + psi_m( surf%z0(m) / surf%ol(m) ) ) … … 1217 1301 !-- Compute u* at downward-facing surfaces. This case, do not consider any stability. 1218 1302 ELSE 1219 !$OMP PARALLEL DO PRIVATE( z_mo )1220 !$ACC PARALLEL LOOP PRIVATE(z_mo)&1303 !$OMP PARALLEL DO 1304 !$ACC PARALLEL LOOP & 1221 1305 !$ACC PRESENT(surf) 1222 1306 DO m = 1, surf%ns 1223 z_mo = surf%z_mo(m)1224 1307 ! 1225 1308 !-- Compute u* at the scalars' grid points 1226 surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m))1309 surf%us(m) = kappa * surf%uvw_abs(m) / surf%ln_z_z0(m) 1227 1310 ENDDO 1228 1311 ENDIF … … 1231 1314 !-- No stability is considered in this case. 1232 1315 ELSE 1233 !$OMP PARALLEL DO PRIVATE( z_mo )1234 !$ACC PARALLEL LOOP PRIVATE(z_mo)&1316 !$OMP PARALLEL DO 1317 !$ACC PARALLEL LOOP & 1235 1318 !$ACC PRESENT(surf) 1236 1319 DO m = 1, surf%ns 1237 z_mo = surf%z_mo(m) 1238 surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) ) 1320 surf%us(m) = kappa * surf%uvw_abs(m) / surf%ln_z_z0(m) 1239 1321 ENDDO 1240 1322 ENDIF … … 1411 1493 z_mo = surf%z_mo(m) 1412 1494 surf%ts(m) = kappa * ( surf%pt1(m) - surf%pt_surface(m) ) & 1413 / ( LOG( z_mo / surf%z0h(m) ) - psi_h( z_mo / surf%ol(m) )&1414 1495 / ( surf%ln_z_z0h(m) - psi_h( z_mo / surf%ol(m) ) & 1496 + psi_h( surf%z0h(m) / surf%ol(m) ) ) 1415 1497 ENDDO 1416 1498 … … 1482 1564 z_mo = surf%z_mo(m) 1483 1565 surf%qs(m) = kappa * ( surf%qv1(m) - surf%q_surface(m) ) & 1484 / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )&1485 1566 / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) ) & 1567 + psi_h( surf%z0q(m) / surf%ol(m) ) ) 1486 1568 ENDDO 1487 1569 ELSE … … 1493 1575 z_mo = surf%z_mo(m) 1494 1576 surf%qs(m) = kappa * ( q(k,j,i) - q(k-1,j,i) ) & 1495 / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )&1496 1577 / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) ) & 1578 + psi_h( surf%z0q(m) / surf%ol(m) ) ) 1497 1579 ENDDO 1498 1580 ENDIF … … 1536 1618 1537 1619 surf%ss(m) = kappa * ( s(k,j,i) - s(k-1,j,i) ) & 1538 / ( LOG( z_mo / surf%z0h(m) ) - psi_h( z_mo / surf%ol(m) )&1539 1620 / ( surf%ln_z_z0h(m) - psi_h( z_mo / surf%ol(m) ) & 1621 + psi_h( surf%z0h(m) / surf%ol(m) ) ) 1540 1622 ENDDO 1541 1623 ENDIF … … 1594 1676 1595 1677 surf%qcs(m) = kappa * ( qc(k,j,i) - qc(k-1,j,i) ) & 1596 / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )&1597 1678 / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) ) & 1679 + psi_h( surf%z0q(m) / surf%ol(m) ) ) 1598 1680 1599 1681 surf%ncs(m) = kappa * ( nc(k,j,i) - nc(k-1,j,i) ) & 1600 / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )&1601 1682 / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) ) & 1683 + psi_h( surf%z0q(m) / surf%ol(m) ) ) 1602 1684 ENDDO 1603 1685 … … 1616 1698 1617 1699 surf%qrs(m) = kappa * ( qr(k,j,i) - qr(k-1,j,i) ) & 1618 / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )&1619 1700 / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) ) & 1701 + psi_h( surf%z0q(m) / surf%ol(m) ) ) 1620 1702 1621 1703 surf%nrs(m) = kappa * ( nr(k,j,i) - nr(k-1,j,i) ) & 1622 / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) )&1623 1704 / ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) ) & 1705 + psi_h( surf%z0q(m) / surf%ol(m) ) ) 1624 1706 ENDDO 1625 1707 … … 1671 1753 1672 1754 surf%usws(m) = kappa * ( u(k,j,i) - u(k-1,j,i) ) & 1673 / ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / surf%ol(m) )&1674 1755 / ( surf%ln_z_z0(m) - psi_m( z_mo / surf%ol(m) ) & 1756 + psi_m( surf%z0(m) / surf%ol(m) ) ) 1675 1757 ! 1676 1758 !-- Please note, the computation of usws is not fully accurate. Actually a further … … 1684 1766 !-- At downward-facing surfaces 1685 1767 ELSE 1686 !$OMP PARALLEL DO PRIVATE( i, j, k , z_mo)1768 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1687 1769 DO m = 1, surf%ns 1688 1770 i = surf%i(m) … … 1690 1772 k = surf%k(m) 1691 1773 1692 z_mo = surf%z_mo(m) 1693 1694 surf%usws(m) = kappa * u(k,j,i) / LOG( z_mo / surf%z0(m) ) 1774 surf%usws(m) = kappa * u(k,j,i) / surf%ln_z_z0(m) 1695 1775 surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k) 1696 1776 ENDDO … … 1713 1793 1714 1794 surf%vsws(m) = kappa * ( v(k,j,i) - v(k-1,j,i) ) & 1715 / ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / surf%ol(m) )&1716 1795 / ( surf%ln_z_z0(m) - psi_m( z_mo / surf%ol(m) ) & 1796 + psi_m( surf%z0(m) / surf%ol(m) ) ) 1717 1797 ! 1718 1798 !-- Please note, the computation of vsws is not fully accurate. Actually a further … … 1726 1806 !-- Downward-facing surfaces 1727 1807 ELSE 1728 !$OMP PARALLEL DO PRIVATE( i, j, k , z_mo)1808 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1729 1809 DO m = 1, surf%ns 1730 1810 i = surf%i(m) … … 1732 1812 k = surf%k(m) 1733 1813 1734 z_mo = surf%z_mo(m) 1735 1736 surf%vsws(m) = kappa * v(k,j,i) / LOG( z_mo / surf%z0(m) ) 1814 surf%vsws(m) = kappa * v(k,j,i) / surf%ln_z_z0(m) 1737 1815 surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k) 1738 1816 ENDDO 1739 1817 ENDIF 1740 1818 ! 1741 !-- Compute the vertical kinematic heat flux 1819 !-- Compute the vertical kinematic heat flux. Note, only upward-facing surfaces are considered, 1820 !-- at downward-facing surfaces the flux is not parametrized with a scaling parameter. 1742 1821 IF ( .NOT. constant_heatflux .AND. ( ( time_since_reference_point <= skip_time_do_lsm & 1743 1822 .AND. simulated_time > 0.0_wp ) .OR. .NOT. land_surface ) .AND. & 1744 1823 .NOT. urban_surface .AND. .NOT. downward ) THEN 1745 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1746 DO m = 1, surf%ns 1747 i = surf%i(m) 1748 j = surf%j(m) 1824 !$OMP PARALLEL DO PRIVATE( k ) 1825 DO m = 1, surf%ns 1749 1826 k = surf%k(m) 1750 1827 surf%shf(m) = -surf%ts(m) * surf%us(m) * rho_air_zw(k-1) … … 1757 1834 .OR. .NOT. land_surface ) .AND. .NOT. urban_surface .AND. .NOT. downward ) & 1758 1835 THEN 1759 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1760 DO m = 1, surf%ns 1761 i = surf%i(m) 1762 j = surf%j(m) 1836 !$OMP PARALLEL DO PRIVATE( k ) 1837 DO m = 1, surf%ns 1763 1838 k = surf%k(m) 1764 1839 surf%qsws(m) = -surf%qs(m) * surf%us(m) * rho_air_zw(k-1) … … 1768 1843 !-- Compute the vertical scalar flux 1769 1844 IF ( .NOT. constant_scalarflux .AND. passive_scalar .AND. .NOT. downward ) THEN 1770 !$OMP PARALLEL DO PRIVATE( i, j ) 1771 DO m = 1, surf%ns 1772 i = surf%i(m) 1773 j = surf%j(m) 1774 surf%ssws(m) = -surf%ss(m) * surf%us(m) 1845 !$OMP PARALLEL DO 1846 DO m = 1, surf%ns 1847 surf%ssws(m) = -surf%ss(m) * surf%us(m) * rho_air_zw(k-1) 1775 1848 ENDDO 1776 1849 ENDIF … … 1779 1852 DO lsp = 1, nvar 1780 1853 IF ( .NOT. constant_csflux(lsp) .AND. air_chemistry .AND. .NOT. downward ) THEN 1781 !$OMP PARALLEL DO PRIVATE( i, j )1854 !$OMP PARALLEL DO 1782 1855 DO m = 1, surf%ns 1783 i = surf%i(m) 1784 j = surf%j(m) 1785 surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m) 1856 surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m) * rho_air_zw(k-1) 1786 1857 ENDDO 1787 1858 ENDIF … … 1791 1862 !-- Compute (turbulent) fluxes of cloud water content and cloud drop conc. 1792 1863 IF ( bulk_cloud_model .AND. microphysics_morrison .AND. .NOT. downward) THEN 1793 !$OMP PARALLEL DO PRIVATE( i, j ) 1794 DO m = 1, surf%ns 1795 i = surf%i(m) 1796 j = surf%j(m) 1797 1798 surf%qcsws(m) = -surf%qcs(m) * surf%us(m) 1799 surf%ncsws(m) = -surf%ncs(m) * surf%us(m) 1864 !$OMP PARALLEL DO 1865 DO m = 1, surf%ns 1866 surf%qcsws(m) = -surf%qcs(m) * surf%us(m) * rho_air_zw(k-1) 1867 surf%ncsws(m) = -surf%ncs(m) * surf%us(m) * rho_air_zw(k-1) 1800 1868 ENDDO 1801 1869 ENDIF … … 1803 1871 !-- Compute (turbulent) fluxes of rain water content and rain drop conc. 1804 1872 IF ( bulk_cloud_model .AND. microphysics_seifert .AND. .NOT. downward) THEN 1805 !$OMP PARALLEL DO PRIVATE( i, j ) 1806 DO m = 1, surf%ns 1807 i = surf%i(m) 1808 j = surf%j(m) 1809 1810 surf%qrsws(m) = -surf%qrs(m) * surf%us(m) 1811 surf%nrsws(m) = -surf%nrs(m) * surf%us(m) 1873 !$OMP PARALLEL DO 1874 DO m = 1, surf%ns 1875 surf%qrsws(m) = -surf%qrs(m) * surf%us(m) * rho_air_zw(k-1) 1876 surf%nrsws(m) = -surf%nrs(m) * surf%us(m) * rho_air_zw(k-1) 1812 1877 ENDDO 1813 1878 ENDIF … … 1832 1897 ! 1833 1898 !-- Calcuate surface fluxes at vertical surfaces. No stability is considered. 1899 !-- Further, no density needs to be considered here. 1834 1900 ELSE 1835 1901 ! … … 1841 1907 flag_u = MERGE( 1.0_wp, 0.0_wp, l == 0 .OR. l == 1 ) 1842 1908 flag_v = MERGE( 1.0_wp, 0.0_wp, l == 2 .OR. l == 3 ) 1843 !$OMP PARALLEL DO PRIVATE( i, j, k , z_mo)1909 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1844 1910 DO m = 1, surf%ns 1845 1911 i = surf%i(m) … … 1847 1913 k = surf%k(m) 1848 1914 1849 z_mo = surf%z_mo(m)1850 1851 1915 surf%mom_flux_uv(m) = kappa * ( flag_u * u(k,j,i) + flag_v * v(k,j,i) ) / & 1852 LOG( z_mo / surf%z0(m))1916 surf%ln_z_z0(m) 1853 1917 1854 1918 surf%mom_flux_uv(m) = - surf%mom_flux_uv(m) * surf%us(m) … … 1858 1922 !-- Compute wsus l={0,1} and wsvs l={2,3} 1859 1923 IF ( mom_w ) THEN 1860 !$OMP PARALLEL DO PRIVATE( i, j, k , z_mo)1924 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1861 1925 DO m = 1, surf%ns 1862 1926 i = surf%i(m) … … 1864 1928 k = surf%k(m) 1865 1929 1866 z_mo = surf%z_mo(m) 1867 1868 surf%mom_flux_w(m) = kappa * w(k,j,i) / LOG( z_mo / surf%z0(m) ) 1930 surf%mom_flux_w(m) = kappa * w(k,j,i) / surf%ln_z_z0(m) 1869 1931 1870 1932 surf%mom_flux_w(m) = - surf%mom_flux_w(m) * surf%us(m) … … 1907 1969 ENDIF 1908 1970 1909 !$OMP PARALLEL DO PRIVATE( i, j, dum , z_mo)1971 !$OMP PARALLEL DO PRIVATE( i, j, dum ) 1910 1972 DO m = 1, surf%ns 1911 1973 i = surf%i(m) 1912 1974 j = surf%j(m) 1913 1975 1914 z_mo = surf%z_mo(m) 1915 1916 dum = kappa / LOG( z_mo / surf%z0(m) ) 1976 dum = kappa / surf%ln_z_z0(m) 1917 1977 ! 1918 1978 !-- usvs (l=0,1) and vsus (l=2,3) … … 1946 2006 CHARACTER(LEN = *), INTENT(IN) :: z_char !< string identifier to identify z level 1947 2007 1948 INTEGER(iwp) :: i !< grid index x-dimension1949 INTEGER(iwp) :: j !< grid index y-dimension1950 INTEGER(iwp) :: k !< grid index z-dimension1951 2008 INTEGER(iwp) :: m !< running index for surface elements 1952 2009 1953 1954 2010 SELECT CASE ( z_char) 1955 2011 … … 1957 2013 1958 2014 DO m = 1, surf%ns 1959 1960 i = surf%i(m)1961 j = surf%j(m)1962 k = surf%k(m)1963 1964 2015 surf%pt_10cm(m) = surf%pt_surface(m) + surf%ts(m) / kappa & 1965 2016 * ( LOG( 0.1_wp / surf%z0h(m) ) - psi_h( 0.1_wp / surf%ol(m) ) &
Note: See TracChangeset
for help on using the changeset viewer.
