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