[4687] | 1 | module lmdz_atke_exchange_coeff |
---|
[4449] | 2 | |
---|
| 3 | implicit none |
---|
| 4 | |
---|
| 5 | contains |
---|
| 6 | |
---|
[4631] | 7 | subroutine atke_compute_km_kh(ngrid,nlay,dtime, & |
---|
[4653] | 8 | wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv, & |
---|
[4478] | 9 | tke,Km_out,Kh_out) |
---|
[4449] | 10 | |
---|
| 11 | !======================================================================== |
---|
| 12 | ! Routine that computes turbulent Km / Kh coefficients with a |
---|
| 13 | ! 1.5 order closure scheme (TKE) with or without stationarity assumption |
---|
| 14 | ! |
---|
| 15 | ! This parameterization has been constructed in the framework of a |
---|
| 16 | ! collective and collaborative workshop, |
---|
| 17 | ! the so-called 'Atelier TKE (ATKE)' with |
---|
[4478] | 18 | ! K. Arjdal, L. Raillard, C. Dehondt, P. Tiengou, A. Spiga, F. Cheruy, T Dubos, |
---|
[4449] | 19 | ! M. Coulon-Decorzens, S. Fromang, G. Riviere, A. Sima, F. Hourdin, E. Vignon |
---|
| 20 | ! |
---|
| 21 | ! Main assumptions of the model : |
---|
[4714] | 22 | ! (1) horizontal homogeneity (Dx=Dy=0.) |
---|
[4449] | 23 | !======================================================================= |
---|
| 24 | |
---|
| 25 | |
---|
| 26 | |
---|
[4687] | 27 | USE lmdz_atke_turbulence_ini, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, atke_ok_virtual |
---|
| 28 | USE lmdz_atke_turbulence_ini, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes,rg, rd, rv, atke_ok_vdiff |
---|
| 29 | USE lmdz_atke_turbulence_ini, ONLY : viscom, viscoh, clmix, clmixshear, iflag_atke_lmix, lmin, smmin |
---|
[4449] | 30 | |
---|
| 31 | implicit none |
---|
| 32 | |
---|
| 33 | |
---|
| 34 | ! Declarations: |
---|
| 35 | !============= |
---|
| 36 | |
---|
| 37 | INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid) |
---|
[4631] | 38 | INTEGER, INTENT(IN) :: nlay ! number of vertical index |
---|
[4449] | 39 | |
---|
[4631] | 40 | REAL, INTENT(IN) :: dtime ! physics time step (s) |
---|
[4449] | 41 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: wind_u ! zonal velocity (m/s) |
---|
| 42 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: wind_v ! meridional velocity (m/s) |
---|
| 43 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: temp ! temperature (K) |
---|
[4653] | 44 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: qvap ! specific humidity (kg/kg) |
---|
[4449] | 45 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: play ! pressure (Pa) |
---|
| 46 | REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: pinterf ! pressure at interfaces(Pa) |
---|
[4644] | 47 | REAL, DIMENSION(ngrid), INTENT(IN) :: cdrag_uv ! surface drag coefficient for momentum |
---|
[4449] | 48 | |
---|
| 49 | REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers |
---|
| 50 | |
---|
[4478] | 51 | REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Km_out ! output: Exchange coefficient for momentum at interface between layers |
---|
| 52 | REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Kh_out ! output: Exchange coefficient for heat flux at interface between layers |
---|
[4449] | 53 | |
---|
| 54 | ! Local variables |
---|
[4478] | 55 | REAL, DIMENSION(ngrid,nlay+1) :: Km ! Exchange coefficient for momentum at interface between layers |
---|
| 56 | REAL, DIMENSION(ngrid,nlay+1) :: Kh ! Exchange coefficient for heat flux at interface between layers |
---|
| 57 | REAL, DIMENSION(ngrid,nlay) :: theta ! Potential temperature |
---|
| 58 | REAL, DIMENSION(ngrid,nlay+1) :: l_exchange ! Length of exchange (at interface) |
---|
| 59 | REAL, DIMENSION(ngrid,nlay+1) :: z_interf ! Altitude at the interface |
---|
| 60 | REAL, DIMENSION(ngrid,nlay) :: z_lay ! Altitude of layers |
---|
| 61 | REAL, DIMENSION(ngrid,nlay) :: dz_interf ! distance between two consecutive interfaces |
---|
| 62 | REAL, DIMENSION(ngrid,nlay) :: dz_lay ! distance between two layer middles (NB: first and last are half layers) |
---|
[4631] | 63 | REAL, DIMENSION(ngrid,nlay+1) :: N2 ! square of Brunt Vaisala pulsation (at interface) |
---|
| 64 | REAL, DIMENSION(ngrid,nlay+1) :: shear2 ! square of wind shear (at interface) |
---|
[4478] | 65 | REAL, DIMENSION(ngrid,nlay+1) :: Ri ! Richardson's number (at interface) |
---|
| 66 | REAL, DIMENSION(ngrid,nlay+1) :: Prandtl ! Turbulent Prandtl's number (at interface) |
---|
| 67 | REAL, DIMENSION(ngrid,nlay+1) :: Sm ! Stability function for momentum (at interface) |
---|
| 68 | REAL, DIMENSION(ngrid,nlay+1) :: Sh ! Stability function for heat (at interface) |
---|
[4449] | 69 | |
---|
| 70 | INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid) |
---|
[4481] | 71 | REAL :: cn,Ri0,Ri1 ! parameter for Sm stability function and Prandlt |
---|
[4478] | 72 | REAL :: preff ! reference pressure for potential temperature calculations |
---|
| 73 | REAL :: thetam ! mean potential temperature at interface |
---|
[4631] | 74 | REAL :: delta ! discriminant of the second order polynomial |
---|
| 75 | REAL :: qq ! tke=qq**2/2 |
---|
| 76 | REAL :: shear ! wind shear |
---|
| 77 | REAL :: lstrat ! mixing length depending on local stratification |
---|
| 78 | REAL :: taustrat ! caracteristic timescale for turbulence in very stable conditions |
---|
[4644] | 79 | REAL :: netloss ! net loss term of tke |
---|
| 80 | REAL :: netsource ! net source term of tke |
---|
| 81 | REAL :: ustar ! friction velocity estimation |
---|
[4653] | 82 | REAL :: invtau |
---|
| 83 | REAL :: rvap |
---|
[4449] | 84 | |
---|
| 85 | ! Initializations: |
---|
| 86 | !================ |
---|
| 87 | |
---|
| 88 | DO igrid=1,ngrid |
---|
[4478] | 89 | dz_interf(igrid,1) = 0.0 |
---|
[4449] | 90 | z_interf(igrid,1) = 0.0 |
---|
| 91 | END DO |
---|
| 92 | |
---|
[4653] | 93 | ! Calculation of potential temperature: (if vapor -> virtual potential temperature) |
---|
[4478] | 94 | !===================================== |
---|
[4449] | 95 | |
---|
[4478] | 96 | preff=100000. |
---|
[4653] | 97 | ! results should not depend on the choice of preff |
---|
[4478] | 98 | DO ilay=1,nlay |
---|
| 99 | DO igrid = 1, ngrid |
---|
| 100 | theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd) |
---|
| 101 | END DO |
---|
| 102 | END DO |
---|
[4449] | 103 | |
---|
[4653] | 104 | ! account for water vapor mass for buoyancy calculation |
---|
| 105 | IF (atke_ok_virtual) THEN |
---|
| 106 | DO ilay=1,nlay |
---|
| 107 | DO igrid = 1, ngrid |
---|
| 108 | rvap=max(0.,qvap(igrid,ilay)/(1.-qvap(igrid,ilay))) |
---|
| 109 | theta(igrid,ilay)=theta(igrid,ilay)*(1.+rvap/(RD/RV))/(1.+rvap) |
---|
| 110 | END DO |
---|
| 111 | END DO |
---|
| 112 | ENDIF |
---|
[4449] | 113 | |
---|
[4478] | 114 | |
---|
| 115 | ! Calculation of altitude of layers' middle and bottom interfaces: |
---|
| 116 | !================================================================= |
---|
| 117 | |
---|
[4449] | 118 | DO ilay=2,nlay+1 |
---|
| 119 | DO igrid=1,ngrid |
---|
[4478] | 120 | dz_interf(igrid,ilay-1) = rd*temp(igrid,ilay-1)/rg/play(igrid,ilay-1)*(pinterf(igrid,ilay-1)-pinterf(igrid,ilay)) |
---|
| 121 | z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1) |
---|
[4449] | 122 | ENDDO |
---|
| 123 | ENDDO |
---|
| 124 | |
---|
[4478] | 125 | DO ilay=1,nlay |
---|
| 126 | DO igrid=1,ngrid |
---|
| 127 | z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay)) |
---|
| 128 | ENDDO |
---|
| 129 | ENDDO |
---|
[4449] | 130 | |
---|
[4478] | 131 | |
---|
[4449] | 132 | ! Computes the gradient Richardson's number and stability functions: |
---|
| 133 | !=================================================================== |
---|
| 134 | |
---|
[4481] | 135 | ! calculation of cn = Sm value at Ri=0 |
---|
| 136 | ! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0 |
---|
| 137 | cn=(1./sqrt(cepsilon))**(2/3) |
---|
| 138 | ! calculation of Ri0 such that continuity in slope of Sm at Ri=0 |
---|
| 139 | Ri0=2./rpi*(cinf - cn)*ric/cn |
---|
| 140 | ! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0 |
---|
[4688] | 141 | Ri1 = -2./rpi * (pr_asym - pr_neut) |
---|
[4481] | 142 | |
---|
| 143 | |
---|
[4478] | 144 | DO ilay=2,nlay |
---|
[4449] | 145 | DO igrid=1,ngrid |
---|
[4478] | 146 | dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1) |
---|
| 147 | thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1)) |
---|
[4631] | 148 | N2(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay) |
---|
| 149 | shear2(igrid,ilay)= (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & |
---|
| 150 | ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 ) |
---|
| 151 | Ri(igrid,ilay) = N2(igrid,ilay) / MAX(shear2(igrid,ilay),1E-10) |
---|
[4481] | 152 | |
---|
| 153 | IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases |
---|
[4478] | 154 | Sm(igrid,ilay) = 2./rpi * (cinf-cn) * atan(-Ri(igrid,ilay)/Ri0) + cn |
---|
| 155 | Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut |
---|
[4481] | 156 | ELSE ! stable cases |
---|
[4644] | 157 | Sm(igrid,ilay) = max(smmin,cn*(1.-Ri(igrid,ilay)/Ric)) |
---|
[4688] | 158 | ! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019 |
---|
| 159 | Prandtl(igrid,ilay) = pr_neut*exp(-pr_slope/pr_neut*Ri(igrid,ilay)+Ri(igrid,ilay)/pr_neut) & |
---|
| 160 | + Ri(igrid,ilay) * pr_slope |
---|
[4481] | 161 | IF (Ri(igrid,ilay) .GE. Prandtl(igrid,ilay)) THEN |
---|
| 162 | call abort_physic("atke_compute_km_kh", & |
---|
| 163 | 'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1) |
---|
| 164 | ENDIF |
---|
[4449] | 165 | END IF |
---|
| 166 | |
---|
| 167 | Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay) |
---|
| 168 | |
---|
| 169 | ENDDO |
---|
| 170 | ENDDO |
---|
| 171 | |
---|
[4631] | 172 | |
---|
| 173 | ! Computing the mixing length: |
---|
| 174 | !============================================================== |
---|
| 175 | |
---|
| 176 | |
---|
| 177 | IF (iflag_atke_lmix .EQ. 1 ) THEN |
---|
| 178 | |
---|
| 179 | DO ilay=2,nlay |
---|
| 180 | DO igrid=1,ngrid |
---|
| 181 | l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) |
---|
| 182 | IF (N2(igrid,ilay) .GT. 0.) THEN |
---|
| 183 | lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)) |
---|
[4663] | 184 | lstrat=max(lstrat,lmin) |
---|
[4688] | 185 | !Inverse interpolation, Van de Wiel et al. 2010 |
---|
[4687] | 186 | l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0) |
---|
[4631] | 187 | ENDIF |
---|
| 188 | ENDDO |
---|
| 189 | ENDDO |
---|
| 190 | |
---|
[4632] | 191 | ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN |
---|
[4663] | 192 | ! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms |
---|
[4632] | 193 | DO ilay=2,nlay |
---|
| 194 | DO igrid=1,ngrid |
---|
| 195 | l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) |
---|
[4663] | 196 | IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN |
---|
| 197 | lstrat=min(clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)), & |
---|
| 198 | clmixshear*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay))) |
---|
| 199 | lstrat=max(lstrat,lmin) |
---|
[4688] | 200 | !Inverse interpolation, Van de Wiel et al. 2010 |
---|
[4687] | 201 | l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0) |
---|
[4632] | 202 | ENDIF |
---|
| 203 | ENDDO |
---|
| 204 | ENDDO |
---|
| 205 | |
---|
[4714] | 206 | ELSE IF (iflag_atke_lmix .EQ. 3 ) THEN |
---|
| 207 | ! add effect of wind shear on lstrat following grisogono 2010, qjrms |
---|
| 208 | DO ilay=2,nlay |
---|
| 209 | DO igrid=1,ngrid |
---|
| 210 | l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) |
---|
| 211 | IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN |
---|
| 212 | lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay))*(1.0+Ri(igrid,ilay)/(2.*Prandtl(igrid,ilay))) |
---|
| 213 | lstrat=max(lstrat,lmin) |
---|
| 214 | !Inverse interpolation, Van de Wiel et al. 2010 |
---|
| 215 | l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0) |
---|
| 216 | ENDIF |
---|
| 217 | ENDDO |
---|
| 218 | ENDDO |
---|
[4632] | 219 | |
---|
| 220 | |
---|
[4714] | 221 | |
---|
[4631] | 222 | ELSE |
---|
[4632] | 223 | ! default: neglect effect of local stratification and shear |
---|
[4631] | 224 | |
---|
| 225 | DO ilay=2,nlay+1 |
---|
| 226 | DO igrid=1,ngrid |
---|
| 227 | l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) |
---|
| 228 | ENDDO |
---|
| 229 | |
---|
| 230 | ENDDO |
---|
| 231 | ENDIF |
---|
| 232 | |
---|
| 233 | |
---|
[4644] | 234 | ! Computing the TKE k>=2: |
---|
| 235 | !======================== |
---|
[4449] | 236 | IF (iflag_atke == 0) THEN |
---|
| 237 | |
---|
[4644] | 238 | ! stationary solution (dtke/dt=0) |
---|
| 239 | |
---|
[4631] | 240 | DO ilay=2,nlay |
---|
[4449] | 241 | DO igrid=1,ngrid |
---|
| 242 | tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * & |
---|
[4631] | 243 | shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)) |
---|
[4449] | 244 | ENDDO |
---|
| 245 | ENDDO |
---|
| 246 | |
---|
[4631] | 247 | ELSE IF (iflag_atke == 1) THEN |
---|
| 248 | |
---|
| 249 | ! full implicit scheme resolved with a second order polynomial equation |
---|
| 250 | |
---|
| 251 | DO ilay=2,nlay |
---|
| 252 | DO igrid=1,ngrid |
---|
[4644] | 253 | qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) |
---|
[4675] | 254 | delta=(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime)**2. & |
---|
| 255 | +4.*(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime*qq + & |
---|
| 256 | 2.*l_exchange(igrid,ilay)*l_exchange(igrid,ilay)*cepsilon*Sm(igrid,ilay) & |
---|
| 257 | *shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))) |
---|
| 258 | qq=(-2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime + sqrt(delta))/2. |
---|
| 259 | qq=max(0.,qq) |
---|
[4631] | 260 | tke(igrid,ilay)=0.5*(qq**2) |
---|
| 261 | ENDDO |
---|
| 262 | ENDDO |
---|
| 263 | |
---|
[4644] | 264 | |
---|
[4631] | 265 | ELSE IF (iflag_atke == 2) THEN |
---|
| 266 | |
---|
[4644] | 267 | ! semi implicit scheme when l does not depend on tke |
---|
| 268 | ! positive-guaranteed if pr slope in stable condition >1 |
---|
| 269 | |
---|
| 270 | DO ilay=2,nlay |
---|
| 271 | DO igrid=1,ngrid |
---|
| 272 | qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) |
---|
| 273 | qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.) & |
---|
| 274 | *shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) & |
---|
| 275 | /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) |
---|
| 276 | tke(igrid,ilay)=0.5*(qq**2) |
---|
| 277 | ENDDO |
---|
| 278 | ENDDO |
---|
| 279 | |
---|
| 280 | |
---|
| 281 | ELSE IF (iflag_atke == 3) THEN |
---|
| 282 | ! numerical resolution adapted from that in MAR (Deleersnijder 1992) |
---|
| 283 | ! positively defined by construction |
---|
| 284 | |
---|
[4631] | 285 | DO ilay=2,nlay |
---|
| 286 | DO igrid=1,ngrid |
---|
[4644] | 287 | qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) |
---|
| 288 | IF (Ri(igrid,ilay) .LT. 0.) THEN |
---|
| 289 | netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)) |
---|
| 290 | netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) |
---|
[4631] | 291 | ELSE |
---|
[4644] | 292 | netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))+ & |
---|
| 293 | l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*N2(igrid,ilay)/Prandtl(igrid,ilay) |
---|
| 294 | netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay) |
---|
[4631] | 295 | ENDIF |
---|
[4644] | 296 | qq=((qq**2)/dtime+qq*netsource)/(qq/dtime+netloss) |
---|
[4631] | 297 | tke(igrid,ilay)=0.5*(qq**2) |
---|
| 298 | ENDDO |
---|
| 299 | ENDDO |
---|
| 300 | |
---|
[4644] | 301 | ELSE IF (iflag_atke == 4) THEN |
---|
| 302 | ! semi implicit scheme from Arpege (V. Masson methodology with |
---|
| 303 | ! Taylor expansion of the dissipation term) |
---|
[4631] | 304 | DO ilay=2,nlay |
---|
| 305 | DO igrid=1,ngrid |
---|
[4644] | 306 | qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) |
---|
| 307 | qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) & |
---|
| 308 | +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) & |
---|
| 309 | /(1.+2.*qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) |
---|
| 310 | qq=max(0.,qq) |
---|
| 311 | tke(igrid,ilay)=0.5*(qq**2) |
---|
[4631] | 312 | ENDDO |
---|
| 313 | ENDDO |
---|
| 314 | |
---|
| 315 | |
---|
| 316 | ELSE |
---|
[4463] | 317 | call abort_physic("atke_compute_km_kh", & |
---|
[4631] | 318 | 'numerical treatment of TKE not possible yet', 1) |
---|
[4449] | 319 | |
---|
| 320 | END IF |
---|
| 321 | |
---|
[4644] | 322 | ! We impose a 0 tke at nlay+1 |
---|
| 323 | !============================== |
---|
[4449] | 324 | |
---|
[4644] | 325 | DO igrid=1,ngrid |
---|
| 326 | tke(igrid,nlay+1)=0. |
---|
| 327 | END DO |
---|
| 328 | |
---|
| 329 | |
---|
| 330 | ! Calculation of surface TKE (k=1) |
---|
| 331 | !================================= |
---|
| 332 | ! surface TKE calculation inspired from what is done in Arpege (see E. Bazile note) |
---|
| 333 | DO igrid=1,ngrid |
---|
| 334 | ustar=sqrt(cdrag_uv(igrid)*(wind_u(igrid,1)**2+wind_v(igrid,1)**2)) |
---|
| 335 | tke(igrid,1)=ctkes*(ustar**2) |
---|
| 336 | END DO |
---|
| 337 | |
---|
| 338 | |
---|
| 339 | ! vertical diffusion of TKE |
---|
| 340 | !========================== |
---|
[4653] | 341 | IF (atke_ok_vdiff) THEN |
---|
[4644] | 342 | CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke) |
---|
| 343 | ENDIF |
---|
| 344 | |
---|
| 345 | |
---|
[4449] | 346 | ! Computing eddy diffusivity coefficients: |
---|
| 347 | !======================================== |
---|
[4478] | 348 | DO ilay=2,nlay ! TODO: also calculate for nlay+1 ? |
---|
[4449] | 349 | DO igrid=1,ngrid |
---|
[4478] | 350 | ! we add the molecular viscosity to Km,h |
---|
| 351 | Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5 |
---|
| 352 | Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5 |
---|
[4449] | 353 | END DO |
---|
| 354 | END DO |
---|
| 355 | |
---|
[4478] | 356 | ! for output: |
---|
| 357 | !=========== |
---|
| 358 | Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay) |
---|
| 359 | Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay) |
---|
[4449] | 360 | |
---|
| 361 | end subroutine atke_compute_km_kh |
---|
| 362 | |
---|
[4644] | 363 | !=============================================================================================== |
---|
| 364 | subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke) |
---|
[4449] | 365 | |
---|
[4644] | 366 | ! routine that computes the vertical diffusion of TKE by the turbulence |
---|
[4653] | 367 | ! using an implicit resolution (See note by Dufresne and Ghattas (2009)) |
---|
[4644] | 368 | ! E Vignon, July 2023 |
---|
| 369 | |
---|
[4687] | 370 | USE lmdz_atke_turbulence_ini, ONLY : rd, cke, viscom |
---|
[4644] | 371 | |
---|
| 372 | |
---|
| 373 | INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid) |
---|
| 374 | INTEGER, INTENT(IN) :: nlay ! number of vertical index |
---|
| 375 | |
---|
| 376 | REAL, INTENT(IN) :: dtime ! physics time step (s) |
---|
| 377 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: z_lay ! altitude of mid-layers (m) |
---|
| 378 | REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: z_interf ! altitude of bottom interfaces (m) |
---|
| 379 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: temp ! temperature (K) |
---|
| 380 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: play ! pressure (Pa) |
---|
| 381 | REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: l_exchange ! mixing length at interfaces between layers |
---|
| 382 | REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: Sm ! stability function for eddy diffusivity for momentum at interface between layers |
---|
| 383 | |
---|
| 384 | REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers |
---|
| 385 | |
---|
| 386 | |
---|
| 387 | |
---|
| 388 | INTEGER :: igrid,ilay |
---|
| 389 | REAL, DIMENSION(ngrid,nlay+1) :: Ke ! eddy diffusivity for TKE |
---|
| 390 | REAL, DIMENSION(ngrid,nlay+1) :: dtke |
---|
| 391 | REAL, DIMENSION(ngrid,nlay+1) :: ak, bk, ck, CCK, DDK |
---|
| 392 | REAL :: gammak,Kem,KKb,KKt |
---|
| 393 | |
---|
| 394 | |
---|
| 395 | ! Few initialisations |
---|
| 396 | CCK(:,:)=0. |
---|
| 397 | DDK(:,:)=0. |
---|
| 398 | dtke(:,:)=0. |
---|
| 399 | |
---|
| 400 | |
---|
| 401 | ! Eddy diffusivity for TKE |
---|
| 402 | |
---|
| 403 | DO ilay=2,nlay |
---|
| 404 | DO igrid=1,ngrid |
---|
| 405 | Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay)))*cke |
---|
| 406 | ENDDO |
---|
| 407 | ENDDO |
---|
| 408 | ! at the top of the atmosphere set to 0 |
---|
| 409 | Ke(:,nlay+1)=0. |
---|
| 410 | ! at the surface, set it equal to that at the first model level |
---|
| 411 | Ke(:,1)=Ke(:,2) |
---|
| 412 | |
---|
| 413 | |
---|
| 414 | ! calculate intermediary variables |
---|
| 415 | |
---|
| 416 | DO ilay=2,nlay |
---|
| 417 | DO igrid=1,ngrid |
---|
| 418 | Kem=0.5*(Ke(igrid,ilay+1)+Ke(igrid,ilay)) |
---|
| 419 | KKt=Kem*play(igrid,ilay)/rd/temp(igrid,ilay)/(z_interf(igrid,ilay+1)-z_interf(igrid,ilay)) |
---|
| 420 | Kem=0.5*(Ke(igrid,ilay)+Ke(igrid,ilay-1)) |
---|
| 421 | KKb=Kem*play(igrid,ilay-1)/rd/temp(igrid,ilay-1)/(z_interf(igrid,ilay)-z_interf(igrid,ilay-1)) |
---|
| 422 | gammak=1./(z_lay(igrid,ilay)-z_lay(igrid,ilay-1)) |
---|
| 423 | ak(igrid,ilay)=-gammak*dtime*KKb |
---|
| 424 | ck(igrid,ilay)=-gammak*dtime*KKt |
---|
| 425 | bk(igrid,ilay)=1.+gammak*dtime*(KKt+KKb) |
---|
| 426 | ENDDO |
---|
| 427 | ENDDO |
---|
| 428 | |
---|
| 429 | ! calculate CCK and DDK coefficients |
---|
| 430 | ! downhill phase |
---|
| 431 | |
---|
| 432 | DO igrid=1,ngrid |
---|
| 433 | CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay) |
---|
| 434 | DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay) |
---|
| 435 | ENDDO |
---|
| 436 | |
---|
| 437 | |
---|
| 438 | DO ilay=nlay-1,2,-1 |
---|
| 439 | DO igrid=1,ngrid |
---|
| 440 | CCK(igrid,ilay)=(tke(igrid,ilay)/bk(igrid,ilay)-ck(igrid,ilay)/bk(igrid,ilay)*CCK(igrid,ilay+1)) & |
---|
| 441 | / (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1)) |
---|
| 442 | DDK(igrid,ilay)=-ak(igrid,ilay)/bk(igrid,ilay)/(1+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1)) |
---|
| 443 | ENDDO |
---|
| 444 | ENDDO |
---|
| 445 | |
---|
| 446 | ! calculate TKE |
---|
| 447 | ! uphill phase |
---|
| 448 | |
---|
| 449 | DO ilay=2,nlay+1 |
---|
| 450 | DO igrid=1,ngrid |
---|
| 451 | dtke(igrid,ilay)=CCK(igrid,ilay)+DDK(igrid,ilay)*tke(igrid,ilay-1)-tke(igrid,ilay) |
---|
| 452 | ENDDO |
---|
| 453 | ENDDO |
---|
| 454 | |
---|
| 455 | ! update TKE |
---|
| 456 | tke(:,:)=tke(:,:)+dtke(:,:) |
---|
| 457 | |
---|
| 458 | |
---|
| 459 | end subroutine atke_vdiff_tke |
---|
| 460 | |
---|
| 461 | |
---|
| 462 | |
---|
[4687] | 463 | end module lmdz_atke_exchange_coeff |
---|