Ignore:
Timestamp:
Jun 19, 2008, 12:24:22 PM (16 years ago)
Author:
lmdzadmin
Message:

Version thermique FH/CRio
Ajout tests cas physiques non pris en comptes et ajout/enleve prints
Nouvelle routine thermcell_flux2.F90
IM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/trunk/libf/phylmd/thermcell_plume.F90

    r938 r972  
    1       SUBROUTINE thermcell_plume(ngrid,klev,ztv,zthl,po,zl,rhobarz,  &
     1      SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
    22     &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,  &
    33     &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
    … …  
    1515#include "iniprint.h"
    1616
     17      INTEGER itap
     18
    1719      INTEGER ngrid,klev
     20      REAL ptimestep
    1821      REAL ztv(ngrid,klev)
    1922      REAL zthl(ngrid,klev)
    … …  
    3235      INTEGER lalim(ngrid)
    3336      integer lev_out                           ! niveau pour les print
     37      real zcon2(ngrid)
    3438
    3539      REAL ztva(ngrid,klev)
    3640      REAL ztla(ngrid,klev)
    3741      REAL zqla(ngrid,klev)
     42      REAL zqla0(ngrid,klev)
    3843      REAL zqta(ngrid,klev)
    3944      REAL zha(ngrid,klev)
    4045
    4146      REAL detr_star(ngrid,klev)
     47      REAL coefc
     48      REAL detr_stara(ngrid,klev)
     49      REAL detr_starb(ngrid,klev)
     50      REAL detr_starc(ngrid,klev)
     51      REAL detr_star0(ngrid,klev)
     52      REAL detr_star1(ngrid,klev)
     53      REAL detr_star2(ngrid,klev)
     54
    4255      REAL entr_star(ngrid,klev)
     56      REAL entr_star1(ngrid,klev)
     57      REAL entr_star2(ngrid,klev)
    4358      REAL detr(ngrid,klev)
    4459      REAL entr(ngrid,klev)
    … …  
    5570      REAL linter(ngrid)
    5671      INTEGER lmix(ngrid)
     72      INTEGER lmix_bis(ngrid)
    5773      REAL    wmaxa(ngrid)
    5874
    … …  
    85101            zqla(ig,k)=0.
    86102            zqta(ig,k)=po(ig,k)
     103!
     104            ztva(ig,k) = ztla(ig,k)*zpspsk(ig,k)+RLvCp*zqla(ig,k)
     105            ztva(ig,k) = ztva(ig,k)/zpspsk(ig,k)
     106            zha(ig,k) = ztva(ig,k)
     107!
    87108         enddo
    88109      enddo
    … …  
    91112           detr_star(ig,k)=0.
    92113           entr_star(ig,k)=0.
     114
     115           detr_stara(ig,k)=0.
     116           detr_starb(ig,k)=0.
     117           detr_starc(ig,k)=0.
     118           detr_star0(ig,k)=0.
     119           zqla0(ig,k)=0.
     120           detr_star1(ig,k)=0.
     121           detr_star2(ig,k)=0.
     122           entr_star1(ig,k)=0.
     123           entr_star2(ig,k)=0.
     124
    93125           detr(ig,k)=0.
    94126           entr(ig,k)=0.
    … …  
    117149      do l=1,klev-1
    118150         do ig=1,ngrid
     151
     152
     153
     154! Calcul dans la premiere couche active du thermique (ce qu'on teste
     155! en disant que la couche est instable et que w2 en bas de la couche
     156! est nulle.
     157
    119158            if (ztv(ig,l).gt.ztv(ig,l+1)  &
    120159     &         .and.alim_star(ig,l).gt.1.e-10  &
    121160     &         .and.zw2(ig,l).lt.1e-10) then
     161
     162
     163! Le panache va prendre au debut les caracteristiques de l'air contenu
     164! dans cette couche.
    122165               ztla(ig,l)=zthl(ig,l)
    123166               zqta(ig,l)=po(ig,l)
    124167               zqla(ig,l)=zl(ig,l)
    125168               f_star(ig,l+1)=alim_star(ig,l)
     169
    126170               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
    127171     &                     *(zlev(ig,l+1)-zlev(ig,l))  &
    … …  
    129173               w_est(ig,l+1)=zw2(ig,l+1)
    130174!
     175
     176
    131177            else if ((zw2(ig,l).ge.1e-10).and.  &
    132178     &         (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
    … …  
    147193!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    148194
     195
     196
     197! Premier calcul de la vitesse verticale a partir de la temperature
     198! potentielle virtuelle
     199
     200! FH CESTQUOI CA ????
     201#define int1d2
     202!#undef int1d2
     203#ifdef int1d2
     204      if (l.ge.2) then
     205#else
    149206      if (l.gt.2) then
     207#endif
     208
     209      if (1.eq.1) then
    150210          w_est(ig,3)=zw2(ig,2)* &
    151211     &      ((f_star(ig,2))**2) &
    … …  
    153213     &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
    154214     &      *(zlev(ig,3)-zlev(ig,2))
     215       endif
    155216
    156217
    … …  
    160221!
    161222!test:estimation de ztva_new_est sans entrainement
     223
    162224               Tbef=ztla(ig,l-1)*zpspsk(ig,l)
    163225               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    … …  
    204266!
    205267!calcul du detrainement
     268!=======================
     269
     270
     271! Modifications du calcul du detrainement.
     272! Dans la version de la these de Catherine, on passe brusquement
     273! de la version seche a la version nuageuse pour le detrainement
     274! ce qui peut occasioner des oscillations.
     275! dans la nouvelle version, on commence par calculer un detrainement sec.
     276! Puis un autre en cas de nuages.
     277! Puis on combine les deux lineairement en fonction de la quantite d'eau.
     278
     279#define int1d3
     280!#undef int1d3
     281#define RIO_TH
     282#ifdef RIO_TH
     283!1. Cas non nuageux
     284! 1.1 on est sous le zmax_sec et w croit
    206285          if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
    207286     &       (zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
     287#ifdef int1d3
     288     &       (zqla_est(ig,l).lt.1.e-10)) then
     289#else
    208290     &       (zqla(ig,l-1).lt.1.e-10)) then
     291#endif
    209292             detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
    210293     &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
    211294     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
    212295     &       /(r_aspect*zmax_sec(ig)))
    213        if (prt_level.ge.20) print*,'coucou calcul detr 1'
     296             detr_stara(ig,l)=detr_star(ig,l)
     297
     298       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l',ig,l
     299
     300! 1.2 on est sous le zmax_sec et w decroit
    214301          else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
     302#ifdef int1d3
     303     &            (zqla_est(ig,l).lt.1.e-10)) then
     304#else
    215305     &            (zqla(ig,l-1).lt.1.e-10)) then
     306#endif
    216307             detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
    217308     &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
    … …  
    222313     &       *((zmax_sec(ig)-zlev(ig,l))/  &
    223314     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
    224         if (prt_level.ge.20) print*,'coucou calcul detr 2'
     315             detr_starb(ig,l)=detr_star(ig,l)
     316
     317        if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l',ig,l
     318
    225319          else
     320
     321! 1.3 dans les autres cas
    226322             detr_star(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
    227323     &                      *(zlev(ig,l+1)-zlev(ig,l))
    228         if (prt_level.ge.20) print*,'coucou calcul detr 3'
     324             detr_starc(ig,l)=detr_star(ig,l)
     325
     326        if (prt_level.ge.20) print*,'coucou calcul detr 3 n: ig, l',ig, l
    229327             
    230328          endif
    231           detr_star(ig,l)=detr_star(ig,l)/f0(ig)
     329
     330#else
     331
     332! 1.1 on est sous le zmax_sec et w croit
     333          if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
     334     &       (zlev(ig,l+1).lt.zmax_sec(ig)) ) then
     335             detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
     336     &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
     337     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
     338     &       /(r_aspect*zmax_sec(ig)))
     339
     340       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
     341
     342! 1.2 on est sous le zmax_sec et w decroit
     343          else if ((zlev(ig,l+1).lt.zmax_sec(ig)) ) then
     344             detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
     345     &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
     346     &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
     347     &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
     348     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
     349     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
     350     &       *((zmax_sec(ig)-zlev(ig,l))/  &
     351     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
     352       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
     353
     354          else
     355             detr_star=0.
     356          endif
     357
     358! 1.3 dans les autres cas
     359          detr_starc(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
     360     &                      *(zlev(ig,l+1)-zlev(ig,l))
     361
     362          coefc=min(zqla(ig,l-1)/1.e-3,1.)
     363          if (zlev(ig,l+1).ge.zmax_sec(ig)) coefc=1.
     364          coefc=1.
     365! il semble qu'il soit important de baser le calcul sur
     366! zqla_est(ig,l-1) plutot que sur zqla_est(ig,l)
     367          detr_star(ig,l)=detr_starc(ig,l)*coefc+detr_star(ig,l)*(1.-coefc)
     368
     369        if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l', ig, l
     370
     371#endif
     372
     373
     374        if (prt_level.ge.20) print*,'coucou calcul detr 444: ig, l', ig, l
     375!IM 730508 beg
     376!        if(itap.GE.7200) THEN
     377!         print*,'th_plume ig,l,itap,zqla_est=',ig,l,itap,zqla_est(ig,l)
     378!        endif
     379!IM 730508 end
     380         
     381         zqla0(ig,l)=zqla_est(ig,l)
     382         detr_star0(ig,l)=detr_star(ig,l)
     383!IM 060508 beg
     384!         if(detr_star(ig,l).GT.1.) THEN
     385!          print*,'th_plumeBEF ig l detr_star detr_starc coefc',ig,l,detr_star(ig,l) &
     386!   &      ,detr_starc(ig,l),coefc
     387!         endif
     388!IM 060508 end
     389!IM 160508 beg
     390!IM 160508       IF (f0(ig).NE.0.) THEN
     391           detr_star(ig,l)=detr_star(ig,l)/f0(ig)
     392!IM 160508       ELSE IF(detr_star(ig,l).EQ.0.) THEN
     393!IM 160508        print*,'WARNING1  : th_plume f0=0, detr_star=0: ig, l, itap',ig,l,itap
     394!IM 160508       ELSE
     395!IM 160508        print*,'WARNING2  : th_plume f0=0, ig, l, itap, detr_star',ig,l,itap,detr_star(ig,l)
     396!IM 160508       ENDIF
     397!IM 160508 end
     398!IM 060508 beg
     399!        if(detr_star(ig,l).GT.1.) THEN
     400!         print*,'th_plumeAFT ig l detr_star f0 1/f0',ig,l,detr_star(ig,l),f0(ig), &
     401!   &     float(1)/f0(ig)
     402!        endif
     403!IM 060508 end
     404        if (prt_level.ge.20) print*,'coucou calcul detr 445: ig, l', ig, l
    232405!
    233406!calcul de entr_star
     407
     408! #undef test2
     409! #ifdef test2
     410! La version test2 destabilise beaucoup le modele.
     411! Il semble donc que ca aide d'avoir un entrainement important sous
     412! le nuage.
     413!         if (zqla_est(ig,l-1).ge.1.e-10.and.l.gt.lalim(ig)) then
     414!          entr_star(ig,l)=0.4*detr_star(ig,l)
     415!         else
     416!          entr_star(ig,l)=0.
     417!         endif
     418! #else
    234419!
    235420! Deplacement du calcul de entr_star pour eviter d'avoir aussi
    … …  
    237422! Redeplacer suite a la transformation du cas detr>f
    238423! FH
     424
     425        if (prt_level.ge.20) print*,'coucou calcul detr 446: ig, l', ig, l
     426#define int1d
     427!FH 070508 #define int1d4
     428!#undef int1d4
     429! L'option int1d4 correspond au choix dans le cas ou le detrainement
     430! devient trop grand.
     431
     432#ifdef int1d
     433
     434#ifdef int1d4
     435#else
     436       detr_star(ig,l)=min(detr_star(ig,l),f_star(ig,l))
     437!FH 070508 plus
     438       detr_star(ig,l)=min(detr_star(ig,l),1.)
     439#endif
     440
     441       entr_star(ig,l)=max(0.4*detr_star(ig,l)-alim_star(ig,l),0.)
     442
     443        if (prt_level.ge.20) print*,'coucou calcul detr 447: ig, l', ig, l
     444#ifdef int1d4
     445! Si le detrainement excede le flux en bas + l'entrainement, le thermique
     446! doit disparaitre.
     447       if (detr_star(ig,l)>f_star(ig,l)+entr_star(ig,l)) then
     448          detr_star(ig,l)=f_star(ig,l)+entr_star(ig,l)
     449          f_star(ig,l+1)=0.
     450          linter(ig)=l+1
     451          zw2(ig,l+1)=-1.e-10
     452       endif
     453#endif
     454
     455
     456#else
     457
     458        if (prt_level.ge.20) print*,'coucou calcul detr 448: ig, l', ig, l
    239459        if(l.gt.lalim(ig)) then
    240460         entr_star(ig,l)=0.4*detr_star(ig,l)
    … …  
    249469! d* non nul) on a une discontinuité de e* ou d* en haut de la couche
    250470! d'alimentation, ce qui n'est pas forcement heureux.
     471
     472        if (prt_level.ge.20) print*,'coucou calcul detr 449: ig, l', ig, l
     473#undef pre_int1c
     474#ifdef pre_int1c
    251475         entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.)
    252476         detr_star(ig,l)=entr_star(ig,l)
     477#else
     478         entr_star(ig,l)=0.
     479#endif
     480
    253481        endif
    254482
    255 !
     483#endif
     484
     485        if (prt_level.ge.20) print*,'coucou calcul detr 440: ig, l', ig, l
     486        entr_star1(ig,l)=entr_star(ig,l)
     487        detr_star1(ig,l)=detr_star(ig,l)
     488!
     489
     490#ifdef int1d
     491#else
    256492        if (detr_star(ig,l).gt.f_star(ig,l)) then
    257493
    … …  
    261497
    262498           detr_star(ig,l)=f_star(ig,l)
     499#ifdef pre_int1c
    263500           if (l.gt.lalim(ig)+1) then
    264501               entr_star(ig,l)=0.
    … …  
    269506               linter(ig)=l+1
    270507            else
    271                entr_star(ig,l)=detr_star(ig,l)
     508               entr_star(ig,l)=0.4*detr_star(ig,l)
    272509            endif
     510#else
     511           entr_star(ig,l)=0.4*detr_star(ig,l)
     512#endif
    273513        endif
    274 
    275       else
     514#endif
     515
     516      else !l > 2
    276517         detr_star(ig,l)=0.
    277518         entr_star(ig,l)=0.
    278519      endif
     520
     521        entr_star2(ig,l)=entr_star(ig,l)
     522        detr_star2(ig,l)=detr_star(ig,l)
     523        if (prt_level.ge.20) print*,'coucou calcul detr 450: ig, l', ig, l
    279524
    280525!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    292537
    293538!test sur le signe de f_star
     539        if (prt_level.ge.20) print*,'coucou calcul detr 451: ig, l', ig, l
    294540       if (f_star(ig,l+1).gt.1.e-10) then
    295541!----------------------------------------------------------------------------
    … …  
    336582              endif
    337583!   
     584        if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
    338585! on ecrit de maniere conservative (sat ou non)
    339586!          T = Tl +Lv/Cp ql
    … …  
    356603            endif
    357604        endif
     605        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
    358606!
    359607!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
    … …  
    383631      enddo
    384632
     633        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
     634#ifdef troisD
     635      call wrgradsfi(1,klev,zqla_est,'ql_es     ','ql_es     ')
     636      call wrgradsfi(1,klev,entr_star1,'e_st1     ','e_st1     ')
     637      call wrgradsfi(1,klev,entr_star2,'e_st2     ','e_st2     ')
     638      call wrgradsfi(1,klev,detr_stara,'d_sta     ','d_sta     ')
     639      call wrgradsfi(1,klev,detr_starb,'d_stb     ','d_stb     ')
     640      call wrgradsfi(1,klev,detr_starc,'d_stc     ','d_stc     ')
     641      call wrgradsfi(1,klev,zqla0     ,'zqla0     ','zqla0     ')
     642      call wrgradsfi(1,klev,detr_star0,'d_st0     ','d_st0     ')
     643      call wrgradsfi(1,klev,detr_star1,'d_st1     ','d_st1     ')
     644      call wrgradsfi(1,klev,detr_star2,'d_st2     ','d_st2     ')
     645#endif
     646
     647!     print*,'thermcell_plume OK'
    385648
    386649      return
Note: See TracChangeset for help on using the changeset viewer.