Changeset 2127 for trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90
- Timestamp:
- Apr 29, 2019, 10:07:47 AM (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90
r2104 r2127 2 2 ! 3 3 ! 4 SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep, & 5 pplay,pplev,pphi,firstcall, & 6 pu,pv,pt,po, & 7 pduadj,pdvadj,pdtadj,pdoadj, & 8 f0,fm0,entr0,detr0, & 9 zqta,zqla,ztv,ztva,ztla,zthl,zqsatth, & 10 zw2,fraca, & 11 lmin,lmix,lalim,lmax, & 12 zpspsk,ratqscth,ratqsdiff, & 13 Ale_bl,Alp_bl,lalim_conv,wght_th, & 14 !!! nrlmd le 10/04/2012 15 pbl_tke,pctsrf,omega,airephy, & 16 zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0, & 17 n2,s2,ale_bl_stat, & 18 therm_tke_max,env_tke_max, & 19 alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke, & 20 alp_bl_conv,alp_bl_stat) 21 !!! fin nrlmd le 10/04/2012 22 23 24 !============================================================================== 4 SUBROUTINE thermcell_main(ngrid,nlay,nq,ptimestep,firstcall, & 5 pplay,pplev,pphi,zpopsk, & 6 pu,pv,pt,pq, & 7 pduadj,pdvadj,pdtadj,pdqadj, & 8 f0,fm0,entr0,detr0,zw2,fraca, & 9 zqta,zqla,ztv,ztva,zhla,zhl,zqsa, & 10 lmin,lmix,lalim,lmax) 11 12 13 !=============================================================================== 25 14 ! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu 26 15 ! Version du 09.02.07 … … 44 33 ! Introduction of an implicit computation of vertical advection in 45 34 ! the environment of thermal plumes in thermcell_dq 46 ! impl = 0 : explicit, 1 : implicit,-1 : old version35 ! impl = 0 : explicit ; impl = 1 : implicit ; impl =-1 : old version 47 36 ! controled by iflag_thermals = 48 37 ! 15, 16 run with impl=-1 : numerical convergence with NPv3 49 38 ! 17, 18 run with impl=1 : more stable 50 39 ! 15 and 17 correspond to the activation of the stratocumulus "bidouille" 51 ! 52 !============================================================================== 40 ! 41 ! Major changes 2018-19 (AB alexandre.boissinot@lmd.jussieu.fr) 42 ! New detr and entre formulae (no longer alimentation) 43 ! lmin can be greater than 1 44 ! Mix every tracer (EN COURS) 45 ! Old version of thermcell_dq is removed 46 ! Alternative version thermcell_dv2 is removed 47 ! 48 !=============================================================================== 53 49 54 50 USE thermcell_mod 51 USE tracer_h, ONLY: igcm_h2o_vap 55 52 USE print_control_mod, ONLY: lunout, prt_level 56 53 … … 58 55 59 56 60 !============================================================================== 57 !=============================================================================== 61 58 ! Declaration 62 !============================================================================== 63 64 ! inputs: 65 ! ------- 66 67 INTEGER itap 68 INTEGER ngrid 69 INTEGER nlay 59 !=============================================================================== 60 61 ! Inputs: 62 ! ------- 63 64 INTEGER ngrid, nlay, nq 70 65 71 66 REAL ptimestep 72 REAL pplay(ngrid,nlay) 73 REAL pplev(ngrid,nlay+1) 74 REAL pphi(ngrid,nlay) 75 76 REAL p t(ngrid,nlay) !77 REAL p u(ngrid,nlay) !78 REAL p v(ngrid,nlay) !79 REAL p o(ngrid,nlay) !67 REAL pplay(ngrid,nlay) ! Layer pressure 68 REAL pplev(ngrid,nlay+1) ! Level pressure 69 REAL pphi(ngrid,nlay) ! Geopotential 70 71 REAL pu(ngrid,nlay) ! Zonal wind 72 REAL pv(ngrid,nlay) ! Meridional wind 73 REAL pt(ngrid,nlay) ! Temperature 74 REAL pq(ngrid,nlay,nq) ! Tracers mass mixing ratio 80 75 81 76 LOGICAL firstcall 82 77 83 ! outputs: 84 ! -------- 85 86 REAL pdtadj(ngrid,nlay) ! t convective variations 78 ! Outputs: 79 ! -------- 80 87 81 REAL pduadj(ngrid,nlay) ! u convective variations 88 82 REAL pdvadj(ngrid,nlay) ! v convective variations 89 REAL pdoadj(ngrid,nlay) ! water convective variations 90 83 REAL pdtadj(ngrid,nlay) ! t convective variations 84 REAL pdqadj(ngrid,nlay,nq) ! q convective variations 85 86 REAL f0(ngrid) ! mass flux norm (after possible time relaxation) 91 87 REAL fm0(ngrid,nlay+1) ! mass flux (after possible time relaxation) 92 88 REAL entr0(ngrid,nlay) ! entrainment (after possible time relaxation) 93 89 REAL detr0(ngrid,nlay) ! detrainment (after possible time relaxation) 94 REAL f0(ngrid) ! mass flux norm (after possible time relaxation) 95 96 ! local: 97 ! ------ 98 99 INTEGER, save :: dvimpl=1 100 !$OMP THREADPRIVATE(dvimpl) 101 102 INTEGER, save :: dqimpl=-1 103 !$OMP THREADPRIVATE(dqimpl) 104 105 INTEGER, SAVE :: igout=1 106 !$OMP THREADPRIVATE(igout) 107 108 INTEGER, SAVE :: lunout1=6 109 !$OMP THREADPRIVATE(lunout1) 110 111 INTEGER, SAVE :: lev_out=10 112 !$OMP THREADPRIVATE(lev_out) 113 114 INTEGER ig,k,l,ll,ierr 115 INTEGER lmix_bis(ngrid) 116 INTEGER lmax(ngrid) ! 117 INTEGER lmin(ngrid) ! 118 INTEGER lalim(ngrid) ! 119 INTEGER lmix(ngrid) ! 120 121 REAL linter(ngrid) 122 REAL zmix(ngrid) 123 REAL zmax(ngrid) 124 REAL zw2(ngrid,nlay+1) 125 REAL zw_est(ngrid,nlay+1) 126 REAL zmax_sec(ngrid) 127 128 REAL zlay(ngrid,nlay) ! layers altitude 129 REAL zlev(ngrid,nlay+1) ! levels altitude 130 REAL rho(ngrid,nlay) ! layers density 131 REAL rhobarz(ngrid,nlay) ! levels density 132 REAL deltaz(ngrid,nlay) ! layers heigth 133 REAL masse(ngrid,nlay) ! layers mass 134 REAL zpspsk(ngrid,nlay) ! Exner function 135 136 REAL zu(ngrid,nlay) ! environment zonal wind 137 REAL zv(ngrid,nlay) ! environment meridional wind 138 REAL zo(ngrid,nlay) ! environment water tracer 139 REAL zva(ngrid,nlay) ! plume zonal wind 140 REAL zua(ngrid,nlay) ! plume meridional wind 141 REAL zoa(ngrid,nlay) ! plume water tracer 142 143 REAL zt(ngrid,nlay) ! T environment 144 REAL zh(ngrid,nlay) ! T,TP environment 145 REAL zthl(ngrid,nlay) ! TP environment 146 REAL ztv(ngrid,nlay) ! TPV environment ? 147 REAL zl(ngrid,nlay) ! ql environment 148 149 REAL zta(ngrid,nlay) ! 150 REAL zha(ngrid,nlay) ! TRPV plume 151 REAL ztla(ngrid,nlay) ! TP plume 152 REAL ztva(ngrid,nlay) ! TRPV plume 153 REAL ztva_est(ngrid,nlay) ! TRPV plume (temporary) 90 91 ! Local: 92 ! ------ 93 94 INTEGER ig, k, l, iq 95 INTEGER lmax(ngrid) ! Highest layer reached by the plume 96 INTEGER lmix(ngrid) ! Layer in which plume vertical speed is maximal 97 INTEGER lmin(ngrid) ! First unstable layer 98 99 REAL zlay(ngrid,nlay) ! Layers altitudes 100 REAL zlev(ngrid,nlay+1) ! Levels altitudes 101 REAL rho(ngrid,nlay) ! Layers densities 102 REAL rhobarz(ngrid,nlay) ! Levels densities 103 REAL masse(ngrid,nlay) ! Layers masses 104 REAL zpopsk(ngrid,nlay) ! Exner function 105 106 REAL zu(ngrid,nlay) ! u environment 107 REAL zv(ngrid,nlay) ! v environment 108 REAL zt(ngrid,nlay) ! TR environment 109 REAL zqt(ngrid,nlay) ! qt environment 110 REAL zql(ngrid,nlay) ! ql environment 111 REAL zhl(ngrid,nlay) ! TP environment 112 REAL ztv(ngrid,nlay) ! TRPV environment 113 REAL zqs(ngrid,nlay) ! qsat environment 114 115 REAL zua(ngrid,nlay) ! u plume 116 REAL zva(ngrid,nlay) ! v plume 117 REAL zta(ngrid,nlay) ! TR plume 154 118 REAL zqla(ngrid,nlay) ! qv plume 155 119 REAL zqta(ngrid,nlay) ! qt plume 156 157 REAL wmax(ngrid) ! maximal vertical speed 158 REAL wmax_tmp(ngrid) ! temporary maximal vertical speed 159 REAL wmax_sec(ngrid) ! maximal vertical speed if dry case 160 161 REAL fraca(ngrid,nlay+1) ! updraft fraction 162 REAL f_star(ngrid,nlay+1) ! normalized mass flux 163 REAL entr_star(ngrid,nlay) ! normalized entrainment 164 REAL detr_star(ngrid,nlay) ! normalized detrainment 165 REAL alim_star_tot(ngrid) ! integrated alimentation 166 REAL alim_star(ngrid,nlay) ! normalized alimentation 167 REAL alim_star_clos(ngrid,nlay) ! closure alimentation 168 169 REAL fm(ngrid,nlay+1) ! mass flux 170 REAL entr(ngrid,nlay) ! entrainment 171 REAL detr(ngrid,nlay) ! detrainment 172 REAL f(ngrid) ! mass flux norm 173 174 REAL zdthladj(ngrid,nlay) ! 175 REAL lambda ! time relaxation coefficent 176 177 REAL zsortie(ngrid,nlay) 178 REAL zsortie1d(ngrid) 179 REAL susqr2pi, Reuler 180 REAL zf 181 REAL zf2 182 REAL thetath2(ngrid,nlay) 183 REAL wth2(ngrid,nlay) 184 REAL wth3(ngrid,nlay) 185 REAL q2(ngrid,nlay) 186 ! FH : probleme de dimensionnement avec l'allocation dynamique 187 ! common/comtherm/thetath2,wth2 188 real wq(ngrid,nlay) 189 real wthl(ngrid,nlay) 190 real wthv(ngrid,nlay) 191 real ratqscth(ngrid,nlay) 192 real var 193 real vardiff 194 real ratqsdiff(ngrid,nlay) 195 ! niveau de condensation 196 integer nivcon(ngrid) 197 real zcon(ngrid) 198 real CHI 199 real zcon2(ngrid) 200 real pcon(ngrid) 201 real zqsat(ngrid,nlay) 202 real zqsatth(ngrid,nlay) 203 real zlevinter(ngrid) 204 real seuil 205 206 !!! nrlmd le 10/04/2012 207 208 !------Entrees 209 real pbl_tke(ngrid,nlay+1,nbsrf) 210 real pctsrf(ngrid,nbsrf) 211 real omega(ngrid,nlay) 212 real airephy(ngrid) 213 !------Sorties 214 real zlcl(ngrid),fraca0(ngrid),w0(ngrid),w_conv(ngrid) 215 real therm_tke_max0(ngrid) 216 real env_tke_max0(ngrid) 217 real n2(ngrid),s2(ngrid) 218 real ale_bl_stat(ngrid) 219 real therm_tke_max(ngrid,nlay) 220 real env_tke_max(ngrid,nlay) 221 real alp_bl_det(ngrid) 222 real alp_bl_fluct_m(ngrid) 223 real alp_bl_fluct_tke(ngrid) 224 real alp_bl_conv(ngrid) 225 real alp_bl_stat(ngrid) 226 !------Local 227 integer nsrf 228 real rhobarz0(ngrid) ! Densite au LCL 229 logical ok_lcl(ngrid) ! Existence du LCL des thermiques 230 integer klcl(ngrid) ! Niveau du LCL 231 real interp(ngrid) ! Coef d'interpolation pour le LCL 232 !------Triggering 233 real,parameter :: Su=4e4 ! Surface unite: celle d'un updraft elementaire 234 real,parameter :: hcoef=1. ! Coefficient directeur pour le calcul de s2 235 real,parameter :: hmincoef=0.3 ! Coefficient directeur pour l'ordonnee a l'origine pour le calcul de s2 (hmincoef=0.3) 236 real,parameter :: eps1=0.3 ! Fraction de surface occupee par la population 1 : eps1=n1*s1/(fraca0*Sd) 237 real hmin(ngrid) ! Ordonnee a l'origine pour le calcul de s2 238 real zmax_moy(ngrid) ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl) 239 real,parameter :: zmax_moy_coef=0.33 240 real depth(ngrid) ! Epaisseur moyenne du cumulus 241 real w_max(ngrid) ! Vitesse max statistique 242 real s_max(ngrid) 243 !------Closure 244 real pbl_tke_max(ngrid,nlay) ! Profil de TKE moyenne 245 real pbl_tke_max0(ngrid) ! TKE moyenne au LCL 246 real w_ls(ngrid,nlay) ! Vitesse verticale grande echelle (m/s) 247 real,parameter :: coef_m=1. ! On considere un rendement pour alp_bl_fluct_m 248 real,parameter :: coef_tke=1. ! On considere un rendement pour alp_bl_fluct_tke 249 250 !!! fin nrlmd le 10/04/2012 251 252 ! Nouvelles variables pour la convection 253 integer lalim_conv(ngrid) 254 integer n_int(ngrid) 255 real Ale_bl(ngrid) 256 real Alp_bl(ngrid) 257 real alp_int(ngrid) 258 real dp_int(ngrid),zdp 259 real ale_int(ngrid) 260 real fm_tot(ngrid) 261 real wght_th(ngrid,nlay) 262 263 CHARACTER*2 str2 264 CHARACTER*10 str10 120 REAL zhla(ngrid,nlay) ! TP plume 121 REAL ztva(ngrid,nlay) ! TRPV plume 122 REAL zqsa(ngrid,nlay) ! qsat plume 123 124 REAL zqa(ngrid,nlay,nq) ! q plume (ql=0, qv=qt) 125 126 REAL linter(ngrid) ! Level (continuous) of maximal vertical speed 127 REAL zmix(ngrid) ! Altitude of maximal vertical speed 128 REAL zmax(ngrid) ! Maximal altitude reached by the plume 129 REAL wmax(ngrid) ! Maximal vertical speed 130 REAL zw2(ngrid,nlay+1) ! Plume vertical speed 131 132 REAL fraca(ngrid,nlay+1) ! Updraft fraction 133 134 REAL f_star(ngrid,nlay+1) ! Normalized mass flux 135 REAL entr_star(ngrid,nlay) ! Normalized entrainment 136 REAL detr_star(ngrid,nlay) ! Normalized detrainment 137 138 REAL f(ngrid) ! Mass flux norm 139 REAL fm(ngrid,nlay+1) ! Mass flux 140 REAL entr(ngrid,nlay) ! Entrainment 141 REAL detr(ngrid,nlay) ! Detrainment 142 143 REAL lambda ! Time relaxation coefficent 144 145 REAL zdthladj(ngrid,nlay) ! Potential temperature variations 146 REAL dummy(ngrid,nlay) ! Dummy argument for thermcell_dq() 265 147 266 148 CHARACTER (len=20) :: modname='thermcell_main' 267 149 CHARACTER (len=80) :: abort_message 268 150 269 EXTERNAL SCOPY 270 271 !============================================================================== 151 ! AB: I remove alimentation, which is in fact included in entr_star. EN COURS 152 INTEGER lalim(ngrid) ! Highest alimentation level 153 REAL alim_star(ngrid,nlay) ! Normalized alimentation 154 REAL alim_star_tot(ngrid) ! Integrated alimentation 155 REAL alim_star_clos(ngrid,nlay) ! Closure alimentation 156 157 !=============================================================================== 272 158 ! Initialization 273 !============================================================================== 274 275 seuil = 0.25 159 !=============================================================================== 276 160 277 161 IF (firstcall) THEN 278 IF (iflag_thermals==15.or.iflag_thermals==16) THEN279 dvimpl = 0280 dqimpl = -1281 ELSE282 dvimpl = 1283 dqimpl = 1284 ENDIF285 286 162 fm0(:,:) = 0. 287 163 entr0(:,:) = 0. … … 289 165 ENDIF 290 166 167 f_star(:,:) = 0. 168 entr_star(:,:) = 0. 169 detr_star(:,:) = 0. 170 171 f(:) = 0. 172 291 173 fm(:,:) = 0. 292 174 entr(:,:) = 0. 293 175 detr(:,:) = 0. 294 f(:) = 0. 295 176 177 lmax(:) = 1 178 lmix(:) = 1 179 lmin(:) = 1 180 181 pduadj(:,:) = 0.0 182 pdvadj(:,:) = 0.0 183 pdtadj(:,:) = 0.0 184 pdqadj(:,:,:) = 0.0 185 186 ! AB: I remove alimentation, which is in fact included in entr_star. EN COURS 187 alim_star(:,:) = 0. 188 alim_star_tot(:) = 0. 189 alim_star_clos(:,:) = 0. 190 lalim(:) = 1 191 192 ! AB: Careful, hard-coded value from Earth tuned version of the thermal plume model! 296 193 DO ig=1,ngrid 297 194 f0(ig) = max(f0(ig), 1.e-2) … … 304 201 ENDIF 305 202 306 wmax_tmp(:) = 0. 307 308 !------------------------------------------------------------------------------ 203 !=============================================================================== 204 ! Environment settings 205 !=============================================================================== 206 207 !------------------------------------------------------------------------------- 309 208 ! Calcul de T,q,ql a partir de Tl et qt dans l environnement 310 !------------------------------------------------------------------------------ 311 312 CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, & 313 & pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out) 314 315 !------------------------------------------------------------------------------ 316 ! 317 ! -------------------- 318 ! 319 ! 320 ! + + + + + + + + + + + 321 ! 322 ! 323 ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz 324 ! wh,wt,wo ... 325 ! 326 ! + + + + + + + + + + + zh,zu,zv,zo,rho 327 ! 328 ! 329 ! -------------------- zlev(1) 330 ! \\\\\\\\\\\\\\\\\\\\ 331 ! 332 !------------------------------------------------------------------------------ 333 ! Calcul des altitudes des couches 334 !------------------------------------------------------------------------------ 209 !------------------------------------------------------------------------------- 210 211 CALL thermcell_env(ngrid,nlay,nq,pq,pt,pu,pv,pplay,pplev, & 212 & zqt,zql,zt,ztv,zhl,zu,zv,zpopsk,zqs) 213 214 !------------------------------------------------------------------------------- 215 ! Levels and layers altitudes 216 !------------------------------------------------------------------------------- 335 217 336 218 DO l=2,nlay … … 345 227 ENDDO 346 228 347 !------------------------------------------------------------------------------ 348 ! Calcul de l'epaisseur des couches 349 !------------------------------------------------------------------------------ 350 351 DO l=1,nlay 352 deltaz(:,l) = zlev(:,l+1)-zlev(:,l) 353 ENDDO 354 355 !------------------------------------------------------------------------------ 356 ! Calcul des densites 357 !------------------------------------------------------------------------------ 358 359 rho(:,:) = pplay(:,:) / (zpspsk(:,:) * RD * ztv(:,:)) 229 !------------------------------------------------------------------------------- 230 ! Levels and layers densities 231 !------------------------------------------------------------------------------- 232 233 rho(:,:) = pplay(:,:) / (zpopsk(:,:) * RD * ztv(:,:)) 360 234 361 235 IF (prt_level.ge.10) THEN … … 369 243 ENDDO 370 244 371 !------------------------------------------------------------------------------ 372 ! Calcul de la masse373 !------------------------------------------------------------------------------ 245 !------------------------------------------------------------------------------- 246 ! Layers masses 247 !------------------------------------------------------------------------------- 374 248 375 249 DO l=1,nlay … … 377 251 ENDDO 378 252 379 IF (prt_level.ge.1) print *, 'thermcell_main apres initialisation' 380 381 !------------------------------------------------------------------------------ 382 ! 383 ! /|\ 384 ! -------- | F_k+1 ------- 385 ! ----> D_k 386 ! /|\ <---- E_k , A_k 387 ! -------- | F_k --------- 388 ! ----> D_k-1 389 ! <---- E_k-1 , A_k-1 390 ! 391 ! 392 ! 393 ! 394 ! 395 ! --------------------------- 396 ! 397 ! ----- F_lmax+1=0 ---------- \ 398 ! lmax (zmax) | 399 ! --------------------------- | 400 ! | 401 ! --------------------------- | 402 ! | 403 ! --------------------------- | 404 ! | 405 ! --------------------------- | 406 ! | 407 ! --------------------------- | 408 ! | E 409 ! --------------------------- | D 410 ! | 411 ! --------------------------- | 412 ! | 413 ! --------------------------- \ | 414 ! lalim | | 415 ! --------------------------- | | 416 ! | | 417 ! --------------------------- | | 418 ! | A | 419 ! --------------------------- | | 420 ! | | 421 ! --------------------------- | | 422 ! lmin | | 423 ! ----- F_lmin=0 ------------ / / 424 ! 425 ! --------------------------- 426 ! //////////////////////////// 427 ! 428 !------------------------------------------------------------------------------ 429 430 !============================================================================== 431 ! Calculs initiaux ne faisant pas intervenir les changements de phase 432 !============================================================================== 433 434 !------------------------------------------------------------------------------ 435 ! 1. alim_star est le profil vertical de l'alimentation a la base du 436 ! panache thermique, calcule a partir de la flotabilite de l'air sec 437 ! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star 438 ! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un 439 ! panache sec conservatif (e=d=0) alimente selon alim_star 440 ! Il s'agit d'un calcul de type CAPE 441 ! zmax_sec est utilise pour determiner la geometrie du thermique. 442 !------------------------------------------------------------------------------ 443 444 entr_star(:,:) = 0. 445 detr_star(:,:) = 0. 446 alim_star(:,:) = 0. 447 448 alim_star_tot(:) = 0. 449 450 lmin(:) = 1 451 452 !============================================================================== 453 ! Calcul du melange et des variables dans le thermique 454 !============================================================================== 455 456 CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl, & 457 & po,zl,rhobarz,zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & 458 & lalim,f0,detr_star,entr_star,f_star,ztva, & 459 & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis, & 460 & lmin,lev_out,lunout1,igout) 461 462 ! CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ') 463 ! CALL test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ') 464 465 !============================================================================== 253 !=============================================================================== 254 ! Explicative schemes 255 !=============================================================================== 256 257 !------------------------------------------------------------------------------- 258 ! Thermal plume variables 259 !------------------------------------------------------------------------------- 260 261 ! top of the model 262 ! =========================== 263 ! 264 ! --------------------------- 265 ! _ 266 ! ----- F_lmax+1=0 ------zmax \ 267 ! lmax | 268 ! ------F_lmax>0------------- | 269 ! | 270 ! --------------------------- | 271 ! | 272 ! --------------------------- | 273 ! | 274 ! ------------------wmax,zmix | 275 ! lmix | 276 ! --------------------------- | 277 ! | 278 ! --------------------------- | 279 ! | E, D 280 ! --------------------------- | 281 ! | 282 ! --------------------------- rhobarz, f_star, fm, fm0, zw2, fraca 283 ! zt, zu, zv, zo, rho | 284 ! --------------------------- | 285 ! | 286 ! --------------------------- | 287 ! | 288 ! --------------------------- | 289 ! | 290 ! ------F_lmin+1>0----------- | 291 ! lmin | 292 ! ----- F_lmin=0 ------------ _/ 293 ! 294 ! --------------------------- 295 ! 296 ! =========================== 297 ! bottom of the model 298 299 !------------------------------------------------------------------------------- 300 ! Zoom on layers k and k-1 301 !------------------------------------------------------------------------------- 302 303 ! | /|\ | | 304 ! |---- | F_k+1 -----------|--------------------------| level k+1 305 ! | | w_k+1 | | 306 ! | --|--> D_k | 307 ! | | | layer k 308 ! | <--|-- E_k | 309 ! | /|\ | | 310 ! |---- | F_k ----------|-----------------------------| level k 311 ! | | w_k | | 312 ! | --|--> D_k-1 | 313 ! | | | layer k-1 314 ! | <--|-- E_k-1 | 315 ! | /|\ | | 316 ! |---- | F_k-1 -----|--------------------------------| level k-1 317 ! | w_k-1 318 ! 0 fraca 1 319 ! \__________________/ \______________________________/ 320 ! plume (fraca) environment (1-fraca) 321 322 !=============================================================================== 323 ! Thermal plumes computation 324 !=============================================================================== 325 326 !------------------------------------------------------------------------------- 327 ! Thermal plumes speeds, fluxes, tracers and temperatures 328 !------------------------------------------------------------------------------- 329 330 CALL thermcell_plume(ngrid,nlay,nq,ptimestep,ztv, & 331 & zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk, & 332 & detr_star,entr_star,f_star, & 333 & ztva,zhla,zqla,zqta,zta, & 334 & zw2,zqsa,lmix,lmin) 335 336 !------------------------------------------------------------------------------- 466 337 ! Thermal plumes characteristics: zmax, zmix, wmax 467 !============================================================================== 468 469 CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, & 470 & zlev,lmax,zmax,zmix,wmax,f_star,lev_out) 471 472 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 473 ! AB : WARNING: zw2 became its square root in thermcell_height! 474 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 475 476 ! CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ') 477 ! CALL test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ') 478 ! CALL test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ') 479 ! CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ') 480 481 !============================================================================== 338 !------------------------------------------------------------------------------- 339 340 ! AB: Careful, zw2 became its square root in thermcell_height! 341 CALL thermcell_height(ngrid,nlay,lmin,linter,lmix,zw2, & 342 & zlev,lmax,zmax,zmix,wmax,f_star) 343 344 !=============================================================================== 482 345 ! Closure and mass fluxes computation 483 !============================================================================== 484 485 CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, & 486 & lalim,lmin,zmax_sec,wmax_sec,lev_out) 487 488 ! CALL test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ') 489 ! CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ') 490 491 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 492 ! Choix de la fonction d'alimentation utilisee pour la fermeture. 493 ! Apparemment sans importance 494 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 495 alim_star_clos(:,:) = alim_star(:,:) 496 alim_star_clos(:,:) = entr_star(:,:) + alim_star(:,:) 497 498 !------------------------------------------------------------------------------ 499 ! Closure (dry if iflag_thermals_closure=1, moist if iflag_thermals_closure=2) 500 !------------------------------------------------------------------------------ 501 502 IF (iflag_thermals_closure.eq.1) THEN 503 CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev, & 504 & lalim,alim_star_clos,f_star, & 505 & zmax_sec,wmax_sec,f,lev_out) 506 ELSEIF (iflag_thermals_closure.eq.2) THEN 507 CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev, & 508 & lalim,alim_star,f_star, & 509 & zmax,wmax,f,lev_out) 510 ELSE 511 print *, 'ERROR: no closure had been selected!' 512 call abort 513 ENDIF 346 !=============================================================================== 347 348 !------------------------------------------------------------------------------- 349 ! Closure 350 !------------------------------------------------------------------------------- 351 352 CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev, & 353 & lmax,entr_star,zmax,wmax,f) 514 354 515 355 IF (tau_thermals>1.) THEN … … 520 360 ENDIF 521 361 522 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 362 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 523 363 ! Test valable seulement en 1D mais pas genant 524 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 364 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 525 365 IF (.not. (f0(1).ge.0.) ) THEN 526 366 abort_message = '.not. (f0(1).ge.0.)' … … 529 369 ENDIF 530 370 531 !------------------------------------------------------------------------------ 371 !------------------------------------------------------------------------------- 532 372 ! Mass fluxes 533 !------------------------------------------------------------------------------ 373 !------------------------------------------------------------------------------- 534 374 535 375 CALL thermcell_flux(ngrid,nlay,ptimestep,masse, & 536 & lalim,lmin,lmax,alim_star,entr_star,detr_star, & 537 & f,rhobarz,zlev,zw2,fm,entr,detr,zqla, & 538 & lev_out,lunout1,igout) 539 540 ! CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ') 541 ! CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ') 542 543 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 376 ! AB: I remove alimentation, which is in fact included in entr_star. EN COURS 377 ! That is not already done for thermcell_flux. 378 & lalim,lmin,lmax,entr_star,detr_star, & 379 & f,rhobarz,zlev,zw2,fm,entr,detr,zqla) 380 381 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 544 382 ! On ne prend pas directement les profils issus des calculs precedents mais on 545 383 ! s'autorise genereusement une relaxation vers ceci avec une constante de temps 546 ! tau_thermals (typiquement 1800s ).547 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 384 ! tau_thermals (typiquement 1800s sur Terre). 385 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 548 386 549 387 IF (tau_thermals>1.) THEN … … 558 396 ENDIF 559 397 560 !------------------------------------------------------------------------------ 398 !------------------------------------------------------------------------------- 561 399 ! Updraft fraction 562 !------------------------------------------------------------------------------ 400 !------------------------------------------------------------------------------- 563 401 564 402 DO ig=1,ngrid … … 577 415 ENDDO 578 416 579 !============================================================================== 417 !=============================================================================== 580 418 ! Transport vertical 581 !============================================================================== 582 583 !------------------------------------------------------------------------------ 584 ! Calcul du transport vertical (de qt et tp) 585 !------------------------------------------------------------------------------ 586 587 CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, & 588 & zthl,zdthladj,zta,lmin,lev_out) 589 590 CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, & 591 & po,pdoadj,zoa,lmin,lev_out) 419 !=============================================================================== 420 421 !------------------------------------------------------------------------------- 422 ! Calcul du transport vertical de la temperature potentielle 423 !------------------------------------------------------------------------------- 424 425 CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & 426 & zhl,zdthladj,dummy,lmin) 592 427 593 428 DO l=1,nlay 594 429 DO ig=1,ngrid 595 pdtadj(ig,l) = zdthladj(ig,l) * zpspsk(ig,l)430 pdtadj(ig,l) = zdthladj(ig,l) * zpopsk(ig,l) 596 431 ENDDO 597 432 ENDDO 598 433 599 !------------------------------------------------------------------------------ 434 !------------------------------------------------------------------------------- 435 ! Calcul du transport vertical des traceurs 436 !------------------------------------------------------------------------------- 437 438 DO iq=1,nq 439 CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & 440 & pq(:,:,iq),pdqadj(:,:,iq),zqa(:,:,iq),lmin) 441 ENDDO 442 443 !------------------------------------------------------------------------------- 600 444 ! Calcul du transport vertical du moment horizontal 601 !------------------------------------------------------------------------------ 602 603 IF (dvimpl.eq.0) THEN 604 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 605 ! Calcul du transport de V tenant compte d'echange par gradient 606 ! de pression horizontal avec l'environnement 607 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 608 CALL thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca, & 609 & zmax,zu,zv,pduadj,pdvadj,zua,zva,lev_out) 610 ELSE 611 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 612 ! Calcul purement conservatif pour le transport de V=(u,v) 613 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 614 CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, & 615 & zu,pduadj,zua,lmin,lev_out) 616 617 CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, & 618 & zv,pdvadj,zva,lmin,lev_out) 619 ENDIF 620 621 !============================================================================== 622 ! Calculs de diagnostiques pour les sorties 623 !============================================================================== 624 625 IF (sorties) THEN 626 627 !------------------------------------------------------------------------------ 628 ! Calcul du niveau de condensation 629 !------------------------------------------------------------------------------ 630 631 DO ig=1,ngrid 632 nivcon(ig) = 0 633 zcon(ig) = 0. 634 ENDDO 635 !nouveau calcul 636 do ig=1,ngrid 637 CHI = zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1)) 638 pcon(ig) = pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**CHI 639 ENDDO 640 641 !IM do k=1,nlay 642 do k=1,nlay-1 643 do ig=1,ngrid 644 IF ((pcon(ig).le.pplay(ig,k)).and.(pcon(ig).gt.pplay(ig,k+1))) then 645 zcon2(ig) = zlay(ig,k) - (pcon(ig) - pplay(ig,k)) & 646 / (RG * rho(ig,k)) / 100. 647 ENDIF 648 ENDDO 649 ENDDO 650 651 ierr = 0 652 653 do ig=1,ngrid 654 IF (pcon(ig).le.pplay(ig,nlay)) then 655 zcon2(ig) = zlay(ig,nlay) - (pcon(ig) - pplay(ig,nlay)) & 656 / (RG * rho(ig,nlay)) / 100. 657 ierr = 1 658 ENDIF 659 ENDDO 660 661 IF (ierr==1) then 662 write(*,*) 'ERROR: thermal plumes rise the model top' 663 CALL abort 664 ENDIF 665 666 IF (prt_level.ge.1) print*,'14b OK convect8' 667 668 do k=nlay,1,-1 669 do ig=1,ngrid 670 IF (zqla(ig,k).gt.1e-10) then 671 nivcon(ig) = k 672 zcon(ig) = zlev(ig,k) 673 ENDIF 674 ENDDO 675 ENDDO 676 677 IF (prt_level.ge.1) print*,'14c OK convect8' 678 679 !------------------------------------------------------------------------------ 680 ! Calcul des moments 681 !------------------------------------------------------------------------------ 682 683 do l=1,nlay 684 do ig=1,ngrid 685 q2(ig,l) = 0. 686 wth2(ig,l) = 0. 687 wth3(ig,l) = 0. 688 ratqscth(ig,l) = 0. 689 ratqsdiff(ig,l) = 0. 690 ENDDO 691 ENDDO 692 693 IF (prt_level.ge.1) print*,'14d OK convect8' 694 695 do l=1,nlay 696 do ig=1,ngrid 697 zf = fraca(ig,l) 698 zf2 = zf/(1.-zf) 699 thetath2(ig,l) = zf2*(ztla(ig,l)-zthl(ig,l))**2 700 701 IF (zw2(ig,l).gt.1.e-10) then 702 wth2(ig,l) = zf2*(zw2(ig,l))**2 703 else 704 wth2(ig,l) = 0. 705 ENDIF 706 707 wth3(ig,l) = zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l)) & 708 & *zw2(ig,l)*zw2(ig,l)*zw2(ig,l) 709 q2(ig,l) = zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2 710 !test: on calcul q2/po=ratqsc 711 ratqscth(ig,l) = sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.)) 712 ENDDO 713 ENDDO 714 715 !------------------------------------------------------------------------------ 716 ! Calcul des flux: q, thetal et thetav 717 !------------------------------------------------------------------------------ 718 719 do l=1,nlay 720 do ig=1,ngrid 721 wq(ig,l)=fraca(ig,l)*zw2(ig,l)*(zqta(ig,l)*1000.-po(ig,l)*1000.) 722 wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l)) 723 wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l)) 724 ENDDO 725 ENDDO 726 727 CALL thermcell_alp(ngrid,nlay,ptimestep, & 728 & pplay,pplev, & 729 & fm0,entr0,lmax, & 730 & Ale_bl,Alp_bl,lalim_conv,wght_th, & 731 & zw2,fraca,pcon, & 732 & rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax, & 733 & pbl_tke,pctsrf,omega,airephy, & 734 & zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0,& 735 & n2,s2,ale_bl_stat, & 736 & therm_tke_max,env_tke_max, & 737 & alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke, & 738 & alp_bl_conv,alp_bl_stat) 739 740 !------------------------------------------------------------------------------ 741 ! Calcul du ratqscdiff 742 !------------------------------------------------------------------------------ 743 744 var = 0. 745 vardiff = 0. 746 ratqsdiff(:,:) = 0. 747 748 DO l=1,nlay 749 DO ig=1,ngrid 750 IF (l<=lalim(ig)) THEN 751 var = var + alim_star(ig,l) * zqta(ig,l) * 1000. 752 ENDIF 753 ENDDO 754 ENDDO 755 756 IF (prt_level.ge.1) print*,'14f OK convect8' 757 758 DO l=1,nlay 759 DO ig=1,ngrid 760 IF (l<=lalim(ig)) THEN 761 zf = fraca(ig,l) 762 zf2 = zf / (1.-zf) 763 vardiff = vardiff + alim_star(ig,l) & 764 * (zqta(ig,l) * 1000. - var)**2 765 ENDIF 766 ENDDO 767 ENDDO 768 769 IF (prt_level.ge.1) print*,'14g OK convect8' 770 771 DO l=1,nlay 772 DO ig=1,ngrid 773 ratqsdiff(ig,l) = sqrt(vardiff) / (po(ig,l) * 1000.) 774 ENDDO 775 ENDDO 776 777 ENDIF 445 !------------------------------------------------------------------------------- 446 447 CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & 448 & zu,pduadj,zua,lmin) 449 450 CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & 451 & zv,pdvadj,zva,lmin) 778 452 779 453 780 454 RETURN 781 455 END 782 783 784 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!785 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!786 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!787 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!788 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!789 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!790 791 792 SUBROUTINE test_ltherm(ngrid,nlay,pplev,pplay,long,seuil,ztv,po, &793 ztva,zqla,f_star,zw2,comment)794 795 796 USE print_control_mod, ONLY: prt_level797 798 IMPLICIT NONE799 800 801 !==============================================================================802 ! Declaration803 !==============================================================================804 805 ! inputs:806 ! -------807 808 INTEGER ngrid809 INTEGER nlay810 INTEGER long(ngrid)811 812 REAL pplev(ngrid,nlay+1)813 REAL pplay(ngrid,nlay)814 REAL ztv(ngrid,nlay)815 REAL po(ngrid,nlay)816 REAL ztva(ngrid,nlay)817 REAL zqla(ngrid,nlay)818 REAL f_star(ngrid,nlay)819 REAL zw2(ngrid,nlay)820 REAL seuil821 822 CHARACTER*21 comment823 824 ! local:825 ! ------826 827 INTEGER i, k828 829 !==============================================================================830 ! Test831 !==============================================================================832 833 IF (prt_level.ge.1) THEN834 write(*,*) 'WARNING: in test, ', comment835 ENDIF836 837 return838 839 ! test sur la hauteur des thermiques ...840 do i=1,ngrid841 !IMtemp IF (pplay(i,long(i)).lt.seuil*pplev(i,1)) then842 IF (prt_level.ge.10) then843 print *, 'WARNING ',comment,' au point ',i,' K= ',long(i)844 print *, ' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F* W2'845 do k=1,nlay846 write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)847 ENDDO848 ENDIF849 ENDDO850 851 852 RETURN853 END854 855 856 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!857 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!858 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!859 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!860 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!861 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!862 863 864 SUBROUTINE thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &865 rg,pplev,therm_tke_max)866 867 868 !==============================================================================869 ! Calcul du transport verticale dans la couche limite en presence870 ! de "thermiques" explicitement representes871 ! calcul du dq/dt une fois qu'on connait les ascendances872 !873 ! Transport de la TKE par le thermique moyen pour la fermeture en ALP874 ! On transporte pbl_tke pour donner therm_tke875 !==============================================================================876 877 USE print_control_mod, ONLY: prt_level878 879 IMPLICIT NONE880 881 882 !==============================================================================883 ! Declarations884 !==============================================================================885 886 integer ngrid,nlay,nsrf887 888 real ptimestep889 real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)890 real entr0(ngrid,nlay),rg891 real therm_tke_max(ngrid,nlay)892 real detr0(ngrid,nlay)893 894 real masse(ngrid,nlay),fm(ngrid,nlay+1)895 real entr(ngrid,nlay)896 real q(ngrid,nlay)897 898 real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)899 900 real zzm901 902 integer ig,k903 integer isrf904 905 !------------------------------------------------------------------------------906 ! Calcul du detrainement907 !------------------------------------------------------------------------------908 909 do k=1,nlay910 detr0(:,k) = fm0(:,k) - fm0(:,k+1) + entr0(:,k)911 masse0(:,k) = (pplev(:,k) - pplev(:,k+1)) / RG912 ENDDO913 914 !------------------------------------------------------------------------------915 ! Decalage vertical des entrainements et detrainements.916 !------------------------------------------------------------------------------917 918 masse(:,1)=0.5*masse0(:,1)919 entr(:,1)=0.5*entr0(:,1)920 detr(:,1)=0.5*detr0(:,1)921 fm(:,1)=0.922 923 do k=1,nlay-1924 masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))925 entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))926 detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))927 fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)928 ENDDO929 930 fm(:,nlay+1)=0.931 932 !!! nrlmd le 16/09/2010933 ! calcul de la valeur dans les ascendances934 ! do ig=1,ngrid935 ! qa(ig,1) = q(ig,1)936 ! ENDDO937 938 q(:,:)=therm_tke_max(:,:)939 940 do ig=1,ngrid941 qa(ig,1)=q(ig,1)942 ENDDO943 944 IF (1==1) then945 do k=2,nlay946 do ig=1,ngrid947 IF ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.1.e-5*masse(ig,k)) then948 qa(ig,k) = (fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k)) &949 & / (fm(ig,k+1)+detr(ig,k))950 else951 qa(ig,k)=q(ig,k)952 ENDIF953 954 ! IF (qa(ig,k).lt.0.) print *, 'WARNING: qa is negative!'955 ! IF (q(ig,k).lt.0.) print *, 'WARNING: q is negative!'956 ENDDO957 ENDDO958 959 !------------------------------------------------------------------------------960 ! Calcul du flux subsident961 !------------------------------------------------------------------------------962 963 do k=2,nlay964 do ig=1,ngrid965 wqd(ig,k) = fm(ig,k) * q(ig,k)966 IF (wqd(ig,k).lt.0.) print*,'WARNING: wqd is negative!'967 ENDDO968 ENDDO969 970 do ig=1,ngrid971 wqd(ig,1) = 0.972 wqd(ig,nlay+1) = 0.973 ENDDO974 975 !------------------------------------------------------------------------------976 ! Calcul des tendances977 !------------------------------------------------------------------------------978 979 do k=1,nlay980 do ig=1,ngrid981 q(ig,k) = q(ig,k) + ptimestep / masse(ig,k) &982 & * (detr(ig,k) * qa(ig,k) - entr(ig,k) * q(ig,k) &983 & - wqd(ig,k) + wqd(ig,k+1))984 ENDDO985 ENDDO986 ENDIF987 988 therm_tke_max(:,:) = q(:,:)989 990 991 RETURN992 END993
Note: See TracChangeset
for help on using the changeset viewer.
