Changeset 4671 for palm/trunk/SOURCE/turbulence_closure_mod.f90
- Timestamp:
- Sep 9, 2020 8:27:58 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/turbulence_closure_mod.f90 (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/turbulence_closure_mod.f90
r4625 r4671 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Implementation of downward facing USM and LSM surfaces 28 ! 29 ! 4625 2020-07-24 12:13:16Z oliver.maas 27 30 ! for turbulence_closure = '1.5-order-dai': copy value of distance_to_wall to 28 31 ! ghostpoints at dirichlet/radiation boundaries … … 268 271 surf_lsm_v, & 269 272 surf_usm_h, & 270 surf_usm_v 273 surf_usm_v, & 274 surf_type 271 275 272 276 IMPLICIT NONE … … 475 479 ! 476 480 !-- Natural surfaces 477 DO m = 1, surf_lsm_h %ns478 i = surf_lsm_h %i(m)479 j = surf_lsm_h %j(m)480 k = surf_lsm_h %k(m)481 e_p(k,j,i) = ( surf_lsm_h %us(m) / c_0 )**2481 DO m = 1, surf_lsm_h(0)%ns 482 i = surf_lsm_h(0)%i(m) 483 j = surf_lsm_h(0)%j(m) 484 k = surf_lsm_h(0)%k(m) 485 e_p(k,j,i) = ( surf_lsm_h(0)%us(m) / c_0 )**2 482 486 ENDDO 483 487 ! 484 488 !-- Urban surfaces 485 DO m = 1, surf_usm_h %ns486 i = surf_usm_h %i(m)487 j = surf_usm_h %j(m)488 k = surf_usm_h %k(m)489 e_p(k,j,i) = ( surf_usm_h %us(m) / c_0 )**2489 DO m = 1, surf_usm_h(0)%ns 490 i = surf_usm_h(0)%i(m) 491 j = surf_usm_h(0)%j(m) 492 k = surf_usm_h(0)%k(m) 493 e_p(k,j,i) = ( surf_usm_h(0)%us(m) / c_0 )**2 490 494 ENDDO 491 495 ! … … 579 583 ! 580 584 !-- Natural surfaces 581 DO m = 1, surf_lsm_h %ns582 i = surf_lsm_h %i(m)583 j = surf_lsm_h %j(m)584 k = surf_lsm_h %k(m)585 diss_p(k,j,i) = surf_lsm_h %us(m)**3 / ( kappa * surf_lsm_h%z_mo(m) )585 DO m = 1, surf_lsm_h(0)%ns 586 i = surf_lsm_h(0)%i(m) 587 j = surf_lsm_h(0)%j(m) 588 k = surf_lsm_h(0)%k(m) 589 diss_p(k,j,i) = surf_lsm_h(0)%us(m)**3 / ( kappa * surf_lsm_h(0)%z_mo(m) ) 586 590 ENDDO 587 591 ! 588 592 !-- Urban surfaces 589 DO m = 1, surf_usm_h %ns590 i = surf_usm_h %i(m)591 j = surf_usm_h %j(m)592 k = surf_usm_h %k(m)593 diss_p(k,j,i) = surf_usm_h %us(m)**3 / ( kappa * surf_usm_h%z_mo(m) )593 DO m = 1, surf_usm_h(0)%ns 594 i = surf_usm_h(0)%i(m) 595 j = surf_usm_h(0)%j(m) 596 k = surf_usm_h(0)%k(m) 597 diss_p(k,j,i) = surf_usm_h(0)%us(m)**3 / ( kappa * surf_usm_h(0)%z_mo(m) ) 594 598 ENDDO 595 599 ! … … 2150 2154 ENDDO 2151 2155 ! 2152 !-- Natural surfaces, upward- facing2156 !-- Natural surfaces, upward- and downward facing 2153 2157 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 2154 2158 !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) & 2155 !$ACC PRESENT(surf_lsm_h , u, v, drho_air_zw, zu)2156 DO m = 1, surf_lsm_h %ns2157 2158 i = surf_lsm_h %i(m)2159 j = surf_lsm_h %j(m)2160 k = surf_lsm_h %k(m)2159 !$ACC PRESENT(surf_lsm_h(0), u, v, drho_air_zw, zu) 2160 DO m = 1, surf_lsm_h(0)%ns 2161 2162 i = surf_lsm_h(0)%i(m) 2163 j = surf_lsm_h(0)%j(m) 2164 k = surf_lsm_h(0)%k(m) 2161 2165 ! 2162 2166 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v and km are not on the same 2163 2167 !-- grid. Actually, a further interpolation of km onto the u/v-grid is necessary. However, the 2164 2168 !-- effect of this error is negligible. 2165 km_sfc = kappa * surf_lsm_h %us(m) * surf_lsm_h%z_mo(m) / &2166 phi_m( surf_lsm_h %z_mo(m) / surf_lsm_h%ol(m) )2167 2168 surf_lsm_h %u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m) * drho_air_zw(k-1) * &2169 km_sfc = kappa * surf_lsm_h(0)%us(m) * surf_lsm_h(0)%z_mo(m) / & 2170 phi_m( surf_lsm_h(0)%z_mo(m) / surf_lsm_h(0)%ol(m) ) 2171 2172 surf_lsm_h(0)%u_0(m) = u(k+1,j,i) + surf_lsm_h(0)%usws(m) * drho_air_zw(k-1) * & 2169 2173 ( zu(k+1) - zu(k-1) ) / ( km_sfc + 1.0E-20_wp ) 2170 surf_lsm_h %v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m) * drho_air_zw(k-1) * &2174 surf_lsm_h(0)%v_0(m) = v(k+1,j,i) + surf_lsm_h(0)%vsws(m) * drho_air_zw(k-1) * & 2171 2175 ( zu(k+1) - zu(k-1)) / ( km_sfc + 1.0E-20_wp ) 2172 2176 2173 IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) ) > ABS( u(k+1,j,i) - u(k-1,j,i) ) ) & 2174 surf_lsm_h%u_0(m) = u(k-1,j,i) 2175 2176 IF ( ABS( v(k+1,j,i) - surf_lsm_h%v_0(m) ) > ABS( v(k+1,j,i) - v(k-1,j,i) ) ) & 2177 surf_lsm_h%v_0(m) = v(k-1,j,i) 2177 IF ( ABS( u(k+1,j,i) - surf_lsm_h(0)%u_0(m) ) > ABS( u(k+1,j,i) - u(k-1,j,i) ) ) & 2178 surf_lsm_h(0)%u_0(m) = u(k-1,j,i) 2179 2180 IF ( ABS( v(k+1,j,i) - surf_lsm_h(0)%v_0(m) ) > ABS( v(k+1,j,i) - v(k-1,j,i) ) ) & 2181 surf_lsm_h(0)%v_0(m) = v(k-1,j,i) 2182 2183 ENDDO 2184 ! 2185 !-- Natural surfaces, downward-facing surfaces 2186 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 2187 !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) & 2188 !$ACC PRESENT(surf_lsm_h(1), u, v, drho_air_zw, zu, km) 2189 DO m = 1, surf_lsm_h(1)%ns 2190 2191 i = surf_lsm_h(1)%i(m) 2192 j = surf_lsm_h(1)%j(m) 2193 k = surf_lsm_h(1)%k(m) 2194 ! 2195 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v and km are not on the same 2196 !-- grid. Actually, a further interpolation of km onto the u/v-grid is necessary. However, the 2197 !-- effect of this error is negligible. 2198 surf_lsm_h(1)%u_0(m) = u(k-1,j,i) - surf_lsm_h(1)%usws(m) * drho_air_zw(k-1) * & 2199 ( zu(k+1) - zu(k-1) ) / ( km(k,j,i) + 1.0E-20_wp ) 2200 surf_lsm_h(1)%v_0(m) = v(k-1,j,i) - surf_lsm_h(1)%vsws(m) * drho_air_zw(k-1) * & 2201 ( zu(k+1) - zu(k-1) ) / (km(k,j,i) + 1.0E-20_wp ) 2202 2203 IF ( ABS( surf_lsm_h(1)%u_0(m) - u(k-1,j,i) ) > ABS( u(k+1,j,i) - u(k-1,j,i) ) ) & 2204 surf_lsm_h(1)%u_0(m) = u(k+1,j,i) 2205 2206 IF ( ABS( surf_lsm_h(1)%v_0(m) - v(k-1,j,i) ) > ABS( v(k+1,j,i) - v(k-1,j,i) ) ) & 2207 surf_lsm_h(1)%v_0(m) = v(k+1,j,i) 2178 2208 2179 2209 ENDDO … … 2182 2212 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 2183 2213 !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) & 2184 !$ACC PRESENT(surf_usm_h , u, v, drho_air_zw, zu)2185 DO m = 1, surf_usm_h %ns2186 2187 i = surf_usm_h %i(m)2188 j = surf_usm_h %j(m)2189 k = surf_usm_h %k(m)2214 !$ACC PRESENT(surf_usm_h(0), u, v, drho_air_zw, zu) 2215 DO m = 1, surf_usm_h(0)%ns 2216 2217 i = surf_usm_h(0)%i(m) 2218 j = surf_usm_h(0)%j(m) 2219 k = surf_usm_h(0)%k(m) 2190 2220 ! 2191 2221 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v and km are not on the same 2192 2222 !-- grid. Actually, a further interpolation of km onto the u/v-grid is necessary. However, the 2193 2223 !-- effect of this error is negligible. 2194 km_sfc = kappa * surf_usm_h %us(m) * surf_usm_h%z_mo(m) / &2195 phi_m( surf_usm_h %z_mo(m) / surf_usm_h%ol(m) )2196 2197 surf_usm_h %u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m) * drho_air_zw(k-1) * &2224 km_sfc = kappa * surf_usm_h(0)%us(m) * surf_usm_h(0)%z_mo(m) / & 2225 phi_m( surf_usm_h(0)%z_mo(m) / surf_usm_h(0)%ol(m) ) 2226 2227 surf_usm_h(0)%u_0(m) = u(k+1,j,i) + surf_usm_h(0)%usws(m) * drho_air_zw(k-1) * & 2198 2228 ( zu(k+1) - zu(k-1) ) / ( km_sfc + 1.0E-20_wp ) 2199 surf_usm_h %v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m) * drho_air_zw(k-1) * &2229 surf_usm_h(0)%v_0(m) = v(k+1,j,i) + surf_usm_h(0)%vsws(m) * drho_air_zw(k-1) * & 2200 2230 ( zu(k+1) - zu(k-1) ) / ( km_sfc + 1.0E-20_wp ) 2201 2231 2202 IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) ) > ABS( u(k+1,j,i) - u(k-1,j,i) ) ) & 2203 surf_usm_h%u_0(m) = u(k-1,j,i) 2204 2205 IF ( ABS( v(k+1,j,i) - surf_usm_h%v_0(m) ) > ABS( v(k+1,j,i) - v(k-1,j,i) ) ) & 2206 surf_usm_h%v_0(m) = v(k-1,j,i) 2232 IF ( ABS( u(k+1,j,i) - surf_usm_h(0)%u_0(m) ) > ABS( u(k+1,j,i) - u(k-1,j,i) ) ) & 2233 surf_usm_h(0)%u_0(m) = u(k-1,j,i) 2234 2235 IF ( ABS( v(k+1,j,i) - surf_usm_h(0)%v_0(m) ) > ABS( v(k+1,j,i) - v(k-1,j,i) ) ) & 2236 surf_usm_h(0)%v_0(m) = v(k-1,j,i) 2237 2238 ENDDO 2239 ! 2240 !-- Urban surfaces, downward-facing surfaces 2241 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 2242 !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) & 2243 !$ACC PRESENT(surf_usm_h(1), u, v, drho_air_zw, zu, km) 2244 DO m = 1, surf_usm_h(1)%ns 2245 2246 i = surf_usm_h(1)%i(m) 2247 j = surf_usm_h(1)%j(m) 2248 k = surf_usm_h(1)%k(m) 2249 ! 2250 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v and km are not on the same 2251 !-- grid. Actually, a further interpolation of km onto the u/v-grid is necessary. However, the 2252 !-- effect of this error is negligible. 2253 surf_usm_h(1)%u_0(m) = u(k-1,j,i) - surf_usm_h(1)%usws(m) * drho_air_zw(k-1) * & 2254 ( zu(k+1) - zu(k-1) ) / ( km(k,j,i) + 1.0E-20_wp ) 2255 surf_usm_h(1)%v_0(m) = v(k-1,j,i) - surf_usm_h(1)%vsws(m) * drho_air_zw(k-1) * & 2256 ( zu(k+1) - zu(k-1) ) / (km(k,j,i) + 1.0E-20_wp ) 2257 2258 IF ( ABS( surf_usm_h(1)%u_0(m) - u(k-1,j,i) ) > ABS( u(k+1,j,i) - u(k-1,j,i) ) ) & 2259 surf_usm_h(1)%u_0(m) = u(k+1,j,i) 2260 2261 IF ( ABS( surf_usm_h(1)%v_0(m) - v(k-1,j,i) ) > ABS( v(k+1,j,i) - v(k-1,j,i) ) ) & 2262 surf_usm_h(1)%v_0(m) = v(k+1,j,i) 2207 2263 2208 2264 ENDDO … … 2963 3019 ! 2964 3020 !-- Natural surfaces 2965 surf_s = surf_lsm_h %start_index(j,i)2966 surf_e = surf_lsm_h %end_index(j,i)3021 surf_s = surf_lsm_h(0)%start_index(j,i) 3022 surf_e = surf_lsm_h(0)%end_index(j,i) 2967 3023 !$ACC LOOP PRIVATE(m, k) 2968 3024 DO m = surf_s, surf_e 2969 k = surf_lsm_h %k(m)2970 2971 dudz(k) = ( u(k+1,j,i) - surf_lsm_h %u_0(m) ) * dd2zu(k)2972 dvdz(k) = ( v(k+1,j,i) - surf_lsm_h %v_0(m) ) * dd2zu(k)3025 k = surf_lsm_h(0)%k(m) 3026 3027 dudz(k) = ( u(k+1,j,i) - surf_lsm_h(0)%u_0(m) ) * dd2zu(k) 3028 dvdz(k) = ( v(k+1,j,i) - surf_lsm_h(0)%v_0(m) ) * dd2zu(k) 2973 3029 2974 3030 ENDDO 2975 3031 ! 2976 3032 !-- Urban surfaces 2977 surf_s = surf_usm_h %start_index(j,i)2978 surf_e = surf_usm_h %end_index(j,i)3033 surf_s = surf_usm_h(0)%start_index(j,i) 3034 surf_e = surf_usm_h(0)%end_index(j,i) 2979 3035 !$ACC LOOP PRIVATE(m, k) 2980 3036 DO m = surf_s, surf_e 2981 k = surf_usm_h %k(m)2982 2983 dudz(k) = ( u(k+1,j,i) - surf_usm_h %u_0(m) ) * dd2zu(k)2984 dvdz(k) = ( v(k+1,j,i) - surf_usm_h %v_0(m) ) * dd2zu(k)2985 2986 ENDDO 2987 ! 2988 !-- Compute gradients at downward-facing walls, only for non-naturaldefault surfaces3037 k = surf_usm_h(0)%k(m) 3038 3039 dudz(k) = ( u(k+1,j,i) - surf_usm_h(0)%u_0(m) ) * dd2zu(k) 3040 dvdz(k) = ( v(k+1,j,i) - surf_usm_h(0)%v_0(m) ) * dd2zu(k) 3041 3042 ENDDO 3043 ! 3044 !-- Compute gradients at downward-facing surfaces, default surfaces 2989 3045 surf_s = surf_def_h(1)%start_index(j,i) 2990 3046 surf_e = surf_def_h(1)%end_index(j,i) … … 2995 3051 dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k) 2996 3052 dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) 3053 3054 ENDDO 3055 ! 3056 !-- Compute gradients at downward-facing surfaces, natural surfaces 3057 surf_s = surf_lsm_h(1)%start_index(j,i) 3058 surf_e = surf_lsm_h(1)%end_index(j,i) 3059 !$ACC LOOP PRIVATE(m, k) 3060 DO m = surf_s, surf_e 3061 k = surf_lsm_h(1)%k(m) 3062 3063 dudz(k) = ( surf_lsm_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k) 3064 dvdz(k) = ( surf_lsm_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) 3065 3066 ENDDO 3067 ! 3068 !-- Compute gradients at downward-facing surfaces, urban surfaces 3069 surf_s = surf_usm_h(1)%start_index(j,i) 3070 surf_e = surf_usm_h(1)%end_index(j,i) 3071 !$ACC LOOP PRIVATE(m, k) 3072 DO m = surf_s, surf_e 3073 k = surf_usm_h(1)%k(m) 3074 3075 dudz(k) = ( surf_usm_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k) 3076 dvdz(k) = ( surf_usm_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) 2997 3077 2998 3078 ENDDO … … 3125 3205 tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m) 3126 3206 ENDDO 3127 ENDDO 3128 ! 3129 !-- Natural surfaces 3130 surf_s = surf_lsm_h%start_index(j,i)3131 surf_e = surf_lsm_h%end_index(j,i)3132 !$ACC LOOP PRIVATE(m, k)3133 DO m = surf_s, surf_e3134 k = surf_lsm_h%k(m)3135 tmp_flux(k) = drho_air_zw(k-1) * surf_lsm_h%shf(m)3136 ENDDO 3137 ! 3138 !-- Urban surfaces 3139 surf_s = surf_usm_h%start_index(j,i)3140 surf_e = surf_usm_h%end_index(j,i)3141 !$ACC LOOP PRIVATE(m, k)3142 DO m = surf_s, surf_e3143 k = surf_usm_h%k(m)3144 tmp_flux(k) = drho_air_zw(k-1) * surf_usm_h%shf(m)3207 ! 3208 !-- Natural surfaces 3209 surf_s = surf_lsm_h(l)%start_index(j,i) 3210 surf_e = surf_lsm_h(l)%end_index(j,i) 3211 !$ACC LOOP PRIVATE(m, k) 3212 DO m = surf_s, surf_e 3213 k = surf_lsm_h(l)%k(m) 3214 tmp_flux(k) = drho_air_zw(k-1) * surf_lsm_h(l)%shf(m) 3215 ENDDO 3216 ! 3217 !-- Urban surfaces 3218 surf_s = surf_usm_h(l)%start_index(j,i) 3219 surf_e = surf_usm_h(l)%end_index(j,i) 3220 !$ACC LOOP PRIVATE(m, k) 3221 DO m = surf_s, surf_e 3222 k = surf_usm_h(l)%k(m) 3223 tmp_flux(k) = drho_air_zw(k-1) * surf_usm_h(l)%shf(m) 3224 ENDDO 3145 3225 ENDDO 3146 3226 ENDIF … … 3261 3341 ! 3262 3342 !-- Treat horizontal natural surfaces 3263 surf_s = surf_lsm_h%start_index(j,i) 3264 surf_e = surf_lsm_h%end_index(j,i) 3265 DO m = surf_s, surf_e 3266 k = surf_lsm_h%k(m) 3267 3268 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN 3269 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3270 k2 = 0.61_wp * pt(k,j,i) 3271 ELSE IF ( bulk_cloud_model ) THEN 3272 IF ( ql(k,j,i) == 0.0_wp ) THEN 3343 DO l = 0, 1 3344 surf_s = surf_lsm_h(l)%start_index(j,i) 3345 surf_e = surf_lsm_h(l)%end_index(j,i) 3346 DO m = surf_s, surf_e 3347 k = surf_lsm_h(l)%k(m) 3348 3349 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN 3273 3350 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3274 3351 k2 = 0.61_wp * pt(k,j,i) 3275 ELSE 3276 theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i) 3277 temp = theta * exner(k) 3278 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) * & 3279 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3280 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3281 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 3282 k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp ) 3352 ELSE IF ( bulk_cloud_model ) THEN 3353 IF ( ql(k,j,i) == 0.0_wp ) THEN 3354 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3355 k2 = 0.61_wp * pt(k,j,i) 3356 ELSE 3357 theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i) 3358 temp = theta * exner(k) 3359 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) * & 3360 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3361 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3362 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 3363 k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp ) 3364 ENDIF 3365 ELSE IF ( cloud_droplets ) THEN 3366 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 3367 k2 = 0.61_wp * pt(k,j,i) 3283 3368 ENDIF 3284 ELSE IF ( cloud_droplets ) THEN 3285 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 3286 k2 = 0.61_wp * pt(k,j,i) 3287 ENDIF 3288 3289 tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) + k2 * surf_lsm_h%qsws(m) ) & 3290 * drho_air_zw(k-1) 3369 3370 tmp_flux(k) = ( k1 * surf_lsm_h(l)%shf(m) + k2 * surf_lsm_h(l)%qsws(m) ) & 3371 * drho_air_zw(k-1) 3372 ENDDO 3291 3373 ENDDO 3292 3374 ! 3293 3375 !-- Treat horizontal urban surfaces 3294 surf_s = surf_usm_h%start_index(j,i) 3295 surf_e = surf_usm_h%end_index(j,i) 3296 DO m = surf_s, surf_e 3297 k = surf_usm_h%k(m) 3298 3299 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN 3300 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3301 k2 = 0.61_wp * pt(k,j,i) 3302 ELSE IF ( bulk_cloud_model ) THEN 3303 IF ( ql(k,j,i) == 0.0_wp ) THEN 3376 DO l = 0, 1 3377 surf_s = surf_usm_h(l)%start_index(j,i) 3378 surf_e = surf_usm_h(l)%end_index(j,i) 3379 DO m = surf_s, surf_e 3380 k = surf_usm_h(l)%k(m) 3381 3382 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN 3304 3383 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3305 3384 k2 = 0.61_wp * pt(k,j,i) 3306 ELSE 3307 theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i) 3308 temp = theta * exner(k) 3309 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) * & 3310 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3311 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3312 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 3313 k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp ) 3385 ELSE IF ( bulk_cloud_model ) THEN 3386 IF ( ql(k,j,i) == 0.0_wp ) THEN 3387 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3388 k2 = 0.61_wp * pt(k,j,i) 3389 ELSE 3390 theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i) 3391 temp = theta * exner(k) 3392 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) * & 3393 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3394 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3395 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 3396 k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp ) 3397 ENDIF 3398 ELSE IF ( cloud_droplets ) THEN 3399 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 3400 k2 = 0.61_wp * pt(k,j,i) 3314 3401 ENDIF 3315 ELSE IF ( cloud_droplets ) THEN 3316 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 3317 k2 = 0.61_wp * pt(k,j,i) 3318 ENDIF 3319 3320 tmp_flux(k) = ( k1 * surf_usm_h%shf(m) + k2 * surf_usm_h%qsws(m) ) & 3321 * drho_air_zw(k-1) 3402 3403 tmp_flux(k) = ( k1 * surf_usm_h(l)%shf(m) + k2 * surf_usm_h(l)%qsws(m) ) & 3404 * drho_air_zw(k-1) 3405 ENDDO 3322 3406 ENDDO 3323 3407 … … 3613 3697 ! 3614 3698 !-- Natural surfaces 3615 surf_s = surf_lsm_h %start_index(j,i)3616 surf_e = surf_lsm_h %end_index(j,i)3699 surf_s = surf_lsm_h(0)%start_index(j,i) 3700 surf_e = surf_lsm_h(0)%end_index(j,i) 3617 3701 DO m = surf_s, surf_e 3618 k = surf_lsm_h %k(m)3619 3620 dudz(k) = ( u(k+1,j,i) - surf_lsm_h %u_0(m) ) * dd2zu(k)3621 dvdz(k) = ( v(k+1,j,i) - surf_lsm_h %v_0(m) ) * dd2zu(k)3702 k = surf_lsm_h(0)%k(m) 3703 3704 dudz(k) = ( u(k+1,j,i) - surf_lsm_h(0)%u_0(m) ) * dd2zu(k) 3705 dvdz(k) = ( v(k+1,j,i) - surf_lsm_h(0)%v_0(m) ) * dd2zu(k) 3622 3706 3623 3707 ENDDO 3624 3708 ! 3625 3709 !-- Urban surfaces 3626 surf_s = surf_usm_h %start_index(j,i)3627 surf_e = surf_usm_h %end_index(j,i)3710 surf_s = surf_usm_h(0)%start_index(j,i) 3711 surf_e = surf_usm_h(0)%end_index(j,i) 3628 3712 DO m = surf_s, surf_e 3629 k = surf_usm_h %k(m)3630 3631 dudz(k) = ( u(k+1,j,i) - surf_usm_h %u_0(m) ) * dd2zu(k)3632 dvdz(k) = ( v(k+1,j,i) - surf_usm_h %v_0(m) ) * dd2zu(k)3713 k = surf_usm_h(0)%k(m) 3714 3715 dudz(k) = ( u(k+1,j,i) - surf_usm_h(0)%u_0(m) ) * dd2zu(k) 3716 dvdz(k) = ( v(k+1,j,i) - surf_usm_h(0)%v_0(m) ) * dd2zu(k) 3633 3717 3634 3718 ENDDO 3635 3719 ! 3636 !-- Compute gradients at downward-facing walls, only for non-naturaldefault surfaces3720 !-- Compute gradients at downward-facing surfaces, default surfaces 3637 3721 surf_s = surf_def_h(1)%start_index(j,i) 3638 3722 surf_e = surf_def_h(1)%end_index(j,i) 3639 3723 DO m = surf_s, surf_e 3640 3724 k = surf_def_h(1)%k(m) 3641 3642 3725 dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k) 3643 3726 dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) 3644 3727 ENDDO 3728 ! 3729 !-- Compute gradients at downward-facing surfaces, natural surfaces 3730 surf_s = surf_lsm_h(1)%start_index(j,i) 3731 surf_e = surf_lsm_h(1)%end_index(j,i) 3732 DO m = surf_s, surf_e 3733 k = surf_lsm_h(1)%k(m) 3734 dudz(k) = ( surf_lsm_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k) 3735 dvdz(k) = ( surf_lsm_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) 3736 ENDDO 3737 ! 3738 !-- Compute gradients at downward-facing surfaces, urban surfaces 3739 surf_s = surf_usm_h(1)%start_index(j,i) 3740 surf_e = surf_usm_h(1)%end_index(j,i) 3741 DO m = surf_s, surf_e 3742 k = surf_usm_h(1)%k(m) 3743 dudz(k) = ( surf_usm_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k) 3744 dvdz(k) = ( surf_usm_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) 3645 3745 ENDDO 3646 3746 … … 3747 3847 tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m) 3748 3848 ENDDO 3749 ENDDO 3750 ! 3751 !-- Natural surfaces 3752 surf_s = surf_lsm_h%start_index(j,i)3753 surf_e = surf_lsm_h%end_index(j,i)3754 DO m = surf_s, surf_e3755 k = surf_lsm_h%k(m)3756 tmp_flux(k) = drho_air_zw(k-1) * surf_lsm_h%shf(m)3757 ENDDO 3758 ! 3759 !-- Urban surfaces 3760 surf_s = surf_usm_h%start_index(j,i)3761 surf_e = surf_usm_h%end_index(j,i)3762 DO m = surf_s, surf_e3763 k = surf_usm_h%k(m)3764 tmp_flux(k) = drho_air_zw(k-1) * surf_usm_h%shf(m)3849 ! 3850 !-- Natural surfaces 3851 surf_s = surf_lsm_h(l)%start_index(j,i) 3852 surf_e = surf_lsm_h(l)%end_index(j,i) 3853 DO m = surf_s, surf_e 3854 k = surf_lsm_h(l)%k(m) 3855 tmp_flux(k) = drho_air_zw(k-1) * surf_lsm_h(l)%shf(m) 3856 ENDDO 3857 ! 3858 !-- Urban surfaces 3859 surf_s = surf_usm_h(l)%start_index(j,i) 3860 surf_e = surf_usm_h(l)%end_index(j,i) 3861 DO m = surf_s, surf_e 3862 k = surf_usm_h(l)%k(m) 3863 tmp_flux(k) = drho_air_zw(k-1) * surf_usm_h(l)%shf(m) 3864 ENDDO 3765 3865 ENDDO 3766 3866 ENDIF … … 3836 3936 IF ( use_surface_fluxes ) THEN 3837 3937 ! 3838 !-- Treat horizontal defaultsurfaces3938 !-- Treat horizontal surfaces 3839 3939 DO l = 0, 1 3940 !-- Default surfaces 3840 3941 surf_s = surf_def_h(l)%start_index(j,i) 3841 3942 surf_e = surf_def_h(l)%end_index(j,i) … … 3867 3968 drho_air_zw(k-1) 3868 3969 ENDDO 3869 ENDDO 3870 ! 3871 !-- Treat horizontal natural surfaces 3872 surf_s = surf_lsm_h%start_index(j,i) 3873 surf_e = surf_lsm_h%end_index(j,i) 3874 DO m = surf_s, surf_e 3875 k = surf_lsm_h%k(m) 3876 3877 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN 3878 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3879 k2 = 0.61_wp * pt(k,j,i) 3880 ELSE IF ( bulk_cloud_model ) THEN 3881 IF ( ql(k,j,i) == 0.0_wp ) THEN 3970 ! 3971 !-- natural surfaces 3972 surf_s = surf_lsm_h(l)%start_index(j,i) 3973 surf_e = surf_lsm_h(l)%end_index(j,i) 3974 DO m = surf_s, surf_e 3975 k = surf_lsm_h(l)%k(m) 3976 3977 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN 3882 3978 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3883 3979 k2 = 0.61_wp * pt(k,j,i) 3884 ELSE 3885 theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i) 3886 temp = theta * exner(k) 3887 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) * & 3888 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3889 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3890 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 3891 k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp ) 3980 ELSE IF ( bulk_cloud_model ) THEN 3981 IF ( ql(k,j,i) == 0.0_wp ) THEN 3982 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3983 k2 = 0.61_wp * pt(k,j,i) 3984 ELSE 3985 theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i) 3986 temp = theta * exner(k) 3987 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) * & 3988 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3989 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3990 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 3991 k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp ) 3992 ENDIF 3993 ELSE IF ( cloud_droplets ) THEN 3994 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 3995 k2 = 0.61_wp * pt(k,j,i) 3892 3996 ENDIF 3893 ELSE IF ( cloud_droplets ) THEN 3894 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 3895 k2 = 0.61_wp * pt(k,j,i) 3896 ENDIF 3897 3898 tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) + k2 * surf_lsm_h%qsws(m) ) * & 3899 drho_air_zw(k-1) 3900 ENDDO 3901 ! 3902 !-- Treat horizontal urban surfaces 3903 surf_s = surf_usm_h%start_index(j,i) 3904 surf_e = surf_usm_h%end_index(j,i) 3905 DO m = surf_s, surf_e 3906 k = surf_usm_h%k(m) 3907 3908 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN 3909 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3910 k2 = 0.61_wp * pt(k,j,i) 3911 ELSE IF ( bulk_cloud_model ) THEN 3912 IF ( ql(k,j,i) == 0.0_wp ) THEN 3997 3998 tmp_flux(k) = ( k1 * surf_lsm_h(l)%shf(m) + k2 * surf_lsm_h(l)%qsws(m) ) * & 3999 drho_air_zw(k-1) 4000 ENDDO 4001 ! 4002 !-- urban surfaces 4003 surf_s = surf_usm_h(l)%start_index(j,i) 4004 surf_e = surf_usm_h(l)%end_index(j,i) 4005 DO m = surf_s, surf_e 4006 k = surf_usm_h(l)%k(m) 4007 4008 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN 3913 4009 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3914 4010 k2 = 0.61_wp * pt(k,j,i) 3915 ELSE 3916 theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i) 3917 temp = theta * exner(k) 3918 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) * & 3919 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3920 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3921 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 3922 k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp ) 4011 ELSE IF ( bulk_cloud_model ) THEN 4012 IF ( ql(k,j,i) == 0.0_wp ) THEN 4013 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 4014 k2 = 0.61_wp * pt(k,j,i) 4015 ELSE 4016 theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i) 4017 temp = theta * exner(k) 4018 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) * & 4019 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 4020 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 4021 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 4022 k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp ) 4023 ENDIF 4024 ELSE IF ( cloud_droplets ) THEN 4025 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 4026 k2 = 0.61_wp * pt(k,j,i) 3923 4027 ENDIF 3924 ELSE IF ( cloud_droplets ) THEN 3925 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 3926 k2 = 0.61_wp * pt(k,j,i) 3927 ENDIF 3928 3929 tmp_flux(k) = ( k1 * surf_usm_h%shf(m) + k2 * surf_usm_h%qsws(m) ) * & 3930 drho_air_zw(k-1) 4028 4029 tmp_flux(k) = ( k1 * surf_usm_h(l)%shf(m) + k2 * surf_usm_h(l)%qsws(m) ) * & 4030 drho_air_zw(k-1) 4031 ENDDO 3931 4032 ENDDO 3932 4033 … … 4610 4711 INTEGER(iwp) :: j !< loop index 4611 4712 INTEGER(iwp) :: k !< loop index 4713 INTEGER(iwp) :: l !< loop index 4612 4714 INTEGER(iwp) :: m !< loop index 4613 INTEGER(iwp) :: n !< loop index4614 4715 4615 4716 REAL(wp) :: var_reference !< reference temperature … … 4647 4748 !-- Upward facing surfaces 4648 4749 !-- Default surfaces 4649 n = 0 4650 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 4651 DO m = 1, surf_def_h(0)%ns 4652 i = surf_def_h(0)%i(m) 4653 j = surf_def_h(0)%j(m) 4654 k = surf_def_h(0)%k(m) 4655 km(k,j,i) = kappa * surf_def_h(0)%us(m) * surf_def_h(0)%z_mo(m) / & 4656 phi_m( surf_def_h(0)%z_mo(m) / surf_def_h(0)%ol(m) ) 4657 kh(k,j,i) = 1.35_wp * km(k,j,i) 4750 DO l = 0, 1 4751 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 4752 DO m = 1, surf_def_h(l)%ns 4753 i = surf_def_h(l)%i(m) 4754 j = surf_def_h(l)%j(m) 4755 k = surf_def_h(l)%k(m) 4756 km(k,j,i) = kappa * surf_def_h(l)%us(m) * surf_def_h(l)%z_mo(m) / & 4757 phi_m( surf_def_h(l)%z_mo(m) / surf_def_h(l)%ol(m) ) 4758 kh(k,j,i) = 1.35_wp * km(k,j,i) 4759 ENDDO 4760 ! 4761 !-- Natural surfaces 4762 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 4763 DO m = 1, surf_lsm_h(l)%ns 4764 i = surf_lsm_h(l)%i(m) 4765 j = surf_lsm_h(l)%j(m) 4766 k = surf_lsm_h(l)%k(m) 4767 km(k,j,i) = kappa * surf_lsm_h(l)%us(m) * surf_lsm_h(l)%z_mo(m) / & 4768 phi_m( surf_lsm_h(l)%z_mo(m) / surf_lsm_h(l)%ol(m) ) 4769 kh(k,j,i) = 1.35_wp * km(k,j,i) 4770 ENDDO 4771 ! 4772 !-- Urban surfaces 4773 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 4774 DO m = 1, surf_usm_h(l)%ns 4775 i = surf_usm_h(l)%i(m) 4776 j = surf_usm_h(l)%j(m) 4777 k = surf_usm_h(l)%k(m) 4778 km(k,j,i) = kappa * surf_usm_h(l)%us(m) * surf_usm_h(l)%z_mo(m) / & 4779 phi_m( surf_usm_h(l)%z_mo(m) / surf_usm_h(l)%ol(m) ) 4780 kh(k,j,i) = 1.35_wp * km(k,j,i) 4781 ENDDO 4658 4782 ENDDO 4659 !4660 !-- Natural surfaces4661 !$OMP PARALLEL DO PRIVATE(i,j,k,m)4662 DO m = 1, surf_lsm_h%ns4663 i = surf_lsm_h%i(m)4664 j = surf_lsm_h%j(m)4665 k = surf_lsm_h%k(m)4666 km(k,j,i) = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) / &4667 phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) )4668 kh(k,j,i) = 1.35_wp * km(k,j,i)4669 ENDDO4670 !4671 !-- Urban surfaces4672 !$OMP PARALLEL DO PRIVATE(i,j,k,m)4673 DO m = 1, surf_usm_h%ns4674 i = surf_usm_h%i(m)4675 j = surf_usm_h%j(m)4676 k = surf_usm_h%k(m)4677 km(k,j,i) = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) / &4678 phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) )4679 kh(k,j,i) = 1.35_wp * km(k,j,i)4680 ENDDO4681 4682 4783 ! 4683 4784 !-- North-, south-, west and eastward facing surfaces 4684 4785 !-- Do not consider stratification at these surfaces. 4685 DO n= 0, 34786 DO l = 0, 3 4686 4787 ! 4687 4788 !-- Default surfaces 4688 4789 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 4689 DO m = 1, surf_def_v( n)%ns4690 i = surf_def_v( n)%i(m)4691 j = surf_def_v( n)%j(m)4692 k = surf_def_v( n)%k(m)4693 km(k,j,i) = kappa * surf_def_v( n)%us(m) * surf_def_v(n)%z_mo(m)4790 DO m = 1, surf_def_v(l)%ns 4791 i = surf_def_v(l)%i(m) 4792 j = surf_def_v(l)%j(m) 4793 k = surf_def_v(l)%k(m) 4794 km(k,j,i) = kappa * surf_def_v(l)%us(m) * surf_def_v(l)%z_mo(m) 4694 4795 kh(k,j,i) = 1.35_wp * km(k,j,i) 4695 4796 ENDDO … … 4697 4798 !-- Natural surfaces 4698 4799 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 4699 DO m = 1, surf_lsm_v( n)%ns4700 i = surf_lsm_v( n)%i(m)4701 j = surf_lsm_v( n)%j(m)4702 k = surf_lsm_v( n)%k(m)4703 km(k,j,i) = kappa * surf_lsm_v( n)%us(m) * surf_lsm_v(n)%z_mo(m)4800 DO m = 1, surf_lsm_v(l)%ns 4801 i = surf_lsm_v(l)%i(m) 4802 j = surf_lsm_v(l)%j(m) 4803 k = surf_lsm_v(l)%k(m) 4804 km(k,j,i) = kappa * surf_lsm_v(l)%us(m) * surf_lsm_v(l)%z_mo(m) 4704 4805 kh(k,j,i) = 1.35_wp * km(k,j,i) 4705 4806 ENDDO … … 4707 4808 !-- Urban surfaces 4708 4809 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 4709 DO m = 1, surf_usm_v( n)%ns4710 i = surf_usm_v( n)%i(m)4711 j = surf_usm_v( n)%j(m)4712 k = surf_usm_v( n)%k(m)4713 km(k,j,i) = kappa * surf_usm_v( n)%us(m) * surf_usm_v(n)%z_mo(m)4810 DO m = 1, surf_usm_v(l)%ns 4811 i = surf_usm_v(l)%i(m) 4812 j = surf_usm_v(l)%j(m) 4813 k = surf_usm_v(l)%k(m) 4814 km(k,j,i) = kappa * surf_usm_v(l)%us(m) * surf_usm_v(l)%z_mo(m) 4714 4815 kh(k,j,i) = 1.35_wp * km(k,j,i) 4715 4816 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
