Changeset 3167
- Timestamp:
- Dec 30, 2023, 6:37:47 PM (11 months ago)
- Location:
- trunk/LMDZ.MARS/libf/phymars
- Files:
-
- 7 edited
-
callkeys.h (modified) (2 diffs)
-
conf_phys.F (modified) (1 diff)
-
pbl_parameters_mod.F90 (modified) (4 diffs)
-
physiq_mod.F (modified) (5 diffs)
-
soilwater.F90 (modified) (1 diff)
-
vdif_cd_mod.F90 (modified) (5 diffs)
-
vdifc_mod.F (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/phymars/callkeys.h
r3062 r3167 18 18 & ,latentheat_surfwater,gwd_convective_source,startphy_file & 19 19 & ,hdo,hdofrac,cst_cap_albedo,temp_dependent_m,refill_watercap & 20 & ,cloud_adapt_ts 20 & ,cloud_adapt_ts, callatke 21 21 !$OMP THREADPRIVATE(/callkeys_l/) 22 22 … … 36 36 & ,callnirco2,callnlte,callthermos,callconduct, & 37 37 & calleuv,callmolvis,callmoldiff,thermochem,thermoswater & 38 & ,calltherm,callrichsl,callslope,tituscap,callyamada4 38 & ,calltherm,callrichsl,callslope,tituscap,callyamada4,callatke 39 39 40 40 COMMON/aeroutput/dustiropacity -
trunk/LMDZ.MARS/libf/phymars/conf_phys.F
r3144 r3167 243 243 call getin_p("callyamada4",callyamada4) 244 244 write(*,*) " callyamada4 = ",callyamada4 245 245 246 write(*,*) "used latest version of ATKE scheme?" 247 callatke =.false. ! default value 248 call getin_p("callatke",callatke) 249 write(*,*) " callatke = ",callatke 250 251 if(callyamada4.and.callatke) then 252 print*,"You can't use both yamada and atke" 253 print*,"Choose either Yamada or ATKE" 254 call abort_physic(modname, 255 & "Can't use both yamada and ATKE",1) 256 endif 257 246 258 if (calltherm .and. .not.callyamada4) then 247 259 print*,'!!!! WARNING WARNING WARNING !!!!' -
trunk/LMDZ.MARS/libf/phymars/pbl_parameters_mod.F90
r3151 r3167 25 25 use write_output_mod, only: write_output 26 26 use turb_mod, only: turb_resolved 27 27 use lmdz_atke_turbulence_ini, only : smmin, cinf, cepsilon, pr_slope, pr_asym, pr_neut, ri0, ri1, cn, rpi 28 28 29 IMPLICIT NONE 29 30 … … 136 137 REAL vhf(ngrid),vvv(ngrid) ! vertical heat flux (W/m^2) and vertical velocity variance (m/s) 137 138 REAL Teta_out(ngrid,n_out) ! Temporary variable to compute interpolated potential temperature (K) 139 REAL sm ! Stability function computed with the ATKE scheme 140 REAL f_ri_cd_min ! minimum of the stability function with the ATKE scheme 138 141 139 142 ! Code: … … 160 163 itemax= 10 161 164 tol_iter = 0.01 162 165 f_ri_cd_min = 0.01 166 163 167 ! this formulation assumes alphah=1., implying betah=betam 164 168 ! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers : … … 194 198 ENDIF 195 199 196 ! Compute the stability functions fm; fh depending on the stability of the surface layer, from D.E. England et al. (95) 197 ! STABLE BOUNDARY LAYER : 198 IF (rib(ig) .gt. 0.) THEN 199 prandtl(ig) = 1. 200 IF(rib(ig) .lt. ric) THEN 201 ! Assuming alphah=1. and bh=bm for stable conditions : 202 fm(ig) = ((ric-rib(ig))/ric)**2 203 fh(ig) = ((ric-rib(ig))/ric)**2 200 ! Compute the stability functions fm; fh depending on the stability of the surface layer 201 202 IF(callatke) THEN 203 ! Computation following parameterizaiton from ATKE 204 IF(rib(ig) .gt. 0) THEN 205 ! STABLE BOUNDARY LAYER : 206 sm = MAX(smmin,cn*(1.-rib(ig)/ric)) 207 ! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019 208 prandtl(ig) = pr_neut*exp(-pr_slope/pr_neut*rib(ig)+rib(ig)/pr_neut) + rib(ig) * pr_slope 209 fm(ig) = max((sm**(3./2.) * sqrt(cepsilon) * (1 - rib(ig) / prandtl(ig))**(1./2.)),f_ri_cd_min) 210 fh(ig) = max((fm(ig) / prandtl(ig)), f_ri_cd_min) 204 211 ELSE 205 ! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface 206 fm(ig) = 1. 207 fh(ig) = 1. 208 ! should be 0 no? 209 ENDIF 210 ! UNSTABLE BOUNDARY LAYER : 211 ELSE 212 fm(ig) = sqrt(1.-lambda*bm*rib(ig)) 213 fh(ig) = (1./alphah)*((1.-lambda*bh*rib(ig))**0.5)*(1.-lambda*bm*rib(ig))**0.25 214 prandtl(ig) = alphah*((1.-lambda*bm*rib(ig))**0.25)/((1.-lambda*bh*rib(ig))**0.5) 215 ENDIF 216 ! Recompute the reynolds and z0t 217 reynolds(ig) = karman*sqrt(fm(ig))*sqrt(zu2(ig))*pz0(ig)/(log(z1z0)*nu) 218 pz0tcomp(ig) = pz0(ig)*exp(-karman*7.3*(reynolds(ig)**0.25)*(prandtl(ig)**0.5)+5*karman) 219 residual = abs(pz0t-pz0tcomp(ig)) 220 ite = ite+1 212 ! UNSTABLE BOUNDARY LAYER : 213 sm = 2./rpi * (cinf - cn) * atan(-rib(ig)/ri0) + cn 214 prandtl(ig) = -2./rpi * (pr_asym - pr_neut) * atan(rib(ig)/ri1) + pr_neut 215 fm(ig) = MAX((sm**(3./2.) * sqrt(cepsilon) * (1 - rib(ig) / prandtl(ig))**(1./2.)),f_ri_cd_min) 216 fh(ig) = MAX((fm(ig) / prandtl(ig)), f_ri_cd_min) 217 ENDIF ! Rib < 0 218 ELSE 219 ! Computation following parameterizaiton from from D.E. England et al. (95) 220 IF (rib(ig) .gt. 0.) THEN 221 ! STABLE BOUNDARY LAYER : 222 prandtl(ig) = 1. 223 IF(rib(ig) .lt. ric) THEN 224 ! Assuming alphah=1. and bh=bm for stable conditions : 225 fm(ig) = ((ric-rib(ig))/ric)**2 226 fh(ig) = ((ric-rib(ig))/ric)**2 227 ELSE 228 ! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface 229 fm(ig) = 1. 230 fh(ig) = 1. 231 ENDIF 232 ELSE 233 ! UNSTABLE BOUNDARY LAYER : 234 fm(ig) = sqrt(1.-lambda*bm*rib(ig)) 235 fh(ig) = (1./alphah)*((1.-lambda*bh*rib(ig))**0.5)*(1.-lambda*bm*rib(ig))**0.25 236 prandtl(ig) = alphah*((1.-lambda*bm*rib(ig))**0.25)/((1.-lambda*bh*rib(ig))**0.5) 237 ENDIF ! Rib < 0 238 ENDIF ! atke 239 ! Recompute the reynolds and z0t 240 reynolds(ig) = karman*sqrt(fm(ig))*sqrt(zu2(ig))*pz0(ig)/(log(z1z0)*nu) 241 pz0tcomp(ig) = pz0(ig)*exp(-karman*7.3*(reynolds(ig)**0.25)*(prandtl(ig)**0.5)+5*karman) 242 residual = abs(pz0t-pz0tcomp(ig)) 243 ite = ite+1 221 244 ENDDO ! of while 222 245 ! Compute the coefficients Cdv; Cdh : neutral coefficient x stability functions -
trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
r3159 r3167 124 124 use write_output_mod, only: write_output 125 125 use pbl_parameters_mod, only: pbl_parameters 126 use lmdz_atke_turbulence_ini, only : atke_ini 126 127 IMPLICIT NONE 127 128 c======================================================================= … … 574 575 REAL :: ztmp1,ztmp2 ! intermediate variables to compute the mean molar mass of the layer 575 576 577 ! Variable for the computation of the TKE with parameterization from ATKE working group 578 REAL :: viscom ! kinematic molecular viscosity for momentum 579 REAL :: viscoh ! kinematic molecular viscosity for heat 580 581 576 582 c======================================================================= 577 583 pdq(:,:,:) = 0. … … 737 743 ENDIF 738 744 icount=1 745 746 739 747 740 748 #ifndef MESOSCALE … … 794 802 c ~~~~~~~~~~~~~~~ 795 803 if (topflows) call topmons_setup(ngrid) 804 805 c Parameterization of the ATKE routine 806 c ~~~~~~~~~~~~~~~ 807 if (callatke) then 808 viscom = 0.001 809 viscoh = 0.001 810 CALL atke_ini(g, r, pi, cpp, 0., viscom, viscoh) 811 endif 796 812 797 813 #endif … … 2487 2503 2488 2504 ENDIF ! refill_watercap 2489 2490 2505 2491 2506 c----------------------------------------------------------------------- -
trunk/LMDZ.MARS/libf/phymars/soilwater.F90
r3128 r3167 1443 1443 1444 1444 ! call write_output("coeff_diffusion_soil",'interlayer diffusion coefficient','m^2/s',D_inter) 1445 1446 1447 1445 1448 1446 !endif 1449 1447 endif -
trunk/LMDZ.MARS/libf/phymars/vdif_cd_mod.F90
r3165 r3167 21 21 use turb_mod, only: turb_resolved 22 22 use watersat_mod, only: watersat 23 23 use lmdz_atke_turbulence_ini, only : smmin, cinf, cepsilon, pr_slope, pr_asym, pr_neut, ri0, ri1, cn, rpi 24 24 25 IMPLICIT NONE 25 26 include "callkeys.h" 26 27 !======================================================================= 27 28 ! … … 63 64 ! ------------- 64 65 65 #include "callkeys.h"66 67 66 68 67 ! Arguments: … … 134 133 REAL qsat(ngrid) ! saturated mixing ratio of water vapor 135 134 135 REAL sm ! Stability function computed with the ATKE scheme 136 REAL f_ri_cd_min ! minimum of the stability function with the ATKE scheme 137 136 138 ! Code: 137 139 ! -------- … … 148 150 pcdv(:,:) = 0. 149 151 pcdh(:,:) = 0. 150 152 f_ri_cd_min = 0.01 151 153 ! this formulation assumes alphah=1., implying betah=betam 152 154 ! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers : … … 215 217 ENDIF 216 218 217 ! Compute the stability functions fm; fh depending on the stability of the surface layer, from D.E. England et al. (95) 218 219 ! STABLE BOUNDARY LAYER : 220 IF (rib(ig) .gt. 0.) THEN 221 prandtl(ig) = 1. 222 IF(rib(ig) .lt. ric) THEN 219 ! Compute the stability functions fm; fh depending on the stability of the surface layer 220 221 IF(callatke) THEN 222 ! Computation following parameterizaiton from ATKE 223 IF(rib(ig) .gt. 0) THEN 224 ! STABLE BOUNDARY LAYER : 225 sm = MAX(smmin,cn*(1.-rib(ig)/ric)) 226 ! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019 227 prandtl(ig) = pr_neut*exp(-pr_slope/pr_neut*rib(ig)+rib(ig)/pr_neut) + rib(ig) * pr_slope 228 fm(ig) = max((sm**(3./2.) * sqrt(cepsilon) * (1 - rib(ig) / prandtl(ig))**(1./2.)),f_ri_cd_min) 229 fh(ig) = max((fm(ig) / prandtl(ig)), f_ri_cd_min) 230 ELSE 231 ! UNSTABLE BOUNDARY LAYER : 232 sm = 2./rpi * (cinf - cn) * atan(-rib(ig)/ri0) + cn 233 prandtl(ig) = -2./rpi * (pr_asym - pr_neut) * atan(rib(ig)/ri1) + pr_neut 234 fm(ig) = MAX((sm**(3./2.) * sqrt(cepsilon) * (1 - rib(ig) / prandtl(ig))**(1./2.)),f_ri_cd_min) 235 fh(ig) = MAX((fm(ig) / prandtl(ig)), f_ri_cd_min) 236 ENDIF ! Rib < 0 237 ELSE 238 ! Computation following parameterizaiton from from D.E. England et al. (95) 239 IF (rib(ig) .gt. 0.) THEN 240 ! STABLE BOUNDARY LAYER : 241 prandtl(ig) = 1. 242 IF(rib(ig) .lt. ric) THEN 223 243 ! Assuming alphah=1. and bh=bm for stable conditions : 224 fm(ig) = ((ric-rib(ig))/ric)**2 225 fh(ig) = ((ric-rib(ig))/ric)**2 244 fm(ig) = ((ric-rib(ig))/ric)**2 245 fh(ig) = ((ric-rib(ig))/ric)**2 246 ELSE 247 ! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface 248 fm(ig) = 1. 249 fh(ig) = 1. 250 ENDIF 226 251 ELSE 227 ! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface 228 fm(ig) = 1. 229 fh(ig) = 1. 230 ENDIF 231 ! UNSTABLE BOUNDARY LAYER : 232 ELSE 233 fm(ig) = sqrt(1.-lambda*bm*rib(ig)) 234 fh(ig) = (1./alphah)*((1.-lambda*bh*rib(ig))**0.5)*(1.-lambda*bm*rib(ig))**0.25 235 prandtl(ig) = alphah*((1.-lambda*bm*rib(ig))**0.25)/((1.-lambda*bh*rib(ig))**0.5) 236 ENDIF 252 ! UNSTABLE BOUNDARY LAYER : 253 fm(ig) = sqrt(1.-lambda*bm*rib(ig)) 254 fh(ig) = (1./alphah)*((1.-lambda*bh*rib(ig))**0.5)*(1.-lambda*bm*rib(ig))**0.25 255 prandtl(ig) = alphah*((1.-lambda*bm*rib(ig))**0.25)/((1.-lambda*bh*rib(ig))**0.5) 256 ENDIF ! Rib < 0 257 ENDIF ! atke 237 258 ! Recompute the reynolds and z0t 238 259 reynolds(ig) = karman*sqrt(fm(ig))*sqrt(zu2(ig))*pz0(ig)/(log(z1z0)*nu) -
trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
r3166 r3167 36 36 use comsoil_h, only: layer, mlayer,adsorption_soil 37 37 use vdif_cd_mod, only: vdif_cd 38 use lmdz_call_atke, only: call_atke 38 39 IMPLICIT NONE 39 40 … … 466 467 c ** schema de diffusion turbulente dans la couche limite 467 468 c ---------------------------------------------------- 468 IF (.not. callyamada4) THEN469 470 CALL vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay471 & ,pu,pv,ph,zcdv_true_tmp472 & ,pq2,zkv,zkh,zq)473 474 ELSE475 476 469 pt(:,:)=ph(:,:)*ppopsk(:,:) 477 CALL yamada4(ngrid,nlay,nq,ptimestep,g,r,pplev,pt 470 if (callyamada4) then 471 call yamada4(ngrid,nlay,nq,ptimestep,g,r,pplev,pt 478 472 s ,pzlev,pzlay,pu,pv,ph,pq,zcdv_true_tmp,pq2,zkv,zkh,zkq,ustar 479 s ,9) 480 ENDIF 481 473 s ,9) 474 475 elseif (callatke) then 476 call call_atke(ptimestep,ngrid,nlay,zcdv_true_tmp, 477 s zcdh_true_tmp,pu(:,1),pv(:,1),ptsrf_tmp, 478 s pu,pv,pt,zq(:,1,igcm_h2o_vap),pplay,pplev, 479 s pzlay,pzlev,pq2,zkv(:,1:nlay),zkh(:,1:nlay)) 480 481 zkv(:,nlay+1) = zkv(:,nlay) 482 zkh(:,nlay+1) = zkh(:,nlay) 483 else 484 call vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay 485 s ,pu,pv,ph,zcdv_true_tmp 486 s ,pq2,zkv,zkh,zq) 487 488 endif 482 489 if ((doubleq).and.(ngrid.eq.1)) then 483 490 kmixmin = 80. !80.! minimum eddy mix coeff in 1D 484 do ilev= 1,nlay491 do ilev=2,nlay 485 492 do ig=1,ngrid 486 493 zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev)) … … 1013 1020 DO tsub=1,nsubtimestep(ig) 1014 1021 if(tsub.eq.nsubtimestep(ig)) writeoutput = .true. 1015 1016 1022 c C'est parti ! 1017 1023 zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
Note: See TracChangeset
for help on using the changeset viewer.
