| 1 | MODULE vdif_cd_mod |
|---|
| 2 | |
|---|
| 3 | !======================================================================================================================! |
|---|
| 4 | ! Subject: |
|---|
| 5 | !--------- |
|---|
| 6 | ! Module used to compute the exchange coefficient in the surface layers Cd; Ch |
|---|
| 7 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 8 | ! Reference: |
|---|
| 9 | !----------- |
|---|
| 10 | ! Colaïtis, A., Spiga, A., Hourdin, F., Rio, C., Forget, F., and Millour, E. (2013), A thermal plume model for the Martian convective boundary layer, J. Geophys. Res. Planets, 118, 1468–1487, doi:10.1002/jgre.20104. |
|---|
| 11 | ! |
|---|
| 12 | !======================================================================================================================! |
|---|
| 13 | |
|---|
| 14 | IMPLICIT NONE |
|---|
| 15 | |
|---|
| 16 | CONTAINS |
|---|
| 17 | |
|---|
| 18 | SUBROUTINE vdif_cd(ngrid,nlay,nslope,pz0,pg,pz,pp,pu,pv,wstar,pts,ph,mumean,pqvap,pqsurf,write_outputs,pcdv,pcdh) |
|---|
| 19 | |
|---|
| 20 | use comcstfi_h |
|---|
| 21 | use turb_mod, only: turb_resolved |
|---|
| 22 | use watersat_mod, only: watersat |
|---|
| 23 | use lmdz_atke_turbulence_ini, only : smmin, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut, ri0, ri1, cn, rpi |
|---|
| 24 | use paleoclimate_mod, only: include_waterbuoyancy |
|---|
| 25 | use write_output_mod, only: write_output |
|---|
| 26 | use comslope_mod, ONLY: iflat |
|---|
| 27 | IMPLICIT NONE |
|---|
| 28 | include "callkeys.h" |
|---|
| 29 | |
|---|
| 30 | |
|---|
| 31 | !======================================================================= |
|---|
| 32 | ! |
|---|
| 33 | ! Subject: computation of the surface drag coefficient using the |
|---|
| 34 | ! ------- approch developed by Loui for ECMWF. |
|---|
| 35 | ! |
|---|
| 36 | ! Author: Frederic Hourdin 15 /10 /93 |
|---|
| 37 | ! Modified by : Arnaud Colaitis 03/08/11; rewritten by LL to switch to F90 |
|---|
| 38 | ! ------- |
|---|
| 39 | ! |
|---|
| 40 | ! Arguments: |
|---|
| 41 | ! ---------- |
|---|
| 42 | ! |
|---|
| 43 | ! inputs: |
|---|
| 44 | ! ------ |
|---|
| 45 | ! ngrid size of the horizontal grid |
|---|
| 46 | ! pg gravity (m s -2) |
|---|
| 47 | ! pz(ngrid,nlay) height of layers |
|---|
| 48 | ! pp(ngrid,nlay) pressure of layers |
|---|
| 49 | ! pu(ngrid,nlay) u component of the wind |
|---|
| 50 | ! pv(ngrid,nlay) v component of the wind |
|---|
| 51 | ! pts(ngrid) surface temperature |
|---|
| 52 | ! ph(ngrid) potential temperature T*(p/ps)^kappa |
|---|
| 53 | ! mumean Molecular mass of the atmosphere (kg/mol) |
|---|
| 54 | ! pqvap(ngrid,nlay) Water vapor mixing ratio (kg/kg) to account for vapor flottability |
|---|
| 55 | ! pqsurf(ngrid) Water ice frost on the surface(kg/m^2) to account for vapor flottability |
|---|
| 56 | ! |
|---|
| 57 | ! outputs: |
|---|
| 58 | ! -------- |
|---|
| 59 | ! pcdv(ngrid) Cd for the wind |
|---|
| 60 | ! pcdh(ngrid) Cd for potential temperature |
|---|
| 61 | ! |
|---|
| 62 | !======================================================================= |
|---|
| 63 | |
|---|
| 64 | |
|---|
| 65 | !----------------------------------------------------------------------- |
|---|
| 66 | ! Declarations: |
|---|
| 67 | ! ------------- |
|---|
| 68 | |
|---|
| 69 | |
|---|
| 70 | ! Arguments: |
|---|
| 71 | |
|---|
| 72 | ! Inputs: |
|---|
| 73 | ! ---------- |
|---|
| 74 | INTEGER, INTENT(IN) :: ngrid,nlay,nslope ! Number of points in the physical grid and atmospheric grid |
|---|
| 75 | REAL, INTENT(IN) :: pz0(ngrid) ! Surface Roughness (m) |
|---|
| 76 | REAL, INTENT(IN) :: pg ! Mars gravity (m/s^2) |
|---|
| 77 | REAL, INTENT(IN) :: pz(ngrid,nlay),pp(ngrid,nlay) ! Layers pseudo altitudes (m) and pressure (Pa) |
|---|
| 78 | REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay) ! Zonal and Meriditionnal winds (m/s) |
|---|
| 79 | REAL, INTENT(IN) :: pts(ngrid,nslope),ph(ngrid,nlay) ! Surface Temperature and atmospheric temperature (K) |
|---|
| 80 | REAL, INTENT(IN) :: wstar(ngrid) ! Vertical velocity due to turbulence (m/s) |
|---|
| 81 | REAL, INTENT(IN) :: mumean(ngrid) ! Molecular mass of the atmosphere (kg/mol) |
|---|
| 82 | REAL, INTENT(IN) :: pqvap(ngrid,nlay) ! Water vapor mixing ratio (kg/kg) to account for vapor flottability |
|---|
| 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. |
|---|
| 85 | |
|---|
| 86 | ! Outputs: |
|---|
| 87 | ! ---------- |
|---|
| 88 | REAL, INTENT(OUT) :: pcdv(ngrid,nslope) ! momentum drag coefficient (1) |
|---|
| 89 | REAL, INTENT(OUT) :: pcdh(ngrid,nslope) ! heat drag coefficient (1) |
|---|
| 90 | |
|---|
| 91 | ! Local: |
|---|
| 92 | ! ------ |
|---|
| 93 | |
|---|
| 94 | INTEGER ig,islope ! Loop variable |
|---|
| 95 | |
|---|
| 96 | REAL karman,nu ! Von Karman constant and fluid kinematic viscosity |
|---|
| 97 | |
|---|
| 98 | LOGICAL firstcal |
|---|
| 99 | DATA karman,nu/.41,0.001/ |
|---|
| 100 | DATA firstcal/.true./ |
|---|
| 101 | SAVE karman,nu |
|---|
| 102 | |
|---|
| 103 | !$OMP THREADPRIVATE(karman,nu) |
|---|
| 104 | |
|---|
| 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) |
|---|
| 109 | |
|---|
| 110 | ! Formalism for stability functions fm= 1/phim^2; fh = 1/(phim*phih) where |
|---|
| 111 | ! phim = 1+betam*zeta or (1-bm*zeta)**am |
|---|
| 112 | ! phih = alphah + betah*zeta or alphah(1.-bh*zeta)**ah |
|---|
| 113 | ! ah and am are assumed to be -0.25 and -0.5 respectively |
|---|
| 114 | ! lambda is an intermediate variable to simplify the expression |
|---|
| 115 | REAL betam, betah, alphah, bm, bh, lambda |
|---|
| 116 | |
|---|
| 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 |
|---|
| 139 | |
|---|
| 140 | REAL sm ! Stability function computed with the ATKE scheme |
|---|
| 141 | REAL f_ri_cd_min ! minimum of the stability function with the ATKE scheme |
|---|
| 142 | |
|---|
| 143 | ! Code: |
|---|
| 144 | ! -------- |
|---|
| 145 | |
|---|
| 146 | ! Initialisation |
|---|
| 147 | itemax= 10 |
|---|
| 148 | tol_iter = 0.01 |
|---|
| 149 | mu_h2o = 18e-3 |
|---|
| 150 | tol_frost = 1e-4 |
|---|
| 151 | reynolds(:,:) = 0. |
|---|
| 152 | pz0t = 0. |
|---|
| 153 | pz0tcomp(:) = 0.1*pz0(:) |
|---|
| 154 | rib(:,:) = 0. |
|---|
| 155 | rib_dry(:,:) = 0. |
|---|
| 156 | pcdv(:,:) = 0. |
|---|
| 157 | pcdh(:,:) = 0. |
|---|
| 158 | z0t(:,:) = 0. |
|---|
| 159 | fm(:,:) = 0. |
|---|
| 160 | fh(:,:) = 0. |
|---|
| 161 | f_ri_cd_min = 0.01 |
|---|
| 162 | ! this formulation assumes alphah=1., implying betah=betam |
|---|
| 163 | ! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers : |
|---|
| 164 | bm = 16. |
|---|
| 165 | bh = 16. |
|---|
| 166 | alphah = 1. |
|---|
| 167 | betam = 5. |
|---|
| 168 | betah = 5. |
|---|
| 169 | lambda = (sqrt(bh/bm))/alphah |
|---|
| 170 | |
|---|
| 171 | ric_colaitis = betah/(betam**2) |
|---|
| 172 | |
|---|
| 173 | IF(include_waterbuoyancy) then |
|---|
| 174 | temp_v(:) = ph(:,1)*(1.+pqvap(:,1)/(mu_h2o/mumean(:)))/(1.+pqvap(:,1)) |
|---|
| 175 | DO islope = 1,nslope |
|---|
| 176 | DO ig = 1,ngrid |
|---|
| 177 | IF(pqsurf(ig,islope).gt.tol_frost) then |
|---|
| 178 | call watersat(1,pts(ig,islope),pp(ig,1),qsat(ig)) |
|---|
| 179 | tsurf_v(ig,islope) = pts(ig,islope)*(1.+qsat(ig)/(mu_h2o/mumean(ig)))/(1.+qsat(ig)) |
|---|
| 180 | ELSE |
|---|
| 181 | tsurf_v(ig,islope) = pts(ig,islope)*(1.+pqvap(ig,1)/(mu_h2o/mumean(ig)))/(1.+pqvap(ig,1)) |
|---|
| 182 | ENDIF |
|---|
| 183 | ENDDO |
|---|
| 184 | ENDDO |
|---|
| 185 | ELSE |
|---|
| 186 | tsurf_v(:,:) = pts(:,:) |
|---|
| 187 | temp_v(:) = ph(:,1) |
|---|
| 188 | ENDIF |
|---|
| 189 | |
|---|
| 190 | IF(.not.callrichsl) then |
|---|
| 191 | ! Original formulation as in LMDZ Earth: Cd = Ch = (kappa/(ln(1+z1/z0)))^2 |
|---|
| 192 | ! NB: In forget et al., 1999, it's Cd = Ch = (kappa/(ln(z1/z0)))^2 |
|---|
| 193 | DO ig = 1,ngrid |
|---|
| 194 | z1 = 1.E+0 + pz(ig,1)/pz0(ig) |
|---|
| 195 | zcd0 = karman/log(z1) |
|---|
| 196 | zcd0 = zcd0*zcd0 |
|---|
| 197 | DO islope = 1,nslope |
|---|
| 198 | pcdv(ig,islope) = zcd0 |
|---|
| 199 | pcdh(ig,islope) = zcd0 |
|---|
| 200 | ENDDO |
|---|
| 201 | ENDDO |
|---|
| 202 | ELSE |
|---|
| 203 | ! We follow the parameterization from Colaitis et al., 2013 (supplementary material) |
|---|
| 204 | DO islope = 1,nslope |
|---|
| 205 | DO ig=1,ngrid |
|---|
| 206 | ite = 0. |
|---|
| 207 | residual = 100*tol_iter*pz0(ig) |
|---|
| 208 | z1z0=pz(ig,1)/pz0(ig) |
|---|
| 209 | cdn(ig)=karman/log(z1z0) |
|---|
| 210 | cdn(ig)=cdn(ig)*cdn(ig) |
|---|
| 211 | DO WHILE((residual .gt. tol_iter*pz0(ig)) .and. (ite .lt. itemax)) |
|---|
| 212 | ! Computations of z0T; iterated until z0T converges |
|---|
| 213 | pz0t = pz0tcomp(ig) |
|---|
| 214 | z1z0t=pz(ig,1)/pz0t |
|---|
| 215 | chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t) |
|---|
| 216 | IF ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then |
|---|
| 217 | zu2(ig) = pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) + (log(1.+0.7*wstar(ig) + 2.3*wstar(ig)**2))**2 ! Near surface winds account for buoyancy |
|---|
| 218 | IF(turb_resolved) zu2(ig)=MAX(zu2(ig),1.) |
|---|
| 219 | ! Richardson number formulation proposed by D.E. England et al. (1995) |
|---|
| 220 | 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) |
|---|
| 221 | 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) |
|---|
| 222 | ELSE |
|---|
| 223 | print*,'warning, infinite Richardson at surface' |
|---|
| 224 | print*,pu(ig,1),pv(ig,1) |
|---|
| 225 | rib(ig,islope) = ric_colaitis |
|---|
| 226 | rib_dry(ig,islope) = ric_colaitis |
|---|
| 227 | ENDIF |
|---|
| 228 | ! Compute the stability functions fm; fh depending on the stability of the surface layer |
|---|
| 229 | |
|---|
| 230 | IF(callatke) THEN |
|---|
| 231 | ! Computation following parameterizaiton from ATKE |
|---|
| 232 | IF(rib(ig,islope) .gt. 0) THEN |
|---|
| 233 | ! STABLE BOUNDARY LAYER :include_waterbuoyancy |
|---|
| 234 | sm = MAX(smmin,cn*(1.-rib(ig,islope)/ric)) |
|---|
| 235 | ! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019 |
|---|
| 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) |
|---|
| 239 | ELSE |
|---|
| 240 | ! UNSTABLE BOUNDARY LAYER : |
|---|
| 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) |
|---|
| 245 | ENDIF ! Rib < 0 |
|---|
| 246 | ELSE |
|---|
| 247 | ! Computation following parameterizaiton from from D.E. England et al. (95) |
|---|
| 248 | IF (rib(ig,islope) .gt. 0.) THEN |
|---|
| 249 | ! STABLE BOUNDARY LAYER : |
|---|
| 250 | prandtl(ig) = 1. |
|---|
| 251 | IF(rib(ig,islope) .lt. ric_colaitis) THEN |
|---|
| 252 | ! Assuming alphah=1. and bh=bm for stable conditions : |
|---|
| 253 | fm(ig,islope) = ((ric_colaitis-rib(ig,islope))/ric_colaitis)**2 |
|---|
| 254 | fh(ig,islope) = ((ric_colaitis-rib(ig,islope))/ric_colaitis)**2 |
|---|
| 255 | ELSE |
|---|
| 256 | ! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface |
|---|
| 257 | fm(ig,islope) = 1. |
|---|
| 258 | fh(ig,islope) = 1. |
|---|
| 259 | ENDIF |
|---|
| 260 | ELSE |
|---|
| 261 | ! UNSTABLE BOUNDARY LAYER : |
|---|
| 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) |
|---|
| 265 | ENDIF ! Rib < 0 |
|---|
| 266 | ENDIF ! atke |
|---|
| 267 | ! Recompute the reynolds and z0t |
|---|
| 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) |
|---|
| 270 | residual = abs(pz0t-pz0tcomp(ig)) |
|---|
| 271 | ite = ite+1 |
|---|
| 272 | ENDDO ! of while |
|---|
| 273 | |
|---|
| 274 | ! Compute the coefficients Cdv; Cdh : neutral coefficient x stability functions |
|---|
| 275 | pcdv(ig,islope)=cdn(ig)*fm(ig,islope) |
|---|
| 276 | pcdh(ig,islope)=chn(ig)*fh(ig,islope) |
|---|
| 277 | z0t(ig,islope) = pz0tcomp(ig) |
|---|
| 278 | pz0t = 0. ! for next grid point |
|---|
| 279 | ENDDO ! of ngrid |
|---|
| 280 | enddo |
|---|
| 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 |
|---|
| 300 | |
|---|
| 301 | RETURN |
|---|
| 302 | |
|---|
| 303 | END SUBROUTINE vdif_cd |
|---|
| 304 | |
|---|
| 305 | END MODULE vdif_cd_mod |
|---|
| 306 | |
|---|
| 307 | |
|---|
| 308 | |
|---|
| 309 | |
|---|
| 310 | |
|---|
| 311 | |
|---|
| 312 | |
|---|
| 313 | |
|---|
| 314 | |
|---|
| 315 | |
|---|
| 316 | |
|---|
| 317 | |
|---|
