| 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 turb_mod, only: turb_resolved |
|---|
| 21 | use watersat_mod, only: watersat |
|---|
| 22 | use lmdz_atke_turbulence_ini, only : smmin, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut, ri0, ri1, cn, rpi |
|---|
| 23 | use paleoclimate_mod, only: include_waterbuoyancy |
|---|
| 24 | use write_output_mod, only: write_output |
|---|
| 25 | use comslope_mod, ONLY: iflat |
|---|
| 26 | use callkeys_mod, only: callrichsl, callatke |
|---|
| 27 | IMPLICIT NONE |
|---|
| 28 | |
|---|
| 29 | !======================================================================= |
|---|
| 30 | ! |
|---|
| 31 | ! Subject: computation of the surface drag coefficient using the |
|---|
| 32 | ! ------- approch developed by Loui for ECMWF. |
|---|
| 33 | ! |
|---|
| 34 | ! Author: Frederic Hourdin 15 /10 /93 |
|---|
| 35 | ! Modified by : Arnaud Colaitis 03/08/11; rewritten by LL to switch to F90 |
|---|
| 36 | ! ------- |
|---|
| 37 | ! |
|---|
| 38 | ! Arguments: |
|---|
| 39 | ! ---------- |
|---|
| 40 | ! |
|---|
| 41 | ! inputs: |
|---|
| 42 | ! ------ |
|---|
| 43 | ! ngrid size of the horizontal grid |
|---|
| 44 | ! pg gravity (m s -2) |
|---|
| 45 | ! pz(ngrid,nlay) height of layers |
|---|
| 46 | ! pp(ngrid,nlay) pressure of layers |
|---|
| 47 | ! pu(ngrid,nlay) u component of the wind |
|---|
| 48 | ! pv(ngrid,nlay) v component of the wind |
|---|
| 49 | ! pts(ngrid) surface temperature |
|---|
| 50 | ! ph(ngrid) potential temperature T*(p/ps)^kappa |
|---|
| 51 | ! mumean Molecular mass of the atmosphere (kg/mol) |
|---|
| 52 | ! pqvap(ngrid,nlay) Water vapor mixing ratio (kg/kg) to account for vapor flottability |
|---|
| 53 | ! pqsurf(ngrid) Water ice frost on the surface(kg/m^2) to account for vapor flottability |
|---|
| 54 | ! |
|---|
| 55 | ! outputs: |
|---|
| 56 | ! -------- |
|---|
| 57 | ! pcdv(ngrid) Cd for the wind |
|---|
| 58 | ! pcdh(ngrid) Cd for potential temperature |
|---|
| 59 | ! |
|---|
| 60 | !======================================================================= |
|---|
| 61 | |
|---|
| 62 | |
|---|
| 63 | !----------------------------------------------------------------------- |
|---|
| 64 | ! Declarations: |
|---|
| 65 | ! ------------- |
|---|
| 66 | |
|---|
| 67 | |
|---|
| 68 | ! Arguments: |
|---|
| 69 | |
|---|
| 70 | ! Inputs: |
|---|
| 71 | ! ---------- |
|---|
| 72 | INTEGER, INTENT(IN) :: ngrid,nlay,nslope ! Number of points in the physical grid and atmospheric grid |
|---|
| 73 | REAL, INTENT(IN) :: pz0(ngrid) ! Surface Roughness (m) |
|---|
| 74 | REAL, INTENT(IN) :: pg ! Mars gravity (m/s^2) |
|---|
| 75 | REAL, INTENT(IN) :: pz(ngrid,nlay),pp(ngrid,nlay) ! Layers pseudo altitudes (m) and pressure (Pa) |
|---|
| 76 | REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay) ! Zonal and Meriditionnal winds (m/s) |
|---|
| 77 | REAL, INTENT(IN) :: pts(ngrid,nslope),ph(ngrid,nlay) ! Surface Temperature and atmospheric temperature (K) |
|---|
| 78 | REAL, INTENT(IN) :: wstar(ngrid) ! Vertical velocity due to turbulence (m/s) |
|---|
| 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 | LOGICAL, INTENT(IN) :: write_outputs ! write_output in xios or not. |
|---|
| 83 | |
|---|
| 84 | ! Outputs: |
|---|
| 85 | ! ---------- |
|---|
| 86 | REAL, INTENT(OUT) :: pcdv(ngrid,nslope) ! momentum drag coefficient (1) |
|---|
| 87 | REAL, INTENT(OUT) :: pcdh(ngrid,nslope) ! heat drag coefficient (1) |
|---|
| 88 | |
|---|
| 89 | ! Local: |
|---|
| 90 | ! ------ |
|---|
| 91 | |
|---|
| 92 | INTEGER ig,islope ! Loop variable |
|---|
| 93 | |
|---|
| 94 | REAL karman,nu ! Von Karman constant and fluid kinematic viscosity |
|---|
| 95 | |
|---|
| 96 | LOGICAL firstcal |
|---|
| 97 | DATA karman,nu/.41,0.001/ |
|---|
| 98 | DATA firstcal/.true./ |
|---|
| 99 | SAVE karman,nu |
|---|
| 100 | |
|---|
| 101 | !$OMP THREADPRIVATE(karman,nu) |
|---|
| 102 | |
|---|
| 103 | REAL rib(ngrid,nslope) ! Bulk Richardson number [virtual or dry] (1) |
|---|
| 104 | REAL rib_dry(ngrid,nslope) ! Dry bulk Richardson number [virtual or dry] (1) |
|---|
| 105 | REAL fm(ngrid,nslope) ! stability function for momentum (1) |
|---|
| 106 | REAL fh(ngrid,nslope) ! stability function for heat (1) |
|---|
| 107 | |
|---|
| 108 | ! Formalism for stability functions fm= 1/phim^2; fh = 1/(phim*phih) where |
|---|
| 109 | ! phim = 1+betam*zeta or (1-bm*zeta)**am |
|---|
| 110 | ! phih = alphah + betah*zeta or alphah(1.-bh*zeta)**ah |
|---|
| 111 | ! ah and am are assumed to be -0.25 and -0.5 respectively |
|---|
| 112 | ! lambda is an intermediate variable to simplify the expression |
|---|
| 113 | REAL betam, betah, alphah, bm, bh, lambda |
|---|
| 114 | |
|---|
| 115 | REAL pz0t ! initial thermal roughness length z0t for the iterative scheme (m) |
|---|
| 116 | REAL ric_colaitis ! critical richardson number (1) |
|---|
| 117 | REAL reynolds(ngrid,nslope) ! Reynolds number (1) |
|---|
| 118 | REAL prandtl(ngrid) ! Prandtl number (1) |
|---|
| 119 | REAL pz0tcomp(ngrid) ! computed z0t during the iterative scheme (m) |
|---|
| 120 | REAL z0t(ngrid,nslope) ! computed z0t at the last step the iterative scheme (m) |
|---|
| 121 | REAL ite ! Number of iteration (1) |
|---|
| 122 | REAL itemax ! Maximum number of iteration (1) |
|---|
| 123 | REAL residual ! Residual between pz0t and pz0tcomp (m) |
|---|
| 124 | REAL tol_iter ! Tolerance for the residual to ensure convergence (1= |
|---|
| 125 | |
|---|
| 126 | REAL zu2(ngrid) ! Near surface winds (m/s) |
|---|
| 127 | |
|---|
| 128 | REAL cdn(ngrid) ! neutral momentum drag coefficient (1) |
|---|
| 129 | REAL chn(ngrid) ! neutral heat drag coefficient (1) |
|---|
| 130 | REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T |
|---|
| 131 | REAL z1,zcd0 ! Neutral roughness (m) and Cd/Ch coefficient when call richls is not called |
|---|
| 132 | REAL tsurf_v(ngrid,nslope) ! Virtual surface temperature, accounting for vapor flottability |
|---|
| 133 | REAL temp_v(ngrid) ! Potential virtual air temperature in the frist layer, accounting for vapor flottability |
|---|
| 134 | REAL mu_h2o ! Molecular mass of water vapor (kg/mol) |
|---|
| 135 | REAL tol_frost ! Tolerance to consider the effect of frost on the surface |
|---|
| 136 | REAL qsat(ngrid) ! saturated mixing ratio of water vapor |
|---|
| 137 | |
|---|
| 138 | REAL sm ! Stability function computed with the ATKE scheme |
|---|
| 139 | REAL f_ri_cd_min ! minimum of the stability function with the ATKE scheme |
|---|
| 140 | |
|---|
| 141 | ! Code: |
|---|
| 142 | ! -------- |
|---|
| 143 | |
|---|
| 144 | ! Initialisation |
|---|
| 145 | itemax= 10 |
|---|
| 146 | tol_iter = 0.01 |
|---|
| 147 | mu_h2o = 18e-3 |
|---|
| 148 | tol_frost = 1e-4 |
|---|
| 149 | reynolds(:,:) = 0. |
|---|
| 150 | pz0t = 0. |
|---|
| 151 | pz0tcomp(:) = 0.1*pz0(:) |
|---|
| 152 | rib(:,:) = 0. |
|---|
| 153 | rib_dry(:,:) = 0. |
|---|
| 154 | pcdv(:,:) = 0. |
|---|
| 155 | pcdh(:,:) = 0. |
|---|
| 156 | z0t(:,:) = 0. |
|---|
| 157 | fm(:,:) = 0. |
|---|
| 158 | fh(:,:) = 0. |
|---|
| 159 | f_ri_cd_min = 0.01 |
|---|
| 160 | ! this formulation assumes alphah=1., implying betah=betam |
|---|
| 161 | ! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers : |
|---|
| 162 | bm = 16. |
|---|
| 163 | bh = 16. |
|---|
| 164 | alphah = 1. |
|---|
| 165 | betam = 5. |
|---|
| 166 | betah = 5. |
|---|
| 167 | lambda = (sqrt(bh/bm))/alphah |
|---|
| 168 | |
|---|
| 169 | ric_colaitis = betah/(betam**2) |
|---|
| 170 | |
|---|
| 171 | IF(include_waterbuoyancy) then |
|---|
| 172 | temp_v(:) = ph(:,1)*(1.+pqvap(:,1)/(mu_h2o/mumean(:)))/(1.+pqvap(:,1)) |
|---|
| 173 | DO islope = 1,nslope |
|---|
| 174 | DO ig = 1,ngrid |
|---|
| 175 | IF(pqsurf(ig,islope).gt.tol_frost) then |
|---|
| 176 | call watersat(1,pts(ig,islope),pp(ig,1),qsat(ig)) |
|---|
| 177 | tsurf_v(ig,islope) = pts(ig,islope)*(1.+qsat(ig)/(mu_h2o/mumean(ig)))/(1.+qsat(ig)) |
|---|
| 178 | ELSE |
|---|
| 179 | tsurf_v(ig,islope) = pts(ig,islope)*(1.+pqvap(ig,1)/(mu_h2o/mumean(ig)))/(1.+pqvap(ig,1)) |
|---|
| 180 | ENDIF |
|---|
| 181 | ENDDO |
|---|
| 182 | ENDDO |
|---|
| 183 | ELSE |
|---|
| 184 | tsurf_v(:,:) = pts(:,:) |
|---|
| 185 | temp_v(:) = ph(:,1) |
|---|
| 186 | ENDIF |
|---|
| 187 | |
|---|
| 188 | IF(.not.callrichsl) then |
|---|
| 189 | ! Original formulation as in LMDZ Earth: Cd = Ch = (kappa/(ln(1+z1/z0)))^2 |
|---|
| 190 | ! NB: In forget et al., 1999, it's Cd = Ch = (kappa/(ln(z1/z0)))^2 |
|---|
| 191 | DO ig = 1,ngrid |
|---|
| 192 | z1 = 1.E+0 + pz(ig,1)/pz0(ig) |
|---|
| 193 | zcd0 = karman/log(z1) |
|---|
| 194 | zcd0 = zcd0*zcd0 |
|---|
| 195 | DO islope = 1,nslope |
|---|
| 196 | pcdv(ig,islope) = zcd0 |
|---|
| 197 | pcdh(ig,islope) = zcd0 |
|---|
| 198 | ENDDO |
|---|
| 199 | ENDDO |
|---|
| 200 | ELSE |
|---|
| 201 | ! We follow the parameterization from Colaitis et al., 2013 (supplementary material) |
|---|
| 202 | DO islope = 1,nslope |
|---|
| 203 | DO ig=1,ngrid |
|---|
| 204 | ite = 0. |
|---|
| 205 | residual = 100*tol_iter*pz0(ig) |
|---|
| 206 | z1z0=pz(ig,1)/pz0(ig) |
|---|
| 207 | cdn(ig)=karman/log(z1z0) |
|---|
| 208 | cdn(ig)=cdn(ig)*cdn(ig) |
|---|
| 209 | DO WHILE((residual .gt. tol_iter*pz0(ig)) .and. (ite .lt. itemax)) |
|---|
| 210 | ! Computations of z0T; iterated until z0T converges |
|---|
| 211 | pz0t = pz0tcomp(ig) |
|---|
| 212 | z1z0t=pz(ig,1)/pz0t |
|---|
| 213 | chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t) |
|---|
| 214 | IF ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then |
|---|
| 215 | 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 |
|---|
| 216 | IF(turb_resolved) zu2(ig)=MAX(zu2(ig),1.) |
|---|
| 217 | ! Richardson number formulation proposed by D.E. England et al. (1995) |
|---|
| 218 | 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) |
|---|
| 219 | 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) |
|---|
| 220 | ELSE |
|---|
| 221 | print*,'warning, infinite Richardson at surface' |
|---|
| 222 | print*,pu(ig,1),pv(ig,1) |
|---|
| 223 | rib(ig,islope) = ric_colaitis |
|---|
| 224 | rib_dry(ig,islope) = ric_colaitis |
|---|
| 225 | ENDIF |
|---|
| 226 | ! Compute the stability functions fm; fh depending on the stability of the surface layer |
|---|
| 227 | |
|---|
| 228 | IF(callatke) THEN |
|---|
| 229 | ! Computation following parameterizaiton from ATKE |
|---|
| 230 | IF(rib(ig,islope) .gt. 0) THEN |
|---|
| 231 | ! STABLE BOUNDARY LAYER :include_waterbuoyancy |
|---|
| 232 | sm = MAX(smmin,cn*(1.-rib(ig,islope)/ric)) |
|---|
| 233 | ! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019 |
|---|
| 234 | prandtl(ig) = pr_neut*exp(-pr_slope/pr_neut*rib(ig,islope)+rib(ig,islope)/pr_neut) + rib(ig,islope) * pr_slope |
|---|
| 235 | fm(ig,islope) = max((sm**(3./2.) * sqrt(cepsilon) * (1 - rib(ig,islope) / prandtl(ig))**(1./2.)),f_ri_cd_min) |
|---|
| 236 | fh(ig,islope) = max((fm(ig,islope) / prandtl(ig)), f_ri_cd_min) |
|---|
| 237 | ELSE |
|---|
| 238 | ! UNSTABLE BOUNDARY LAYER : |
|---|
| 239 | sm = 2./rpi * (cinf - cn) * atan(-rib(ig,islope)/ri0) + cn |
|---|
| 240 | prandtl(ig) = -2./rpi * (pr_asym - pr_neut) * atan(rib(ig,islope)/ri1) + pr_neut |
|---|
| 241 | fm(ig,islope) = MAX((sm**(3./2.) * sqrt(cepsilon) * (1 - rib(ig,islope) / prandtl(ig))**(1./2.)),f_ri_cd_min) |
|---|
| 242 | fh(ig,islope) = MAX((fm(ig,islope) / prandtl(ig)), f_ri_cd_min) |
|---|
| 243 | ENDIF ! Rib < 0 |
|---|
| 244 | ELSE |
|---|
| 245 | ! Computation following parameterizaiton from from D.E. England et al. (95) |
|---|
| 246 | IF (rib(ig,islope) .gt. 0.) THEN |
|---|
| 247 | ! STABLE BOUNDARY LAYER : |
|---|
| 248 | prandtl(ig) = 1. |
|---|
| 249 | IF(rib(ig,islope) .lt. ric_colaitis) THEN |
|---|
| 250 | ! Assuming alphah=1. and bh=bm for stable conditions : |
|---|
| 251 | fm(ig,islope) = ((ric_colaitis-rib(ig,islope))/ric_colaitis)**2 |
|---|
| 252 | fh(ig,islope) = ((ric_colaitis-rib(ig,islope))/ric_colaitis)**2 |
|---|
| 253 | ELSE |
|---|
| 254 | ! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface |
|---|
| 255 | fm(ig,islope) = 1. |
|---|
| 256 | fh(ig,islope) = 1. |
|---|
| 257 | ENDIF |
|---|
| 258 | ELSE |
|---|
| 259 | ! UNSTABLE BOUNDARY LAYER : |
|---|
| 260 | fm(ig,islope) = sqrt(1.-lambda*bm*rib(ig,islope)) |
|---|
| 261 | fh(ig,islope) = (1./alphah)*((1.-lambda*bh*rib(ig,islope))**0.5)*(1.-lambda*bm*rib(ig,islope))**0.25 |
|---|
| 262 | prandtl(ig) = alphah*((1.-lambda*bm*rib(ig,islope))**0.25)/((1.-lambda*bh*rib(ig,islope))**0.5) |
|---|
| 263 | ENDIF ! Rib < 0 |
|---|
| 264 | ENDIF ! atke |
|---|
| 265 | ! Recompute the reynolds and z0t |
|---|
| 266 | reynolds(ig,islope) = karman*sqrt(fm(ig,islope))*sqrt(zu2(ig))*pz0(ig)/(log(z1z0)*nu) |
|---|
| 267 | pz0tcomp(ig) = pz0(ig)*exp(-karman*7.3*(reynolds(ig,islope)**0.25)*(prandtl(ig)**0.5)+5*karman) |
|---|
| 268 | residual = abs(pz0t-pz0tcomp(ig)) |
|---|
| 269 | ite = ite+1 |
|---|
| 270 | ENDDO ! of while |
|---|
| 271 | |
|---|
| 272 | ! Compute the coefficients Cdv; Cdh : neutral coefficient x stability functions |
|---|
| 273 | pcdv(ig,islope)=cdn(ig)*fm(ig,islope) |
|---|
| 274 | pcdh(ig,islope)=chn(ig)*fh(ig,islope) |
|---|
| 275 | z0t(ig,islope) = pz0tcomp(ig) |
|---|
| 276 | pz0t = 0. ! for next grid point |
|---|
| 277 | ENDDO ! of ngrid |
|---|
| 278 | enddo |
|---|
| 279 | ENDIF !of if call richardson surface layer |
|---|
| 280 | |
|---|
| 281 | IF(write_outputs) then |
|---|
| 282 | |
|---|
| 283 | call write_output("rib_dry_vdif_cd","Dry Richardson number in vdif_cd","1",rib_dry(:,iflat)) |
|---|
| 284 | |
|---|
| 285 | call write_output("rib_vdif_cd","Richardson number in vdif_cd","1",rib(:,iflat)) |
|---|
| 286 | |
|---|
| 287 | call write_output("fm_vdif_cd","Momentum stability function in vdif_cd","1",fm(:,iflat)) |
|---|
| 288 | |
|---|
| 289 | call write_output("fh_vdif_cd","Heat stability function in vdif_cd","1",fh(:,iflat)) |
|---|
| 290 | |
|---|
| 291 | call write_output("z0t_vdif_cd","Thermal roughness length in vdif_cd","m",z0t(:,iflat)) |
|---|
| 292 | |
|---|
| 293 | call write_output("z0_vdif_cd","Momentum roughness length in vdif_cd","m",pz0(:)) |
|---|
| 294 | |
|---|
| 295 | call write_output("Reynolds_vdif_cd","Reynolds number in vdif_cd","1",reynolds(:,iflat)) |
|---|
| 296 | |
|---|
| 297 | ENDIF |
|---|
| 298 | |
|---|
| 299 | END SUBROUTINE vdif_cd |
|---|
| 300 | |
|---|
| 301 | END MODULE vdif_cd_mod |
|---|
| 302 | |
|---|
| 303 | |
|---|
| 304 | |
|---|
| 305 | |
|---|
| 306 | |
|---|
| 307 | |
|---|
| 308 | |
|---|
| 309 | |
|---|
| 310 | |
|---|
| 311 | |
|---|
| 312 | |
|---|
| 313 | |
|---|
