Ignore:
Timestamp:
Jun 11, 2014, 3:46:46 PM (10 years ago)
Author:
Laurent Fairhead
Message:

Merged trunk changes r1997:2055 into testing branch

Location:
LMDZ5/branches/testing
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/testing

  • LMDZ5/branches/testing/libf/phylmd/thermcell_plume.F90

    r1999 r2056  
    77     &           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
    88     &           ,lev_out,lunout1,igout)
    9 
    109!--------------------------------------------------------------------------
    11 ! thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
    12 ! Last modified : Arnaud Jam 2014/02/11
    13 !                 Better representation of stratocumulus
     10!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
    1411!--------------------------------------------------------------------------
    1512
    … …  
    8380      real zbuoyjam(ngrid,klev)
    8481      real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis
     82      real fact_shell
     83      real ztv1,ztv2,factinv,zinv,zlmel
     84      real ztv_est1,ztv_est2
    8585      real zcor,zdelta,zcvm5,qlbef
    8686      real betalpha,zbetalpha
    … …  
    9898      fact_epsilon=0.002
    9999      betalpha=0.9
    100       afact=2./3.           
     100      afact=2./3.   
     101      fact_shell=0.85       
    101102
    102103      zbetalpha=betalpha/(1.+betalpha)
    … …  
    164165               lalim(ig)=l+1
    165166               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
    166 !              print*,'alim2',l,ztv(ig,l),ztv(ig,l+1),alim_star(ig,l)
     167!               print*,'alim2',l,ztv(ig,l),ztv(ig,l+1),alim_star(ig,l)
    167168            endif
    168169         enddo
    … …  
    234235   call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
    235236    do ig=1,ngrid
    236 !      print*,'active',active(ig),ig,l
     237!       print*,'active',active(ig),ig,l
    237238        if(active(ig)) then
    238239        zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
    … …  
    260261!    &                     Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
    261262!              w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact))
     263
     264!--------------------------------------------------
     265!AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
     266!--------------------------------------------------
     267         if (iflag_thermals_ed==8) then
     268! Ancienne version
    262269             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
    263270    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
    264271    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
    265              if (w_est(ig,l+1).lt.0.) then
     272 
     273! Nouvelle version Arnaud
     274         else
     275              w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     276    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     277    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
     278         endif
     279
     280
     281         if (iflag_thermals_ed<6) then
     282             zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
     283!              fact_epsilon=0.0005/(zalpha+0.025)**0.5
     284!              fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
     285              fact_epsilon=0.0002/(zalpha+0.1)
     286              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
     287              zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
     288              zdw2=afact*zbuoy(ig,l)/fact_epsilon
     289              zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
     290              w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     291    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     292    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
     293
     294         endif
     295!--------------------------------------------------
     296!AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
     297!on fait max(0.0001,.....)
     298!--------------------------------------------------     
     299
     300!             if (w_est(ig,l+1).lt.0.) then
    266301!               w_est(ig,l+1)=zw2(ig,l)
    267                 w_est(ig,l+1)=0.0001
    268              endif
     302!                w_est(ig,l+1)=0.0001
     303!             endif
     304
    269305       endif
    270306    enddo
    … …  
    297333!Modif AJAM
    298334         
    299         lmel=0.1*zlev(ig,l)
     335        lmel=fact_thermals_ed_dz*zlev(ig,l)
     336        zlmel=zlev(ig,l)+lmel
    300337!        lmel=2.5*(zlev(ig,l)-zlev(ig,l-1))
    301338        lt=l+1
    302          do it=1,klev-(l+1)
    303           zdz2=zlev(ig,lt)-zlev(ig,l)
    304           if (zdz2.gt.lmel) then
    305           zdz3=zlev(ig,lt)-zlev(ig,lt-1)
     339        zdz2=zlev(ig,lt)-zlev(ig,l)
     340!--------------------------------------------------
     341!AJ052014: J'ai remplac? la boucle do par un do while
     342! afin de faire moins de calcul dans la boucle
     343!--------------------------------------------------
     344         do while (zdz2.lt.lmel)
     345         lt=lt+1
     346         zdz2=zlev(ig,lt)-zlev(ig,l)
     347         end do
     348
     349         zdz3=zlev(ig,lt)-zlev(ig,lt-1)
     350
     351!--------------------------------------------------
     352!AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
     353! on cherche o? se trouve l'altitude d'inversion
     354! en calculant ztv1 (interpolation de la valeur de
     355! theta au niveau lt en utilisant les niveaux lt-1 et
     356! lt-2) et ztv2 (interpolation avec les niveaux lt+1
     357! et lt+2). Si theta r?ellement calcul?e au niveau lt
     358! comprise entre ztv1 et ztv2, alors il y a inversion
     359! et on calcule son altitude zinv en supposant que ztv(lt)
     360! est une combinaison lin?aire de ztv1 et ztv2.
     361! Ensuite, on calcule la flottabilit? en comparant
     362! la temp?rature de la couche l ? celle de l'air situ?
     363! l+lmel plus haut, ce qui n?cessite de savoir quel fraction
     364! de cet air est au-dessus ou en-dessous de l'inversion   
     365!--------------------------------------------------
     366
     367
     368         if (iflag_thermals_ed.lt.8) then
     369
     370          ztv1=(ztv(ig,lt-1)-ztv(ig,lt-2))*zlev(ig,lt)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
     371    &          +(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
     372    &          /(zlev(ig,lt-1)-zlev(ig,lt-2))
     373
     374          ztv2=(ztv(ig,lt+2)-ztv(ig,lt+1))*zlev(ig,lt)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
     375    &          +(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
     376    &          /(zlev(ig,lt+2)-zlev(ig,lt+1))
     377
     378          if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then 
     379
     380          factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1)
     381          zinv=zlev(ig,lt-1)+factinv*(zlev(ig,lt)-zlev(ig,lt-1))
     382
     383          if (zlmel+0.5*zdz.ge.zinv) then
     384           if (zlmel-0.5*zdz.ge.zinv) then
     385
     386          ztv_est(ig,l)=(ztv(ig,lt+2)-ztv(ig,lt+1))*(zlmel-0.*zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
     387    &          +(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
     388    &          /(zlev(ig,lt+2)-zlev(ig,lt+1))
     389
     390          zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)+(1.-fact_shell)*zbuoy(ig,l)
     391
     392          else
     393
     394          ztv_est1=(ztv(ig,lt+2)-ztv(ig,lt+1))*0.5*(zlmel+zinv+0.5*zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
     395    &          +(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
     396    &          /(zlev(ig,lt+2)-zlev(ig,lt+1))
     397          ztv_est2=(ztv(ig,lt-1)-ztv(ig,lt-2))*0.5*(zinv+zlmel-0.5*zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
     398    &          +(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
     399    &          /(zlev(ig,lt-1)-zlev(ig,lt-2))
     400          zbuoyjam(ig,l)=fact_shell*RG*(((zlmel+0.5*zdz-zinv)/zdz)*(ztva_est(ig,l)- &
     401    &          ztv_est1)/ztv_est1+((zinv-zlmel+0.5*zdz)/zdz)*(ztva_est(ig,l)- &
     402    &          ztv_est2)/ztv_est2)+(1.-fact_shell)*zbuoy(ig,l)
     403
     404           endif
     405
     406          else
     407
     408          ztv_est(ig,l)=(ztv(ig,lt-1)-ztv(ig,lt-2))*(zlmel-0.*zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
     409    &          +(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
     410    &          /(zlev(ig,lt-1)-zlev(ig,lt-2))
     411
     412          zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)+(1.-fact_shell)*zbuoy(ig,l)
     413!         ztv_est1=(ztv(ig,lt+2)-ztv(ig,lt+1))*0.5*(zlmel+zinv+0.5*zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
     414!    &          +(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
     415!    &          /(zlev(ig,lt+2)-zlev(ig,lt+1))
     416!          ztv_est2=(ztv(ig,lt-1)-ztv(ig,lt-2))*0.5*(zinv+zlmel-0.5*zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
     417!    &          +(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
     418!    &          /(zlev(ig,lt-1)-zlev(ig,lt-2))
     419!          zbuoyjam(ig,l)=fact_shell*RG*(((zlmel+0.5*zdz-zinv)/zdz)*(ztva_est(ig,l)- &
     420!    &          ztv_est1)/ztv_est1+((zinv-zlmel+0.5*zdz)/zdz)*(ztva_est(ig,l)- &
     421!    &          ztv_est2)/ztv_est2)+(1.-fact_shell)*zbuoy(ig,l)
     422
     423
     424
     425
     426!          print*,'on est pass? par l?',l,lt,zbuoyjam(ig,l),zbuoy(ig,l)
     427          endif
     428
     429
     430          else
     431         
    306432!          ztv_est(ig,l)=(lmel/zdz2)*(ztv(ig,lt)-ztv(ig,l))+ztv(ig,l)
    307433!          zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)
     434
     435         zbuoyjam(ig,l)=fact_shell*RG*(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- &
     436    &          ztv(ig,lt))/ztv(ig,lt)+((zdz2-lmel)/zdz3)*(ztva_est(ig,l)- &
     437    &          ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
     438
     439!         zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
     440!    &          po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
     441!     &          po(ig,lt-1))/po(ig,lt-1))
     442
     443          endif
     444
     445          else
    308446
    309447         zbuoyjam(ig,l)=1.*RG*(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- &
    … …  
    311449    &          ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
    312450
    313 !         zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
    314 !    &          po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
    315 !    &          po(ig,lt-1))/po(ig,lt-1))
     451          endif
    316452         
    317           else
    318           lt=lt+1
    319           endif
    320           enddo
    321453
    322454!          zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    323455
     456!          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
     457!    &     afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
     458
     459!          entrbis=entr_star(ig,l)
     460
     461          if (iflag_thermals_ed.lt.6) then
     462          fact_epsilon=0.0002/(zalpha+0.1)
     463          endif
     464
     465          detr_star(ig,l)=f_star(ig,l)*zdz             &
     466    &     *MAX(1.e-4, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m          &
     467    &     + 0.012*(zdqt(ig,l)/zw2m)**0.5)
     468
     469          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     470
    324471          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
    325     &     afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
    326 
    327           entrbis=entr_star(ig,l)
    328 
    329 
    330           detr_star(ig,l)=f_star(ig,l)*zdz                        &
    331     &     *MAX(1.e-4, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m          &
    332     &     + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
    333 
    334 
    335 !           zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    336 !
    337 !           entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha*  &     
    338 !     &     afact*zbuoy(ig,l)/zw2m &
    339 !     &     - 1.*fact_epsilon)
     472    &     afact*zbuoy(ig,l)/zw2m - fact_epsilon)
     473
     474!          entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha*  &     
     475!    &     afact*zbuoy(ig,l)/zw2m &
     476!    &     - 1.*fact_epsilon)
    340477
    341478         
    … …  
    350487!        endif
    351488
    352 !print*,'alim0',l,lalim(ig),alim_star(ig,l),entrbis,f_star(ig,l)
     489!       print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
    353490! Calcul du flux montant normalise
    354491      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
    … …  
    393530           zdzbis=zlev(ig,l+1)-zlev(ig,l-1)
    394531           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
    395 
     532           fact_epsilon=0.002
    396533            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
    397534            zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
    … …  
    402539    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
    403540    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
     541
     542           if (iflag_thermals_ed.lt.6) then
     543           zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l)
     544!           fact_epsilon=0.0005/(zalpha+0.025)**0.5
     545!           fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
     546           fact_epsilon=0.0002/(zalpha+0.1)**1
     547            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
     548            zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
     549            zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
     550            zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
     551
     552            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     553    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     554    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
     555
     556           endif
     557
     558
    404559      endif
    405560   enddo
    … …  
    425580     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
    426581           zw2(ig,l+1)=0.
     582!+CR:04/05/12:correction calcul linter pour calcul de zmax continu
     583        elseif (f_star(ig,l+1).lt.0.) then
     584           linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l))  &
     585     &               -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
     586           zw2(ig,l+1)=0.
     587!fin CR:04/05/12
    427588        endif
    428589
    … …  
    462623        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
    463624
     625#undef wrgrads_thermcell
     626#ifdef wrgrads_thermcell
     627         call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta      ','esta      ')
     628         call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta      ','dsta      ')
     629         call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy      ','buoy      ')
     630         call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt      ','dqt      ')
     631         call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est     ','w_est     ')
     632         call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2     ','w_es2     ')
     633         call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A      ','zw2A      ')
     634#endif
     635
     636
    464637     return
    465638     end
     639
     640
     641
     642
     643
     644
     645
     646
     647
     648
     649
     650
     651
     652
     653
     654
     655
     656
     657
     658
     659
     660
     661
     662
     663
     664
     665
     666
     667
     668
     669
     670
     671
    466672
    467673
    … …  
    536742      REAL zqsatth(ngrid,klev)
    537743      REAL zta_est(ngrid,klev)
     744      REAL zbuoyjam(ngrid,klev)
    538745      REAL ztemp(ngrid),zqsat(ngrid)
    539746      REAL zdw2
    … …  
    572779
    573780! Initialisations des variables reeles
    574 if (1==0) then
     781if (1==1) then
    575782      ztva(:,:)=ztv(:,:)
    576783      ztva_est(:,:)=ztva(:,:)
    … …  
    598805      zw2(:,:)=0.
    599806      zbuoy(:,:)=0.
     807      zbuoyjam(:,:)=0.
    600808      gamma(:,:)=0.
    601809      zeps(:,:)=0.
    … …  
    8621070        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
    8631071
     1072#undef wrgrads_thermcell
     1073#ifdef wrgrads_thermcell
     1074         call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta      ','esta      ')
     1075         call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta      ','dsta      ')
     1076         call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy      ','buoy      ')
     1077         call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt      ','dqt      ')
     1078         call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est     ','w_est     ')
     1079         call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2     ','w_es2     ')
     1080         call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A      ','zw2A      ')
     1081#endif
     1082
     1083
    8641084     return
    8651085     end
Note: See TracChangeset for help on using the changeset viewer.