| 1 | ! $Header$ |
|---|
| 2 | |
|---|
| 3 | SUBROUTINE yamada_c(ngrid, timestep, plev, play & |
|---|
| 4 | , pu, pv, pt, d_u, d_v, d_t, cd, q2, km, kn, kq, d_t_diss, ustar & |
|---|
| 5 | , iflag_pbl) |
|---|
| 6 | USE dimphy, ONLY: klon, klev |
|---|
| 7 | USE lmdz_print_control, ONLY: prt_level |
|---|
| 8 | USE lmdz_ioipsl_getin_p, ONLY: getin_p |
|---|
| 9 | USE lmdz_yomcst |
|---|
| 10 | |
|---|
| 11 | IMPLICIT NONE |
|---|
| 12 | |
|---|
| 13 | ! timestep : pas de temps |
|---|
| 14 | ! g : g |
|---|
| 15 | ! zlev : altitude a chaque niveau (interface inferieure de la couche |
|---|
| 16 | ! de meme indice) |
|---|
| 17 | ! zlay : altitude au centre de chaque couche |
|---|
| 18 | ! u,v : vitesse au centre de chaque couche |
|---|
| 19 | ! (en entree : la valeur au debut du pas de temps) |
|---|
| 20 | ! teta : temperature potentielle au centre de chaque couche |
|---|
| 21 | ! (en entree : la valeur au debut du pas de temps) |
|---|
| 22 | ! cd : cdrag |
|---|
| 23 | ! (en entree : la valeur au debut du pas de temps) |
|---|
| 24 | ! q2 : $q^2$ au bas de chaque couche |
|---|
| 25 | ! (en entree : la valeur au debut du pas de temps) |
|---|
| 26 | ! (en sortie : la valeur a la fin du pas de temps) |
|---|
| 27 | ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque |
|---|
| 28 | ! couche) |
|---|
| 29 | ! (en sortie : la valeur a la fin du pas de temps) |
|---|
| 30 | ! kn : diffusivite turbulente des scalaires (au bas de chaque couche) |
|---|
| 31 | ! (en sortie : la valeur a la fin du pas de temps) |
|---|
| 32 | |
|---|
| 33 | ! iflag_pbl doit valoir entre 6 et 9 |
|---|
| 34 | ! l=6, on prend systematiquement une longueur d'equilibre |
|---|
| 35 | ! iflag_pbl=6 : MY 2.0 |
|---|
| 36 | ! iflag_pbl=7 : MY 2.0.Fournier |
|---|
| 37 | ! iflag_pbl=8/9 : MY 2.5 |
|---|
| 38 | ! iflag_pbl=8 with special obsolete treatments for convergence |
|---|
| 39 | ! with Cmpi5 NPv3.1 simulations |
|---|
| 40 | ! iflag_pbl=10/11 : New scheme M2 and N2 explicit and dissiptation exact |
|---|
| 41 | ! iflag_pbl=12 = 11 with vertical diffusion off q2 |
|---|
| 42 | |
|---|
| 43 | ! 2013/04/01 (FH hourdin@lmd.jussieu.fr) |
|---|
| 44 | ! Correction for very stable PBLs (iflag_pbl=10 and 11) |
|---|
| 45 | ! iflag_pbl=8 converges numerically with NPv3.1 |
|---|
| 46 | ! iflag_pbl=11 -> the model starts with NP from start files created by ce0l |
|---|
| 47 | ! -> the model can run with longer time-steps. |
|---|
| 48 | !....................................................................... |
|---|
| 49 | |
|---|
| 50 | REAL, DIMENSION(klon, klev) :: d_u, d_v, d_t |
|---|
| 51 | REAL, DIMENSION(klon, klev) :: pu, pv, pt |
|---|
| 52 | REAL, DIMENSION(klon, klev) :: d_t_diss |
|---|
| 53 | |
|---|
| 54 | REAL timestep |
|---|
| 55 | REAL plev(klon, klev + 1) |
|---|
| 56 | REAL play(klon, klev) |
|---|
| 57 | REAL ustar(klon) |
|---|
| 58 | REAL kmin, qmin, pblhmin(klon), coriol(klon) |
|---|
| 59 | REAL zlev(klon, klev + 1) |
|---|
| 60 | REAL zlay(klon, klev) |
|---|
| 61 | REAL zu(klon, klev) |
|---|
| 62 | REAL zv(klon, klev) |
|---|
| 63 | REAL zt(klon, klev) |
|---|
| 64 | REAL teta(klon, klev) |
|---|
| 65 | REAL cd(klon) |
|---|
| 66 | REAL q2(klon, klev + 1), qpre |
|---|
| 67 | REAL unsdz(klon, klev) |
|---|
| 68 | REAL unsdzdec(klon, klev + 1) |
|---|
| 69 | |
|---|
| 70 | REAL km(klon, klev) |
|---|
| 71 | REAL kmpre(klon, klev + 1), tmp2 |
|---|
| 72 | REAL mpre(klon, klev + 1) |
|---|
| 73 | REAL kn(klon, klev) |
|---|
| 74 | REAL kq(klon, klev) |
|---|
| 75 | REAL ff(klon, klev + 1), delta(klon, klev + 1) |
|---|
| 76 | REAL aa(klon, klev + 1), aa0, aa1 |
|---|
| 77 | INTEGER iflag_pbl, ngrid |
|---|
| 78 | INTEGER nlay, nlev |
|---|
| 79 | |
|---|
| 80 | LOGICAL first |
|---|
| 81 | INTEGER ipas |
|---|
| 82 | save first, ipas |
|---|
| 83 | !FH/IM data first,ipas/.TRUE.,0/ |
|---|
| 84 | data first, ipas/.FALSE., 0/ |
|---|
| 85 | !$OMP THREADPRIVATE( first,ipas) |
|---|
| 86 | INTEGER, SAVE :: iflag_tke_diff = 0 |
|---|
| 87 | !$OMP THREADPRIVATE(iflag_tke_diff) |
|---|
| 88 | |
|---|
| 89 | INTEGER ig, k |
|---|
| 90 | |
|---|
| 91 | REAL ri, zrif, zalpha, zsm, zsn |
|---|
| 92 | REAL rif(klon, klev + 1), sm(klon, klev + 1), alpha(klon, klev) |
|---|
| 93 | |
|---|
| 94 | REAL m2(klon, klev + 1), dz(klon, klev + 1), zq, n2(klon, klev + 1) |
|---|
| 95 | REAL, DIMENSION(klon, klev + 1) :: km2, kn2, sqrtq |
|---|
| 96 | REAL dtetadz(klon, klev + 1) |
|---|
| 97 | REAL m2cstat, mcstat, kmcstat |
|---|
| 98 | REAL l(klon, klev + 1) |
|---|
| 99 | REAL leff(klon, klev + 1) |
|---|
| 100 | REAL, ALLOCATABLE, save :: l0(:) |
|---|
| 101 | !$OMP THREADPRIVATE(l0) |
|---|
| 102 | REAL sq(klon), sqz(klon), zz(klon, klev + 1) |
|---|
| 103 | INTEGER iter |
|---|
| 104 | |
|---|
| 105 | REAL ric, rifc, b1, kap |
|---|
| 106 | save ric, rifc, b1, kap |
|---|
| 107 | data ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.4/ |
|---|
| 108 | !$OMP THREADPRIVATE(ric,rifc,b1,kap) |
|---|
| 109 | REAL frif, falpha, fsm |
|---|
| 110 | REAL fl, zzz, zl0, zq2, zn2 |
|---|
| 111 | |
|---|
| 112 | REAL rino(klon, klev + 1), smyam(klon, klev), styam(klon, klev) |
|---|
| 113 | REAL lyam(klon, klev), knyam(klon, klev) |
|---|
| 114 | REAL w2yam(klon, klev), t2yam(klon, klev) |
|---|
| 115 | logical, save :: firstcall = .TRUE. |
|---|
| 116 | !$OMP THREADPRIVATE(firstcall) |
|---|
| 117 | CHARACTER(len = 20), PARAMETER :: modname = "yamada_c" |
|---|
| 118 | REAL, DIMENSION(klon, klev + 1) :: fluxu, fluxv, fluxt |
|---|
| 119 | REAL, DIMENSION(klon, klev + 1) :: dddu, dddv, dddt |
|---|
| 120 | REAL, DIMENSION(klon, klev) :: exner, masse |
|---|
| 121 | REAL, DIMENSION(klon, klev + 1) :: masseb, q2old, q2neg |
|---|
| 122 | LOGICAL okiophys |
|---|
| 123 | |
|---|
| 124 | frif(ri) = 0.6588 * (ri + 0.1776 - sqrt(ri * ri - 0.3221 * ri + 0.03156)) |
|---|
| 125 | falpha(ri) = 1.318 * (0.2231 - ri) / (0.2341 - ri) |
|---|
| 126 | fsm(ri) = 1.96 * (0.1912 - ri) * (0.2341 - ri) / ((1. - ri) * (0.2231 - ri)) |
|---|
| 127 | fl(zzz, zl0, zq2, zn2) = & |
|---|
| 128 | max(min(l0(ig) * kap * zlev(ig, k) / (kap * zlev(ig, k) + l0(ig)) & |
|---|
| 129 | , 0.5 * sqrt(q2(ig, k)) / sqrt(max(n2(ig, k), 1.e-10))), 1.) |
|---|
| 130 | |
|---|
| 131 | okiophys = klon==1 |
|---|
| 132 | IF (firstcall) THEN |
|---|
| 133 | CALL getin_p('iflag_tke_diff', iflag_tke_diff) |
|---|
| 134 | allocate(l0(klon)) |
|---|
| 135 | #define IOPHYS |
|---|
| 136 | #ifdef IOPHYS |
|---|
| 137 | ! CALL iophys_ini(timestep) |
|---|
| 138 | #endif |
|---|
| 139 | firstcall = .FALSE. |
|---|
| 140 | endif |
|---|
| 141 | |
|---|
| 142 | IF (ngrid<=0) RETURN ! Bizarre : on n a pas ce probeleme pour coef_diff_turb |
|---|
| 143 | |
|---|
| 144 | #ifdef IOPHYS |
|---|
| 145 | IF (okiophys) THEN |
|---|
| 146 | CALL iophys_ecrit('q2i', klev, 'q2 debut my', 'm2/s2', q2(:, 1:klev)) |
|---|
| 147 | CALL iophys_ecrit('kmi', klev, 'Kz debut my', 'm/s2', km(:, 1:klev)) |
|---|
| 148 | END IF |
|---|
| 149 | #endif |
|---|
| 150 | |
|---|
| 151 | nlay = klev |
|---|
| 152 | nlev = klev + 1 |
|---|
| 153 | |
|---|
| 154 | |
|---|
| 155 | !------------------------------------------------------------------------- |
|---|
| 156 | ! Computation of conservative source terms from the turbulent tendencies |
|---|
| 157 | !------------------------------------------------------------------------- |
|---|
| 158 | |
|---|
| 159 | zalpha = 0.5 ! Anciennement 0.5. Essayer de voir pourquoi ? |
|---|
| 160 | zu(:, :) = pu(:, :) + zalpha * d_u(:, :) |
|---|
| 161 | zv(:, :) = pv(:, :) + zalpha * d_v(:, :) |
|---|
| 162 | zt(:, :) = pt(:, :) + zalpha * d_t(:, :) |
|---|
| 163 | |
|---|
| 164 | DO k = 1, klev |
|---|
| 165 | exner(:, k) = (play(:, k) / plev(:, 1))**RKAPPA |
|---|
| 166 | masse(:, k) = (plev(:, k) - plev(:, k + 1)) / RG |
|---|
| 167 | teta(:, k) = zt(:, k) / exner(:, k) |
|---|
| 168 | enddo |
|---|
| 169 | |
|---|
| 170 | ! Atmospheric mass at layer interfaces, where the TKE is computed |
|---|
| 171 | masseb(:, :) = 0. |
|---|
| 172 | DO k = 1, klev |
|---|
| 173 | masseb(:, k) = masseb(:, k) + masse(:, k) |
|---|
| 174 | masseb(:, k + 1) = masseb(:, k + 1) + masse(:, k) |
|---|
| 175 | enddo |
|---|
| 176 | masseb(:, :) = 0.5 * masseb(:, :) |
|---|
| 177 | |
|---|
| 178 | zlev(:, 1) = 0. |
|---|
| 179 | zlay(:, 1) = RCPD * teta(:, 1) * (1. - exner(:, 1)) |
|---|
| 180 | DO k = 1, klev - 1 |
|---|
| 181 | zlay(:, k + 1) = zlay(:, k) + 0.5 * RCPD * (teta(:, k) + teta(:, k + 1)) * (exner(:, k) - exner(:, k + 1)) / RG |
|---|
| 182 | zlev(:, k) = 0.5 * (zlay(:, k) + zlay(:, k + 1)) ! PASBO |
|---|
| 183 | enddo |
|---|
| 184 | |
|---|
| 185 | fluxu(:, klev + 1) = 0. |
|---|
| 186 | fluxv(:, klev + 1) = 0. |
|---|
| 187 | fluxt(:, klev + 1) = 0. |
|---|
| 188 | |
|---|
| 189 | DO k = klev, 1, -1 |
|---|
| 190 | fluxu(:, k) = fluxu(:, k + 1) + masse(:, k) * d_u(:, k) |
|---|
| 191 | fluxv(:, k) = fluxv(:, k + 1) + masse(:, k) * d_v(:, k) |
|---|
| 192 | fluxt(:, k) = fluxt(:, k + 1) + masse(:, k) * d_t(:, k) / exner(:, k) ! Flux de theta |
|---|
| 193 | enddo |
|---|
| 194 | |
|---|
| 195 | dddu(:, 1) = 2 * zu(:, 1) * fluxu(:, 1) |
|---|
| 196 | dddv(:, 1) = 2 * zv(:, 1) * fluxv(:, 1) |
|---|
| 197 | dddt(:, 1) = (exner(:, 1) - 1.) * fluxt(:, 1) |
|---|
| 198 | |
|---|
| 199 | DO k = 2, klev |
|---|
| 200 | dddu(:, k) = (zu(:, k) - zu(:, k - 1)) * fluxu(:, k) |
|---|
| 201 | dddv(:, k) = (zv(:, k) - zv(:, k - 1)) * fluxv(:, k) |
|---|
| 202 | dddt(:, k) = (exner(:, k) - exner(:, k - 1)) * fluxt(:, k) |
|---|
| 203 | enddo |
|---|
| 204 | dddu(:, klev + 1) = 0. |
|---|
| 205 | dddv(:, klev + 1) = 0. |
|---|
| 206 | dddt(:, klev + 1) = 0. |
|---|
| 207 | |
|---|
| 208 | #ifdef IOPHYS |
|---|
| 209 | IF (okiophys) THEN |
|---|
| 210 | CALL iophys_ecrit('zlay', klev, 'Geop', 'm', zlay) |
|---|
| 211 | CALL iophys_ecrit('teta', klev, 'teta', 'K', teta) |
|---|
| 212 | CALL iophys_ecrit('temp', klev, 'temp', 'K', zt) |
|---|
| 213 | CALL iophys_ecrit('pt', klev, 'temp', 'K', pt) |
|---|
| 214 | CALL iophys_ecrit('pu', klev, 'u', 'm/s', pu) |
|---|
| 215 | CALL iophys_ecrit('pv', klev, 'v', 'm/s', pv) |
|---|
| 216 | CALL iophys_ecrit('d_u', klev, 'd_u', 'm/s2', d_u) |
|---|
| 217 | CALL iophys_ecrit('d_v', klev, 'd_v', 'm/s2', d_v) |
|---|
| 218 | CALL iophys_ecrit('d_t', klev, 'd_t', 'K/s', d_t) |
|---|
| 219 | CALL iophys_ecrit('exner', klev, 'exner', '', exner) |
|---|
| 220 | CALL iophys_ecrit('masse', klev, 'masse', '', masse) |
|---|
| 221 | CALL iophys_ecrit('masseb', klev, 'masseb', '', masseb) |
|---|
| 222 | END IF |
|---|
| 223 | #endif |
|---|
| 224 | |
|---|
| 225 | |
|---|
| 226 | |
|---|
| 227 | ipas = ipas + 1 |
|---|
| 228 | |
|---|
| 229 | |
|---|
| 230 | !....................................................................... |
|---|
| 231 | ! les increments verticaux |
|---|
| 232 | !....................................................................... |
|---|
| 233 | |
|---|
| 234 | !!!!!! allerte !!!!!c |
|---|
| 235 | !!!!!! zlev n'est pas declare a nlev !!!!!c |
|---|
| 236 | !!!!!! ----> |
|---|
| 237 | DO ig = 1, ngrid |
|---|
| 238 | zlev(ig, nlev) = zlay(ig, nlay) & |
|---|
| 239 | + (zlay(ig, nlay) - zlev(ig, nlev - 1)) |
|---|
| 240 | ENDDO |
|---|
| 241 | !!!!!! <---- |
|---|
| 242 | !!!!!! allerte !!!!!c |
|---|
| 243 | |
|---|
| 244 | DO k = 1, nlay |
|---|
| 245 | DO ig = 1, ngrid |
|---|
| 246 | unsdz(ig, k) = 1.E+0 / (zlev(ig, k + 1) - zlev(ig, k)) |
|---|
| 247 | ENDDO |
|---|
| 248 | ENDDO |
|---|
| 249 | DO ig = 1, ngrid |
|---|
| 250 | unsdzdec(ig, 1) = 1.E+0 / (zlay(ig, 1) - zlev(ig, 1)) |
|---|
| 251 | ENDDO |
|---|
| 252 | DO k = 2, nlay |
|---|
| 253 | DO ig = 1, ngrid |
|---|
| 254 | unsdzdec(ig, k) = 1.E+0 / (zlay(ig, k) - zlay(ig, k - 1)) |
|---|
| 255 | ENDDO |
|---|
| 256 | ENDDO |
|---|
| 257 | DO ig = 1, ngrid |
|---|
| 258 | unsdzdec(ig, nlay + 1) = 1.E+0 / (zlev(ig, nlay + 1) - zlay(ig, nlay)) |
|---|
| 259 | ENDDO |
|---|
| 260 | |
|---|
| 261 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 262 | ! Computing M^2, N^2, Richardson numbers, stability functions |
|---|
| 263 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 264 | |
|---|
| 265 | DO k = 2, klev |
|---|
| 266 | DO ig = 1, ngrid |
|---|
| 267 | dz(ig, k) = zlay(ig, k) - zlay(ig, k - 1) |
|---|
| 268 | m2(ig, k) = ((zu(ig, k) - zu(ig, k - 1))**2 + (zv(ig, k) - zv(ig, k - 1))**2) / (dz(ig, k) * dz(ig, k)) |
|---|
| 269 | dtetadz(ig, k) = (teta(ig, k) - teta(ig, k - 1)) / dz(ig, k) |
|---|
| 270 | n2(ig, k) = RG * 2. * dtetadz(ig, k) / (teta(ig, k - 1) + teta(ig, k)) |
|---|
| 271 | ! n2(ig,k)=0. |
|---|
| 272 | ri = n2(ig, k) / max(m2(ig, k), 1.e-10) |
|---|
| 273 | IF (ri<ric) THEN |
|---|
| 274 | rif(ig, k) = frif(ri) |
|---|
| 275 | else |
|---|
| 276 | rif(ig, k) = rifc |
|---|
| 277 | endif |
|---|
| 278 | IF(rif(ig, k)<0.16) THEN |
|---|
| 279 | alpha(ig, k) = falpha(rif(ig, k)) |
|---|
| 280 | sm(ig, k) = fsm(rif(ig, k)) |
|---|
| 281 | else |
|---|
| 282 | alpha(ig, k) = 1.12 |
|---|
| 283 | sm(ig, k) = 0.085 |
|---|
| 284 | endif |
|---|
| 285 | zz(ig, k) = b1 * m2(ig, k) * (1. - rif(ig, k)) * sm(ig, k) |
|---|
| 286 | enddo |
|---|
| 287 | enddo |
|---|
| 288 | |
|---|
| 289 | |
|---|
| 290 | |
|---|
| 291 | !==================================================================== |
|---|
| 292 | ! Computing the mixing length |
|---|
| 293 | !==================================================================== |
|---|
| 294 | |
|---|
| 295 | ! Mise a jour de l0 |
|---|
| 296 | IF (iflag_pbl==8.OR.iflag_pbl==10) THEN |
|---|
| 297 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 298 | ! Iterative computation of l0 |
|---|
| 299 | ! This version is kept for iflag_pbl only for convergence |
|---|
| 300 | ! with NPv3.1 Cmip5 simulations |
|---|
| 301 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 302 | |
|---|
| 303 | DO ig = 1, ngrid |
|---|
| 304 | sq(ig) = 1.e-10 |
|---|
| 305 | sqz(ig) = 1.e-10 |
|---|
| 306 | enddo |
|---|
| 307 | DO k = 2, klev - 1 |
|---|
| 308 | DO ig = 1, ngrid |
|---|
| 309 | zq = sqrt(q2(ig, k)) |
|---|
| 310 | sqz(ig) = sqz(ig) + zq * zlev(ig, k) * (zlay(ig, k) - zlay(ig, k - 1)) |
|---|
| 311 | sq(ig) = sq(ig) + zq * (zlay(ig, k) - zlay(ig, k - 1)) |
|---|
| 312 | enddo |
|---|
| 313 | enddo |
|---|
| 314 | DO ig = 1, ngrid |
|---|
| 315 | l0(ig) = 0.2 * sqz(ig) / sq(ig) |
|---|
| 316 | enddo |
|---|
| 317 | DO k = 2, klev |
|---|
| 318 | DO ig = 1, ngrid |
|---|
| 319 | l(ig, k) = fl(zlev(ig, k), l0(ig), q2(ig, k), n2(ig, k)) |
|---|
| 320 | enddo |
|---|
| 321 | enddo |
|---|
| 322 | ! PRINT*,'L0 cas 8 ou 10 ',l0 |
|---|
| 323 | |
|---|
| 324 | else |
|---|
| 325 | |
|---|
| 326 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 327 | ! In all other case, the assymptotic mixing length l0 is imposed (100m) |
|---|
| 328 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 329 | |
|---|
| 330 | l0(:) = 150. |
|---|
| 331 | DO k = 2, klev |
|---|
| 332 | DO ig = 1, ngrid |
|---|
| 333 | l(ig, k) = fl(zlev(ig, k), l0(ig), q2(ig, k), n2(ig, k)) |
|---|
| 334 | enddo |
|---|
| 335 | enddo |
|---|
| 336 | ! PRINT*,'L0 cas autres ',l0 |
|---|
| 337 | |
|---|
| 338 | endif |
|---|
| 339 | |
|---|
| 340 | |
|---|
| 341 | #ifdef IOPHYS |
|---|
| 342 | IF (okiophys) THEN |
|---|
| 343 | CALL iophys_ecrit('rif', klev, 'Flux Richardson', 'm', rif(:, 1:klev)) |
|---|
| 344 | CALL iophys_ecrit('m2', klev, 'm2 ', 'm/s', m2(:, 1:klev)) |
|---|
| 345 | CALL iophys_ecrit('Km2app', klev, 'm2 conserv', 'm/s', km(:, 1:klev) * m2(:, 1:klev)) |
|---|
| 346 | CALL iophys_ecrit('Km', klev, 'Km', 'm2/s', km(:, 1:klev)) |
|---|
| 347 | END IF |
|---|
| 348 | #endif |
|---|
| 349 | |
|---|
| 350 | |
|---|
| 351 | IF (iflag_pbl<20) THEN |
|---|
| 352 | ! For diagnostics only |
|---|
| 353 | RETURN |
|---|
| 354 | |
|---|
| 355 | ELSE |
|---|
| 356 | |
|---|
| 357 | ! PRINT*,'OK1' |
|---|
| 358 | |
|---|
| 359 | ! Evolution of TKE under source terms K M2 and K N2 |
|---|
| 360 | leff(:, :) = max(l(:, :), 1.) |
|---|
| 361 | |
|---|
| 362 | !################################################################## |
|---|
| 363 | !# IF (iflag_pbl==29) THEN |
|---|
| 364 | !# stop 'Ne pas utiliser iflag_pbl=29' |
|---|
| 365 | !# km2(:,:)=km(:,:)*m2(:,:) |
|---|
| 366 | !# kn2(:,:)=kn2(:,:)*rif(:,:) |
|---|
| 367 | !# ELSEIF (iflag_pbl==25) THEN |
|---|
| 368 | ! VERSION AVEC LA TKE EN MILIEU DE COUCHE |
|---|
| 369 | !# stop 'Ne pas utiliser iflag_pbl=25' |
|---|
| 370 | !# DO k=1,klev |
|---|
| 371 | !# km2(:,k)=-0.5*(dddu(:,k)+dddv(:,k)+dddu(:,k+1)+dddv(:,k+1)) & |
|---|
| 372 | !# & /(masse(:,k)*timestep) |
|---|
| 373 | !# kn2(:,k)=rcpd*0.5*(dddt(:,k)+dddt(:,k+1))/(masse(:,k)*timestep) |
|---|
| 374 | !# leff(:,k)=0.5*(leff(:,k)+leff(:,k+1)) |
|---|
| 375 | !# ENDDO |
|---|
| 376 | !# km2(:,klev+1)=0. ; kn2(:,klev+1)=0. |
|---|
| 377 | !# ELSE |
|---|
| 378 | !################################################################# |
|---|
| 379 | |
|---|
| 380 | km2(:, :) = -(dddu(:, :) + dddv(:, :)) / (masseb(:, :) * timestep) |
|---|
| 381 | kn2(:, :) = rcpd * dddt(:, :) / (masseb(:, :) * timestep) |
|---|
| 382 | ! ENDIF |
|---|
| 383 | q2neg(:, :) = q2(:, :) + timestep * (km2(:, :) - kn2(:, :)) |
|---|
| 384 | q2(:, :) = min(max(q2neg(:, :), 1.e-10), 1.e4) |
|---|
| 385 | |
|---|
| 386 | |
|---|
| 387 | #ifdef IOPHYS |
|---|
| 388 | IF (okiophys) THEN |
|---|
| 389 | CALL iophys_ecrit('km2', klev, 'm2 conserv', 'm/s', km2(:, 1:klev)) |
|---|
| 390 | CALL iophys_ecrit('kn2', klev, 'n2 conserv', 'm/s', kn2(:, 1:klev)) |
|---|
| 391 | END IF |
|---|
| 392 | #endif |
|---|
| 393 | |
|---|
| 394 | ! Dissipation of TKE |
|---|
| 395 | q2old(:, :) = q2(:, :) |
|---|
| 396 | q2(:, :) = 1. / (1. / sqrt(q2(:, :)) + timestep / (2 * leff(:, :) * b1)) |
|---|
| 397 | q2(:, :) = q2(:, :) * q2(:, :) |
|---|
| 398 | ! IF (iflag_pbl<=24) THEN |
|---|
| 399 | DO k = 1, klev |
|---|
| 400 | d_t_diss(:, k) = (masseb(:, k) * (q2neg(:, k) - q2(:, k)) + masseb(:, k + 1) * (q2neg(:, k + 1) - q2(:, k + 1))) / (2. * rcpd * masse(:, k)) |
|---|
| 401 | ENDDO |
|---|
| 402 | |
|---|
| 403 | !################################################################### |
|---|
| 404 | ! ELSE IF (iflag_pbl<=27) THEN |
|---|
| 405 | ! DO k=1,klev |
|---|
| 406 | ! d_t_diss(:,k)=(q2neg(:,k)-q2(:,k))/rcpd |
|---|
| 407 | ! ENDDO |
|---|
| 408 | ! ENDIF |
|---|
| 409 | ! PRINT*,'iflag_pbl ',d_t_diss |
|---|
| 410 | !################################################################### |
|---|
| 411 | |
|---|
| 412 | |
|---|
| 413 | ! Compuation of stability functions |
|---|
| 414 | ! IF (iflag_pbl/=29) THEN |
|---|
| 415 | DO k = 1, klev |
|---|
| 416 | DO ig = 1, ngrid |
|---|
| 417 | IF (ABS(km2(ig, k))<=1.e-20) THEN |
|---|
| 418 | rif(ig, k) = 0. |
|---|
| 419 | ELSE |
|---|
| 420 | rif(ig, k) = min(kn2(ig, k) / km2(ig, k), rifc) |
|---|
| 421 | ENDIF |
|---|
| 422 | IF (rif(ig, k)<0.16) THEN |
|---|
| 423 | alpha(ig, k) = falpha(rif(ig, k)) |
|---|
| 424 | sm(ig, k) = fsm(rif(ig, k)) |
|---|
| 425 | else |
|---|
| 426 | alpha(ig, k) = 1.12 |
|---|
| 427 | sm(ig, k) = 0.085 |
|---|
| 428 | endif |
|---|
| 429 | ENDDO |
|---|
| 430 | ENDDO |
|---|
| 431 | ! ENDIF |
|---|
| 432 | |
|---|
| 433 | ! Computation of turbulent diffusivities |
|---|
| 434 | ! IF (25<=iflag_pbl.AND.iflag_pbl<=28) THEN |
|---|
| 435 | ! DO k=2,klev |
|---|
| 436 | ! sqrtq(:,k)=sqrt(0.5*(q2(:,k)+q2(:,k-1))) |
|---|
| 437 | ! ENDDO |
|---|
| 438 | ! ELSE |
|---|
| 439 | kq(:, :) = 0. |
|---|
| 440 | DO k = 1, klev |
|---|
| 441 | ! Coefficient au milieu des couches pour diffuser la TKE |
|---|
| 442 | kq(:, k) = 0.5 * leff(:, k) * sqrt(q2(:, k)) * 0.2 |
|---|
| 443 | ENDDO |
|---|
| 444 | |
|---|
| 445 | #ifdef IOPHYS |
|---|
| 446 | IF (okiophys) THEN |
|---|
| 447 | CALL iophys_ecrit('q2b', klev, 'KTE inter', 'm2/s', q2(:, 1:klev)) |
|---|
| 448 | END IF |
|---|
| 449 | #endif |
|---|
| 450 | |
|---|
| 451 | IF (iflag_tke_diff==1) THEN |
|---|
| 452 | CALL vdif_q2(timestep, RG, RD, ngrid, plev, pt, kq, q2) |
|---|
| 453 | ENDIF |
|---|
| 454 | |
|---|
| 455 | km(:, :) = 0. |
|---|
| 456 | kn(:, :) = 0. |
|---|
| 457 | DO k = 1, klev |
|---|
| 458 | km(:, k) = leff(:, k) * sqrt(q2(:, k)) * sm(:, k) |
|---|
| 459 | kn(:, k) = km(:, k) * alpha(:, k) |
|---|
| 460 | ENDDO |
|---|
| 461 | |
|---|
| 462 | |
|---|
| 463 | #ifdef IOPHYS |
|---|
| 464 | IF (okiophys) THEN |
|---|
| 465 | CALL iophys_ecrit('mixingl', klev, 'Mixing length', 'm', leff(:, 1:klev)) |
|---|
| 466 | CALL iophys_ecrit('rife', klev, 'Flux Richardson', 'm', rif(:, 1:klev)) |
|---|
| 467 | CALL iophys_ecrit('q2f', klev, 'KTE finale', 'm2/s', q2(:, 1:klev)) |
|---|
| 468 | CALL iophys_ecrit('q2neg', klev, 'KTE non bornee', 'm2/s', q2neg(:, 1:klev)) |
|---|
| 469 | CALL iophys_ecrit('alpha', klev, 'alpha', '', alpha(:, 1:klev)) |
|---|
| 470 | CALL iophys_ecrit('sm', klev, 'sm', '', sm(:, 1:klev)) |
|---|
| 471 | CALL iophys_ecrit('q2f', klev, 'KTE finale', 'm2/s', q2(:, 1:klev)) |
|---|
| 472 | CALL iophys_ecrit('kmf', klev, 'Kz final', 'm2/s', km(:, 1:klev)) |
|---|
| 473 | CALL iophys_ecrit('knf', klev, 'Kz final', 'm2/s', kn(:, 1:klev)) |
|---|
| 474 | CALL iophys_ecrit('kqf', klev, 'Kz final', 'm2/s', kq(:, 1:klev)) |
|---|
| 475 | END IF |
|---|
| 476 | #endif |
|---|
| 477 | |
|---|
| 478 | ENDIF |
|---|
| 479 | |
|---|
| 480 | |
|---|
| 481 | ! PRINT*,'OK2' |
|---|
| 482 | RETURN |
|---|
| 483 | END |
|---|
