Changeset 1028 for trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90
- Timestamp:
- Sep 5, 2013, 12:29:13 PM (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90
r682 r1028 1 !======================================================================= 2 ! THERMCELL_MAIN_MARS 3 ! Author : A Colaitis 4 ! 5 ! This routine is called by calltherm_interface and is inside a sub-timestep 6 ! loop. It computes thermals properties from parametrized entrainment and 7 ! detrainment rate as well as the source profile. 8 ! Mass flux are then computed and temperature and CO2 MMR are transported. 9 ! 10 ! This routine is based on the terrestrial version thermcell_main of the 11 ! LMDZ model, written by C. Rio and F. Hourdin 12 !======================================================================= 1 13 ! 2 14 ! … … 4 16 & ,pplay,pplev,pphi,zlev,zlay & 5 17 & ,pu,pv,pt,pq,pq2 & 6 & ,pd uadj,pdvadj,pdtadj,pdqadj,pdq2adj &18 & ,pdtadj,pdqadj & 7 19 & ,fm,entr,detr,lmax,zmax & 8 20 & ,r_aspect & 9 21 & ,zw2,fraca & 10 & ,zpopsk, ztla,heatFlux,heatFlux_down &22 & ,zpopsk,heatFlux,heatFlux_down & 11 23 & ,buoyancyOut, buoyancyEst) 12 24 … … 14 26 15 27 !======================================================================= 16 ! Mars version of thermcell_main. Author : A Colaitis 17 !======================================================================= 18 19 #include "dimensions.h" 20 #include "dimphys.h" 21 #include "comcstfi.h" 22 #include "tracer.h" 23 #include "callkeys.h" 24 28 29 #include "callkeys.h" !contains logical values determining which subroutines and which options are used. 30 ! includes logical "tracer", which is .true. if tracers are present and to be transported 31 ! includes logical "lwrite", which is .true. if one wants more verbose outputs as GCM runs 32 #include "dimensions.h" !contains global GCM grid dimension informations 33 #include "dimphys.h" !similar to dimensions.h, and based on it; includes 34 ! "ngridmx": number of horizontal grid points 35 ! "nlayermx": number of atmospheric layers 36 #include "comcstfi.h" !contains physical constant values such as 37 ! "g" : gravitational acceleration (m.s-2) 38 ! "r" : recuced gas constant (J.K-1.mol-1) 39 #include "tracer.h" !contains tracer information such as 40 ! "nqmx": number of tracers 41 ! "noms()": tracer name 25 42 ! ============== INPUTS ============== 26 43 27 REAL, INTENT(IN) :: ptimestep,r_aspect 28 REAL, INTENT(IN) :: pt(ngridmx,nlayermx) 29 REAL, INTENT(IN) :: pu(ngridmx,nlayermx) 30 REAL, INTENT(IN) :: pv(ngridmx,nlayermx) 31 REAL, INTENT(IN) :: pq(ngridmx,nlayermx,nqmx) 32 REAL, INTENT(IN) :: pq2(ngridmx,nlayermx) 33 REAL, INTENT(IN) :: pplay(ngridmx,nlayermx) 34 REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1) 35 REAL, INTENT(IN) :: pphi(ngridmx,nlayermx) 36 REAL, INTENT(IN) :: zlay(ngridmx,nlayermx) 37 REAL, INTENT(IN) :: zlev(ngridmx,nlayermx+1) 44 REAL, INTENT(IN) :: ptimestep !subtimestep (s) 45 REAL, INTENT(IN) :: r_aspect !aspect ratio (see paragraph 45 of paper and appendix S4) 46 REAL, INTENT(IN) :: pt(ngridmx,nlayermx) !temperature (K) 47 REAL, INTENT(IN) :: pu(ngridmx,nlayermx) !u component of the wind (ms-1) 48 REAL, INTENT(IN) :: pv(ngridmx,nlayermx) !v component of the wind (ms-1) 49 REAL, INTENT(IN) :: pq(ngridmx,nlayermx,nqmx) !tracer concentration (kg/kg) 50 REAL, INTENT(IN) :: pq2(ngridmx,nlayermx) ! Turbulent Kinetic Energy 51 REAL, INTENT(IN) :: pplay(ngridmx,nlayermx) !Pressure at the middle of the layers (Pa) 52 REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1) !intermediate pressure levels (Pa) 53 REAL, INTENT(IN) :: pphi(ngridmx,nlayermx) !Geopotential at the middle of the layers (m2s-2) 54 REAL, INTENT(IN) :: zlay(ngridmx,nlayermx) ! altitude at the middle of the layers 55 REAL, INTENT(IN) :: zlev(ngridmx,nlayermx+1) ! altitude at layer boundaries 38 56 39 57 ! ============== OUTPUTS ============== 40 58 41 REAL, INTENT(OUT) :: pdtadj(ngridmx,nlayermx) 42 REAL :: pduadj(ngridmx,nlayermx) 43 REAL :: pdvadj(ngridmx,nlayermx) 44 REAL :: pdqadj(ngridmx,nlayermx,nqmx) 45 ! REAL, INTENT(OUT) :: pdq2adj(ngridmx,nlayermx) 46 REAL :: pdq2adj(ngridmx,nlayermx) 47 REAL, INTENT(OUT) :: zw2(ngridmx,nlayermx+1) 48 49 ! Diagnostics 59 ! TEMPERATURE 60 REAL, INTENT(OUT) :: pdtadj(ngridmx,nlayermx) !temperature change from thermals dT/dt (K/s) 61 62 ! DIAGNOSTICS 63 REAL, INTENT(OUT) :: zw2(ngridmx,nlayermx+1) ! vertical velocity (m/s) 50 64 REAL, INTENT(OUT) :: heatFlux(ngridmx,nlayermx) ! interface heatflux 51 REAL, INTENT(OUT) :: heatFlux_down(ngridmx,nlayermx) ! interface heat flux from downdraft 52 ! REAL, INTENT(OUT) :: buoyancyOut(ngridmx,nlayermx) ! interlayer buoyancy term 53 ! REAL, INTENT(OUT) :: buoyancyEst(ngridmx,nlayermx) ! interlayer estimated buoyancy term 65 REAL, INTENT(OUT) :: heatFlux_down(ngridmx,nlayermx) ! interface heat flux from downdraft 66 67 ! ============== LOCAL ================ 68 REAL :: pdqadj(ngridmx,nlayermx,nqmx) !tracer change from thermals dq/dt, only for CO2 (the rest can be advected outside of the loop) 54 69 55 70 ! dummy variables when output not needed : 56 71 57 ! REAL :: heatFlux(ngridmx,nlayermx) ! interface heatflux58 ! REAL :: heatFlux_down(ngridmx,nlayermx) ! interface heat flux from downdraft59 72 REAL :: buoyancyOut(ngridmx,nlayermx) ! interlayer buoyancy term 60 73 REAL :: buoyancyEst(ngridmx,nlayermx) ! interlayer estimated buoyancy term 61 74 62 63 75 ! ============== LOCAL ================ 64 76 65 77 INTEGER ig,k,l,ll,iq 66 78 INTEGER lmax(ngridmx),lmin(ngridmx),lalim(ngridmx) 67 REAL linter(ngridmx)68 79 REAL zmax(ngridmx) 69 80 REAL ztva(ngridmx,nlayermx),zw_est(ngridmx,nlayermx+1),ztva_est(ngridmx,nlayermx) 70 REAL zmax_sec(ngridmx)71 81 REAL zh(ngridmx,nlayermx) 72 82 REAL zdthladj(ngridmx,nlayermx) … … 85 95 86 96 REAL wmax(ngridmx) 87 REAL wmax_sec(ngridmx)88 97 REAL fm(ngridmx,nlayermx+1),entr(ngridmx,nlayermx),detr(ngridmx,nlayermx) 89 98 … … 115 124 REAL zdz,zbuoy(ngridmx,nlayermx),zw2m 116 125 LOGICAL activecell(ngridmx),activetmp(ngridmx) 117 REAL a1,b1,ae,be,ad,bd,fdfu,b1inv,a1inv,omega ,adalim126 REAL a1,b1,ae,be,ad,bd,fdfu,b1inv,a1inv,omega 118 127 INTEGER tic 119 128 … … 145 154 REAL f_old,ddd0,eee0,ddd,eee,zzz 146 155 REAL fomass_max,alphamax 147 148 ! =========================================149 150 ! ============= DTETA VARIABLES ===========151 152 ! rien : on prends la divergence du flux turbulent153 154 ! =========================================155 156 ! ============= DV2 VARIABLES =============157 ! not used for now158 159 REAL qa(ngridmx,nlayermx),detr_dv2(ngridmx,nlayermx),zf,zf2160 REAL wvd(ngridmx,nlayermx+1),wud(ngridmx,nlayermx+1)161 REAL gamma0(ngridmx,nlayermx+1),gamma(ngridmx,nlayermx+1)162 REAL ue(ngridmx,nlayermx),ve(ngridmx,nlayermx)163 LOGICAL ltherm(ngridmx,nlayermx)164 REAL dua(ngridmx,nlayermx),dva(ngridmx,nlayermx)165 INTEGER iter166 INTEGER nlarga0167 156 168 157 ! ========================================= … … 183 172 184 173 !----------------------------------------------------------------------- 185 ! initiali sation:174 ! initialization: 186 175 ! --------------- 187 176 188 entr(:,:)=0. 189 detr(:,:)=0. 190 fm(:,:)=0. 191 ! zu(:,:)=pu(:,:) 192 ! zv(:,:)=pv(:,:) 193 zhc(:,:)=pt(:,:)/zpopsk(:,:) 177 entr(:,:)=0. ! entrainment mass flux 178 detr(:,:)=0. ! detrainment mass flux 179 fm(:,:)=0. ! upward mass flux 180 zhc(:,:)=pt(:,:)/zpopsk(:,:) ! potential temperature 194 181 ndt=1 195 182 196 183 ! ********************************************************************** 197 ! Taking into account vertical molar mass gradients 184 ! In order to take into account the effect of vertical molar mass 185 ! gradient on convection, we define modified theta that depends 186 ! on the mass mixing ratio of Co2 in the cell. 187 ! See for details: 188 ! 189 ! Forget, F. and Millour, E. et al. "Non condensable gas enrichment and depletion 190 ! in the martian polar regions", third international workshop on the Mars Atmosphere: 191 ! Modeling and Observations, 1447, 9106. year: 2008 192 ! 193 ! This is especially important for modelling polar convection. 198 194 ! ********************************************************************** 195 196 197 !....................................................................... 198 ! Special treatment for co2 at firstcall: compute coefficients and 199 ! get index of co2 tracer 200 !....................................................................... 199 201 200 202 if(firstcall) then … … 235 237 end if 236 238 237 239 !------------------------------------------------------------------------ 240 ! where are the different quantities defined ? 238 241 !------------------------------------------------------------------------ 239 242 ! -------------------- … … 255 258 256 259 !----------------------------------------------------------------------- 257 ! Calcul des altitudes des couches260 ! Densities at layer and layer interface (see above), mass: 258 261 !----------------------------------------------------------------------- 259 262 260 ! do l=2,nlayermx261 ! zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/g262 ! enddo263 ! zlev(:,1)=0.264 ! zlev(:,nlayermx+1)=(2.*pphi(:,nlayermx)-pphi(:,nlayermx-1))/g265 266 ! zlay(:,:)=pphi(:,:)/g267 !-----------------------------------------------------------------------268 ! Calcul des densites269 !-----------------------------------------------------------------------270 271 263 rho(:,:)=pplay(:,:)/(r*pt(:,:)) 272 264 … … 274 266 275 267 do l=2,nlayermx 276 ! rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))277 268 rhobarz(:,l)=pplev(:,l)/(r*0.5*(pt(:,l)+pt(:,l-1))) 278 269 enddo 279 270 280 ! calcul de la masse271 ! mass computation 281 272 do l=1,nlayermx 282 273 masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g … … 284 275 285 276 277 !----------------------------------------------------------------- 278 ! Schematic representation of an updraft: 286 279 !------------------------------------------------------------------ 287 280 ! … … 330 323 331 324 !============================================================================= 332 ! Calculs initiaux ne faisant pas intervenir les changements de phase 325 ! Mars version: no phase change is considered, we use a "dry" definition 326 ! for the potential temperature. 333 327 !============================================================================= 334 328 335 329 !------------------------------------------------------------------ 336 ! 1. alim_star est le profil vertical de l'alimentation a la base du337 ! panache thermique, calcule a partir de la flotabilite de l'air sec338 ! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star330 ! 1. alim_star is the source layer vertical profile in the lowest layers 331 ! of the thermal plume. Computed from the air buoyancy 332 ! 2. lmin and lalim are the indices of begining and end of source profile 339 333 !------------------------------------------------------------------ 340 334 ! … … 343 337 344 338 !----------------------------------------------------------------------------- 345 ! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un 346 ! panache sec conservatif (e=d=0) alimente selon alim_star 347 ! Il s'agit d'un calcul de type CAPE 348 ! zmax_sec est utilise pour determiner la geometrie du thermique. 349 !------------------------------------------------------------------------------ 350 !--------------------------------------------------------------------------------- 351 !calcul du melange et des variables dans le thermique 352 !-------------------------------------------------------------------------------- 339 ! 3. wmax and zmax are maximum vertical velocity and altitude of a 340 ! conservative plume (entrainment = detrainment = 0) using only 341 ! the source layer. This is a CAPE computation used for determining 342 ! the closure mass flux. 343 !----------------------------------------------------------------------------- 353 344 354 345 ! =========================================================================== … … 356 347 ! =========================================================================== 357 348 358 ! Initialisations des variables reeles 359 ztva(:,:)=ztv(:,:) 360 ztva_est(:,:)=ztva(:,:) 361 ztla(:,:)=0. 362 zdz=0. 363 zbuoy(:,:)=0. 364 w_est(:,:)=0. 365 f_star(:,:)=0. 366 wa_moy(:,:)=0. 367 linter(:)=1. 349 ! Initialization 350 ztva(:,:)=ztv(:,:) ! temperature in the updraft = temperature of the env. 351 ztva_est(:,:)=ztva(:,:) ! estimated temp. in the updraft 352 ztla(:,:)=0. !intermediary variable 353 zdz=0. !layer thickness 354 zbuoy(:,:)=0. !buoyancy 355 w_est(:,:)=0. !estimated vertical velocity 356 f_star(:,:)=0. !non-dimensional upward mass flux f* 357 wa_moy(:,:)=0. !vertical velocity 368 358 369 359 ! -------------------------------------------------------------------------- 370 360 ! -------------- MAIN PARAMETERS FOR THERMALS MODEL ------------------------ 371 ! -------------- see thermiques.pro and getfit.py ------------------------- 372 373 ! a1=2.5 ; b1=0.0015 ; ae=0.045 ; be = 0.6 ! svn baseline 374 375 ! Using broad downdraft selection 376 ! a1=1.60226 ; b1=0.0006 ; ae=0.0454 ; be = 0.57 377 ! ad = 0.0005114 ; bd = -0.662 378 ! fdfu = -1.9 379 380 ! Using conditional sampling downdraft selection 381 a1=1.4716 ; b1=0.0005698 ; ae=0.03683 ; be = 0.57421 382 ad = 0.00048088 ; bd = -0.6697 383 ! fdfu = -1.3 384 fdfu=-0.8 385 a1inv=a1 386 b1inv=b1 387 omega=0. 388 adalim=0. 389 390 ! One good config for 34/35 levels 391 ! a1inv=a1 392 ! b1inv=b1 393 ! be=1.1*be 394 395 ! Best configuration for 222 levels: 396 397 ! omega=0.06 398 ! b1=0. 399 ! a1=1. 400 ! a1inv=0.25*a1 401 ! b1inv=0.0002 402 !! 403 !! 404 !! ae=0.9*ae 405 406 ! Best config for norad 222 levels: 407 ! and more specifically to C.large : 408 409 ! omega=0.06 410 ! omega=0. 411 omega=-0.03 412 ! omega=0. 413 a1=1. 414 ! b1=0. 415 b1=0.0001 416 a1inv=a1 417 be=1.1*be 418 ad = 0.0004 419 ! ad=0.0003 420 ! adalim=0. 421 422 ! b1inv=0.00025 423 b1inv=b1 424 425 ! b1inv = 0.0003 426 427 ! omega=0.06 428 ! Trying stuff : 429 430 ! ad=0.00035 431 ! ae=0.95*ae 432 ! b1=0.00055 433 ! omega=0.04 434 ! 435 ! ad = 0.0003 436 ! ae=0.9*ae 437 438 ! omega=0.04 439 !! b1=0. 440 ! a1=1. 441 ! a1inv=a1 442 ! b1inv=0.0005689 443 !! be=1.1*be 444 !! ae=0.96*ae 445 446 447 ! omega=0.06 448 ! a1=1. 449 ! b1=0. 450 ! be=be 451 ! a1inv=0.25*a1 452 ! b1inv=0.0002 453 ! ad=1.1*ad 454 ! ae=1.*ae 361 362 ! Detrainment 363 ad = 0.0004 !D_2 in paper, see paragraph 44 364 bd = -0.6697 !D_1 in paper, see paragraph 44 365 366 ! Entrainment 367 ae = 0.03683 !E_1 in paper, see paragraph 43 368 be = 0.631631 !E_2 in paper, see paragraph 43 369 370 ! Downdraft 371 fdfu=-0.8 !downdraft to updraft mass flux ratio, see paper paragraph 48 372 omega=-0.03 !see paper paragraph 48 373 374 ! Vertical velocity equation 375 a1=1. !a in paper, see paragraph 41 376 b1=0.0001 !b in paper, see paragraph 41 377 378 ! Inversion layer 379 a1inv=a1 !a1 in inversion layer 380 b1inv=b1 !b1 in inversion layer 455 381 ! -------------------------------------------------------------------------- 456 382 ! -------------------------------------------------------------------------- 457 383 ! -------------------------------------------------------------------------- 458 384 459 ! Initialisation des variables entieres385 ! Some more initializations 460 386 wmaxa(:)=0. 461 387 lalim(:)=1 462 388 463 389 !------------------------------------------------------------------------- 464 ! On ne considere comme actif que les colonnes dont les deux premieres 465 ! couches sont instables. 390 ! We consider as an activecell columns where the two first layers are 391 ! convectively unstable 392 ! When it is the case, we compute the source layer profile (alim_star) 393 ! see paper appendix 4.1 for details on the source layer 466 394 !------------------------------------------------------------------------- 395 467 396 activecell(:)=ztv(:,1)>ztv(:,2) 468 397 do ig=1,ngridmx 469 398 if (ztv(ig,1)>=(ztv(ig,2))) then 470 399 alim_star(ig,1)=MAX((ztv(ig,1)-ztv(ig,2)),0.) & 471 ! & *log(1.+zlev(ig,2))472 400 & *sqrt(zlev(ig,2)) 473 ! & *sqrt(sqrt(zlev(ig,2)))474 ! & /sqrt(zlev(ig,2))475 ! & *zlev(ig,2)476 ! & *exp(-zlev(ig,2)/1000.)477 401 lalim(ig)=2 478 402 alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,1) … … 481 405 482 406 do l=2,nlayermx-1 483 ! do l=2,4 484 do ig=1,ngridmx485 if (ztv(ig,l)>(ztv(ig,l+1)) .and. ztv(ig,1)>=ztv(ig,l) .and. (alim_star(ig,l-1).ne. 0.)) then ! .and. (zlev(ig,l+1) .lt. 1000.)) then407 do ig=1,ngridmx 408 if (ztv(ig,l)>(ztv(ig,l+1)) .and. ztv(ig,1)>=ztv(ig,l) & 409 & .and. (alim_star(ig,l-1).ne. 0.)) then 486 410 alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) & 487 ! & *log(1.+zlev(ig,l+1))488 411 & *sqrt(zlev(ig,l+1)) 489 ! & *sqrt(sqrt(zlev(ig,l+1)))490 ! & /sqrt(zlev(ig,l+1))491 ! & *zlev(ig,l+1)492 ! & *exp(-zlev(ig,l+1)/1000.)493 412 lalim(ig)=l+1 494 413 alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) … … 505 424 506 425 alim_star_tot(:)=1. 507 ! if(alim_star(1,1) .ne. 0.) then 508 ! print*, 'lalim star' ,lalim(1) 509 ! endif 510 511 !------------------------------------------------------------------------------ 512 ! Calcul dans la premiere couche 513 ! On decide dans cette version que le thermique n'est actif que si la premiere 514 ! couche est instable. 515 ! Pourrait etre change si on veut que le thermiques puisse se déclencher 516 ! dans une couche l>1 517 !------------------------------------------------------------------------------ 518 519 do ig=1,ngridmx 520 ! Le panache va prendre au debut les caracteristiques de l'air contenu 521 ! dans cette couche. 426 427 ! We compute the initial squared velocity (zw2) and non-dimensional upward mass flux 428 ! (f_star) in the first and second layer from the source profile. 429 430 do ig=1,ngridmx 522 431 if (activecell(ig)) then 523 432 ztla(ig,1)=ztv(ig,1) 524 !cr: attention, prise en compte de f*(1)=1 => AC : what ? f*(1) =0. ! (d'ou f*(2)=a*(1)525 ! dans un panache conservatif526 433 f_star(ig,1)=0. 527 528 434 f_star(ig,2)=alim_star(ig,1) 529 530 435 zw2(ig,2)=2.*g*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) & 531 436 & *(zlev(ig,2)-zlev(ig,1)) & 532 & *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1)) 437 & *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1)) !0.4=von Karman constant 533 438 w_est(ig,2)=zw2(ig,2) 534 535 439 endif 536 440 enddo 537 441 538 539 442 !============================================================================== 540 !boucle de calcul de la vitesse verticale dans le thermique 443 !============================================================================== 444 !============================================================================== 445 ! LOOP ON VERTICAL LEVELS 541 446 !============================================================================== 542 447 do l=2,nlayermx-1 543 448 !============================================================================== 544 545 546 ! On decide si le thermique est encore actif ou non 547 ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test 449 !============================================================================== 450 !============================================================================== 451 452 453 ! is the thermal plume still active ? 548 454 do ig=1,ngridmx 549 455 activecell(ig)=activecell(ig) & … … 553 459 554 460 !--------------------------------------------------------------------------- 555 ! calcul des proprietes thermodynamiques et de la vitesse de la couche l 556 ! sans tenir compte du detrainement et de l'entrainement dans cette 557 ! couche 558 ! C'est a dire qu'on suppose 559 ! ztla(l)=ztla(l-1) 560 ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer 561 ! avant) a l'alimentation pour avoir un calcul plus propre 461 ! 462 ! .I. INITIALIZATION 463 ! 464 ! Computations of the temperature and buoyancy properties in layer l, 465 ! without accounting for entrainment and detrainment. We are therefore 466 ! assuming constant temperature in the updraft 467 ! 468 ! This computation yields an estimation of the buoyancy (zbuoy) and thereforce 469 ! an estimation of the velocity squared (w_est) 562 470 !--------------------------------------------------------------------------- 563 471 564 472 do ig=1,ngridmx 565 473 if(activecell(ig)) then 566 ! if(l .lt. lalim(ig)) then567 ! ztva_est(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &568 ! & alim_star(ig,l)*ztv(ig,l)) &569 ! & /(f_star(ig,l)+alim_star(ig,l))570 ! else571 474 ztva_est(ig,l)=ztla(ig,l-1) 572 ! endif573 475 574 476 zdz=zlev(ig,l+1)-zlev(ig,l) 575 477 zbuoy(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) 478 479 ! Estimated vertical velocity squared 480 ! (discretized version of equation 12 in paragraph 40 of paper) 576 481 577 482 if (((a1*zbuoy(ig,l)/w_est(ig,l)-b1) .gt. 0.) .and. (w_est(ig,l) .ne. 0.)) then … … 588 493 589 494 !------------------------------------------------- 590 ! calcul des taux d'entrainement et de detrainement495 ! Compute corresponding non-dimensional (ND) entrainment and detrainment rates 591 496 !------------------------------------------------- 592 497 … … 597 502 598 503 if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then 504 505 ! ND entrainment rate, see equation 16 of paper (paragraph 43) 506 599 507 entr_star(ig,l)=f_star(ig,l)*zdz* & 600 508 & MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be) 601 ! & MAX(0.,log(1. + 0.03*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))602 509 else 603 510 entr_star(ig,l)=0. … … 607 514 if(l .lt. lalim(ig)) then 608 515 609 ! detr_star(ig,l)=0. 610 detr_star(ig,l) = f_star(ig,l)*zdz* & 611 & adalim 516 detr_star(ig,l)=0. 612 517 else 613 518 614 ! detr_star(ig,l) = f_star(ig,l)*zdz* & 615 ! & 0.0105*((zbuoy(ig,l)/zw2m)/0.048)**(1./1.7) 616 ! detr_star(ig,l) = f_star(ig,l)*zdz* & 617 ! & 0.0085*((zbuoy(ig,l)/zw2m)/0.05)**(1./1.55) 618 619 ! last baseline from direct les 620 ! detr_star(ig,l) = f_star(ig,l)*zdz* & 621 ! & 0.065*(2.5*(zbuoy(ig,l)/zw2m))**0.75 622 623 ! new param from continuity eq with a fit on dfdz 519 ! ND detrainment rate, see paragraph 44 of paper 520 624 521 detr_star(ig,l) = f_star(ig,l)*zdz* & 625 522 & ad 626 ! & Max(0., 0.0005 - 0.55*zbuoy(ig,l)/zw2m)627 628 ! & MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005) !svn baseline629 ! & MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0008)630 631 ! & 0.014*((zbuoy(ig,l)/zw2m)/0.05)**(1./1.35)632 ! detr_star(ig,l) = f_star(ig,l)*zdz* &633 ! & ((zbuoy(ig,l)/zw2m)/2.222)! + 0.0002)634 523 635 524 endif … … 637 526 detr_star(ig,l)=f_star(ig,l)*zdz* & 638 527 & MAX(ad,bd*zbuoy(ig,l)/zw2m) 639 ! & Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)640 ! & Max(0., Min(0.001,0.0005 - 0.55*zbuoy(ig,l)/zw2m))641 642 643 ! & MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005) !svn baseline644 645 ! & *5.*(-afact*zbetalpha*zbuoy(ig,l)/zw2m)646 ! & *5.*(-afact*zbuoy(ig,l)/zw2m)647 648 ! last baseline from direct les649 ! & 0.065*(-2.5*(zbuoy(ig,l)/zw2m))**0.75650 651 ! new param from continuity eq with a fit on dfdz652 653 528 654 529 endif 655 530 656 ! En dessous de lalim, on prend le max de alim_star et entr_star pour657 ! alim_star et 0 sinon531 ! If we are still in the source layer, we define the source layer entr. rate (alim_star) as the 532 ! maximum between the source entrainment rate and the estimated entrainment rate. 658 533 659 534 if (l.lt.lalim(ig)) then … … 662 537 endif 663 538 664 ! Calcul du flux montant normalise 539 ! Compute the non-dimensional upward mass flux at layer l+1 540 ! using equation 11 of appendix 4.2 in paper 665 541 666 542 f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & … … 670 546 enddo 671 547 672 673 !---------------------------------------------------------------------------- 674 !calcul de la vitesse verticale en melangeant Tl et qt du thermique 675 !--------------------------------------------------------------------------- 676 548 ! ----------------------------------------------------------------------------------- 549 ! 550 ! .II. CONVERGENCE LOOP 551 ! 552 ! We have estimated a vertical velocity profile and refined the source layer profile 553 ! We now conduct iterations to compute: 554 ! 555 ! - the temperature inside the updraft from the estimated entrainment/source, detrainment, 556 ! and upward mass flux. 557 ! - the buoyancy from the new temperature inside the updraft 558 ! - the vertical velocity from the new buoyancy 559 ! - the entr., detr. and upward mass flux from the new buoyancy and vertical velocity 560 ! 561 ! This loop (tic) converges quickly. We have hardcoded 6 iterations from empirical observations. 562 ! Convergence occurs in 1 or 2 iterations in most cases. 563 ! ----------------------------------------------------------------------------------- 564 565 ! ----------------------------------------------------------------------------------- 566 ! ----------------------------------------------------------------------------------- 677 567 DO tic=0,5 ! internal convergence loop 568 ! ----------------------------------------------------------------------------------- 569 ! ----------------------------------------------------------------------------------- 570 571 ! Is the cell still active ? 678 572 activetmp(:)=activecell(:) .and. f_star(:,l+1)>1.e-10 573 574 ! If the cell is active, compute temperature inside updraft 679 575 do ig=1,ngridmx 680 576 if (activetmp(ig)) then … … 687 583 enddo 688 584 585 ! Is the cell still active with respect to temperature variations ? 689 586 activetmp(:)=activetmp(:).and.(abs(ztla(:,l)-ztva(:,l)).gt.0.01) 690 587 588 ! Compute new buoyancu and vertical velocity 691 589 do ig=1,ngridmx 692 590 if (activetmp(ig)) then … … 694 592 zbuoy(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) 695 593 594 ! (discretized version of equation 12 in paragraph 40 of paper) 696 595 if (((a1*zbuoy(ig,l)/zw2(ig,l)-b1) .gt. 0.) .and. (zw2(ig,l) .ne. 0.) ) then 697 596 zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*a1*zbuoy(ig,l)- & … … 704 603 705 604 ! ================ RECOMPUTE ENTR, DETR, and F FROM NEW W2 =================== 605 ! ND entrainment rate, see equation 16 of paper (paragraph 43) 606 ! ND detrainment rate, see paragraph 44 of paper 706 607 707 608 do ig=1,ngridmx … … 713 614 entr_star(ig,l)=f_star(ig,l)*zdz* & 714 615 & MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be) 715 ! & MAX(0.,log(1. + 0.03*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))716 616 else 717 617 entr_star(ig,l)=0. … … 721 621 if(l .lt. lalim(ig)) then 722 622 723 ! detr_star(ig,l)=0. 724 detr_star(ig,l) = f_star(ig,l)*zdz* & 725 & adalim 623 detr_star(ig,l)=0. 726 624 727 625 else 728 626 detr_star(ig,l) = f_star(ig,l)*zdz* & 729 627 & ad 730 ! & Max(0., 0.0005 - 0.55*zbuoy(ig,l)/zw2m)731 628 732 629 endif … … 734 631 detr_star(ig,l)=f_star(ig,l)*zdz* & 735 632 & MAX(ad,bd*zbuoy(ig,l)/zw2m) 736 ! & Max(0.,Min(0.001,0.0005 - 0.55*zbuoy(ig,l)/zw2m))737 633 738 634 endif … … 742 638 endif 743 639 744 ! En dessous de lalim, on prend le max de alim_star et entr_star pour745 ! alim_star et 0 sinon640 ! If we are still in the source layer, we define the source layer entr. rate (alim_star) as the 641 ! maximum between the source entrainment rate and the estimated entrainment rate. 746 642 747 643 if (l.lt.lalim(ig)) then … … 750 646 endif 751 647 752 ! Calcul du flux montant normalise 648 ! Compute the non-dimensional upward mass flux at layer l+1 649 ! using equation 11 of appendix 4.2 in paper 753 650 754 651 f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & … … 757 654 endif 758 655 enddo 759 760 ENDDO ! of tic 656 ! ----------------------------------------------------------------------------------- 657 ! ----------------------------------------------------------------------------------- 658 ENDDO ! of internal convergence loop 659 ! ----------------------------------------------------------------------------------- 660 ! ----------------------------------------------------------------------------------- 761 661 762 662 !--------------------------------------------------------------------------- 763 ! initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max663 ! Miscellaneous computations for height 764 664 !--------------------------------------------------------------------------- 765 665 … … 767 667 if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then 768 668 IF (lwrite) THEN 769 print*,'On tombe sur le cas particulier de thermcell_plume'669 print*,'thermcell_plume, particular case in velocity profile' 770 670 ENDIF 771 671 zw2(ig,l+1)=0. 772 linter(ig)=l+1773 672 endif 774 673 775 674 if (zw2(ig,l+1).lt.0.) then 776 linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &777 & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))778 675 zw2(ig,l+1)=0. 779 676 endif … … 781 678 782 679 if (wa_moy(ig,l+1).gt.wmaxa(ig)) then 783 ! lmix est le niveau de la couche ou w (wa_moy) est maximum784 680 wmaxa(ig)=wa_moy(ig,l+1) 785 681 endif … … 787 683 788 684 !========================================================================= 789 ! FIN DE LA BOUCLE VERTICALE790 enddo791 685 !========================================================================= 792 793 !on recalcule alim_star_tot 686 !========================================================================= 687 ! END OF THE LOOP ON VERTICAL LEVELS 688 enddo 689 !========================================================================= 690 !========================================================================= 691 !========================================================================= 692 693 ! Recompute the source layer total entrainment alim_star_tot 694 ! as alim_star may have been modified in the above loop. Renormalization of 695 ! alim_star. 696 794 697 do ig=1,ngridmx 795 698 alim_star_tot(ig)=0. … … 817 720 ! =========================================================================== 818 721 819 ! Attention, w2 est transforme en sa racine carree dans cette routine722 ! WARNING, W2 (squared velocity) IS TRANSFORMED IN ITS SQUARE ROOT HERE 820 723 821 724 !------------------------------------------------------------------------------- 822 ! C alcul des caracteristiques du thermique:zmax,wmax725 ! Computations of the thermal height zmax and maximum vertical velocity wmax 823 726 !------------------------------------------------------------------------------- 824 727 825 ! calcul de la hauteur max du thermique728 ! Index of the thermal plume height 826 729 do ig=1,ngridmx 827 730 lmax(ig)=lalim(ig) … … 835 738 enddo 836 739 837 ! On traite le cas particulier qu'il faudrait éviter ou le thermique 838 ! atteind le haut du modele ... 740 ! Particular case when the thermal reached the model top, which is not a good sign 839 741 do ig=1,ngridmx 840 742 if ( zw2(ig,nlayermx) > 1.e-10 ) then 841 print*,'WARNING !!!!! W2 thermiques non nul derniere couche'743 print*,'WARNING !!!!! W2 non-zero in last layer' 842 744 lmax(ig)=nlayermx 843 745 endif 844 746 enddo 845 747 846 ! pas de thermique si couche 1 stable 847 ! do ig=1,ngridmx 848 ! if (lmin(ig).gt.1) then 849 ! lmax(ig)=1 850 ! lmin(ig)=1 851 ! lalim(ig)=1 852 ! endif 853 ! enddo 854 ! 855 ! Determination de zw2 max 748 ! Maximum vertical velocity zw2 856 749 do ig=1,ngridmx 857 750 wmax(ig)=0. … … 872 765 enddo 873 766 enddo 874 ! Longueur caracteristique correspondant a la hauteur des thermiques. 767 768 ! Height of the thermal plume, defined as the following: 769 ! zmax=Integral[z*w(z)*dz]/Integral[w(z)*dz] 770 ! 875 771 do ig=1,ngridmx 876 772 zmax(ig)=0. … … 892 788 enddo 893 789 894 ! Attention, w2 est transforme en sa racine carree dans cette routine895 896 790 ! =========================================================================== 897 791 ! ================= FIN HEIGHT ============================================== … … 904 798 endif 905 799 906 ! Choix de la fonction d'alimentation utilisee pour la fermeture. 800 ! alim_star_clos is the source profile used for closure. It consists of the 801 ! modified source profile in the source layers, and the entrainment profile 802 ! above it. 907 803 908 804 alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:) … … 913 809 914 810 !------------------------------------------------------------------------------- 915 ! Fermeture,determination de f811 ! Closure, determination of the upward mass flux 916 812 !------------------------------------------------------------------------------- 917 ! Appel avec la version seche813 ! Init. 918 814 919 815 alim_star2(:)=0. … … 921 817 f(:)=0. 922 818 923 ! Indice vertical max (max de lalim) atteint par les thermiques sur le domaine819 ! llmax is the index of the heighest thermal in the simulation domain 924 820 llmax=1 925 821 do ig=1,ngridmx … … 927 823 enddo 928 824 929 930 ! Calcul des integrales sur la verticale de alim_star et de 931 ! alim_star^2/(rho dz) 825 ! Integral of a**2/(rho* Delta z), see equation 13 of appendix 4.2 in paper 826 932 827 do k=1,llmax-1 933 828 do ig=1,ngridmx … … 940 835 enddo 941 836 942 ! WARNING : MARS MODIF : we have added 2. : ratio of wmax/vmoy 943 ! True ratio is 3.5 but wetake into account the vmoy is the one alimentating 944 ! the thermal, so there are vs=0 into the vmoy... the true vmoy is lower. (a la louche) 945 ! And r_aspect has been changed from 2 to 1.5 from observations 837 ! Closure mass flux, equation 13 of appendix 4.2 in paper 838 946 839 do ig=1,ngridmx 947 840 if (alim_star2(ig)>1.e-10) then 948 ! f(ig)=wmax_sec(ig)*alim_star_tot_clos(ig)/ &949 ! & (max(500.,zmax_sec(ig))*r_aspect*alim_star2(ig))950 841 f(ig)=wmax(ig)*alim_star_tot_clos(ig)/ & 951 842 & (max(500.,zmax(ig))*r_aspect*alim_star2(ig)) … … 964 855 965 856 !------------------------------------------------------------------------------- 966 !deduction des flux 857 ! With the closure mass flux, we can compute the entrainment, detrainment and 858 ! upward mass flux from the non-dimensional ones. 967 859 !------------------------------------------------------------------------------- 968 860 969 fomass_max=0.8 970 alphamax=0.5 971 861 fomass_max=0.8 !maximum mass fraction of a cell that can go upward in an 862 ! updraft 863 alphamax=0.5 !maximum updraft coverage in a cell 864 865 866 ! these variables allow to follow corrections made to the mass flux when lwrite=.true. 972 867 ncorecfm1=0 973 868 ncorecfm2=0 … … 981 876 982 877 !------------------------------------------------------------------------- 983 ! Multipl ication par le flux de masse issu de la femreture878 ! Multiply by the closure mass flux 984 879 !------------------------------------------------------------------------- 985 880 … … 988 883 detr(:,l)=f(:)*detr_star(:,l) 989 884 enddo 885 886 ! Reconstruct the updraft mass flux everywhere 990 887 991 888 do l=1,zlmax … … 1002 899 enddo 1003 900 1004 ! Test provisoire : pour comprendre pourquoi on corrige plein de fois1005 ! le cas fm6, on commence par regarder une premiere fois avant les1006 ! autres corrections.1007 1008 ! do l=1,nlayermx1009 ! do ig=1,ngridmx1010 ! if (detr(ig,l).gt.fm(ig,l)) then1011 ! ncorecfm8=ncorecfm8+11012 ! endif1013 ! enddo1014 ! enddo1015 901 1016 902 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1017 ! FH Version en cours de test; 1018 ! par rapport a thermcell_flux, on fait une grande boucle sur "l" 1019 ! et on modifie le flux avec tous les contr�les appliques d'affilee 1020 ! pour la meme couche 1021 ! Momentanement, on duplique le calcule du flux pour pouvoir comparer 1022 ! les flux avant et apres modif 903 ! 904 ! Now we will reconstruct once again the upward 905 ! mass flux, but we will apply corrections 906 ! in some cases. We can compare to the 907 ! previously computed mass flux (above) 908 ! 909 ! This verification is done level by level 910 ! 1023 911 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1024 912 1025 do l=1,zlmax 913 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 914 915 do l=1,zlmax !loop on the levels 916 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 917 918 ! Upward mass flux at level l+1 1026 919 1027 920 do ig=1,ngridmx … … 1038 931 1039 932 !------------------------------------------------------------------------- 1040 ! Verification de la positivite des flux de masse933 ! Upward mass flux should be positive 1041 934 !------------------------------------------------------------------------- 1042 935 … … 1061 954 enddo 1062 955 1063 ! Les "optimisations" du flux sont desactivecelles : moins de bidouilles1064 ! je considere que celles ci ne sont pas justifiees ou trop delicates1065 ! pour MARS, d'apres des observations LES.1066 956 !------------------------------------------------------------------------- 1067 !Test sur fraca croissant 1068 !------------------------------------------------------------------------- 1069 ! if (iflag_thermals_optflux==0) then 1070 ! do ig=1,ngridmx 1071 ! if (l.ge.lalim(ig).and.l.le.lmax(ig) & 1072 ! & .and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then 1073 !! zzz est le flux en l+1 a frac constant 1074 ! zzz=fm(ig,l)*rhobarz(ig,l+1)*zw2(ig,l+1) & 1075 ! & /(rhobarz(ig,l)*zw2(ig,l)) 1076 ! if (fm(ig,l+1).gt.zzz) then 1077 ! detr(ig,l)=detr(ig,l)+fm(ig,l+1)-zzz 1078 ! fm(ig,l+1)=zzz 1079 ! ncorecfm4=ncorecfm4+1 1080 ! endif 1081 ! endif 1082 ! enddo 1083 ! endif 1084 ! 1085 !------------------------------------------------------------------------- 1086 !test sur flux de masse croissant 1087 !------------------------------------------------------------------------- 1088 ! if (iflag_thermals_optflux==0) then 1089 ! do ig=1,ngridmx 1090 ! if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then 1091 ! f_old=fm(ig,l+1) 1092 ! fm(ig,l+1)=fm(ig,l) 1093 ! detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1) 1094 ! ncorecfm5=ncorecfm5+1 1095 ! endif 1096 ! enddo 1097 ! endif 1098 ! 1099 !------------------------------------------------------------------------- 1100 !detr ne peut pas etre superieur a fm 957 ! Detrainment should be lower than upward mass flux 1101 958 !------------------------------------------------------------------------- 1102 959 … … 1107 964 entr(ig,l)=fm(ig,l+1) 1108 965 1109 ! Dans le cas ou on est au dessus de la couche d'alimentation et que le 1110 ! detrainement est plus fort que le flux de masse, on stope le thermique. 1111 ! endif 966 ! When detrainment is stronger than upward mass flux, and we are above the 967 ! thermal last level, the plume is stopped 1112 968 1113 969 if(l.gt.lmax(ig)) then 1114 ! if(l.gt.lalim(ig)) then1115 970 detr(ig,l)=0. 1116 971 fm(ig,l+1)=0. … … 1123 978 1124 979 !------------------------------------------------------------------------- 1125 ! fm ne peut pas etre negatif980 ! Check again for mass flux positivity 1126 981 !------------------------------------------------------------------------- 1127 982 … … 1135 990 1136 991 !----------------------------------------------------------------------- 1137 ! la fraction couverte ne peut pas etre superieure a1992 ! Fractional coverage should be less than 1 1138 993 !----------------------------------------------------------------------- 1139 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1140 ! FH Partie a revisiter.1141 ! Il semble qu'etaient codees ici deux optiques dans le cas1142 ! F/ (rho *w) > 11143 ! soit limiter la hauteur du thermique en considerant que c'est1144 ! la derniere chouche, soit limiter F a rho w.1145 ! Dans le second cas, il faut en fait limiter a un peu moins1146 ! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin1147 ! dans thermcell_main et qu'il semble de toutes facons deraisonable1148 ! d'avoir des fractions de 1..1149 ! Ci dessous, et dans l'etat actuel, le premier des deux if est1150 ! sans doute inutile.1151 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1152 994 1153 995 do ig=1,ngridmx … … 1164 1006 enddo 1165 1007 1166 ! Fin de la grande boucle sur les niveaux verticaux 1167 enddo 1008 enddo ! on vertical levels 1168 1009 1169 1010 !----------------------------------------------------------------------- 1170 ! On fait en sorte que la quantite totale d'air entraine dans le 1171 ! panache ne soit pas trop grande comparee a la masse de la maille 1011 ! 1012 ! We limit the total mass going from one level to the next, compared to the 1013 ! initial total mass fo the cell 1014 ! 1172 1015 !----------------------------------------------------------------------- 1173 1016 … … 1182 1025 entr(ig,l)=entr(ig,l)-eee 1183 1026 if ( ddd.gt.0.) then 1184 ! l'entrainement est trop fort mais l'exces peut etre compense par une 1185 ! diminution du detrainement) 1027 ! The entrainment is too strong but we can compensate the excess by a detrainment decrease 1186 1028 detr(ig,l)=ddd 1187 1029 else 1188 ! l'entrainement est trop fort mais l'exces doit etre compense en partie1189 ! par un entrainement plus fort dans la couche superieure1030 ! The entrainment is too strong and we compensate the excess by a stronger entrainment 1031 ! in the layer above 1190 1032 if(l.eq.lmax(ig)) then 1191 1033 detr(ig,l)=fm(ig,l)+entr(ig,l) … … 1200 1042 enddo 1201 1043 enddo 1202 ! 1203 ! ddd=detr(ig)-entre 1204 !on s assure que tout s annule bien en zmax 1044 1045 ! Check again that everything cancels at zmax 1205 1046 do ig=1,ngridmx 1206 1047 fm(ig,lmax(ig)+1)=0. … … 1210 1051 1211 1052 !----------------------------------------------------------------------- 1212 ! Impression du nombre de bidouilles qui ont ete necessaires 1053 ! Summary of the number of modifications that were necessary (if lwrite=.true. 1054 ! and only if there were a lot of them) 1213 1055 !----------------------------------------------------------------------- 1214 1056 … … 1237 1079 1238 1080 !------------------------------------------------------------------ 1239 ! calcul du transport vertical1081 ! vertical transport computation 1240 1082 !------------------------------------------------------------------ 1241 1083 1242 1084 ! ------------------------------------------------------------------ 1243 ! Transport de teta dans l'updraft : (remplace thermcell_dq)1085 ! IN THE UPDRAFT 1244 1086 ! ------------------------------------------------------------------ 1245 1087 1246 1088 zdthladj(:,:)=0. 1247 1248 if (1 .eq. 0) then 1249 ! call thermcell_dqup(ngridmx,nlayermx,ptimestep & 1250 ! & ,fm,entr, & 1251 ! & masse,ztv,zdthladj) 1252 else 1253 1089 ! Based on equation 14 in appendix 4.2 1254 1090 1255 1091 do ig=1,ngridmx … … 1270 1106 enddo 1271 1107 1272 endif1273 1274 1108 ! ------------------------------------------------------------------ 1275 ! Prescription des proprietes du downdraft1109 ! DOWNDRAFT PARAMETERIZATION 1276 1110 ! ------------------------------------------------------------------ 1277 1111 … … 1281 1115 if (lmax(ig) .gt. 1) then 1282 1116 do l=1,lmax(ig) 1283 ! if(zlay(ig,l) .le. 0.8*zmax(ig)) then1284 1117 if(zlay(ig,l) .le. zmax(ig)) then 1118 1119 ! see equation 18 of paragraph 48 in paper 1285 1120 fm_down(ig,l) =fm(ig,l)* & 1286 1121 & max(fdfu,-4*max(0.,(zlay(ig,l)/zmax(ig)))-0.6) 1287 1122 endif 1288 1123 1289 ! if(zlay(ig,l) .le. 0.06*zmax(ig)) then1290 ! ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(sqrt((zlay(ig,l)/zmax(ig))/0.122449) - 1.)*(ztva(ig,l)/ztv(ig,l) - 1.)))1291 ! elseif(zlay(ig,l) .le. 0.4*zmax(ig)) then1292 ! ztvd(ig,l)=ztv(ig,l)*max(0.,1.-0.25*(ztva(ig,l)/ztv(ig,l) - 1.))1293 ! elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then1294 ! ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(((zlay(ig,l)/zmax(ig))-0.7)/1.)*(ztva(ig,l)/ztv(ig,l) - 1.)))1295 ! else1296 ! ztvd(ig,l)=ztv(ig,l)1297 ! endif1298 1299 ! if(zlay(ig,l) .le. 0.6*zmax(ig)) then1300 ! ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/179.848 + 0.997832)1301 ! elseif(zlay(ig,l) .le. 0.8*zmax(ig)) then1302 ! ztvd(ig,l)=-ztv(ig,l)*(((zlay(ig,l)/zmax(ig))-171.74)/170.94)1303 ! else1304 ! ztvd(ig,l)=ztv(ig,l)1305 ! endif1306 1307 1308 ! if(zlay(ig,l) .le. 0.8*zmax(ig)) then1309 ! ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/224.81 + 0.997832)1310 ! elseif(zlay(ig,l) .le. zmax(ig)) then1311 ! ztvd(ig,l)=-ztv(ig,l)*(((zlay(ig,l)/zmax(ig))-144.685)/143.885)1312 ! else1313 ! ztvd(ig,l)=ztv(ig,l)1314 ! endif1315 1316 1317 ! if (zbuoy(ig,l) .gt. 0.) then1318 ! ztvd(ig,l)=ztva(ig,l)*0.99981319 !! ztvd(ig,l)=ztv(ig,l)*0.9978321320 !! else1321 !! if(zlay(ig,l) .le. zmax(ig)) then1322 !! ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/299.7 + 0.997832)1323 !! endif1324 ! endif1325 1326 1124 if(zlay(ig,l) .le. zmax(ig)) then 1125 ! see equation 19 of paragraph 49 in paper 1327 1126 ztvd(ig,l)=min(ztv(ig,l),ztv(ig,l)*((zlay(ig,l)/zmax(ig))/400. + 0.997832)) 1328 ! ztvd(ig,l)=min(ztv(ig,l),ztv(ig,l)*((zlay(ig,l)/zmax(ig))/299.7 + 0.997832))1329 1127 else 1330 1128 ztvd(ig,l)=ztv(ig,l) … … 1336 1134 1337 1135 ! ------------------------------------------------------------------ 1338 ! T ransport de teta dans le downdraft : (remplace thermcell_dq)1136 ! TRANSPORT IN DOWNDRAFT 1339 1137 ! ------------------------------------------------------------------ 1340 1138 … … 1344 1142 if(lmax(ig) .gt. 1) then 1345 1143 ! No downdraft in the very-near surface layer, we begin at k=3 1144 ! Based on equation 14 in appendix 4.2 1346 1145 1347 1146 do k=3,lmax(ig) … … 1361 1160 1362 1161 !------------------------------------------------------------------ 1363 ! Calcul de la fraction de l'ascendance1162 ! Final fraction coverage with the clean upward mass flux, computed at interfaces 1364 1163 !------------------------------------------------------------------ 1365 1164 fraca(:,:)=0. … … 1374 1173 enddo 1375 1174 1376 1377 1378 ! ===========================================================================1379 ! ============= DV2 =========================================================1380 ! ===========================================================================1381 ! ROUTINE OVERIDE : ne prends pas en compte le downdraft1382 ! de plus, le gradient de pression horizontal semble tout deregler... A VOIR1383 1384 if (0 .eq. 1) then1385 1386 1175 !------------------------------------------------------------------ 1387 ! calcul du transport vertical du moment horizontal 1388 !------------------------------------------------------------------ 1389 1390 ! Calcul du transport de V tenant compte d'echange par gradient 1391 ! de pression horizontal avec l'environnement 1392 1393 ! calcul du detrainement 1394 !--------------------------- 1395 1396 nlarga0=0. 1397 1398 do k=1,nlayermx 1399 do ig=1,ngridmx 1400 detr_dv2(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k) 1401 enddo 1402 enddo 1403 1404 ! calcul de la valeur dans les ascendances 1405 do ig=1,ngridmx 1406 zua(ig,1)=zu(ig,1) 1407 zva(ig,1)=zv(ig,1) 1408 ue(ig,1)=zu(ig,1) 1409 ve(ig,1)=zv(ig,1) 1410 enddo 1411 1412 gamma(1:ngridmx,1)=0. 1413 do k=2,nlayermx 1414 do ig=1,ngridmx 1415 ltherm(ig,k)=(fm(ig,k+1)+detr_dv2(ig,k))*ptimestep > 1.e-5*masse(ig,k) 1416 if(ltherm(ig,k).and.zmax(ig)>0.) then 1417 gamma0(ig,k)=masse(ig,k) & 1418 & *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) ) & 1419 & *0.5/zmax(ig) & 1420 & *1. 1421 else 1422 gamma0(ig,k)=0. 1423 endif 1424 if (ltherm(ig,k).and.zmax(ig)<=0.) nlarga0=nlarga0+1 1425 enddo 1426 enddo 1427 1428 gamma(:,:)=0. 1429 1430 do k=2,nlayermx 1431 1432 do ig=1,ngridmx 1433 1434 if (ltherm(ig,k)) then 1435 dua(ig,k)=zua(ig,k-1)-zu(ig,k-1) 1436 dva(ig,k)=zva(ig,k-1)-zv(ig,k-1) 1437 else 1438 zua(ig,k)=zu(ig,k) 1439 zva(ig,k)=zv(ig,k) 1440 ue(ig,k)=zu(ig,k) 1441 ve(ig,k)=zv(ig,k) 1442 endif 1443 enddo 1444 1445 1446 ! Debut des iterations 1447 !---------------------- 1448 1449 ! AC WARNING : SALE ! 1450 1451 do iter=1,5 1452 do ig=1,ngridmx 1453 ! Pour memoire : calcul prenant en compte la fraction reelle 1454 ! zf=0.5*(fraca(ig,k)+fraca(ig,k+1)) 1455 ! zf2=1./(1.-zf) 1456 ! Calcul avec fraction infiniement petite 1457 zf=0. 1458 zf2=1. 1459 1460 ! la première fois on multiplie le coefficient de freinage 1461 ! par le module du vent dans la couche en dessous. 1462 ! Mais pourquoi donc ??? 1463 if (ltherm(ig,k)) then 1464 ! On choisit une relaxation lineaire. 1465 ! gamma(ig,k)=gamma0(ig,k) 1466 ! On choisit une relaxation quadratique. 1467 gamma(ig,k)=gamma0(ig,k)*sqrt(dua(ig,k)**2+dva(ig,k)**2) 1468 zua(ig,k)=(fm(ig,k)*zua(ig,k-1) & 1469 & +(zf2*entr(ig,k)+gamma(ig,k))*zu(ig,k)) & 1470 & /(fm(ig,k+1)+detr_dv2(ig,k)+entr(ig,k)*zf*zf2 & 1471 & +gamma(ig,k)) 1472 zva(ig,k)=(fm(ig,k)*zva(ig,k-1) & 1473 & +(zf2*entr(ig,k)+gamma(ig,k))*zv(ig,k)) & 1474 & /(fm(ig,k+1)+detr_dv2(ig,k)+entr(ig,k)*zf*zf2 & 1475 & +gamma(ig,k)) 1476 1477 ! print*,' OUTPUT DV2 !!!!!!!!!!!!!!!!!!!!',k,zua(ig,k),zva(ig,k),zu(ig,k),zv(ig,k),dua(ig,k),dva(ig,k) 1478 dua(ig,k)=zua(ig,k)-zu(ig,k) 1479 dva(ig,k)=zva(ig,k)-zv(ig,k) 1480 ue(ig,k)=(zu(ig,k)-zf*zua(ig,k))*zf2 1481 ve(ig,k)=(zv(ig,k)-zf*zva(ig,k))*zf2 1482 endif 1483 enddo 1484 ! Fin des iterations 1485 !-------------------- 1486 enddo 1487 1488 enddo ! k=2,nlayermx 1489 1490 ! Calcul du flux vertical de moment dans l'environnement. 1491 !--------------------------------------------------------- 1492 do k=2,nlayermx 1493 do ig=1,ngridmx 1494 wud(ig,k)=fm(ig,k)*ue(ig,k) 1495 wvd(ig,k)=fm(ig,k)*ve(ig,k) 1496 enddo 1497 enddo 1498 do ig=1,ngridmx 1499 wud(ig,1)=0. 1500 wud(ig,nlayermx+1)=0. 1501 wvd(ig,1)=0. 1502 wvd(ig,nlayermx+1)=0. 1503 enddo 1504 1505 ! calcul des tendances. 1506 !----------------------- 1507 do k=1,nlayermx 1508 do ig=1,ngridmx 1509 pduadj(ig,k)=((detr_dv2(ig,k)+gamma(ig,k))*zua(ig,k) & 1510 & -(entr(ig,k)+gamma(ig,k))*ue(ig,k) & 1511 & -wud(ig,k)+wud(ig,k+1)) & 1512 & /masse(ig,k) 1513 pdvadj(ig,k)=((detr_dv2(ig,k)+gamma(ig,k))*zva(ig,k) & 1514 & -(entr(ig,k)+gamma(ig,k))*ve(ig,k) & 1515 & -wvd(ig,k)+wvd(ig,k+1)) & 1516 & /masse(ig,k) 1517 enddo 1518 enddo 1519 1520 1521 ! Sorties eventuelles. 1522 !---------------------- 1523 1524 ! if (nlarga0>0) then 1525 ! print*,'WARNING !!!!!! DANS THERMCELL_DV2 ' 1526 ! print*,nlarga0,' points pour lesquels laraga=0. dans un thermique' 1527 ! print*,'Il faudrait decortiquer ces points' 1528 ! endif 1529 1530 ! =========================================================================== 1531 ! ============= FIN DV2 ===================================================== 1532 ! =========================================================================== 1533 1534 else 1535 ! detrmod(:,:)=0. 1536 ! do k=1,zlmax 1537 ! do ig=1,ngridmx 1538 ! detrmod(ig,k)=fm(ig,k)-fm(ig,k+1) & 1539 ! & +entr(ig,k) 1540 ! if (detrmod(ig,k).lt.0.) then 1541 ! entr(ig,k)=entr(ig,k)-detrmod(ig,k) 1542 ! detrmod(ig,k)=0. 1543 ! endif 1544 ! enddo 1545 ! enddo 1546 ! 1547 ! 1548 ! call thermcell_dqup(ngridmx,nlayermx,ptimestep & 1549 ! & ,fm,entr,detrmod, & 1550 ! & masse,zu,pduadj,ndt,zlmax) 1551 ! 1552 ! call thermcell_dqup(ngridmx,nlayermx,ptimestep & 1553 ! & ,fm,entr,detrmod, & 1554 ! & masse,zv,pdvadj,ndt,zlmax) 1555 1556 endif 1557 1558 !------------------------------------------------------------------ 1559 ! calcul du transport vertical de traceurs 1176 ! Transport of C02 Tracer 1560 1177 !------------------------------------------------------------------ 1561 1178 … … 1600 1217 do ig=1,ngridmx 1601 1218 pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)*ratiom(ig,l) 1602 ! pdtadj(ig,l)=zdthladj(ig,l)*zpopsk(ig,l)*ratiom(ig,l) 1603 enddo 1604 enddo 1219 enddo 1220 enddo 1221 1222 ! =========================================================================== 1223 ! ============= FIN TRANSPORT =============================================== 1224 ! =========================================================================== 1225 1605 1226 1606 1227 !------------------------------------------------------------------ 1607 ! calcul du transport vertical de la tke1228 ! Diagnostics for outputs 1608 1229 !------------------------------------------------------------------ 1609 1610 ! modname='tke'1611 ! call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr, &1612 ! & masse,pq2,pdq2adj,ztvd,fm_down,ztv,modname,lmax)1613 1614 ! ===========================================================================1615 ! ============= FIN TRANSPORT ===============================================1616 ! ===========================================================================1617 1618 1619 !------------------------------------------------------------------1620 ! Calculs de diagnostiques pour les sorties1621 !------------------------------------------------------------------1622 ! DIAGNOSTIQUE1623 1230 ! We compute interface values for teta env and th. The last interface 1624 1231 ! value does not matter, as the mass flux is 0 there.
Note: See TracChangeset
for help on using the changeset viewer.
