Changeset 1338 for LMDZ4/branches


Ignore:
Timestamp:
Apr 6, 2010, 2:49:00 PM (15 years ago)
Author:
idelkadi
Message:

Nouvelle version des thermiques

Location:
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_main.F90

    r1311 r1338  
    130130
    131131      real wmax(klon)
     132      real wmax_tmp(klon)
    132133      real wmax_sec(klon)
    133134      real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
    … …  
    191192
    192193
    193 ! #define wrgrads_thermcell
    194 #undef wrgrads_thermcell
     194#define wrgrads_thermcell
    195195#ifdef wrgrads_thermcell
    196196! Initialisation des sorties grads pour les thermiques.
    … …  
    208208
    209209      fm=0. ; entr=0. ; detr=0.
     210
     211      print*,'THERMCELL MAIN OPT7'
    210212
    211213      icount=icount+1
    … …  
    395397
    396398      elseif (iflag_thermals_ed<=19) then
    397 ! Version d'Arnaud Jam
    398 !         print*,'THERM RIO et al 2010'
     399!        print*,'THERM RIO et al 2010, version d Arnaud'
    399400         CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
    400401     &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
    … …  
    426427      CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,  &
    427428     &           zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
     429! Attention, w2 est transforme en sa racine carree dans cette routine
     430! Le probleme vient du fait que linter et lmix sont souvent égaux à 1.
     431      wmax_tmp=0.
     432      do  l=1,nlay
     433         wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
     434      enddo
     435!     print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax
    428436
    429437
    … …  
    665673         enddo
    666674      enddo
    667 !calcul de ale_bl et alp_bl
    668 !pour le calcul d'une valeur integree entre la surface et lmax
    669       do ig=1,ngrid
    670       alp_int(ig)=0.
    671       ale_int(ig)=0.
    672       n_int(ig)=0
    673       enddo
    674 !
    675       do l=1,nlay
    676       do ig=1,ngrid
    677        if(l.LE.lmax(ig)) THEN
    678         alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)
    679         ale_int(ig)=ale_int(ig)+0.5*zw2(ig,l)**2
    680         n_int(ig)=n_int(ig)+1
    681        endif
    682       enddo
    683       enddo
     675!
    684676!      print*,'avant calcul ale et alp'
    685677!calcul de ALE et ALP pour la convection
     678      Alp_bl(:)=0.
     679      Ale_bl(:)=0.
     680!          print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
     681      do l=1,nlay
    686682      do ig=1,ngrid
    687 !      Alp_bl(ig)=0.5*rhobarz(ig,lmix_bis(ig))*wth3(ig,lmix(ig))
    688 !          Alp_bl(ig)=0.5*rhobarz(ig,nivcon(ig))*wth3(ig,nivcon(ig))
    689 !      Alp_bl(ig)=0.5*rhobarz(ig,lmix(ig))*wth3(ig,lmix(ig))
    690 !     &           *0.1
    691 !valeur integree de alp_bl * 0.5:
    692        if (n_int(ig).gt.0) then
    693        Alp_bl(ig)=0.5*alp_int(ig)/n_int(ig)
    694 !       if (Alp_bl(ig).lt.0.) then
    695 !       Alp_bl(ig)=0.
    696        endif
    697 !       endif
    698 !         write(18,*),'rhobarz,wth3,Alp',rhobarz(ig,nivcon(ig)),
    699 !     s               wth3(ig,nivcon(ig)),Alp_bl(ig)
    700 !            write(18,*),'ALP_BL',Alp_bl(ig),lmix(ig)
    701 !      Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2
    702 !      if (nivcon(ig).eq.1) then
    703 !       Ale_bl(ig)=0.
    704 !       else
    705 !valeur max de ale_bl:
    706        Ale_bl(ig)=0.5*zw2(ig,lmix(ig))**2
    707 !     & /2.
    708 !     & *0.1
    709 !        Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2
    710 !       if (n_int(ig).gt.0) then
    711 !       Ale_bl(ig)=ale_int(ig)/n_int(ig)
    712 !        Ale_bl(ig)=4.
    713 !       endif
    714 !       endif
    715 !            Ale_bl(ig)=0.5*wth2(ig,lmix_bis(ig))
    716 !          Ale_bl(ig)=wth2(ig,nivcon(ig))
    717 !          write(19,*),'wth2,ALE_BL',wth2(ig,nivcon(ig)),Ale_bl(ig)
    718       enddo
     683           Alp_bl(ig)=max(Alp_bl(ig),0.5*rhobarz(ig,l)*wth3(ig,l) )
     684           Ale_bl(ig)=max(Ale_bl(ig),0.5*zw2(ig,l)**2)
     685!          print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
     686      enddo
     687      enddo
     688
     689!     print*,'AAAAAAA ',Alp_bl,Ale_bl,lmix
     690
     691
     692! TEST. IL FAUT REECRIRE LES ALE et ALP
     693!     Ale_bl(:)=0.5*wmax(:)*wmax(:)
     694!     Alp_bl(:)=0.1*wmax(:)*wmax(:)*wmax(:)
     695
    719696!test:calcul de la ponderation des couches pour KE
    720697!initialisations

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90

    r1311 r1338  
    6464      REAL zqsatth(ngrid,klev)
    6565      REAL zta_est(ngrid,klev)
     66      REAL ztemp(ngrid),zqsat(ngrid)
    6667      REAL zdw2
    6768      REAL zw2modif
    … …  
    7778      real zdz,zfact,zbuoy,zalpha
    7879      real zcor,zdelta,zcvm5,qlbef
    79       real Tbef,qsatbef
    80       real dqsat_dT,DT,num,denom
    8180      REAL REPS,RLvCp,DDT0
    8281      PARAMETER (DDT0=.01)
    … …  
    116115!      write(lunout,*)'THERM 31H '
    117116
     117      print*,'THERMCELL_PLUME OPTIMISE V0 '
    118118! Initialisations des variables reeles
    119119if (1==0) then
    … …  
    239239! sans tenir compte du detrainement et de l'entrainement dans cette
    240240! couche
     241! C'est a dire qu'on suppose
     242! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
    241243! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
    242244! avant) a l'alimentation pour avoir un calcul plus propre
    243245!---------------------------------------------------------------------------
    244246
    245      call thermcell_condens(ngrid,active, &
    246 &          zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
     247   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
     248   call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat)
     249
    247250
    248251    do ig=1,ngrid
    249252        if(active(ig)) then
     253        zqla_est(ig,l) = max(0.,zqta(ig,l-1)-zqsat(ig))
    250254        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
    251255        zta_est(ig,l)=ztva_est(ig,l)
    … …  
    366370    enddo
    367371
    368    call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
    369 
     372   ztemp(:)=zpspsk(:,l)*ztla(:,l)
     373   call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
    370374
    371375   do ig=1,ngrid
    … …  
    373377        if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
    374378! on ecrit de maniere conservative (sat ou non)
     379           zqla(ig,l) = max(0.,zqta(ig,l)-zqsatth(ig,l))
    375380!          T = Tl +Lv/Cp ql
    376381           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
    … …  
    382387
    383388!on ecrit zqsat
    384            zqsatth(ig,l)=qsatbef 
    385389
    386390!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    554558      REAL zqsatth(ngrid,klev)
    555559      REAL zta_est(ngrid,klev)
     560      REAL ztemp(ngrid),zqsat(ngrid)
    556561      REAL zdw2
    557562      REAL zw2modif
    … …  
    570575      real zcor,zdelta,zcvm5,qlbef,zdz2
    571576      real betalpha,zbetalpha
    572       real Tbef,qsatbef,b1,eps, afact
    573       real dqsat_dT,DT,num,denom,m
     577      real eps, afact
    574578      REAL REPS,RLvCp,DDT0
    575579      PARAMETER (DDT0=.01)
    … …  
    589593
    590594!      print*,'THERM 31B'
     595      print*,'THERMCELL_PLUME OPTIMISE V1 CCC '
    591596
    592597! Initialisations des variables reeles
    … …  
    624629      linter(:)=1.
    625630!     linter(:)=1.
    626       b1=2.
    627631! Initialisation des variables entieres
    628632      lmix(:)=1
    … …  
    631635      lalim(:)=1
    632636
    633 !      print*,'THERMCELL PLUME ARNAUD DEDANS'
     637   print*,'THERMCELL PLUME QSAT2 NDDDDN'
    634638
    635639!-------------------------------------------------------------------------
    … …  
    708712! sans tenir compte du detrainement et de l'entrainement dans cette
    709713! couche
     714! C'est a dire qu'on suppose
     715! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
    710716! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
    711717! avant) a l'alimentation pour avoir un calcul plus propre
    712718!---------------------------------------------------------------------------
    713719
    714      call thermcell_condens(ngrid,active, &
    715 &          zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
    716 
     720   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
     721   call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
    717722
    718723    do ig=1,ngrid
    719724!       print*,'active',active(ig),ig,l
    720725        if(active(ig)) then
     726        zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
    721727        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
    722728        zta_est(ig,l)=ztva_est(ig,l)
    … …  
    747753!-------------------------------------------------
    748754
     755     print*,'THERM V1 SANS DQ'
    749756     do ig=1,ngrid
    750757        if (active(ig)) then
    … …  
    758765          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
    759766          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
     767          zdqt(ig,l)=0.
    760768
    761769         
    … …  
    800808    enddo
    801809
    802    call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
    803 
     810   ztemp(:)=zpspsk(:,l)*ztla(:,l)
     811   call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
    804812
    805813   do ig=1,ngrid
    … …  
    808816! on ecrit de maniere conservative (sat ou non)
    809817!          T = Tl +Lv/Cp ql
     818           zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
    810819           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
    811820           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
    … …  
    814823           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
    815824     &              -zqla(ig,l))-zqla(ig,l))
    816 
    817 !on ecrit zqsat
    818            zqsatth(ig,l)=qsatbef 
    819 
    820825           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
    821826           zdz=zlev(ig,l+1)-zlev(ig,l)

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_qsat.F90

    r1336 r1338  
    1 subroutine thermcell_qsat(klon,mask,pplev,ptemp,pqt,pqsat)
     1subroutine thermcell_qsat(klon,active,pplev,ztemp,zqta,zqsat)
    22implicit none
    33
    … …  
    55#include "YOETHF.h"
    66#include "FCTTRE.h"
     7
    78
    89!====================================================================
    … …  
    1213! Arguments
    1314INTEGER klon
    14 REAL pplev(klon)
    15 REAL ptemp(klon),pqt(klon),pqsat(klon)
    16 LOGICAL mask(klon)
     15REAL zpspsk(klon),pplev(klon)
     16REAL ztemp(klon),zqta(klon),zqsat(klon)
     17LOGICAL active(klon)
    1718
    1819! Variables locales
    1920INTEGER ig,iter
    2021REAL Tbef(klon),DT(klon)
    21 REAL tdelta,zcor,qlbef,zdelta,zcvm5,dqsat,num,denom,dqsat_dT
     22REAL tdelta,qsatbef,zcor,qlbef,zdelta,zcvm5,dqsat,num,denom,dqsat_dT
    2223logical Zsat
    2324REAL RLvCp
    … …  
    3132RLvCp = RLVTT/RCPD
    3233tout_converge=.false.
    33 afaire(:)=mask(:)
     34afaire(:)=.false.
    3435DT(:)=0.
    3536
    3637
    3738!====================================================================
    38 ! Routine a vectoriser en copiant mask dans converge et en mettant
     39! Routine a vectoriser en copiant active dans converge et en mettant
    3940! la boucle sur les iterations a l'exterieur est en mettant
    4041! converge= false des que la convergence est atteinte.
    41 ! Pr Tprec=Tl calcul de qsat
    42 ! Si qsat>qT T=Tl, q=qT
    43 ! Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt)
    44 ! On cherche DDT < DDT0
    4542!====================================================================
    4643
    4744do ig=1,klon
    48    if (mask(ig)) then
    49                Tbef(ig)=ptemp(ig)
     45   if (active(ig)) then
     46               Tbef(ig)=ztemp(ig)
    5047               zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
    51                pqsat(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
    52                pqsat(ig)=MIN(0.5,pqsat(ig))
    53                zcor=1./(1.-retv*pqsat(ig))
    54                pqsat(ig)=pqsat(ig)*zcor
    55                qlbef=max(0.,pqt(ig)-pqsat(ig))
    56                afaire(ig)=qlbef>1.e-10
     48               qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
     49               qsatbef=MIN(0.5,qsatbef)
     50               zcor=1./(1.-retv*qsatbef)
     51               qsatbef=qsatbef*zcor
     52               qlbef=max(0.,zqta(ig)-qsatbef)
    5753               DT(ig) = 0.5*RLvCp*qlbef
     54               zqsat(ig)=qsatbef
    5855     endif
    5956enddo
    6057
     58! Traitement du cas ou il y a condensation mais faible
     59! On ne condense pas mais on dit que le qsat est le qta
     60do ig=1,klon
     61   if (active(ig)) then
     62      if (0.<abs(DT(ig)).and.abs(DT(ig))<=DDT0) then
     63         zqsat(ig)=zqta(ig)
     64      endif
     65   endif
     66enddo
    6167
    6268do iter=1,10
    63     afaire(:)=(abs(DT(:))>DDT0).and.afaire(:)
     69    afaire(:)=abs(DT(:)).gt.DDT0
    6470    do ig=1,klon
    6571               if (afaire(ig)) then
    6672                 Tbef(ig)=Tbef(ig)+DT(ig)
    6773                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
    68                  pqsat(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
    69                  pqsat(ig)=MIN(0.5,pqsat(ig))
    70                  zcor=1./(1.-retv*pqsat(ig))
    71                  pqsat(ig)=pqsat(ig)*zcor
    72                  qlbef=pqt(ig)-pqsat(ig)
     74                 qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
     75                 qsatbef=MIN(0.5,qsatbef)
     76                 zcor=1./(1.-retv*qsatbef)
     77                 qsatbef=qsatbef*zcor
     78                 qlbef=zqta(ig)-qsatbef
    7379                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
    7480                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
    75                  zcor=1./(1.-retv*pqsat(ig))
    76                  dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,pqsat(ig),zcor)
    77                  num=-Tbef(ig)+ptemp(ig)+RLvCp*qlbef
     81                 zcor=1./(1.-retv*qsatbef)
     82                 dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef,zcor)
     83                 num=-Tbef(ig)+ztemp(ig)+RLvCp*qlbef
    7884                 denom=1.+RLvCp*dqsat_dT
     85                 zqsat(ig) = qsatbef
    7986                 DT(ig)=num/denom
    8087               endif
Note: See TracChangeset for help on using the changeset viewer.