[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, & |
---|
[4449] | 8 | wind_u,wind_v,temp,play,pinterf, & |
---|
[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 | |
---|
[4481] | 28 | USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd |
---|
[4478] | 29 | USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, rg, rd |
---|
[4631] | 30 | USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, iflag_atke_lmix, lmin |
---|
[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) |
---|
| 45 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: play ! pressure (Pa) |
---|
| 46 | REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: pinterf ! pressure at interfaces(Pa) |
---|
| 47 | |
---|
| 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) |
---|
[4631] | 69 | LOGICAL, DIMENSION(ngrid,nlay+1) :: switch_num ! switch of numerical integration in very stable cases |
---|
[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 |
---|
[4449] | 80 | |
---|
| 81 | ! Initializations: |
---|
| 82 | !================ |
---|
| 83 | |
---|
| 84 | DO igrid=1,ngrid |
---|
[4478] | 85 | dz_interf(igrid,1) = 0.0 |
---|
[4449] | 86 | z_interf(igrid,1) = 0.0 |
---|
| 87 | END DO |
---|
| 88 | |
---|
[4478] | 89 | ! Calculation of potential temperature: (if vapor -> todo virtual potential temperature) |
---|
| 90 | !===================================== |
---|
[4449] | 91 | |
---|
[4478] | 92 | preff=100000. |
---|
| 93 | ! The result should not depend on the choice of preff |
---|
| 94 | DO ilay=1,nlay |
---|
| 95 | DO igrid = 1, ngrid |
---|
| 96 | theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd) |
---|
| 97 | END DO |
---|
| 98 | END DO |
---|
[4449] | 99 | |
---|
| 100 | |
---|
[4478] | 101 | |
---|
| 102 | ! Calculation of altitude of layers' middle and bottom interfaces: |
---|
| 103 | !================================================================= |
---|
| 104 | |
---|
[4449] | 105 | DO ilay=2,nlay+1 |
---|
| 106 | DO igrid=1,ngrid |
---|
[4478] | 107 | dz_interf(igrid,ilay-1) = rd*temp(igrid,ilay-1)/rg/play(igrid,ilay-1)*(pinterf(igrid,ilay-1)-pinterf(igrid,ilay)) |
---|
| 108 | z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1) |
---|
[4449] | 109 | ENDDO |
---|
| 110 | ENDDO |
---|
| 111 | |
---|
[4478] | 112 | DO ilay=1,nlay |
---|
| 113 | DO igrid=1,ngrid |
---|
| 114 | z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay)) |
---|
| 115 | ENDDO |
---|
| 116 | ENDDO |
---|
[4449] | 117 | |
---|
[4478] | 118 | |
---|
[4449] | 119 | ! Computes the gradient Richardson's number and stability functions: |
---|
| 120 | !=================================================================== |
---|
| 121 | |
---|
[4481] | 122 | ! calculation of cn = Sm value at Ri=0 |
---|
| 123 | ! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0 |
---|
| 124 | cn=(1./sqrt(cepsilon))**(2/3) |
---|
| 125 | ! calculation of Ri0 such that continuity in slope of Sm at Ri=0 |
---|
| 126 | Ri0=2./rpi*(cinf - cn)*ric/cn |
---|
| 127 | ! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0 |
---|
| 128 | Ri1 = -2./rpi * (pr_asym - pr_neut) / pr_slope |
---|
| 129 | |
---|
| 130 | |
---|
[4478] | 131 | DO ilay=2,nlay |
---|
[4449] | 132 | DO igrid=1,ngrid |
---|
[4478] | 133 | dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1) |
---|
| 134 | thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1)) |
---|
[4631] | 135 | N2(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay) |
---|
| 136 | shear2(igrid,ilay)= (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & |
---|
| 137 | ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 ) |
---|
| 138 | Ri(igrid,ilay) = N2(igrid,ilay) / MAX(shear2(igrid,ilay),1E-10) |
---|
[4481] | 139 | |
---|
| 140 | IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases |
---|
[4478] | 141 | Sm(igrid,ilay) = 2./rpi * (cinf-cn) * atan(-Ri(igrid,ilay)/Ri0) + cn |
---|
| 142 | Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut |
---|
[4481] | 143 | ELSE ! stable cases |
---|
[4449] | 144 | Sm(igrid,ilay) = max(0.,cn*(1.-Ri(igrid,ilay)/Ric)) |
---|
[4478] | 145 | Prandtl(igrid,ilay) = pr_neut + Ri(igrid,ilay) * pr_slope |
---|
[4481] | 146 | IF (Ri(igrid,ilay) .GE. Prandtl(igrid,ilay)) THEN |
---|
| 147 | call abort_physic("atke_compute_km_kh", & |
---|
| 148 | 'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1) |
---|
| 149 | ENDIF |
---|
[4449] | 150 | END IF |
---|
| 151 | |
---|
| 152 | Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay) |
---|
| 153 | |
---|
| 154 | ENDDO |
---|
| 155 | ENDDO |
---|
| 156 | |
---|
[4631] | 157 | |
---|
| 158 | ! Computing the mixing length: |
---|
| 159 | !============================================================== |
---|
| 160 | |
---|
| 161 | switch_num(:,:)=.false. |
---|
| 162 | |
---|
| 163 | IF (iflag_atke_lmix .EQ. 1 ) THEN |
---|
| 164 | |
---|
| 165 | DO ilay=2,nlay |
---|
| 166 | DO igrid=1,ngrid |
---|
| 167 | l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) |
---|
| 168 | IF (N2(igrid,ilay) .GT. 0.) THEN |
---|
| 169 | lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)) |
---|
| 170 | IF (lstrat .LT. l_exchange(igrid,ilay)) THEN |
---|
| 171 | l_exchange(igrid,ilay)=max(lstrat,lmin) |
---|
| 172 | switch_num(igrid,ilay)=.true. |
---|
| 173 | ENDIF |
---|
| 174 | ENDIF |
---|
| 175 | ENDDO |
---|
| 176 | ENDDO |
---|
| 177 | |
---|
[4632] | 178 | ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN |
---|
| 179 | ! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms, eq 2 |
---|
| 180 | DO ilay=2,nlay |
---|
| 181 | DO igrid=1,ngrid |
---|
| 182 | l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) |
---|
| 183 | IF (N2(igrid,ilay) .GT. 0.) THEN |
---|
| 184 | lstrat=clmix*sqrt(tke(igrid,ilay))/(2.*max(sqrt(shear2(igrid,ilay)),1E-10)*(1.+sqrt(Ri(igrid,ilay))/2.)) |
---|
| 185 | IF (lstrat .LT. l_exchange(igrid,ilay)) THEN |
---|
| 186 | l_exchange(igrid,ilay)=max(lstrat,lmin) |
---|
| 187 | ENDIF |
---|
| 188 | ENDIF |
---|
| 189 | ENDDO |
---|
| 190 | ENDDO |
---|
| 191 | |
---|
| 192 | |
---|
| 193 | |
---|
[4631] | 194 | ELSE |
---|
[4632] | 195 | ! default: neglect effect of local stratification and shear |
---|
[4631] | 196 | |
---|
| 197 | DO ilay=2,nlay+1 |
---|
| 198 | DO igrid=1,ngrid |
---|
| 199 | l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) |
---|
| 200 | ENDDO |
---|
| 201 | |
---|
| 202 | ENDDO |
---|
| 203 | ENDIF |
---|
| 204 | |
---|
| 205 | |
---|
[4449] | 206 | ! Computing the TKE: |
---|
| 207 | !=================== |
---|
| 208 | IF (iflag_atke == 0) THEN |
---|
| 209 | |
---|
[4631] | 210 | DO ilay=2,nlay |
---|
[4449] | 211 | DO igrid=1,ngrid |
---|
| 212 | tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * & |
---|
[4631] | 213 | shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)) |
---|
[4449] | 214 | ENDDO |
---|
| 215 | ENDDO |
---|
| 216 | |
---|
[4631] | 217 | ELSE IF (iflag_atke == 1) THEN |
---|
| 218 | |
---|
| 219 | ! full implicit scheme resolved with a second order polynomial equation |
---|
| 220 | |
---|
| 221 | DO ilay=2,nlay |
---|
| 222 | DO igrid=1,ngrid |
---|
| 223 | qq=sqrt(2.*tke(igrid,ilay)) |
---|
| 224 | delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * & |
---|
| 225 | (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*shear2(igrid,ilay) & |
---|
| 226 | *(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) |
---|
| 227 | qq=(-1. + sqrt(delta))/dtime*cepsilon*sqrt(2.)*l_exchange(igrid,ilay) |
---|
| 228 | tke(igrid,ilay)=0.5*(qq**2) |
---|
| 229 | ENDDO |
---|
| 230 | ENDDO |
---|
| 231 | |
---|
| 232 | ELSE IF (iflag_atke == 2) THEN |
---|
| 233 | |
---|
| 234 | ! full implicit scheme resolved with a second order polynomial equation |
---|
[4632] | 235 | ! + exact resolution for very stable cases (iflag_atke_lmix=1) |
---|
[4631] | 236 | DO ilay=2,nlay |
---|
| 237 | DO igrid=1,ngrid |
---|
| 238 | qq=sqrt(2.*tke(igrid,ilay)) |
---|
| 239 | IF (switch_num(igrid,ilay)) THEN |
---|
| 240 | taustrat=clmix/2./sqrt(N2(igrid,ilay))*Sm(igrid,ilay)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) & |
---|
| 241 | - sqrt(N2(igrid,ilay))/2./cepsilon/clmix |
---|
| 242 | taustrat=-1./taustrat |
---|
| 243 | qq=qq*exp(-dtime/taustrat) |
---|
| 244 | ELSE |
---|
| 245 | delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * & |
---|
| 246 | (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay) * shear2(igrid,ilay) & |
---|
| 247 | *(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) |
---|
| 248 | qq=(-1. + sqrt(delta))/dtime*cepsilon*sqrt(2.)*l_exchange(igrid,ilay) |
---|
| 249 | ENDIF |
---|
| 250 | tke(igrid,ilay)=0.5*(qq**2) |
---|
| 251 | ENDDO |
---|
| 252 | ENDDO |
---|
| 253 | |
---|
| 254 | |
---|
| 255 | |
---|
| 256 | ELSE IF (iflag_atke == 3) THEN |
---|
| 257 | |
---|
| 258 | ! semi implicit scheme when l does not depend on tke |
---|
| 259 | ! positive-guaranteed if pr slope in stable condition >1 |
---|
| 260 | DO ilay=2,nlay |
---|
| 261 | DO igrid=1,ngrid |
---|
| 262 | qq=sqrt(2.*tke(igrid,ilay)) |
---|
| 263 | qq=(qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) & |
---|
| 264 | /(1.+dtime*qq/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2))) |
---|
| 265 | tke(igrid,ilay)=0.5*(qq**2) |
---|
| 266 | ENDDO |
---|
| 267 | ENDDO |
---|
| 268 | |
---|
| 269 | |
---|
| 270 | ELSE |
---|
[4463] | 271 | call abort_physic("atke_compute_km_kh", & |
---|
[4631] | 272 | 'numerical treatment of TKE not possible yet', 1) |
---|
[4449] | 273 | |
---|
| 274 | END IF |
---|
| 275 | |
---|
| 276 | |
---|
| 277 | ! Computing eddy diffusivity coefficients: |
---|
| 278 | !======================================== |
---|
[4478] | 279 | DO ilay=2,nlay ! TODO: also calculate for nlay+1 ? |
---|
[4449] | 280 | DO igrid=1,ngrid |
---|
[4478] | 281 | ! we add the molecular viscosity to Km,h |
---|
| 282 | Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5 |
---|
| 283 | Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5 |
---|
[4449] | 284 | END DO |
---|
| 285 | END DO |
---|
| 286 | |
---|
[4478] | 287 | ! for output: |
---|
| 288 | !=========== |
---|
| 289 | Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay) |
---|
| 290 | Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay) |
---|
[4449] | 291 | |
---|
| 292 | end subroutine atke_compute_km_kh |
---|
| 293 | |
---|
| 294 | |
---|
| 295 | end module atke_exchange_coeff_mod |
---|