Changeset 4804
- Timestamp:
- Feb 7, 2024, 4:07:16 PM (11 months ago)
- Location:
- LMDZ6/trunk/libf/phylmd
- Files:
-
- 2 edited
-
lmdz_atke_exchange_coeff.F90 (modified) (12 diffs)
-
lmdz_atke_turbulence_ini.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90
r4780 r4804 23 23 !======================================================================= 24 24 25 26 27 25 USE lmdz_atke_turbulence_ini, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, atke_ok_virtual, ri0, ri1 28 26 USE lmdz_atke_turbulence_ini, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes, rg, rd, rv, atke_ok_vdiff 29 27 USE lmdz_atke_turbulence_ini, ONLY : viscom, viscoh, clmix, clmixshear, iflag_atke_lmix, lmin, smmin, cn 28 29 !!------------------------------------------------------------------------------------------------------------- 30 ! integer :: iflag_atke ! flag that controls options in atke_compute_km_kh 31 ! integer :: iflag_atke_lmix ! flag that controls the calculation of mixing length in atke 32 ! integer :: iflag_num_atke ! flag that controls the numerical treatment of diffusion coeffiient calculation 33 34 ! logical :: atke_ok_vdiff ! activate vertical diffusion of TKE or not 35 ! logical :: atke_ok_virtual ! account for vapor for flottability 36 37 ! real :: kappa = 0.4 ! Von Karman constant 38 ! real :: cn ! Sm value at Ri=0 39 ! real :: cinf ! Sm value at Ri=-Inf 40 ! real :: ri0 ! Richardson number near zero to guarantee continuity in slope of Sm (stability function) at Ri=0 41 ! real :: ri1 ! Richardson number near zero to guarantee continuity in slope of Pr (Prandlt's number) at Ri=0 42 ! real :: lmin ! minimum mixing length corresponding to the Kolmogov dissipation scale in planetary atmospheres (Chen et al 2016, JGR Atmos) 43 ! real :: ctkes ! coefficient for surface TKE 44 ! real :: clmixshear ! coefficient for mixing length depending on local wind shear 45 46 ! Tunable parameters for the ATKE scheme and their range of values 47 !!------------------------------------------------------------------------------------------------------------- 48 ! real :: cepsilon ! controls the value of the dissipation length scale, range [1.2 - 10] 49 ! real :: cke ! controls the value of the diffusion coefficient of TKE, range [1 - 5] 50 ! real :: l0 ! asymptotic mixing length far from the ground [m] (Sun et al 2011, JAMC), range [15 - 75] 51 ! real :: clmix ! controls the value of the mixing length in stratified conditions, range [0.1 - 2] 52 ! real :: ric ! critical Richardson number controlling the slope of Sm in stable conditions, range [0.19 - 0.25] 53 ! real :: smmin ! minimum value of Sm in very stable conditions (defined here as minsqrt(Ez/Ek)) at large Ri, range [0.025 - 0.1] 54 ! real :: pr_neut ! neutral value of the Prandtl number (Ri=0), range [0.7 - 1] 55 ! real :: pr_slope ! linear slope of Pr with Ri in the very stable regime, range [3 - 5] 56 ! real :: cinffac ! ratio between cinf and cn controlling the convective limit of Sm, range [1.2 - 5.0] 57 ! real :: pr_asym ! value of Prandlt in the convective limit(Ri=-Inf), range [0.3 - 0.5] 58 !!------------------------------------------------------------------------------------------------------------- 30 59 31 60 implicit none … … 96 125 ! results should not depend on the choice of preff 97 126 DO ilay=1,nlay 98 127 DO igrid = 1, ngrid 99 128 theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd) 100 129 END DO 101 130 END DO 102 131 103 132 ! account for water vapor mass for buoyancy calculation 104 133 IF (atke_ok_virtual) THEN 105 DO ilay=1,nlay106 DO igrid = 1, ngrid107 rvap=max(0.,qvap(igrid,ilay)/(1.-qvap(igrid,ilay)))108 theta(igrid,ilay)=theta(igrid,ilay)*(1.+rvap/(RD/RV))/(1.+rvap)109 END DO110 END DO134 DO ilay=1,nlay 135 DO igrid = 1, ngrid 136 rvap=max(0.,qvap(igrid,ilay)/(1.-qvap(igrid,ilay))) 137 theta(igrid,ilay)=theta(igrid,ilay)*(1.+rvap/(RD/RV))/(1.+rvap) 138 END DO 139 END DO 111 140 ENDIF 112 141 … … 123 152 124 153 DO ilay=1,nlay 125 DO igrid=1,ngrid126 z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay))127 ENDDO154 DO igrid=1,ngrid 155 z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay)) 156 ENDDO 128 157 ENDDO 129 158 … … 150 179 + Ri(igrid,ilay) * pr_slope 151 180 IF (Ri(igrid,ilay) .GE. Prandtl(igrid,ilay)) THEN 152 call abort_physic("atke_compute_km_kh", &153 'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1)181 call abort_physic("atke_compute_km_kh", & 182 'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1) 154 183 ENDIF 155 184 END IF … … 166 195 167 196 IF (iflag_atke_lmix .EQ. 1 ) THEN 168 ! Blackadar formulation (~kappa l) + buoyancy length scale (Deardoff 1980) for very stable conditions169 DO ilay=2,nlay170 DO igrid=1,ngrid171 l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)172 IF (N2(igrid,ilay) .GT. 0.) THEN173 lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay))174 lstrat=max(lstrat,lmin)175 !Inverse interpolation, Van de Wiel et al. 2010176 l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)177 ENDIF178 ENDDO179 ENDDO180 181 ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN182 ! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms183 ! implies 2 tuning coefficients clmix and clmixshear184 DO ilay=2,nlay185 DO igrid=1,ngrid186 l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)187 IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN188 lstrat=min(clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)), &197 ! Blackadar formulation (~kappa l) + buoyancy length scale (Deardoff 1980) for very stable conditions 198 DO ilay=2,nlay 199 DO igrid=1,ngrid 200 l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) 201 IF (N2(igrid,ilay) .GT. 0.) THEN 202 lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)) 203 lstrat=max(lstrat,lmin) 204 !Inverse interpolation, Van de Wiel et al. 2010 205 l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0) 206 ENDIF 207 ENDDO 208 ENDDO 209 210 ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN 211 ! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms 212 ! implies 2 tuning coefficients clmix and clmixshear 213 DO ilay=2,nlay 214 DO igrid=1,ngrid 215 l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) 216 IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN 217 lstrat=min(clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)), & 189 218 clmixshear*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay))) 190 lstrat=max(lstrat,lmin) 191 !Inverse interpolation, Van de Wiel et al. 2010 192 l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0) 193 ENDIF 194 ENDDO 195 ENDDO 196 197 ELSE IF (iflag_atke_lmix .EQ. 3 ) THEN 198 ! add effect of wind shear on lstrat following grisogono 2010, qjrms 199 ! keeping a single tuning coefficient clmix 200 DO ilay=2,nlay 201 DO igrid=1,ngrid 202 l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) 203 IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN 204 lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay))*(1.0+Ri(igrid,ilay)/(2.*Prandtl(igrid,ilay))) 205 lstrat=max(lstrat,lmin) 206 !Inverse interpolation, Van de Wiel et al. 2010 207 l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0) 208 ENDIF 209 ENDDO 210 ENDDO 211 212 219 lstrat=max(lstrat,lmin) 220 !Inverse interpolation, Van de Wiel et al. 2010 221 l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0) 222 ENDIF 223 ENDDO 224 ENDDO 225 226 ELSE IF (iflag_atke_lmix .EQ. 3 ) THEN 227 ! add effect of wind shear on lstrat following grisogono 2010, qjrms 228 ! keeping a single tuning coefficient clmix 229 DO ilay=2,nlay 230 DO igrid=1,ngrid 231 l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) 232 IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN 233 lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay))*(1.0+Ri(igrid,ilay)/(2.*Prandtl(igrid,ilay))) 234 lstrat=max(lstrat,lmin) 235 !Inverse interpolation, Van de Wiel et al. 2010 236 l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0) 237 ENDIF 238 ENDDO 239 ENDDO 213 240 214 241 ELSE 215 ! default Blackadar formulation: neglect effect of local stratification and shear 216 217 DO ilay=2,nlay+1 218 DO igrid=1,ngrid 219 l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) 220 ENDDO 221 222 ENDDO 242 ! default Blackadar formulation: neglect effect of local stratification and shear 243 DO ilay=2,nlay+1 244 DO igrid=1,ngrid 245 l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) 246 ENDDO 247 ENDDO 223 248 ENDIF 224 249 … … 228 253 IF (iflag_atke == 0) THEN 229 254 230 ! stationary solution (dtke/dt=0)231 232 DO ilay=2,nlay233 DO igrid=1,ngrid234 tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &235 shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))236 ENDDO237 ENDDO255 ! stationary solution (dtke/dt=0) 256 257 DO ilay=2,nlay 258 DO igrid=1,ngrid 259 tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * & 260 shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)) 261 ENDDO 262 ENDDO 238 263 239 264 ELSE IF (iflag_atke == 1) THEN 240 265 241 ! full implicit scheme resolved with a second order polynomial equation 242 ! default solution which shows fair convergence properties 266 ! full implicit scheme resolved with a second order polynomial equation 267 ! default solution which shows fair convergence properties 268 DO ilay=2,nlay 269 DO igrid=1,ngrid 270 qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) 271 delta=(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime)**2. & 272 +4.*(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime*qq + & 273 2.*l_exchange(igrid,ilay)*l_exchange(igrid,ilay)*cepsilon*Sm(igrid,ilay) & 274 *shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))) 275 qq=(-2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime + sqrt(delta))/2. 276 qq=max(0.,qq) 277 tke(igrid,ilay)=0.5*(qq**2) 278 ENDDO 279 ENDDO 280 281 282 ELSE IF (iflag_atke == 2) THEN 283 284 ! semi implicit scheme when l does not depend on tke 285 ! positive-guaranteed if pr slope in stable condition >1 286 243 287 DO ilay=2,nlay 244 DO igrid=1,ngrid 245 qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) 246 delta=(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime)**2. & 247 +4.*(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime*qq + & 248 2.*l_exchange(igrid,ilay)*l_exchange(igrid,ilay)*cepsilon*Sm(igrid,ilay) & 249 *shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))) 250 qq=(-2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime + sqrt(delta))/2. 251 qq=max(0.,qq) 252 tke(igrid,ilay)=0.5*(qq**2) 253 ENDDO 254 ENDDO 255 256 257 ELSE IF (iflag_atke == 2) THEN 258 259 ! semi implicit scheme when l does not depend on tke 260 ! positive-guaranteed if pr slope in stable condition >1 261 262 DO ilay=2,nlay 263 DO igrid=1,ngrid 264 qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) 265 qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.) & 266 *shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) & 267 /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) 268 tke(igrid,ilay)=0.5*(qq**2) 269 ENDDO 270 ENDDO 288 DO igrid=1,ngrid 289 qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) 290 qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.) & 291 *shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) & 292 /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) 293 tke(igrid,ilay)=0.5*(qq**2) 294 ENDDO 295 ENDDO 271 296 272 297 273 298 ELSE IF (iflag_atke == 3) THEN 274 ! numerical resolution adapted from that in MAR (Deleersnijder 1992)275 ! positively defined by construction276 277 DO ilay=2,nlay278 DO igrid=1,ngrid279 qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)280 IF (Ri(igrid,ilay) .LT. 0.) THEN281 netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))282 netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))283 ELSE284 netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))+ &285 l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*N2(igrid,ilay)/Prandtl(igrid,ilay)286 netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)287 ENDIF288 qq=((qq**2)/dtime+qq*netsource)/(qq/dtime+netloss)289 tke(igrid,ilay)=0.5*(qq**2)290 ENDDO291 ENDDO299 ! numerical resolution adapted from that in MAR (Deleersnijder 1992) 300 ! positively defined by construction 301 302 DO ilay=2,nlay 303 DO igrid=1,ngrid 304 qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) 305 IF (Ri(igrid,ilay) .LT. 0.) THEN 306 netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)) 307 netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) 308 ELSE 309 netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))+ & 310 l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*N2(igrid,ilay)/Prandtl(igrid,ilay) 311 netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay) 312 ENDIF 313 qq=((qq**2)/dtime+qq*netsource)/(qq/dtime+netloss) 314 tke(igrid,ilay)=0.5*(qq**2) 315 ENDDO 316 ENDDO 292 317 293 318 ELSE IF (iflag_atke == 4) THEN 294 ! semi implicit scheme from Arpege (V. Masson methodology with295 ! Taylor expansion of the dissipation term)296 DO ilay=2,nlay297 DO igrid=1,ngrid298 qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)299 qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &300 +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) &301 /(1.+2.*qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))302 qq=max(0.,qq)303 tke(igrid,ilay)=0.5*(qq**2)304 ENDDO305 ENDDO319 ! semi implicit scheme from Arpege (V. Masson methodology with 320 ! Taylor expansion of the dissipation term) 321 DO ilay=2,nlay 322 DO igrid=1,ngrid 323 qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) 324 qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) & 325 +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) & 326 /(1.+2.*qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) 327 qq=max(0.,qq) 328 tke(igrid,ilay)=0.5*(qq**2) 329 ENDDO 330 ENDDO 306 331 307 332 308 333 ELSE 309 call abort_physic("atke_compute_km_kh", &334 call abort_physic("atke_compute_km_kh", & 310 335 'numerical treatment of TKE not possible yet', 1) 311 336 … … 316 341 317 342 DO igrid=1,ngrid 318 tke(igrid,nlay+1)=0.343 tke(igrid,nlay+1)=0. 319 344 END DO 320 345 … … 324 349 ! surface TKE calculation inspired from what is done in Arpege (see E. Bazile note) 325 350 DO igrid=1,ngrid 326 ustar=sqrt(cdrag_uv(igrid)*(wind_u(igrid,1)**2+wind_v(igrid,1)**2))327 tke(igrid,1)=ctkes*(ustar**2)351 ustar=sqrt(cdrag_uv(igrid)*(wind_u(igrid,1)**2+wind_v(igrid,1)**2)) 352 tke(igrid,1)=ctkes*(ustar**2) 328 353 END DO 329 354 … … 332 357 !========================== 333 358 IF (atke_ok_vdiff) THEN 334 CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)359 CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke) 335 360 ENDIF 336 361 … … 395 420 DO ilay=2,nlay 396 421 DO igrid=1,ngrid 397 Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay)))*cke422 Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay)))*cke 398 423 ENDDO 399 424 ENDDO … … 423 448 424 449 DO igrid=1,ngrid 425 CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay)426 DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay)450 CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay) 451 DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay) 427 452 ENDDO 428 453 … … 431 456 DO igrid=1,ngrid 432 457 CCK(igrid,ilay)=(tke(igrid,ilay)/bk(igrid,ilay)-ck(igrid,ilay)/bk(igrid,ilay)*CCK(igrid,ilay+1)) & 433 / (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))458 / (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1)) 434 459 DDK(igrid,ilay)=-ak(igrid,ilay)/bk(igrid,ilay)/(1+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1)) 435 460 ENDDO -
LMDZ6/trunk/libf/phylmd/lmdz_atke_turbulence_ini.F90
r4781 r4804 4 4 5 5 ! declaration of constants and parameters 6 save7 6 8 integer :: iflag_atke, iflag_num_atke, iflag_atke_lmix 9 !$OMP THREADPRIVATE(iflag_atke, iflag_num_atke, iflag_atke_lmix) 10 real :: kappa = 0.4 ! Von Karman constant 11 !$OMP THREADPRIVATE(kappa) 12 real :: cinffac 13 !$OMP THREADPRIVATE(cinffac) 14 real :: l0,ric,ri0,ri1,cinf,cn,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke 15 !$OMP THREADPRIVATE(l0,ri0,ri1,ric,cinf,cn,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke) 16 integer :: lunout,prt_level 17 !$OMP THREADPRIVATE(lunout,prt_level) 18 real :: rg, rd, rpi, rcpd, rv 19 !$OMP THREADPRIVATE(rg, rd, rpi, rcpd, rv) 20 real :: viscom, viscoh 21 !$OMP THREADPRIVATE(viscom,viscoh) 22 real :: lmin=0.01 ! minimum mixing length corresponding to the Kolmogov dissipation scale 23 ! in planetary atmospheres (Chen et al 2016, JGR Atmos) 24 !$OMP THREADPRIVATE(lmin) 25 logical :: atke_ok_vdiff, atke_ok_virtual 26 !$OMP THREADPRIVATE(atke_ok_vdiff,atke_ok_virtual) 7 real, save, protected :: rg, rd, rpi, rcpd, rv, viscom, viscoh 8 !$OMP THREADPRIVATE(rg, rd, rpi, rcpd, rv, viscom, viscoh) 9 10 integer, save, protected :: iflag_atke ! flag that controls options in atke_compute_km_kh 11 integer, save, protected :: iflag_atke_lmix ! flag that controls the calculation of mixing length in atke 12 integer, save, protected :: iflag_num_atke ! flag that controls the numerical treatment of diffusion coeffiient calculation 13 !$OMP THREADPRIVATE(iflag_atke, iflag_atke_lmix, iflag_num_atke) 14 15 logical, save, protected :: atke_ok_vdiff ! activate vertical diffusion of TKE or not 16 logical, save, protected :: atke_ok_virtual ! account for vapor for flottability 17 !$OMP THREADPRIVATE(atke_ok_vdiff, atke_ok_virtual) 18 19 real, save, protected :: kappa = 0.4 ! Von Karman constant 20 real, save, protected :: cn ! Sm value at Ri=0 21 real, save, protected :: cinf ! Sm value at Ri=-Inf 22 real, save, protected :: ri0 ! Richardson number near zero to guarantee continuity in slope of Sm (stability function) at Ri=0 23 real, save, protected :: ri1 ! Richardson number near zero to guarantee continuity in slope of Pr (Prandlt's number) at Ri=0 24 real, save, protected :: lmin = 0.01 ! minimum mixing length corresponding to the Kolmogov dissipation scale in planetary atmospheres (Chen et al 2016, JGR Atmos) 25 real, save, protected :: ctkes ! coefficient for surface TKE 26 real, save, protected :: clmixshear ! coefficient for mixing length depending on local wind shear 27 !$OMP THREADPRIVATE(kappa, cn, cinf, ri0, ri1, lmin, ctkes, clmixshear) 28 29 30 ! Tunable parameters for the ATKE scheme and their range of values 31 !!------------------------------------------------------------------------------------------------------------- 32 real, save, protected :: cepsilon ! controls the value of the dissipation length scale, range [1.2 - 10] 33 real, save, protected :: cke ! controls the value of the diffusion coefficient of TKE, range [1 - 5] 34 real, save, protected :: l0 ! asymptotic mixing length far from the ground [m] (Sun et al 2011, JAMC), range [15 - 75] 35 real, save, protected :: clmix ! controls the value of the mixing length in stratified conditions, range [0.1 - 2] 36 real, save, protected :: ric ! critical Richardson number controlling the slope of Sm in stable conditions, range [0.19 - 0.25] 37 real, save, protected :: smmin ! minimum value of Sm in very stable conditions (defined here as minsqrt(Ez/Ek)) at large Ri, range [0.025 - 0.1] 38 real, save, protected :: pr_neut ! neutral value of the Prandtl number (Ri=0), range [0.7 - 1] 39 real, save, protected :: pr_slope ! linear slope of Pr with Ri in the very stable regime, range [3 - 5] 40 real, save, protected :: cinffac ! ratio between cinf and cn controlling the convective limit of Sm, range [1.2 - 5.0] 41 real, save, protected :: pr_asym ! value of Prandlt in the convective limit(Ri=-Inf), range [0.3 - 0.5] 42 !$OMP THREADPRIVATE(cepsilon, cke, l0, clmix, ric, smmin, pr_neut, pr_slope, cinffac, pr_asym) 43 !!------------------------------------------------------------------------------------------------------------- 27 44 28 45 CONTAINS 29 46 30 47 SUBROUTINE atke_ini(rg_in, rd_in, rpi_in, rcpd_in, rv_in, viscom_in, viscoh_in) 48 !!---------------------------------------------------------------------- 49 !! *** ROUTINE atke_ini *** 50 !! 51 !! ** Purpose : Initialization of the atke module and choice of some constants 52 !! 53 !!---------------------------------------------------------------------- 31 54 32 USE ioipsl_getin_p_mod, ONLY : getin_p55 USE ioipsl_getin_p_mod, ONLY : getin_p 33 56 34 real, intent(in) :: rg_in, rd_in, rpi_in, rcpd_in, rv_in, viscom_in, viscoh_in 57 ! input arguments (universal constants for planet) 58 !------------------------------------------------- 59 real, intent(in) :: rg_in, rd_in, rpi_in, rcpd_in, rv_in, viscom_in, viscoh_in 60 !!---------------------------------------------------------------------- 61 62 rg=rg_in ! gravity acceleration 63 rd=rd_in ! dry gas constant (R/M, R=perfect gas constant and M is the molar mass of the fluid) 64 rpi=rpi_in ! Pi number 65 rcpd=rcpd_in ! cp per unit mass of the gas 66 rv=rv_in ! water vapor constant (for simulations in Earth's atmosphere) 67 viscom=viscom_in ! kinematic molecular viscosity for momentum 68 viscoh=viscoh_in ! kinematic molecular viscosity for heat 35 69 36 70 37 ! input arguments (universal constants for planet) 38 !------------------------------------------------- 39 40 ! gravity acceleration 41 rg=rg_in 42 ! dry gas constant (R/M, R=perfect gas constant and M is the molar mass of the fluid) 43 rd=rd_in 44 ! Pi number 45 rpi=rpi_in 46 ! cp per unit mass of the gas 47 rcpd=rcpd_in 48 ! water vapor constant (for simulations in Earth's atmosphere) 49 rv=rv_in 50 ! kinematic molecular viscosity for momentum 51 viscom=viscom_in 52 ! kinematic molecular viscosity for heat 53 viscoh=viscoh_in 71 ! Read flag values in .def files 72 !------------------------------- 73 74 ! flag that controls options in atke_compute_km_kh 75 iflag_atke=0 76 CALL getin_p('iflag_atke',iflag_atke) 77 78 ! flag that controls the calculation of mixing length in atke 79 iflag_atke_lmix=0 80 CALL getin_p('iflag_atke_lmix',iflag_atke_lmix) 81 82 if (iflag_atke .eq. 0 .and. iflag_atke_lmix>0) then 83 call abort_physic("atke_turbulence_ini", & 84 'stationary scheme must use mixing length formulation not depending on tke', 1) 85 endif 86 87 ! activate vertical diffusion of TKE or not 88 atke_ok_vdiff=.false. 89 CALL getin_p('atke_ok_vdiff',atke_ok_vdiff) 90 91 ! account for vapor for flottability 92 atke_ok_virtual=.true. 93 CALL getin_p('atke_ok_virtual',atke_ok_virtual) 54 94 55 95 56 !viscom=1.46E-5 57 !viscoh=2.06E-5 96 ! flag that controls the numerical treatment of diffusion coeffiient calculation 97 iflag_num_atke=0 98 CALL getin_p('iflag_num_atke',iflag_num_atke) 58 99 100 ! asymptotic mixing length in neutral conditions [m] 101 ! Sun et al 2011, JAMC 102 ! between 10 and 40 103 l0=15.0 104 CALL getin_p('atke_l0',l0) 59 105 60 ! Read flag values in .def files 61 !------------------------------- 106 ! critical Richardson number 107 ric=0.25 108 CALL getin_p('atke_ric',ric) 62 109 110 ! constant for tke dissipation calculation 111 cepsilon=5.87 ! default value as in yamada4 112 CALL getin_p('atke_cepsilon',cepsilon) 63 113 64 ! flag that controls options in atke_compute_km_kh65 iflag_atke=066 CALL getin_p('iflag_atke',iflag_atke)114 ! calculation of cn = Sm value at Ri=0 115 ! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0 116 cn=(1./sqrt(cepsilon))**(2/3) 67 117 68 ! flag that controls the calculation of mixing length in atke 69 iflag_atke_lmix=0 70 CALL getin_p('iflag_atke_lmix',iflag_atke_lmix) 118 ! asymptotic value of Sm for Ri=-Inf 119 cinffac=2.0 120 CALL getin_p('atke_cinffac',cinffac) 121 cinf=cinffac*cn 122 if (cinf .le. cn) then 123 call abort_physic("atke_turbulence_ini", & 124 'cinf cannot be lower than cn', 1) 125 endif 71 126 72 if (iflag_atke .eq. 0 .and. iflag_atke_lmix>0) then73 call abort_physic("atke_turbulence_ini", &74 'stationary scheme must use mixing length formulation not depending on tke', 1)75 endif127 ! coefficient for surface TKE 128 ! following Lenderink & Holtslag 2004, ctkes=(cepsilon)**(2/3) 129 ! (provided by limit condition in neutral conditions) 130 ctkes=(cepsilon)**(2./3.) 76 131 77 ! activate vertical diffusion of TKE or not 78 atke_ok_vdiff=.false. 79 CALL getin_p('atke_ok_vdiff',atke_ok_vdiff) 132 ! slope of Pr=f(Ri) for stable conditions 133 pr_slope=5.0 ! default value from Zilitinkevich et al. 2005 134 CALL getin_p('atke_pr_slope',pr_slope) 135 if (pr_slope .le. 1) then 136 call abort_physic("atke_turbulence_ini", & 137 'pr_slope has to be greater than 1 for consistency of the tke scheme', 1) 138 endif 80 139 140 ! value of turbulent prandtl number in neutral conditions (Ri=0) 141 pr_neut=0.8 142 CALL getin_p('atke_pr_neut',pr_neut) 81 143 82 ! account for vapor for flottability 83 atke_ok_virtual=.true. 84 CALL getin_p('atke_ok_virtual',atke_ok_virtual) 144 ! asymptotic turbulent prandt number value for Ri=-Inf 145 pr_asym=0.4 146 CALL getin_p('atke_pr_asym',pr_asym) 147 if (pr_asym .ge. pr_neut) then 148 call abort_physic("atke_turbulence_ini", & 149 'pr_asym must be be lower than pr_neut', 1) 150 endif 85 151 152 ! coefficient for mixing length depending on local stratification 153 clmix=0.5 154 CALL getin_p('atke_clmix',clmix) 86 155 87 ! flag that controls the numerical treatment of diffusion coeffiient calculation88 iflag_num_atke=089 CALL getin_p('iflag_num_atke',iflag_num_atke)156 ! coefficient for mixing length depending on local wind shear 157 clmixshear=0.5 158 CALL getin_p('atke_clmixshear',clmixshear) 90 159 91 ! asymptotic mixing length in neutral conditions [m] 92 ! Sun et al 2011, JAMC 93 ! between 10 and 40 94 l0=15.0 95 CALL getin_p('atke_l0',l0) 160 ! minimum anisotropy coefficient (defined here as minsqrt(Ez/Ek)) at large Ri. 161 ! From Zilitinkevich et al. 2013, it equals sqrt(0.03)~0.17 162 smmin=0.17 163 CALL getin_p('atke_smmin',smmin) 96 164 97 ! critical Richardson number 98 ric=0.25 99 CALL getin_p('atke_ric',ric) 165 ! ratio between the eddy diffusivity coeff for tke wrt that for momentum 166 ! default value from Lenderink et al. 2004 167 cke=2. 168 CALL getin_p('atke_cke',cke) 100 169 170 ! calculation of Ri0 such that continuity in slope of Sm at Ri=0 171 ri0=2./rpi*(cinf - cn)*ric/cn 101 172 102 ! constant for tke dissipation calculation 103 cepsilon=5.87 ! default value as in yamada4 104 CALL getin_p('atke_cepsilon',cepsilon) 173 ! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0 174 ri1 = -2./rpi * (pr_asym - pr_neut) 105 175 106 107 ! calculation of cn = Sm value at Ri=0 108 ! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0 109 cn=(1./sqrt(cepsilon))**(2/3) 110 111 ! asymptotic value of Sm for Ri=-Inf 112 cinffac=2.0 113 CALL getin_p('atke_cinffac',cinffac) 114 cinf=cinffac*cn 115 if (cinf .le. cn) then 116 call abort_physic("atke_turbulence_ini", & 117 'cinf cannot be lower than cn', 1) 118 endif 119 120 121 ! coefficient for surface TKE 122 ! following Lenderink & Holtslag 2004, ctkes=(cepsilon)**(2/3) 123 ! (provided by limit condition in neutral conditions) 124 ctkes=(cepsilon)**(2./3.) 125 126 ! slope of Pr=f(Ri) for stable conditions 127 pr_slope=5.0 ! default value from Zilitinkevich et al. 2005 128 CALL getin_p('atke_pr_slope',pr_slope) 129 if (pr_slope .le. 1) then 130 call abort_physic("atke_turbulence_ini", & 131 'pr_slope has to be greater than 1 for consistency of the tke scheme', 1) 132 endif 133 134 ! value of turbulent prandtl number in neutral conditions (Ri=0) 135 pr_neut=0.8 136 CALL getin_p('atke_pr_neut',pr_neut) 137 138 139 ! asymptotic turbulent prandt number value for Ri=-Inf 140 pr_asym=0.4 141 CALL getin_p('atke_pr_asym',pr_asym) 142 if (pr_asym .ge. pr_neut) then 143 call abort_physic("atke_turbulence_ini", & 144 'pr_asym must be be lower than pr_neut', 1) 145 endif 146 147 148 149 ! coefficient for mixing length depending on local stratification 150 clmix=0.5 151 CALL getin_p('atke_clmix',clmix) 152 153 ! coefficient for mixing length depending on local wind shear 154 clmixshear=0.5 155 CALL getin_p('atke_clmixshear',clmixshear) 156 157 158 ! minimum anisotropy coefficient (defined here as minsqrt(Ez/Ek)) at large Ri. 159 ! From Zilitinkevich et al. 2013, it equals sqrt(0.03)~0.17 160 smmin=0.17 161 CALL getin_p('atke_smmin',smmin) 162 163 164 ! ratio between the eddy diffusivity coeff for tke wrt that for momentum 165 ! default value from Lenderink et al. 2004 166 cke=2. 167 CALL getin_p('atke_cke',cke) 168 169 170 ! calculation of Ri0 such that continuity in slope of Sm at Ri=0 171 ri0=2./rpi*(cinf - cn)*ric/cn 172 173 ! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0 174 ri1 = -2./rpi * (pr_asym - pr_neut) 175 176 177 RETURN 176 RETURN 178 177 179 178 END SUBROUTINE atke_ini
Note: See TracChangeset
for help on using the changeset viewer.
