Changeset 2594 for LMDZ5/branches/testing/libf/phylmd/yamada4.F90
- Timestamp:
- Jul 18, 2016, 9:41:10 PM (8 years ago)
- Location:
- LMDZ5/branches/testing
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
libf/phylmd/yamada4.F90 (modified) (30 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/branches/testing
- Property svn:mergeinfo changed
/LMDZ5/trunk merged: 2547-2567,2569,2571-2574,2576-2589
- Property svn:mergeinfo changed
-
LMDZ5/branches/testing/libf/phylmd/yamada4.F90
r2471 r2594 1 2 ! $Header$ 3 4 SUBROUTINE yamada4(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, & 1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 3 SUBROUTINE yamada4(ni, nsrf, ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, & 5 4 cd, q2, km, kn, kq, ustar, iflag_pbl) 5 6 6 USE dimphy 7 USE print_control_mod, ONLY: prt_level8 7 USE ioipsl_getin_p_mod, ONLY : getin_p 9 8 10 9 IMPLICIT NONE 11 10 include "iniprint.h" 11 ! ....................................................................... 12 ! ym#include "dimensions.h" 13 ! ym#include "dimphy.h" 14 ! ************************************************************************************************ 15 ! 16 ! yamada4: subroutine qui calcule le transfert turbulent avec une fermeture d'ordre 2 ou 2.5 17 ! 18 ! Reference: Simulation of nocturnal drainage flows by a q2l Turbulence Closure Model 19 ! Yamada T. 20 ! J. Atmos. Sci, 40, 91-106, 1983 21 ! 22 !************************************************************************************************ 23 ! Input : 24 !'====== 25 ! ni: indice horizontal sur la grille de base, non restreinte 26 ! nsrf: type de surface 27 ! ngrid: nombre de mailles concern??es sur l'horizontal 12 28 ! dt : pas de temps 13 29 ! g : g 30 ! rconst: constante de l'air sec 14 31 ! zlev : altitude a chaque niveau (interface inferieure de la couche 15 32 ! de meme indice) … … 17 34 ! u,v : vitesse au centre de chaque couche 18 35 ! (en entree : la valeur au debut du pas de temps) 19 ! teta : temperature potentielle au centre de chaque couche36 ! teta : temperature potentielle virtuelle au centre de chaque couche 20 37 ! (en entree : la valeur au debut du pas de temps) 21 ! cd : cdrag 38 ! cd : cdrag pour la quantit?? de mouvement 22 39 ! (en entree : la valeur au debut du pas de temps) 40 ! ustar: vitesse de friction calcul??e par une formule diagnostique 41 ! iflag_pbl: flag pour choisir des options du sch??ma de turbulence 42 ! 43 ! iflag_pbl doit valoir entre 6 et 9 44 ! l=6, on prend systematiquement une longueur d'equilibre 45 ! iflag_pbl=6 : MY 2.0 46 ! iflag_pbl=7 : MY 2.0.Fournier 47 ! iflag_pbl=8/9 : MY 2.5 48 ! iflag_pbl=8 with special obsolete treatments for convergence 49 ! with Cmpi5 NPv3.1 simulations 50 ! iflag_pbl=10/11 : New scheme M2 and N2 explicit and dissiptation exact 51 ! iflag_pbl=12 = 11 with vertical diffusion off q2 52 ! 53 ! 2013/04/01 (FH hourdin@lmd.jussieu.fr) 54 ! Correction for very stable PBLs (iflag_pbl=10 and 11) 55 ! iflag_pbl=8 converges numerically with NPv3.1 56 ! iflag_pbl=11 -> the model starts with NP from start files created by ce0l 57 ! -> the model can run with longer time-steps. 58 ! 59 ! Inpout/Output : 60 !============== 23 61 ! q2 : $q^2$ au bas de chaque couche 24 62 ! (en entree : la valeur au debut du pas de temps) 25 63 ! (en sortie : la valeur a la fin du pas de temps) 64 65 ! Outputs: 66 !========== 26 67 ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque 27 68 ! couche) … … 29 70 ! kn : diffusivite turbulente des scalaires (au bas de chaque couche) 30 71 ! (en sortie : la valeur a la fin du pas de temps) 31 32 ! iflag_pbl doit valoir entre 6 et 9 33 ! l=6, on prend systematiquement une longueur d'equilibre 34 ! iflag_pbl=6 : MY 2.0 35 ! iflag_pbl=7 : MY 2.0.Fournier 36 ! iflag_pbl=8/9 : MY 2.5 37 ! iflag_pbl=8 with special obsolete treatments for convergence 38 ! with Cmpi5 NPv3.1 simulations 39 ! iflag_pbl=10/11 : New scheme M2 and N2 explicit and dissiptation exact 40 ! iflag_pbl=12 = 11 with vertical diffusion off q2 41 42 ! 2013/04/01 (FH hourdin@lmd.jussieu.fr) 43 ! Correction for very stable PBLs (iflag_pbl=10 and 11) 44 ! iflag_pbl=8 converges numerically with NPv3.1 45 ! iflag_pbl=11 -> the model starts with NP from start files created by ce0l 46 ! -> the model can run with longer time-steps. 47 ! ....................................................................... 48 72 ! 73 !....................................................................... 74 75 !======================================================================= 76 ! Declarations: 77 !======================================================================= 78 79 80 ! Inputs/Outputs 81 !---------------- 49 82 REAL dt, g, rconst 50 83 REAL plev(klon, klev+1), temp(klon, klev) … … 57 90 REAL teta(klon, klev) 58 91 REAL cd(klon) 59 REAL q2(klon, klev+1) , qpre92 REAL q2(klon, klev+1) 60 93 REAL unsdz(klon, klev) 61 94 REAL unsdzdec(klon, klev+1) 62 95 REAL kn(klon, klev+1) 63 96 REAL km(klon, klev+1) 64 REAL kmpre(klon, klev+1), tmp2 97 INTEGER iflag_pbl, ngrid, nsrf 98 INTEGER ni(klon) 99 100 101 ! Local 102 !------- 103 104 INCLUDE "clesphys.h" 105 106 107 REAL kmpre(klon, klev+1), tmp2, qpre 65 108 REAL mpre(klon, klev+1) 66 REAL kn(klon, klev+1)67 109 REAL kq(klon, klev+1) 68 110 REAL ff(klon, klev+1), delta(klon, klev+1) 69 111 REAL aa(klon, klev+1), aa0, aa1 70 INTEGER iflag_pbl, ngrid71 112 INTEGER nlay, nlev 72 73 113 LOGICAL first 74 114 INTEGER ipas … … 77 117 DATA first, ipas/.FALSE., 0/ 78 118 !$OMP THREADPRIVATE( first,ipas) 79 REAL,SAVE :: lmixmin=1. 80 !$OMP THREADPRIVATE(lmixmin) 81 82 119 LOGICAL hboville 83 120 INTEGER ig, k 84 85 86 121 REAL ri, zrif, zalpha, zsm, zsn 87 122 REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev) 88 89 123 REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1) 90 124 REAL dtetadz(klon, klev+1) 91 125 REAL m2cstat, mcstat, kmcstat 92 126 REAL l(klon, klev+1) 93 REAL, ALLOCATABLE, SAVE :: l0(:) 94 !$OMP THREADPRIVATE(l0) 95 REAL sq(klon), sqz(klon), zz(klon, klev+1) 127 REAL zz(klon, klev+1) 96 128 INTEGER iter 97 98 129 REAL ric, rifc, b1, kap 99 130 SAVE ric, rifc, b1, kap 100 131 DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.4/ 101 132 !$OMP THREADPRIVATE(ric,rifc,b1,kap) 133 REAL seuilsm, seuilalpha 134 REAL,SAVE :: lmixmin 135 !$OMP THREADPRIVATE(lmixmin) 102 136 REAL frif, falpha, fsm 103 REAL fl, zzz, zl0, zq2, zn2104 105 137 REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), & 106 138 lyam(klon, klev), knyam(klon, klev), w2yam(klon, klev), t2yam(klon, klev) 107 139 LOGICAL, SAVE :: firstcall = .TRUE. 108 140 !$OMP THREADPRIVATE(firstcall) 141 142 143 ! Fonctions utilis??es 144 !-------------------- 145 109 146 frif(ri) = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156)) 110 147 falpha(ri) = 1.318*(0.2231-ri)/(0.2341-ri) 111 148 fsm(ri) = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri)) 112 fl(zzz, zl0, zq2, zn2) = max(min(l0(ig)*kap*zlev(ig, & 113 k)/(kap*zlev(ig,k)+l0(ig)),0.5*sqrt(q2(ig,k))/sqrt( & 114 max(n2(ig,k),1.E-10))), lmixmin) 115 116 117 nlay = klev 118 nlev = klev + 1 149 119 150 120 151 IF (firstcall) THEN 121 ALLOCATE (l0(klon))122 152 firstcall = .FALSE. 153 lmixmin=1. 123 154 CALL getin_p('lmixmin',lmixmin) 124 155 END IF 156 157 158 159 !=============================================================================== 160 ! Flags, tests et ??valuations de constantes 161 !=============================================================================== 162 163 ! On utilise ou non la routine de Holstalg Boville pour les cas tres stables 164 hboville=.TRUE. 125 165 126 166 … … 129 169 END IF 130 170 171 ! Seuil dans le code de turbulence 172 seuilalpha=1.12 173 seuilsm=0.085 174 175 nlay = klev 176 nlev = klev + 1 131 177 ipas = ipas + 1 132 178 133 179 134 ! ....................................................................... 135 ! les increments verticaux 136 ! ....................................................................... 137 138 ! !!!!! allerte !!!!!c 139 ! !!!!! zlev n'est pas declare a nlev !!!!!c 140 ! !!!!! ----> 180 !======================================================================== 181 ! Calcul des increments verticaux 182 !========================================================================= 183 184 185 ! Attention: zlev n'est pas declare a nlev 141 186 DO ig = 1, ngrid 142 187 zlev(ig, nlev) = zlay(ig, nlay) + (zlay(ig,nlay)-zlev(ig,nlev-1)) 143 188 END DO 144 ! !!!!! <---- 145 ! !!!!! allerte !!!!!c 189 146 190 147 191 DO k = 1, nlay … … 162 206 END DO 163 207 164 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 165 ! Computing M^2, N^2, Richardson numbers, stability functions 166 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 167 ! initialize arrays: 168 m2(:, :) = 0.0 169 sm(:, :) = 0.0 170 rif(:, :) = 0.0 171 208 !========================================================================= 209 ! Richardson number and stability functions 210 !========================================================================= 211 212 ! initialize arrays: 213 214 m2(1:ngrid, :) = 0.0 215 sm(1:ngrid, :) = 0.0 216 rif(1:ngrid, :) = 0.0 217 218 !------------------------------------------------------------ 172 219 DO k = 2, klev 220 173 221 DO ig = 1, ngrid 174 222 dz(ig, k) = zlay(ig, k) - zlay(ig, k-1) … … 177 225 dtetadz(ig, k) = (teta(ig,k)-teta(ig,k-1))/dz(ig, k) 178 226 n2(ig, k) = g*2.*dtetadz(ig, k)/(teta(ig,k-1)+teta(ig,k)) 179 ! n2(ig,k)=0.180 227 ri = n2(ig, k)/max(m2(ig,k), 1.E-10) 181 228 IF (ri<ric) THEN … … 184 231 rif(ig, k) = rifc 185 232 END IF 233 186 234 IF (rif(ig,k)<0.16) THEN 187 235 alpha(ig, k) = falpha(rif(ig,k)) 188 236 sm(ig, k) = fsm(rif(ig,k)) 189 237 ELSE 190 alpha(ig, k) = 1.12191 sm(ig, k) = 0.085238 alpha(ig, k) = seuilalpha 239 sm(ig, k) = seuilsm 192 240 END IF 193 241 zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig,k))*sm(ig, k) … … 196 244 197 245 198 ! ==================================================================== 199 ! Computing the mixing length 200 ! ==================================================================== 201 202 ! Mise a jour de l0 203 IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN 204 205 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 206 ! Iterative computation of l0 207 ! This version is kept for iflag_pbl only for convergence 208 ! with NPv3.1 Cmip5 simulations 209 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 210 211 DO ig = 1, ngrid 212 sq(ig) = 1.E-10 213 sqz(ig) = 1.E-10 214 END DO 215 DO k = 2, klev - 1 216 DO ig = 1, ngrid 217 zq = sqrt(q2(ig,k)) 218 sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1)) 219 sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1)) 220 END DO 221 END DO 222 DO ig = 1, ngrid 223 l0(ig) = 0.2*sqz(ig)/sq(ig) 224 END DO 246 247 ! ==================================================================== 248 ! Computing the mixing length 249 ! ==================================================================== 250 251 252 CALL mixinglength(ni,nsrf,ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, l) 253 254 255 256 !======================================================================= 257 ! DIFFERENT TYPES DE SCHEMA de YAMADA 258 !======================================================================= 259 260 !-------------- 261 ! Yamada 2.0 262 !-------------- 263 IF (iflag_pbl==6) THEN 264 265 DO k = 2, klev 266 q2(1:ngrid, k) = l(1:ngrid, k)**2*zz(1:ngrid, k) 267 END DO 268 269 270 !------------------ 271 ! Yamada 2.Fournier 272 !------------------ 273 274 ELSE IF (iflag_pbl==7) THEN 275 276 277 ! Calcul de l, km, au pas precedent 278 !.................................... 225 279 DO k = 2, klev 226 280 DO ig = 1, ngrid 227 l(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))228 END DO229 END DO230 ! print*,'L0 cas 8 ou 10 ',l0231 232 ELSE233 234 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!235 ! In all other case, the assymptotic mixing length l0 is imposed (100m)236 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!237 238 l0(:) = 150.239 DO k = 2, klev240 DO ig = 1, ngrid241 l(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))242 END DO243 END DO244 ! print*,'L0 cas autres ',l0245 246 END IF247 248 249 ! ====================================================================250 ! Yamada 2.0251 ! ====================================================================252 IF (iflag_pbl==6) THEN253 254 DO k = 2, klev255 q2(:, k) = l(:, k)**2*zz(:, k)256 END DO257 258 259 ELSE IF (iflag_pbl==7) THEN260 ! ====================================================================261 ! Yamada 2.Fournier262 ! ====================================================================263 264 ! Calcul de l, km, au pas precedent265 DO k = 2, klev266 DO ig = 1, ngrid267 ! print*,'SMML=',sm(ig,k),l(ig,k)268 281 delta(ig, k) = q2(ig, k)/(l(ig,k)**2*sm(ig,k)) 269 282 kmpre(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k) 270 283 mpre(ig, k) = sqrt(m2(ig,k)) 271 ! print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)272 284 END DO 273 285 END DO … … 278 290 mcstat = sqrt(m2cstat) 279 291 280 ! print*,'M2 L=',k,mpre(ig,k),mcstat 281 282 ! -----{puis on ecrit la valeur de q qui annule l'equation de m 283 ! supposee en q3} 292 ! Puis on ecrit la valeur de q qui annule l'equation de m supposee en q3 293 !......................................................................... 284 294 285 295 IF (k==2) THEN … … 293 303 (unsdz(ig,k)+unsdz(ig,k-1)) 294 304 END IF 295 ! print*,'T2 L=',k,tmp2 305 296 306 tmp2 = kmcstat/(sm(ig,k)/q2(ig,k))/l(ig, k) 297 307 q2(ig, k) = max(tmp2, 1.E-12)**(2./3.) 298 ! print*,'Q2 L=',k,q2(ig,k)299 308 300 309 END DO 301 310 END DO 302 311 312 313 ! ------------------------ 314 ! Yamada 2.5 a la Didi 315 !------------------------- 316 303 317 ELSE IF (iflag_pbl==8 .OR. iflag_pbl==9) THEN 304 ! ==================================================================== 305 ! Yamada 2.5 a la Didi 306 ! ==================================================================== 307 308 309 ! Calcul de l, km, au pas precedent 318 319 ! Calcul de l, km, au pas precedent 320 !.................................... 310 321 DO k = 2, klev 311 322 DO ig = 1, ngrid 312 ! print*,'SMML=',sm(ig,k),l(ig,k)313 323 delta(ig, k) = q2(ig, k)/(l(ig,k)**2*sm(ig,k)) 314 324 IF (delta(ig,k)<1.E-20) THEN 315 ! print*,'ATTENTION L=',k,' Delta=',delta(ig,k)316 325 delta(ig, k) = 1.E-20 317 326 END IF … … 319 328 aa0 = (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1) 320 329 aa1 = (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1) 321 ! abder print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)322 330 aa(ig, k) = aa1*dt/(delta(ig,k)*l(ig,k)) 323 ! print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)324 331 qpre = sqrt(q2(ig,k)) 325 ! if (iflag_pbl.eq.8 ) then326 332 IF (aa(ig,k)>0.) THEN 327 333 q2(ig, k) = (qpre+aa(ig,k)*qpre*qpre)**2 … … 337 343 ! endif 338 344 q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4) 339 ! print*,'Q2 L=',k,q2(ig,k),qpre*qpre340 345 END DO 341 346 END DO … … 343 348 ELSE IF (iflag_pbl>=10) THEN 344 349 345 ! print*,'Schema mixte D'346 ! print*,'Longueur ',l(:,:)347 350 DO k = 2, klev - 1 348 DO ig = 1, ngrid 349 l(ig, k) = max(l(ig,k), lmixmin) 350 km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k) 351 q2(ig, k) = q2(ig, k) + dt*km(ig, k)*m2(ig, k)*(1.-rif(ig,k)) 352 q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4) 353 q2(ig, k) = 1./(1./sqrt(q2(ig,k))+dt/(2*l(ig,k)*b1)) 354 q2(ig, k) = q2(ig, k)*q2(ig, k) 355 END DO 351 km(1:ngrid, k) = l(1:ngrid, k)*sqrt(q2(1:ngrid,k))*sm(1:ngrid, k) 352 q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k)) 353 q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4) 354 q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1)) 355 q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k) 356 356 END DO 357 357 … … 362 362 END IF ! Fin du cas 8 363 363 364 ! print*,'OK8'365 364 366 365 ! ==================================================================== 367 ! Calcul des coefficients de m �ange366 ! Calcul des coefficients de melange 368 367 ! ==================================================================== 368 369 369 DO k = 2, klev 370 ! print*,'k=',k371 370 DO ig = 1, ngrid 372 ! abde print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)373 371 zq = sqrt(q2(ig,k)) 374 km(ig, k) = l(ig, k)*zq*sm(ig, k) 375 kn(ig, k) = km(ig, k)*alpha(ig, k) 376 kq(ig, k) = l(ig, k)*zq*0.2 377 ! print*,'KML=',km(ig,k),kn(ig,k) 378 END DO 379 END DO 372 km(ig, k) = l(ig, k)*zq*sm(ig, k) ! For momentum 373 kn(ig, k) = km(ig, k)*alpha(ig, k) ! For scalars 374 kq(ig, k) = l(ig, k)*zq*0.2 ! For TKE 375 END DO 376 END DO 377 378 379 !==================================================================== 380 ! Transport diffusif vertical de la TKE par la TKE 381 !==================================================================== 382 383 380 384 ! initialize near-surface and top-layer mixing coefficients 381 kq(1:ngrid, 1) = kq(1:ngrid, 2) ! constant (ie no gradient) near the surface 382 kq(1:ngrid, klev+1) = 0 ! zero at the top 383 384 ! Transport diffusif vertical de la TKE. 385 !........................................................... 386 387 kq(1:ngrid, 1) = kq(1:ngrid, 2) ! constant (ie no gradient) near the surface 388 kq(1:ngrid, klev+1) = 0 ! zero at the top 389 390 ! Transport diffusif vertical de la TKE. 391 !....................................... 392 385 393 IF (iflag_pbl>=12) THEN 386 ! print*,'YAMADA VDIF' 387 q2(:, 1) = q2(:, 2) 394 q2(1:ngrid, 1) = q2(1:ngrid, 2) 388 395 CALL vdif_q2(dt, g, rconst, ngrid, plev, temp, kq, q2) 389 396 END IF 390 397 391 ! Traitement des cas noctrunes avec l'introduction d'une longueur 392 ! minilale. 393 394 ! ==================================================================== 395 ! Traitement particulier pour les cas tres stables. 396 ! D'apres Holtslag Boville. 398 399 !==================================================================== 400 ! Traitement particulier pour les cas tres stables, introduction d'une 401 ! longueur de m??lange minimale 402 !==================================================================== 403 ! 404 ! Reference: Local versus Nonlocal boundary-layer diffusion in a global climate model 405 ! Holtslag A.A.M. and Boville B.A. 406 ! J. Clim., 6, 1825-1842, 1993 407 408 409 IF (hboville) THEN 410 397 411 398 412 IF (prt_level>1) THEN 399 413 PRINT *, 'YAMADA4 0' 400 END IF !(prt_level>1) THEN 414 END IF 415 401 416 DO ig = 1, ngrid 402 417 coriol(ig) = 1.E-4 … … 404 419 END DO 405 420 406 ! print*,'pblhmin ',pblhmin407 ! Test a remettre 21 11 02408 ! test abd 13 05 02 if(0.eq.1) then409 421 IF (1==1) THEN 410 422 IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN … … 419 431 END IF 420 432 IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN 421 ! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k) 422 ! s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k) 433 423 434 kn(ig, k) = kmin 424 435 km(ig, k) = kmin 425 436 kq(ig, k) = kmin 426 ! la longueur de melange est suposee etre l= kap z 427 ! K=l q Sm d'ou q2=(K/l Sm)**2 437 438 ! la longueur de melange est suposee etre l= kap z 439 ! K=l q Sm d'ou q2=(K/l Sm)**2 440 428 441 q2(ig, k) = (qmin/sm(ig,k))**2 429 442 END IF … … 432 445 433 446 ELSE 434 435 447 DO k = 2, klev 436 448 DO ig = 1, ngrid … … 442 454 END IF 443 455 IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN 444 ! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)445 ! s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)446 456 kn(ig, k) = kmin 447 457 km(ig, k) = kmin 448 458 kq(ig, k) = kmin 449 450 459 ! la longueur de melange est suposee etre l= kap z 460 ! K=l q Sm d'ou q2=(K/l Sm)**2 451 461 sm(ig, k) = 1. 452 462 alpha(ig, k) = 1. … … 463 473 END IF 464 474 475 END IF ! hboville 476 465 477 IF (prt_level>1) THEN 466 478 PRINT *, 'YAMADA4 1' 467 479 END IF !(prt_level>1) THEN 480 481 482 !====================================================== 483 ! Estimations de w'2 et T'2 d'apres Abdela et McFarlane 484 !====================================================== 485 ! 486 ! Reference: A New Second-Order Turbulence Closure Scheme for the Planetary Boundary Layer 487 ! Abdella K and McFarlane N 488 ! J. Atmos. Sci., 54, 1850-1867, 1997 489 468 490 ! Diagnostique pour stokage 491 !.......................... 469 492 470 493 IF (1==0) THEN … … 482 505 knyam(1:ngrid, 2:klev) = kn(1:ngrid, 2:klev) 483 506 484 ! Estimations de w'2 et T'2 d'apres Abdela et McFarlane 507 508 ! Calcul de w'2 et T'2 509 !....................... 485 510 486 511 w2yam(1:ngrid, 2:klev) = q2(1:ngrid, 2:klev)*0.24 + & … … 493 518 END IF 494 519 495 ! print*,'OKFIN' 520 !============================================================================ 521 496 522 first = .FALSE. 497 523 RETURN 524 525 498 526 END SUBROUTINE yamada4 527 528 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 529 530 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 499 531 SUBROUTINE vdif_q2(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2) 532 500 533 USE dimphy 501 534 IMPLICIT NONE 502 503 ! dt : pas de temps 535 536 include "dimensions.h" 537 538 ! vdif_q2: subroutine qui calcule la diffusion de la TKE par la TKE 539 ! avec un schema implicite en temps avec 540 ! inversion d'un syst??me tridiagonal 541 ! 542 ! Reference: Description of the interface with the surface and 543 ! the computation of the turbulet diffusion in LMDZ 544 ! Technical note on LMDZ 545 ! Dufresne, J-L, Ghattas, J. and Grandpeix, J-Y 546 ! 547 !============================================================================ 548 ! Declarations 549 !============================================================================ 504 550 505 551 REAL plev(klon, klev+1) … … 516 562 INTEGER i, k 517 563 518 ! print*,'RD=',rconst 564 565 !========================================================================= 566 ! Calcul 567 !========================================================================= 568 519 569 DO k = 1, klev 520 570 DO i = 1, ngrid 521 ! test522 ! print*,'i,k',i,k523 ! print*,'temp(i,k)=',temp(i,k)524 ! print*,'(plev(i,k)-plev(i,k+1))=',plev(i,k),plev(i,k+1)525 571 zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k)) 526 572 kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ & … … 546 592 denom(i, k) = deltap(i, k) + (1.-beta(i,k+1))*kstar(i, k) + & 547 593 kstar(i, k-1) 548 ! correction d'un bug 10 01 2001549 594 alpha(i, k) = (q2(i,k)*deltap(i,k)+kstar(i,k)*alpha(i,k+1))/denom(i, k) 550 595 beta(i, k) = kstar(i, k-1)/denom(i, k) … … 553 598 554 599 ! Si on recalcule q2(1) 600 !....................... 555 601 IF (1==0) THEN 556 602 DO i = 1, ngrid … … 559 605 END DO 560 606 END IF 561 ! sinon, on peut sauter cette boucle... 607 562 608 563 609 DO k = 2, klev + 1 … … 569 615 RETURN 570 616 END SUBROUTINE vdif_q2 571 SUBROUTINE vdif_q2e(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2) 572 USE dimphy 617 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 618 619 620 621 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 622 SUBROUTINE vdif_q2e(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2) 623 624 USE dimphy 573 625 IMPLICIT NONE 574 626 575 ! dt : pas de temps 627 include "dimensions.h" 628 ! 629 ! vdif_q2e: subroutine qui calcule la diffusion de TKE par la TKE 630 ! avec un schema explicite en temps 631 632 633 !==================================================== 634 ! Declarations 635 !==================================================== 576 636 577 637 REAL plev(klon, klev+1) … … 585 645 REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1) 586 646 INTEGER ngrid 587 588 647 INTEGER i, k 648 649 650 !================================================== 651 ! Calcul 652 !================================================== 589 653 590 654 DO k = 1, klev … … 621 685 RETURN 622 686 END SUBROUTINE vdif_q2e 687 688 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 689 690 691 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 692 693 SUBROUTINE mixinglength(ni, nsrf, ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, lmix) 694 695 696 697 USE dimphy 698 USE phys_state_var_mod, only: zstd, zsig, zmea 699 USE phys_local_var_mod, only: l_mixmin, l_mix 700 701 ! zstd: ecart type de la'altitud e sous-maille 702 ! zmea: altitude moyenne sous maille 703 ! zsig: pente moyenne de le maille 704 705 USE geometry_mod, only: cell_area 706 ! aire_cell: aire de la maille 707 708 IMPLICIT NONE 709 !************************************************************************* 710 ! Subrourine qui calcule la longueur de m??lange dans le sch??ma de turbulence 711 ! avec la formule de Blackadar 712 ! Calcul d'un minimum en fonction de l'orographie sous-maille: 713 ! L'id??e est la suivante: plus il y a de relief, plus il y a du m??lange 714 ! induit par les circulations meso et submeso ??chelles. 715 ! 716 ! References: * The vertical distribution of wind and turbulent exchange in a neutral atmosphere 717 ! Blackadar A.K. 718 ! J. Geophys. Res., 64, No 8, 1962 719 ! 720 ! * An evaluation of neutral and convective planetary boundary-layer parametrisations relative 721 ! to large eddy simulations 722 ! Ayotte K et al 723 ! Boundary Layer Meteorology, 79, 131-175, 1996 724 ! 725 ! 726 ! * Local Similarity in the Stable Boundary Layer and Mixing length Approaches: consistency of concepts 727 ! Van de Wiel B.J.H et al 728 ! Boundary-Lay Meteorol, 128, 103-166, 2008 729 ! 730 ! 731 ! Histoire: 732 !---------- 733 ! * premi??re r??daction, Etienne et Frederic, 09/06/2016 734 ! 735 ! *********************************************************************** 736 737 !================================================================== 738 ! Declarations 739 !================================================================== 740 741 ! Inputs 742 !------- 743 INTEGER ni(klon) ! indice sur la grille original (non restreinte) 744 INTEGER nsrf ! Type de surface 745 INTEGER ngrid ! Nombre de points concern??s sur l'horizontal 746 INTEGER iflag_pbl ! Choix du sch??ma de turbulence 747 REAL pbl_lmixmin_alpha ! on active ou non le calcul de la longueur de melange minimum 748 REAL lmixmin ! Minimum absolu de la longueur de m??lange 749 REAL zlay(klon, klev) ! altitude du centre de la couche 750 REAL zlev(klon, klev+1) ! atitude de l'interface inf??rieure de la couche 751 REAL u(klon, klev) ! vitesse du vent zonal 752 REAL v(klon, klev) ! vitesse du vent meridional 753 REAL q2(klon, klev+1) ! energie cin??tique turbulente 754 REAL n2(klon, klev+1) ! frequence de Brunt-Vaisala 755 756 !In/out 757 !------- 758 759 LOGICAL, SAVE :: firstcall = .TRUE. 760 !$OMP THREADPRIVATE(firstcall) 761 762 ! Outputs 763 !--------- 764 765 REAL lmix(klon, klev+1) ! Longueur de melange 766 767 768 ! Local 769 !------- 770 771 INTEGER ig,jg, k 772 REAL h_oro(klon) 773 REAL hlim(klon) 774 REAL, SAVE :: kap=0.4,kapb=0.4 775 REAL zq 776 REAL sq(klon), sqz(klon) 777 REAL, ALLOCATABLE, SAVE :: l0(:) 778 !$OMP THREADPRIVATE(l0) 779 REAL fl, zzz, zl0, zq2, zn2 780 REAL famorti, zzzz, zh_oro, zhlim 781 REAL l1(klon, klev+1), l2(klon,klev+1) 782 REAL winds(klon, klev) 783 REAL xcell 784 REAL zstdslope(klon) 785 REAL lmax 786 REAL l2strat, l2neutre, extent 787 REAL l2limit(klon) 788 !=============================================================== 789 ! Fonctions utiles 790 !=============================================================== 791 792 ! Calcul de l suivant la formule de Blackadar 1962 adapt??e par Ayotte 1996 793 !.......................................................................... 794 795 fl(zzz, zl0, zq2, zn2) = max(min(l0(ig)*kap*zlev(ig, & 796 k)/(kap*zlev(ig,k)+l0(ig)),0.5*sqrt(q2(ig,k))/sqrt( & 797 max(n2(ig,k),1.E-10))), 1.E-5) 798 799 ! Fonction d'amortissement de la turbulence au dessus de la montagne 800 ! On augmente l'amortissement en diminuant la valeur de hlim (extent) dans le code 801 !..................................................................... 802 803 famorti(zzzz, zh_oro, zhlim)=(-1.)*ATAN((zzzz-zh_oro)/(zhlim-zh_oro))*2./3.1416+1. 804 805 IF (ngrid==0) RETURN 806 807 IF (firstcall) THEN 808 ALLOCATE (l0(klon)) 809 firstcall = .FALSE. 810 END IF 811 812 813 !===================================================================== 814 ! CALCUL de la LONGUEUR de m??lange suivant BLACKADAR: l1 815 !===================================================================== 816 817 818 IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN 819 820 821 ! Iterative computation of l0 822 ! This version is kept for iflag_pbl only for convergence 823 ! with NPv3.1 Cmip5 simulations 824 !................................................................... 825 826 DO ig = 1, ngrid 827 sq(ig) = 1.E-10 828 sqz(ig) = 1.E-10 829 END DO 830 DO k = 2, klev - 1 831 DO ig = 1, ngrid 832 zq = sqrt(q2(ig,k)) 833 sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1)) 834 sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1)) 835 END DO 836 END DO 837 DO ig = 1, ngrid 838 l0(ig) = 0.2*sqz(ig)/sq(ig) 839 END DO 840 DO k = 2, klev 841 DO ig = 1, ngrid 842 l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k)) 843 END DO 844 END DO 845 846 ELSE 847 848 849 ! In all other case, the assymptotic mixing length l0 is imposed (150m) 850 !...................................................................... 851 852 l0(1:ngrid) = 150. 853 DO k = 2, klev 854 DO ig = 1, ngrid 855 l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k)) 856 END DO 857 END DO 858 859 END IF 860 861 !================================================================================= 862 ! CALCUL d'une longueur de melange en fonctions de la topographie sous maille: l2 863 ! si plb_lmixmin_alpha=TRUE et si on se trouve sur de la terre ( pas actif sur les 864 ! glacier, la glace de mer et les oc??ans) 865 !================================================================================= 866 867 l2(1:ngrid,:)=0.0 868 l_mixmin(1:ngrid,:,nsrf)=0. 869 l_mix(1:ngrid,:,nsrf)=0. 870 871 IF (nsrf .EQ. 1) THEN 872 873 ! coefficients 874 !-------------- 875 876 extent=2. ! On ??tend l'impact du relief jusqu'?? extent*h, extent >1. 877 lmax=150. ! Longueur de m??lange max dans l'absolu 878 879 ! calculs 880 !--------- 881 882 DO ig=1,ngrid 883 884 ! On calcule la hauteur du relief 885 !................................. 886 ! On ne peut pas prendre zstd seulement pour caracteriser le relief sous maille 887 ! car sur un terrain pentu mais sans relief, zstd est non nul (comme en Antarctique, C. Genthon) 888 ! On corrige donc zstd avec l'ecart type de l'altitude dans la maille sans relief 889 ! (en gros, une maille de taille xcell avec une pente constante zstdslope) 890 jg=ni(ig) 891 ! IF (zsig(jg) .EQ. 0.) THEN 892 ! zstdslope(ig)=0. 893 ! ELSE 894 ! xcell=sqrt(cell_area(jg)) 895 ! zstdslope(ig)=max((xcell*zsig(jg)-zmea(jg))**3 /(3.*zsig(jg)),0.) 896 ! zstdslope(ig)=sqrt(zstdslope(ig)) 897 ! END IF 898 899 ! h_oro(ig)=max(zstd(jg)-zstdslope(ig),0.) ! Hauteur du relief 900 h_oro(ig)=zstd(jg) 901 hlim(ig)=extent*h_oro(ig) 902 ENDDO 903 904 l2limit(1:ngrid)=0. 905 906 DO k=2,klev 907 DO ig=1,ngrid 908 winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2) 909 IF (zlev(ig,k) .LE. h_oro(ig)) THEN ! sous l'orographie 910 l2strat= kapb*pbl_lmixmin_alpha*winds(ig,k)/sqrt(max(n2(ig,k),1.E-10)) ! si stratifi??, amplitude d'oscillation * kappab (voir Van de Wiel et al 2008) 911 l2neutre=kap*zlev(ig,k)*h_oro(ig)/(kap*zlev(ig,k)+h_oro(ig)) ! Dans le cas neutre, formule de blackadar. tend asymptotiquement vers h 912 l2neutre=MIN(l2neutre,lmax) ! On majore par lmax 913 l2limit(ig)=MIN(l2neutre,l2strat) ! Calcule de l2 (minimum de la longueur en cas neutre et celle en situation stratifi??e) 914 l2(ig,k)=l2limit(ig) 915 916 ELSE IF (zlev(ig,k) .LE. hlim(ig)) THEN ! Si on est au dessus des montagnes, mais affect?? encore par elles 917 918 ! Au dessus des montagnes, on prend la l2limit au sommet des montagnes 919 ! (la derni??re calcul??e dans la boucle k, vu que k est un indice croissant avec z) 920 ! et on multiplie l2limit par une fonction qui d??croit entre h et hlim 921 l2(ig,k)=l2limit(ig)*famorti(zlev(ig,k),h_oro(ig), hlim(ig)) 922 ELSE ! Au dessus de extent*h, on prend l2=l0 923 l2(ig,k)=0. 924 END IF 925 ENDDO 926 ENDDO 927 ENDIF ! pbl_lmixmin_alpha 928 929 !================================================================================== 930 ! On prend le max entre la longueur de melange de blackadar et celle calcul??e 931 ! en fonction de la topographie 932 !=================================================================================== 933 934 935 DO k=2,klev 936 DO ig=1,ngrid 937 lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin) 938 ENDDO 939 ENDDO 940 941 ! Diagnostics 942 943 DO k=2,klev 944 DO ig=1,ngrid 945 jg=ni(ig) 946 l_mix(jg,k,nsrf)=lmix(ig,k) 947 l_mixmin(jg,k,nsrf)=l2(ig,k) 948 ENDDO 949 ENDDO 950 DO ig=1,ngrid 951 jg=ni(ig) 952 l_mix(jg,1,nsrf)=hlim(ig) 953 ENDDO 954 955 956 957 END SUBROUTINE mixinglength
Note: See TracChangeset
for help on using the changeset viewer.
