Changeset 3325 for trunk/LMDZ.MARS
- Timestamp:
- May 13, 2024, 11:13:14 AM (6 months ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 8 edited
-
changelog.txt (modified) (1 diff)
-
deftank/field_def_physics_mars.xml (modified) (1 diff)
-
libf/phymars/conf_phys.F (modified) (3 diffs)
-
libf/phymars/paleoclimate_mod.F90 (modified) (1 diff)
-
libf/phymars/pbl_parameters_mod.F90 (modified) (7 diffs)
-
libf/phymars/physiq_mod.F (modified) (1 diff)
-
libf/phymars/vdif_cd_mod.F90 (modified) (11 diffs)
-
libf/phymars/vdifc_mod.F (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/changelog.txt
r3316 r3325 4626 4626 == 25/04/2024 == JBC 4627 4627 Reversion of r3305 and r3307. 4628 4629 == 13/05/2024 == LL 4630 Some modifications in vdifc and pbl_parameters to include the effect of water buoyoncy on the sublimation of water ice. 4631 Note: it only works with paleoclimate = .true. (since the model is not tuned with that ...). -
trunk/LMDZ.MARS/deftank/field_def_physics_mars.xml
r3305 r3325 542 542 <field id="zdqsdif_ssi_tot" 543 543 long_name="Total flux echanged with subsurface ice" 544 unit="kg.m-2.s-1" /> 544 unit="kg.m-2.s-1" /> 545 546 547 <!-- Water ice sublimation parameters --> 548 <field id="rib_dry_vdif_cd" 549 long_name="Dry Richardson number in vdif_cd" 550 unit="1" /> 551 <field id="rib_vdif_cd" 552 long_name="Richardson number in vdif_cd" 553 unit="1" /> 554 <field id="fm_vdif_cd" 555 long_name="Momentum stability function in vdif_cd" 556 unit="1" /> 557 <field id="fh_vdif_cd" 558 long_name="Heat stability function in vdif_cd" 559 unit="1" /> 560 <field id="z0t_vdif_cd" 561 long_name="Thermal roughness length in vdif_cd" 562 unit="1" /> 563 <field id="z0_vdif_cd" 564 long_name="Momentum roughness length in vdif_cd" 565 unit="1" /> 566 <field id="Reynolds_vdif_cd" 567 long_name="Reynolds number in vdif_cd" 568 unit="1" /> 545 569 546 570 <!-- CO2 cycle --> -
trunk/LMDZ.MARS/libf/phymars/conf_phys.F
r3272 r3325 38 38 USE mod_phys_lmdz_transfert_para, ONLY: bcast 39 39 USE paleoclimate_mod,ONLY: paleoclimate,albedo_perennialco2, 40 & lag_layer 40 & lag_layer, include_waterbuoyancy 41 41 use microphys_h, only: mteta 42 42 use comsoil_h, only: adsorption_soil, choice_ads,ads_const_D, … … 364 364 lag_layer=.false. 365 365 call getin_p("lag_layer",lag_layer) 366 write(*,*) " 367 368 write(*,*)"Is it paleoclimate run?"366 write(*,*) "lag_layer = ", lag_layer 367 368 write(*,*)"Is it paleoclimate run?" 369 369 paleoclimate=.false. ! default value 370 370 call getin_p("paleoclimate",paleoclimate) 371 write(*,*)" paleoclimate = ",paleoclimate 372 371 write(*,*)"paleoclimate = ",paleoclimate 373 372 374 373 write(*,*)"Albedo for perennial CO2 ice?" … … 376 375 call getin_p("albedo_perennialco2",albedo_perennialco2) 377 376 write(*,*)"albedo_perennialco2 = ",albedo_perennialco2 377 378 write(*,*)"Using lag layer??" 379 lag_layer=.false. 380 call getin_p("include_waterbuoyancy",include_waterbuoyancy) 381 write(*,*) "include_waterbuoyancy = ", include_waterbuoyancy 382 383 if (include_waterbuoyancy.and.(.not.paleoclimate)) then 384 print*,'You should not include the water buoyancy effect 385 & on current climate (model not tunned with this effect' 386 call abort_physic(modname, 387 & "include_waterbuoyancy must be used with paleoclimate",1) 388 endif 378 389 379 390 ! TRACERS: -
trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90
r3200 r3325 15 15 16 16 !$OMP THREADPRIVATE(paleoclimate) 17 real, save, allocatable, dimension(:,:) :: h2o_ice_depth ! Thickness of the lag before H2O ice [m] 18 real, save, allocatable, dimension(:,:) :: lag_co2_ice ! Thickness of the lag before CO2 ice [m] 19 real, save, allocatable, dimension(:,:) :: d_coef ! Diffusion coeficent 20 real, save :: albedo_perennialco2 ! Albedo for perennial co2 ice [1] 21 logical, save :: lag_layer ! Does lag layer is present? 22 !$OMP THREADPRIVATE(h2o_ice_depth,d_coef,lag_co2_ice,albedo_perennialco2) 17 real, save, allocatable, dimension(:,:) :: h2o_ice_depth ! Thickness of the lag before H2O ice [m] 18 real, save, allocatable, dimension(:,:) :: lag_co2_ice ! Thickness of the lag before CO2 ice [m] 19 real, save, allocatable, dimension(:,:) :: d_coef ! Diffusion coeficent 20 real, save :: albedo_perennialco2 ! Albedo for perennial co2 ice [1] 21 logical, save :: lag_layer ! Does lag layer is present? 22 logical, save :: include_waterbuoyancy ! Include the effect of water buoyancy when computing the sublimation of water ice ? 23 !$OMP THREADPRIVATE(h2o_ice_depth,d_coef,lag_co2_ice,albedo_perennialco2,include_waterbuoyancy) 23 24 24 25 !======================================================================= -
trunk/LMDZ.MARS/libf/phymars/pbl_parameters_mod.F90
r3219 r3325 19 19 CONTAINS 20 20 21 SUBROUTINE pbl_parameters(ngrid,nlay,ps,pplay,pz0,pg,zzlay,zzlev,pu,pv,wstar_in,hfmax,zmax,tke,pts,ph, z_out,n_out, &21 SUBROUTINE pbl_parameters(ngrid,nlay,ps,pplay,pz0,pg,zzlay,zzlev,pu,pv,wstar_in,hfmax,zmax,tke,pts,ph,pqvap,pqsurf,mumean,z_out,n_out, & 22 22 T_out,u_out,ustar,tstar,vhf,vvv) 23 23 … … 26 26 use turb_mod, only: turb_resolved 27 27 use lmdz_atke_turbulence_ini, only : smmin, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut, ri0, ri1, cn, rpi 28 use watersat_mod, only: watersat 29 use paleoclimate_mod, only: include_waterbuoyancy 28 30 29 31 IMPLICIT NONE … … 92 94 REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay) ! surface temperature, potentiel temperature 93 95 REAL, INTENT(IN) :: z_out(n_out) ! altitudes of the interpolation in the surface layer 96 REAL, INTENT(IN) :: mumean(ngrid) ! Molecular mass of the atmosphere (kg/mol) 97 REAL, INTENT(IN) :: pqvap(ngrid,nlay) ! Water vapor mixing ratio (kg/kg) to account for vapor flottability 98 REAL, INTENT(IN) :: pqsurf(ngrid) ! Water ice frost on the surface (kg/m^2) to account for vapor flottability 94 99 95 100 ! Outputs: … … 142 147 REAL f_ri_cd_min ! minimum of the stability function with the ATKE scheme 143 148 REAL ric_4interp ! critical richardson number which is either the one from Colaitis et al. (2013) or ATKE (1) 149 150 151 REAL tsurf_v(ngrid) ! Virtual surface temperature, accounting for vapor flottability 152 REAL temp_v(ngrid) ! Potential virtual air temperature in the frist layer, accounting for vapor flottability 153 REAL mu_h2o ! Molecular mass of water vapor (kg/mol) 154 REAL tol_frost ! Tolerance to consider the effect of frost on the surface 155 REAL qsat(ngrid) ! saturated mixing ratio of water vapor 156 157 144 158 ! Code: 145 159 ! -------- … … 166 180 tol_iter = 0.01 167 181 f_ri_cd_min = 0.01 182 mu_h2o = 18e-3 183 tol_frost = 1e-4 184 tsurf_v(:) = 0. 185 temp_v(:) = 0. 168 186 169 187 ! this formulation assumes alphah=1., implying betah=betam … … 183 201 endif 184 202 185 203 204 205 IF(include_waterbuoyancy) then 206 temp_v(:) = ph(:,1)*(1.+pqvap(:,1)/(mu_h2o/mumean(:)))/(1.+pqvap(:,1)) 207 DO ig = 1,ngrid 208 IF(pqsurf(ig).gt.tol_frost) then 209 call watersat(1,pts(ig),pplay(ig,1),qsat(ig)) 210 tsurf_v(ig) = pts(ig)*(1.+qsat(ig)/(mu_h2o/mumean(ig)))/(1.+qsat(ig)) 211 ELSE 212 tsurf_v(ig) = pts(ig)*(1.+pqvap(ig,1)/(mu_h2o/mumean(ig)))/(1.+pqvap(ig,1)) 213 ENDIF 214 ENDDO 215 ELSE 216 tsurf_v(:) = pts(:) 217 temp_v(:) = ph(:,1) 218 ENDIF 219 186 220 DO ig=1,ngrid 187 221 ite = 0. … … 199 233 IF ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then 200 234 ! Richardson number formulation proposed by D.E. England et al. (1995) 201 rib(ig) = (pg/ pts(ig))*sqrt(zzlay(ig,1)*pz0(ig))*(((log(zzlay(ig,1)/pz0(ig)))**2)/(log(zzlay(ig,1)/pz0t)))*(ph(ig,1)-pts(ig))/zu2(ig)235 rib(ig) = (pg/tsurf_v(ig))*sqrt(zzlay(ig,1)*pz0(ig))*(((log(zzlay(ig,1)/pz0(ig)))**2)/(log(zzlay(ig,1)/pz0t)))*(temp_v(ig)-tsurf_v(ig))/zu2(ig) 202 236 ELSE 203 237 print*,'warning, infinite Richardson at surface' -
trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
r3304 r3325 2584 2584 call pbl_parameters(ngrid,nlayer,ps,zplay,z0, 2585 2585 & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,q2,tsurf(:,iflat), 2586 & zh,z_out,n_out,T_out,u_out,ustar,tstar,vhf,vvv) 2586 & zh,zq(:,1,igcm_h2o_vap),qsurf(:,igcm_h2o_ice,iflat),mmean(:,1), 2587 & z_out,n_out,T_out,u_out,ustar,tstar,vhf,vvv) 2587 2588 ! pourquoi ustar recalcule ici? fait dans vdifc. 2588 2589 -
trunk/LMDZ.MARS/libf/phymars/vdif_cd_mod.F90
r3219 r3325 16 16 CONTAINS 17 17 18 SUBROUTINE vdif_cd(ngrid,nlay,nslope,pz0,pg,pz,pp,pu,pv,wstar,pts,ph, virtual,mumean,pqvap,pqsurf,pcdv,pcdh)18 SUBROUTINE vdif_cd(ngrid,nlay,nslope,pz0,pg,pz,pp,pu,pv,wstar,pts,ph,mumean,pqvap,pqsurf,write_outputs,pcdv,pcdh) 19 19 20 20 use comcstfi_h … … 22 22 use watersat_mod, only: watersat 23 23 use lmdz_atke_turbulence_ini, only : smmin, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut, ri0, ri1, cn, rpi 24 24 use paleoclimate_mod, only: include_waterbuoyancy 25 use write_output_mod, only: write_output 26 use comslope_mod, ONLY: iflat 25 27 IMPLICIT NONE 26 28 include "callkeys.h" 29 30 27 31 !======================================================================= 28 32 ! … … 47 51 ! pts(ngrid) surface temperature 48 52 ! ph(ngrid) potential temperature T*(p/ps)^kappa 49 ! virtual Boolean to account for vapor flottability50 53 ! mumean Molecular mass of the atmosphere (kg/mol) 51 54 ! pqvap(ngrid,nlay) Water vapor mixing ratio (kg/kg) to account for vapor flottability … … 76 79 REAL, INTENT(IN) :: pts(ngrid,nslope),ph(ngrid,nlay) ! Surface Temperature and atmospheric temperature (K) 77 80 REAL, INTENT(IN) :: wstar(ngrid) ! Vertical velocity due to turbulence (m/s) 78 LOGICAL, INTENT(IN) :: virtual ! Boolean to account for vapor flottability79 81 REAL, INTENT(IN) :: mumean(ngrid) ! Molecular mass of the atmosphere (kg/mol) 80 82 REAL, INTENT(IN) :: pqvap(ngrid,nlay) ! Water vapor mixing ratio (kg/kg) to account for vapor flottability 81 83 REAL, INTENT(IN) :: pqsurf(ngrid,nslope) ! Water ice frost on the surface (kg/m^2) to account for vapor flottability 84 LOGICAL, INTENT(IN) :: write_outputs ! write_output in xios or not. 82 85 83 86 ! Outputs: … … 100 103 !$OMP THREADPRIVATE(karman,nu) 101 104 102 REAL rib(ngrid) ! Bulk Richardson number (1) 103 REAL fm(ngrid) ! stability function for momentum (1) 104 REAL fh(ngrid) ! stability function for heat (1) 105 REAL rib(ngrid,nslope) ! Bulk Richardson number [virtual or dry] (1) 106 REAL rib_dry(ngrid,nslope) ! Dry bulk Richardson number [virtual or dry] (1) 107 REAL fm(ngrid,nslope) ! stability function for momentum (1) 108 REAL fh(ngrid,nslope) ! stability function for heat (1) 105 109 106 110 ! Formalism for stability functions fm= 1/phim^2; fh = 1/(phim*phih) where … … 111 115 REAL betam, betah, alphah, bm, bh, lambda 112 116 113 REAL pz0t ! initial thermal roughness length z0t for the iterative scheme (m) 114 REAL ric_colaitis ! critical richardson number (1) 115 REAL reynolds(ngrid) ! Reynolds number (1) 116 REAL prandtl(ngrid) ! Prandtl number (1) 117 REAL pz0tcomp(ngrid) ! computed z0t during the iterative scheme (m) 118 REAL ite ! Number of iteration (1) 119 REAL itemax ! Maximum number of iteration (1) 120 REAL residual ! Residual between pz0t and pz0tcomp (m) 121 REAL tol_iter ! Tolerance for the residual to ensure convergence (1= 122 123 REAL zu2(ngrid) ! Near surface winds (m/s) 124 125 REAL cdn(ngrid) ! neutral momentum drag coefficient (1) 126 REAL chn(ngrid) ! neutral heat drag coefficient (1) 127 REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T 128 REAL z1,zcd0 ! Neutral roughness (m) and Cd/Ch coefficient when call richls is not called 129 REAL tsurf_v(ngrid,nslope) ! Virtual surface temperature, accounting for vapor flottability 130 REAL temp_v(ngrid) ! Potential virtual air temperature in the frist layer, accounting for vapor flottability 131 REAL mu_h2o ! Molecular mass of water vapor (kg/mol) 132 REAL tol_frost ! Tolerance to consider the effect of frost on the surface 133 REAL qsat(ngrid) ! saturated mixing ratio of water vapor 117 REAL pz0t ! initial thermal roughness length z0t for the iterative scheme (m) 118 REAL ric_colaitis ! critical richardson number (1) 119 REAL reynolds(ngrid,nslope) ! Reynolds number (1) 120 REAL prandtl(ngrid) ! Prandtl number (1) 121 REAL pz0tcomp(ngrid) ! computed z0t during the iterative scheme (m) 122 REAL z0t(ngrid,nslope) ! computed z0t at the last step the iterative scheme (m) 123 REAL ite ! Number of iteration (1) 124 REAL itemax ! Maximum number of iteration (1) 125 REAL residual ! Residual between pz0t and pz0tcomp (m) 126 REAL tol_iter ! Tolerance for the residual to ensure convergence (1= 127 128 REAL zu2(ngrid) ! Near surface winds (m/s) 129 130 REAL cdn(ngrid) ! neutral momentum drag coefficient (1) 131 REAL chn(ngrid) ! neutral heat drag coefficient (1) 132 REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T 133 REAL z1,zcd0 ! Neutral roughness (m) and Cd/Ch coefficient when call richls is not called 134 REAL tsurf_v(ngrid,nslope) ! Virtual surface temperature, accounting for vapor flottability 135 REAL temp_v(ngrid) ! Potential virtual air temperature in the frist layer, accounting for vapor flottability 136 REAL mu_h2o ! Molecular mass of water vapor (kg/mol) 137 REAL tol_frost ! Tolerance to consider the effect of frost on the surface 138 REAL qsat(ngrid) ! saturated mixing ratio of water vapor 134 139 135 REAL sm ! Stability function computed with the ATKE scheme136 REAL f_ri_cd_min ! minimum of the stability function with the ATKE scheme140 REAL sm ! Stability function computed with the ATKE scheme 141 REAL f_ri_cd_min ! minimum of the stability function with the ATKE scheme 137 142 138 143 ! Code: … … 144 149 mu_h2o = 18e-3 145 150 tol_frost = 1e-4 146 reynolds(: ) = 0.151 reynolds(:,:) = 0. 147 152 pz0t = 0. 148 153 pz0tcomp(:) = 0.1*pz0(:) 149 rib(:) = 0. 154 rib(:,:) = 0. 155 rib_dry(:,:) = 0. 150 156 pcdv(:,:) = 0. 151 157 pcdh(:,:) = 0. 158 z0t(:,:) = 0. 152 159 f_ri_cd_min = 0.01 153 160 ! this formulation assumes alphah=1., implying betah=betam … … 162 169 ric_colaitis = betah/(betam**2) 163 170 164 IF( virtual) then171 IF(include_waterbuoyancy) then 165 172 temp_v(:) = ph(:,1)*(1.+pqvap(:,1)/(mu_h2o/mumean(:)))/(1.+pqvap(:,1)) 166 173 DO islope = 1,nslope … … 210 217 IF(turb_resolved) zu2(ig)=MAX(zu2(ig),1.) 211 218 ! Richardson number formulation proposed by D.E. England et al. (1995) 212 rib(ig) = (pg/tsurf_v(ig,islope))*sqrt(pz(ig,1)*pz0(ig))*(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))*(temp_v(ig)-tsurf_v(ig,islope))/zu2(ig) 219 rib(ig,islope) = (pg/tsurf_v(ig,islope))*sqrt(pz(ig,1)*pz0(ig))*(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))*(temp_v(ig)-tsurf_v(ig,islope))/zu2(ig) 220 rib_dry(ig,islope) = (pg/pts(ig,islope))*sqrt(pz(ig,1)*pz0(ig))*(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))*(ph(ig,1)-pts(ig,islope))/zu2(ig) 213 221 ELSE 214 222 print*,'warning, infinite Richardson at surface' 215 223 print*,pu(ig,1),pv(ig,1) 216 rib(ig) = ric_colaitis 224 rib(ig,islope) = ric_colaitis 225 rib_dry(ig,islope) = ric_colaitis 217 226 ENDIF 218 227 … … 221 230 IF(callatke) THEN 222 231 ! Computation following parameterizaiton from ATKE 223 IF(rib(ig ) .gt. 0) THEN224 ! STABLE BOUNDARY LAYER : 225 sm = MAX(smmin,cn*(1.-rib(ig )/ric))232 IF(rib(ig,islope) .gt. 0) THEN 233 ! STABLE BOUNDARY LAYER :include_waterbuoyancy 234 sm = MAX(smmin,cn*(1.-rib(ig,islope)/ric)) 226 235 ! 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_slope228 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)236 prandtl(ig) = pr_neut*exp(-pr_slope/pr_neut*rib(ig,islope)+rib(ig,islope)/pr_neut) + rib(ig,islope) * pr_slope 237 fm(ig,islope) = max((sm**(3./2.) * sqrt(cepsilon) * (1 - rib(ig,islope) / prandtl(ig))**(1./2.)),f_ri_cd_min) 238 fh(ig,islope) = max((fm(ig,islope) / prandtl(ig)), f_ri_cd_min) 230 239 ELSE 231 240 ! UNSTABLE BOUNDARY LAYER : 232 sm = 2./rpi * (cinf - cn) * atan(-rib(ig )/ri0) + cn233 prandtl(ig) = -2./rpi * (pr_asym - pr_neut) * atan(rib(ig )/ri1) + pr_neut234 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)241 sm = 2./rpi * (cinf - cn) * atan(-rib(ig,islope)/ri0) + cn 242 prandtl(ig) = -2./rpi * (pr_asym - pr_neut) * atan(rib(ig,islope)/ri1) + pr_neut 243 fm(ig,islope) = MAX((sm**(3./2.) * sqrt(cepsilon) * (1 - rib(ig,islope) / prandtl(ig))**(1./2.)),f_ri_cd_min) 244 fh(ig,islope) = MAX((fm(ig,islope) / prandtl(ig)), f_ri_cd_min) 236 245 ENDIF ! Rib < 0 237 246 ELSE 238 247 ! Computation following parameterizaiton from from D.E. England et al. (95) 239 IF (rib(ig ) .gt. 0.) THEN248 IF (rib(ig,islope) .gt. 0.) THEN 240 249 ! STABLE BOUNDARY LAYER : 241 250 prandtl(ig) = 1. 242 IF(rib(ig ) .lt. ric_colaitis) THEN251 IF(rib(ig,islope) .lt. ric_colaitis) THEN 243 252 ! Assuming alphah=1. and bh=bm for stable conditions : 244 fm(ig ) = ((ric_colaitis-rib(ig))/ric_colaitis)**2245 fh(ig ) = ((ric_colaitis-rib(ig))/ric_colaitis)**2253 fm(ig,islope) = ((ric_colaitis-rib(ig,islope))/ric_colaitis)**2 254 fh(ig,islope) = ((ric_colaitis-rib(ig,islope))/ric_colaitis)**2 246 255 ELSE 247 256 ! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface 248 fm(ig ) = 1.249 fh(ig ) = 1.257 fm(ig,islope) = 1. 258 fh(ig,islope) = 1. 250 259 ENDIF 251 260 ELSE 252 261 ! 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.25255 prandtl(ig) = alphah*((1.-lambda*bm*rib(ig ))**0.25)/((1.-lambda*bh*rib(ig))**0.5)262 fm(ig,islope) = sqrt(1.-lambda*bm*rib(ig,islope)) 263 fh(ig,islope) = (1./alphah)*((1.-lambda*bh*rib(ig,islope))**0.5)*(1.-lambda*bm*rib(ig,islope))**0.25 264 prandtl(ig) = alphah*((1.-lambda*bm*rib(ig,islope))**0.25)/((1.-lambda*bh*rib(ig,islope))**0.5) 256 265 ENDIF ! Rib < 0 257 266 ENDIF ! atke 258 267 ! Recompute the reynolds and z0t 259 reynolds(ig ) = karman*sqrt(fm(ig))*sqrt(zu2(ig))*pz0(ig)/(log(z1z0)*nu)260 pz0tcomp(ig) = pz0(ig)*exp(-karman*7.3*(reynolds(ig )**0.25)*(prandtl(ig)**0.5)+5*karman)268 reynolds(ig,islope) = karman*sqrt(fm(ig,islope))*sqrt(zu2(ig))*pz0(ig)/(log(z1z0)*nu) 269 pz0tcomp(ig) = pz0(ig)*exp(-karman*7.3*(reynolds(ig,islope)**0.25)*(prandtl(ig)**0.5)+5*karman) 261 270 residual = abs(pz0t-pz0tcomp(ig)) 262 271 ite = ite+1 … … 264 273 265 274 ! Compute the coefficients Cdv; Cdh : neutral coefficient x stability functions 266 pcdv(ig,islope)=cdn(ig)*fm(ig) 267 pcdh(ig,islope)=chn(ig)*fh(ig) 275 pcdv(ig,islope)=cdn(ig)*fm(ig,islope) 276 pcdh(ig,islope)=chn(ig)*fh(ig,islope) 277 z0t(ig,islope) = pz0tcomp(ig) 268 278 pz0t = 0. ! for next grid point 269 279 ENDDO ! of ngrid 270 280 enddo 271 281 ENDIF !of if call richardson surface layer 282 283 IF(write_outputs) then 284 285 call write_output("rib_dry_vdif_cd","Dry Richardson number in vdif_cd","1",rib_dry(:,iflat)) 286 287 call write_output("rib_vdif_cd","Richardson number in vdif_cd","1",rib(:,iflat)) 288 289 call write_output("fm_vdif_cd","Momentum stability function in vdif_cd","1",fm(:,iflat)) 290 291 call write_output("fh_vdif_cd","Heat stability function in vdif_cd","1",fh(:,iflat)) 292 293 call write_output("z0t_vdif_cd","Thermal roughness length in vdif_cd","m",z0t(:,iflat)) 294 295 call write_output("z0_vdif_cd","Momentum roughness length in vdif_cd","m",pz0(:)) 296 297 call write_output("Reynolds_vdif_cd","Reynolds number in vdif_cd","1",reynolds(:,iflat)) 298 299 ENDIF 272 300 273 301 RETURN -
trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
r3295 r3325 235 235 LOGICAL :: writeoutput ! boolean to say to soilexchange.F if we are at the last iteration and thus if he can write in the diagsoil 236 236 REAL, INTENT(out) :: pore_icefraction(ngrid,nsoil,nslope) ! ice filling fraction in the pores 237 !! Water buyoncy238 LOGICAL :: virtual239 237 240 238 !! Subsurface ice interactions … … 309 307 zq1temp_regolith(1:ngrid)=0 310 308 zdqsdif_tot(1:ngrid)=0 311 virtual = .false.312 309 pore_icefraction(:,:,:) = 0. 313 310 zdqsdif_ssi_atm(:,:) = 0. … … 434 431 435 432 CALL vdif_cd(ngrid,nlay,nslope,pz0,g,pzlay,pplay,pu,pv,wstar, 436 & ptsrf,ph, virtual,mmean(:,1),zq(:,:,igcm_h2o_vap),437 & pqsurf(:,igcm_h2o_ice,:), 433 & ptsrf,ph,mmean(:,1),zq(:,:,igcm_h2o_vap), 434 & pqsurf(:,igcm_h2o_ice,:), .false., 438 435 & zcdv_true,zcdh_true) 439 436 … … 612 609 613 610 CALL vdif_cd(ngrid,nlay,nslope,pz0,g,pzlay,pplay,zu,zv,wstar, 614 & ptsrf,ph, virtual,mmean(:,1),zq(:,:,igcm_h2o_vap),615 & pqsurf(:,igcm_h2o_ice,:), 611 & ptsrf,ph,mmean(:,1),zq(:,:,igcm_h2o_vap), 612 & pqsurf(:,igcm_h2o_ice,:), .true., 616 613 & zcdv_true,zcdh_true) 617 614
Note: See TracChangeset
for help on using the changeset viewer.
