| 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,virtual,mumean,pqvap,pqsurf,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, cinf, cepsilon, pr_slope, pr_asym, pr_neut, ri0, ri1, cn, rpi |
|---|
| 24 | |
|---|
| 25 | IMPLICIT NONE |
|---|
| 26 | include "callkeys.h" |
|---|
| 27 | !======================================================================= |
|---|
| 28 | ! |
|---|
| 29 | ! Subject: computation of the surface drag coefficient using the |
|---|
| 30 | ! ------- approch developed by Loui for ECMWF. |
|---|
| 31 | ! |
|---|
| 32 | ! Author: Frederic Hourdin 15 /10 /93 |
|---|
| 33 | ! Modified by : Arnaud Colaitis 03/08/11; rewritten by LL to switch to F90 |
|---|
| 34 | ! ------- |
|---|
| 35 | ! |
|---|
| 36 | ! Arguments: |
|---|
| 37 | ! ---------- |
|---|
| 38 | ! |
|---|
| 39 | ! inputs: |
|---|
| 40 | ! ------ |
|---|
| 41 | ! ngrid size of the horizontal grid |
|---|
| 42 | ! pg gravity (m s -2) |
|---|
| 43 | ! pz(ngrid,nlay) height of layers |
|---|
| 44 | ! pp(ngrid,nlay) pressure of layers |
|---|
| 45 | ! pu(ngrid,nlay) u component of the wind |
|---|
| 46 | ! pv(ngrid,nlay) v component of the wind |
|---|
| 47 | ! pts(ngrid) surface temperature |
|---|
| 48 | ! ph(ngrid) potential temperature T*(p/ps)^kappa |
|---|
| 49 | ! virtual Boolean to account for vapor flottability |
|---|
| 50 | ! mumean Molecular mass of the atmosphere (kg/mol) |
|---|
| 51 | ! pqvap(ngrid,nlay) Water vapor mixing ratio (kg/kg) to account for vapor flottability |
|---|
| 52 | ! pqsurf(ngrid) Water ice frost on the surface(kg/m^2) to account for vapor flottability |
|---|
| 53 | ! |
|---|
| 54 | ! outputs: |
|---|
| 55 | ! -------- |
|---|
| 56 | ! pcdv(ngrid) Cd for the wind |
|---|
| 57 | ! pcdh(ngrid) Cd for potential temperature |
|---|
| 58 | ! |
|---|
| 59 | !======================================================================= |
|---|
| 60 | |
|---|
| 61 | |
|---|
| 62 | !----------------------------------------------------------------------- |
|---|
| 63 | ! Declarations: |
|---|
| 64 | ! ------------- |
|---|
| 65 | |
|---|
| 66 | |
|---|
| 67 | ! Arguments: |
|---|
| 68 | |
|---|
| 69 | ! Inputs: |
|---|
| 70 | ! ---------- |
|---|
| 71 | INTEGER, INTENT(IN) :: ngrid,nlay,nslope ! Number of points in the physical grid and atmospheric grid |
|---|
| 72 | REAL, INTENT(IN) :: pz0(ngrid) ! Surface Roughness (m) |
|---|
| 73 | REAL, INTENT(IN) :: pg ! Mars gravity (m/s^2) |
|---|
| 74 | REAL, INTENT(IN) :: pz(ngrid,nlay),pp(ngrid,nlay) ! Layers pseudo altitudes (m) and pressure (Pa) |
|---|
| 75 | REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay) ! Zonal and Meriditionnal winds (m/s) |
|---|
| 76 | REAL, INTENT(IN) :: pts(ngrid,nslope),ph(ngrid,nlay) ! Surface Temperature and atmospheric temperature (K) |
|---|
| 77 | REAL, INTENT(IN) :: wstar(ngrid) ! Vertical velocity due to turbulence (m/s) |
|---|
| 78 | LOGICAL, INTENT(IN) :: virtual ! Boolean to account for vapor flottability |
|---|
| 79 | REAL, INTENT(IN) :: mumean(ngrid) ! Molecular mass of the atmosphere (kg/mol) |
|---|
| 80 | REAL, INTENT(IN) :: pqvap(ngrid,nlay) ! Water vapor mixing ratio (kg/kg) to account for vapor flottability |
|---|
| 81 | REAL, INTENT(IN) :: pqsurf(ngrid,nslope) ! Water ice frost on the surface (kg/m^2) to account for vapor flottability |
|---|
| 82 | |
|---|
| 83 | ! Outputs: |
|---|
| 84 | ! ---------- |
|---|
| 85 | REAL, INTENT(OUT) :: pcdv(ngrid,nslope) ! momentum drag coefficient (1) |
|---|
| 86 | REAL, INTENT(OUT) :: pcdh(ngrid,nslope) ! heat drag coefficient (1) |
|---|
| 87 | |
|---|
| 88 | ! Local: |
|---|
| 89 | ! ------ |
|---|
| 90 | |
|---|
| 91 | INTEGER ig,islope ! Loop variable |
|---|
| 92 | |
|---|
| 93 | REAL karman,nu ! Von Karman constant and fluid kinematic viscosity |
|---|
| 94 | |
|---|
| 95 | LOGICAL firstcal |
|---|
| 96 | DATA karman,nu/.41,0.001/ |
|---|
| 97 | DATA firstcal/.true./ |
|---|
| 98 | SAVE karman,nu |
|---|
| 99 | |
|---|
| 100 | !$OMP THREADPRIVATE(karman,nu) |
|---|
| 101 | |
|---|
| 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 | |
|---|
| 106 | ! Formalism for stability functions fm= 1/phim^2; fh = 1/(phim*phih) where |
|---|
| 107 | ! phim = 1+betam*zeta or (1-bm*zeta)**am |
|---|
| 108 | ! phih = alphah + betah*zeta or alphah(1.-bh*zeta)**ah |
|---|
| 109 | ! ah and am are assumed to be -0.25 and -0.5 respectively |
|---|
| 110 | ! lambda is an intermediate variable to simplify the expression |
|---|
| 111 | REAL betam, betah, alphah, bm, bh, lambda |
|---|
| 112 | |
|---|
| 113 | REAL pz0t ! initial thermal roughness length z0t for the iterative scheme (m) |
|---|
| 114 | REAL ric ! 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 |
|---|
| 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 | |
|---|
| 138 | ! Code: |
|---|
| 139 | ! -------- |
|---|
| 140 | |
|---|
| 141 | ! Initialisation |
|---|
| 142 | itemax= 10 |
|---|
| 143 | tol_iter = 0.01 |
|---|
| 144 | mu_h2o = 18e-3 |
|---|
| 145 | tol_frost = 1e-4 |
|---|
| 146 | reynolds(:) = 0. |
|---|
| 147 | pz0t = 0. |
|---|
| 148 | pz0tcomp(:) = 0.1*pz0(:) |
|---|
| 149 | rib(:) = 0. |
|---|
| 150 | pcdv(:,:) = 0. |
|---|
| 151 | pcdh(:,:) = 0. |
|---|
| 152 | f_ri_cd_min = 0.01 |
|---|
| 153 | ! this formulation assumes alphah=1., implying betah=betam |
|---|
| 154 | ! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers : |
|---|
| 155 | bm = 16. |
|---|
| 156 | bh = 16. |
|---|
| 157 | alphah = 1. |
|---|
| 158 | betam = 5. |
|---|
| 159 | betah = 5. |
|---|
| 160 | lambda = (sqrt(bh/bm))/alphah |
|---|
| 161 | |
|---|
| 162 | ric = betah/(betam**2) |
|---|
| 163 | |
|---|
| 164 | IF(virtual) then |
|---|
| 165 | temp_v(:) = ph(:,1)*(1.+pqvap(:,1)/(mu_h2o/mumean(:)))/(1.+pqvap(:,1)) |
|---|
| 166 | DO islope = 1,nslope |
|---|
| 167 | DO ig = 1,ngrid |
|---|
| 168 | IF(pqsurf(ig,islope).gt.tol_frost) then |
|---|
| 169 | call watersat(1,pts(ig,islope),pp(ig,1),qsat(ig)) |
|---|
| 170 | tsurf_v(ig,islope) = pts(ig,islope)*(1.+qsat(ig)/(mu_h2o/mumean(ig)))/(1.+qsat(ig)) |
|---|
| 171 | ELSE |
|---|
| 172 | tsurf_v(ig,islope) = pts(ig,islope)*(1.+pqvap(ig,1)/(mu_h2o/mumean(ig)))/(1.+pqvap(ig,1)) |
|---|
| 173 | ENDIF |
|---|
| 174 | ENDDO |
|---|
| 175 | ENDDO |
|---|
| 176 | ELSE |
|---|
| 177 | tsurf_v(:,:) = pts(:,:) |
|---|
| 178 | temp_v(:) = ph(:,1) |
|---|
| 179 | ENDIF |
|---|
| 180 | |
|---|
| 181 | IF(.not.callrichsl) then |
|---|
| 182 | ! Original formulation as in LMDZ Earth: Cd = Ch = (kappa/(ln(1+z1/z0)))^2 |
|---|
| 183 | ! NB: In forget et al., 1999, it's Cd = Ch = (kappa/(ln(z1/z0)))^2 |
|---|
| 184 | DO ig = 1,ngrid |
|---|
| 185 | z1 = 1.E+0 + pz(ig,1)/pz0(ig) |
|---|
| 186 | zcd0 = karman/log(z1) |
|---|
| 187 | zcd0 = zcd0*zcd0 |
|---|
| 188 | DO islope = 1,nslope |
|---|
| 189 | pcdv(ig,islope) = zcd0 |
|---|
| 190 | pcdh(ig,islope) = zcd0 |
|---|
| 191 | ENDDO |
|---|
| 192 | ENDDO |
|---|
| 193 | ELSE |
|---|
| 194 | ! We follow the parameterization from Colaitis et al., 2013 (supplementary material) |
|---|
| 195 | DO islope = 1,nslope |
|---|
| 196 | DO ig=1,ngrid |
|---|
| 197 | ite = 0. |
|---|
| 198 | residual = abs(pz0tcomp(ig)-pz0t) |
|---|
| 199 | z1z0=pz(ig,1)/pz0(ig) |
|---|
| 200 | cdn(ig)=karman/log(z1z0) |
|---|
| 201 | cdn(ig)=cdn(ig)*cdn(ig) |
|---|
| 202 | |
|---|
| 203 | DO WHILE((residual .gt. tol_iter*pz0(ig)) .and. (ite .lt. itemax)) |
|---|
| 204 | ! Computations of z0T; iterated until z0T converges |
|---|
| 205 | pz0t = pz0tcomp(ig) |
|---|
| 206 | z1z0t=pz(ig,1)/pz0t |
|---|
| 207 | chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t) |
|---|
| 208 | IF ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then |
|---|
| 209 | 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 |
|---|
| 210 | IF(turb_resolved) zu2(ig)=MAX(zu2(ig),1.) |
|---|
| 211 | ! 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) |
|---|
| 213 | ELSE |
|---|
| 214 | print*,'warning, infinite Richardson at surface' |
|---|
| 215 | print*,pu(ig,1),pv(ig,1) |
|---|
| 216 | rib(ig) = ric |
|---|
| 217 | ENDIF |
|---|
| 218 | |
|---|
| 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 |
|---|
| 243 | ! Assuming alphah=1. and bh=bm for stable conditions : |
|---|
| 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 |
|---|
| 251 | ELSE |
|---|
| 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 |
|---|
| 258 | ! 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) |
|---|
| 261 | residual = abs(pz0t-pz0tcomp(ig)) |
|---|
| 262 | ite = ite+1 |
|---|
| 263 | ENDDO ! of while |
|---|
| 264 | |
|---|
| 265 | ! 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) |
|---|
| 268 | pz0t = 0. ! for next grid point |
|---|
| 269 | ENDDO ! of ngrid |
|---|
| 270 | enddo |
|---|
| 271 | ENDIF !of if call richardson surface layer |
|---|
| 272 | |
|---|
| 273 | RETURN |
|---|
| 274 | |
|---|
| 275 | END SUBROUTINE vdif_cd |
|---|
| 276 | |
|---|
| 277 | END MODULE vdif_cd_mod |
|---|
| 278 | |
|---|
| 279 | |
|---|
| 280 | |
|---|
| 281 | |
|---|
| 282 | |
|---|
| 283 | |
|---|
| 284 | |
|---|
| 285 | |
|---|
| 286 | |
|---|
| 287 | |
|---|
