Changeset 1330 for LMDZ4


Ignore:
Timestamp:
Mar 24, 2010, 1:41:35 PM (15 years ago)
Author:
idelkadi
Message:

Optimisation des thermiques

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

Legend:

Unmodified
Added
Removed

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

    r1299 r1330  
    225225!Initialisation
    226226!
    227 !    IF (1.eq.0) THEN
    228 !     do ig=1,klon     
    229 !FH/IM 130308     if ((debut).or.((.not.debut).and.(f0(ig).lt.1.e-10))) then
    230 !     if ((.not.debut).and.(f0(ig).lt.1.e-10)) then
    231 !           f0(ig)=1.e-5
    232 !           zmax0(ig)=40.
    233 !v1d        therm=.false.
    234 !     endif
    235 !     enddo
    236 !    ENDIF !(1.eq.0) THEN
    237227     if (prt_level.ge.10)write(lunout,*)                                &
    238228    &     'WARNING thermcell_main f0=max(f0,1.e-2)'
    239229     do ig=1,klon
    240       if (prt_level.ge.20) then
    241        print*,'th_main ig f0',ig,f0(ig)
    242       endif
    243230         f0(ig)=max(f0(ig),1.e-2)
    244 !IMmarche pas ?!       if (f0(ig)<1.e-2) f0(ig)=1.e-2
    245231     enddo
    246232
    … …  
    637623            zf2=zf/(1.-zf)
    638624!
    639       if (prt_level.ge.10) print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
    640 !
    641       if (prt_level.ge.10) print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
    642625            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))**2
    643626            if(zw2(ig,l).gt.1.e-10) then
    … …  
    649632            wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))  &
    650633     &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
    651       if (prt_level.ge.10) print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
    652634            q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
    653635!test: on calcul q2/po=ratqsc
    … …  
    655637         enddo
    656638      enddo
    657 !calcul de ale_bl et alp_bl
    658 !pour le calcul d'une valeur intégrée entre la surface et lmax
     639
     640      if (prt_level.ge.10) then
     641          print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
     642          ig=igout
     643          do l=1,nlay
     644             print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
     645          enddo
     646          do l=1,nlay
     647             print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
     648          enddo
     649      endif
     650
    659651      do ig=1,ngrid
    660652      alp_int(ig)=0.

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

    r970 r1330  
    11      SUBROUTINE thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
    2      &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
     2     &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,pqsat,lev_out)
    33
    44!--------------------------------------------------------------
    … …  
    3131      REAL zu(ngrid,nlay)
    3232      REAL zv(ngrid,nlay)
    33       REAL zqsat(ngrid,nlay)
     33      REAL pqsat(ngrid,nlay)
    3434
    35       INTEGER ig,l,ll
     35      INTEGER ig,ll
    3636
    37       real zcor,zdelta,zcvm5,qlbef
    38       real Tbef,qsatbef
    39       real dqsat_dT,DT,num,denom
    40       REAL RLvCp,DDT0
    41       PARAMETER (DDT0=.01)
    42       LOGICAL Zsat
     37      real dqsat_dT
     38      real RLvCp
    4339
    44       Zsat=.false.
    45       RLvCp = RLVTT/RCPD
     40logical mask(ngrid,nlay)
    4641
    47 !
    48 ! Pr Tprec=Tl calcul de qsat
    49 ! Si qsat>qT T=Tl, q=qT
    50 ! Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt)
    51 ! On cherche DDT < DDT0
     42
     43!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     44! Initialisations :
     45!------------------
     46
     47print*,'THERMCELL ENV OPTIMISE '
     48mask(:,:)=.true.
     49RLvCp = RLVTT/RCPD
     50
    5251!
    5352! calcul des caracteristiques de l environnement
    … …  
    5756            zl(ig,ll)=0.
    5857            zh(ig,ll)=pt(ig,ll)
    59             zqsat(ig,ll)=0.
    6058         EndDO
    6159       EndDO
    6260!
    6361!
    64 !recherche de saturation dans l environnement
    65        DO ll=1,nlay
    66 ! les points insatures sont definitifs
    67          DO ig=1,ngrid
    68             Tbef=pt(ig,ll)
    69             zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    70             qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,ll)
    71             qsatbef=MIN(0.5,qsatbef)
    72             zcor=1./(1.-retv*qsatbef)
    73             qsatbef=qsatbef*zcor
    74             Zsat = (max(0.,po(ig,ll)-qsatbef) .gt. 1.e-10)
    75             if (Zsat) then
    76             qlbef=max(0.,po(ig,ll)-qsatbef)
    77 ! si sature: ql est surestime, d'ou la sous-relax
    78             DT = 0.5*RLvCp*qlbef
    79 ! on pourra enchainer 2 ou 3 calculs sans Do while
    80             do while (abs(DT).gt.DDT0)
    81 ! il faut verifier si c,a conserve quand on repasse en insature ...
    82               Tbef=Tbef+DT
    83               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    84               qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,ll)
    85               qsatbef=MIN(0.5,qsatbef)
    86               zcor=1./(1.-retv*qsatbef)
    87               qsatbef=qsatbef*zcor
    88 ! on veut le signe de qlbef
    89               qlbef=po(ig,ll)-qsatbef
    90 !          dqsat_dT
    91               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    92               zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
    93               zcor=1./(1.-retv*qsatbef)
    94               dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
    95               num=-Tbef+pt(ig,ll)+RLvCp*qlbef
    96               denom=1.+RLvCp*dqsat_dT
    97               if (denom.lt.1.e-10) then
    98                   print*,'pb denom'
    99               endif
    100               DT=num/denom
    101             enddo
    102 ! on ecrit de maniere conservative (sat ou non)
    103             zl(ig,ll) = max(0.,qlbef)
    104 !          T = Tl +Lv/Cp ql
    105             zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll)
    106             zo(ig,ll) = po(ig,ll)-zl(ig,ll)
    107            endif
    108 !on ecrit zqsat
    109             zqsat(ig,ll)=qsatbef     
    110          EndDO
    111        EndDO
     62! Condensation :
     63!---------------
     64! Calcul de l'humidite a saturation et de la condensation
     65
     66call thermcell_qsat(ngrid*nlay,mask,pplev,pt,po,pqsat)
     67DO ll=1,nlay
     68   DO ig=1,ngrid
     69      zl(ig,ll) = max(0.,po(ig,ll)-pqsat(ig,ll))
     70      zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll)         !   T = Tl + Lv/Cp ql
     71      zo(ig,ll) = po(ig,ll)-zl(ig,ll)
     72   ENDDO
     73ENDDO
    11274!
    11375!
    11476!-----------------------------------------------------------------------
    115 !   incrementation eventuelle de tendances precedentes:
    116 !   ---------------------------------------------------
    11777
    11878      if (prt_level.ge.1) print*,'0 OK convect8'
    11979
    120       DO 1010 l=1,nlay
    121          DO 1015 ig=1,ngrid
    122              zpspsk(ig,l)=(pplay(ig,l)/100000.)**RKAPPA
    123              zu(ig,l)=pu(ig,l)
    124              zv(ig,l)=pv(ig,l)
     80      DO ll=1,nlay
     81         DO ig=1,ngrid
     82             zpspsk(ig,ll)=(pplay(ig,ll)/100000.)**RKAPPA
     83             zu(ig,ll)=pu(ig,ll)
     84             zv(ig,ll)=pv(ig,ll)
    12585!attention zh est maintenant le profil de T et plus le profil de theta !
     86! Quelle horreur ! A eviter.
    12687!
    12788!   T-> Theta
    128             ztv(ig,l)=zh(ig,l)/zpspsk(ig,l)
     89            ztv(ig,ll)=zh(ig,ll)/zpspsk(ig,ll)
    12990!Theta_v
    130             ztv(ig,l)=ztv(ig,l)*(1.+RETV*(zo(ig,l))  &
    131      &           -zl(ig,l))
     91            ztv(ig,ll)=ztv(ig,ll)*(1.+RETV*(zo(ig,ll))-zl(ig,ll))
    13292!Thetal
    133             zthl(ig,l)=pt(ig,l)/zpspsk(ig,l)
     93            zthl(ig,ll)=pt(ig,ll)/zpspsk(ig,ll)
    13494!           
    135 1015     CONTINUE
    136 1010  CONTINUE
     95         ENDDO
     96      ENDDO
    13797 
    13898      RETURN

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

    r1299 r1330  
    4141      REAL zfm
    4242
    43       integer igout
     43      integer igout,lout
    4444      integer lev_out
    4545      integer lunout1
    … …  
    4949      REAL fomass_max,alphamax
    5050      save fomass_max,alphamax
     51
     52      logical check_debug,labort_gcm
    5153
    5254      character (len=20) :: modname='thermcell_flux2'
    … …  
    8486! Verification de la nullite des entrainement et detrainement au dessus
    8587! de lmax(ig)
    86 !-------------------------------------------------------------------------
    87 
     88! Active uniquement si check_debug=.true. ou prt_level>=10
     89!-------------------------------------------------------------------------
     90
     91      check_debug=.false..or.prt_level>=10
     92
     93      if (check_debug) then
    8894      do l=1,klev
    8995         do ig=1,ngrid
    … …  
    102108                    print*,'detr_star(ig,l)',detr_star(ig,l)
    103109                    abort_message = ''
     110                    labort_gcm=.true.
    104111                    CALL abort_gcm (modname,abort_message,1)
    105112               endif
    … …  
    107114         enddo
    108115      enddo
     116      endif
    109117
    110118!-------------------------------------------------------------------------
    … …  
    259267
    260268!     do l=1,klev
     269
     270
     271
     272         labort_gcm=.false.
    261273         do ig=1,ngrid
    262274            if (entr(ig,l)<0.) then
    263                print*,'N1 ig,l,entr',ig,l,entr(ig,l)
    264                abort_message = 'entr negatif'
    265                CALL abort_gcm (modname,abort_message,1)
    266             endif
     275               labort_gcm=.true.
     276               igout=ig
     277               lout=l
     278            endif
     279         enddo
     280
     281         if (labort_gcm) then
     282            print*,'N1 ig,l,entr',igout,lout,entr(igout,lout)
     283            abort_message = 'entr negatif'
     284            CALL abort_gcm (modname,abort_message,1)
     285         endif
     286
     287         do ig=1,ngrid
    267288            if (detr(ig,l).gt.fm(ig,l)) then
    268289               ncorecfm6=ncorecfm6+1
    … …  
    287308               entr(ig,l)=0.
    288309            endif
    289 
     310         enddo
     311
     312         labort_gcm=.false.
     313         do ig=1,ngrid
    290314            if (entr(ig,l).lt.0.) then
    291                print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
    292                print*,'entr(ig,l)',entr(ig,l)
    293                print*,'fm(ig,l)',fm(ig,l)
    294                abort_message = 'probleme dans thermcell flux'
    295                CALL abort_gcm (modname,abort_message,1)
    296             endif
    297          enddo
     315               labort_gcm=.true.
     316               igout=ig
     317            endif
     318         enddo
     319         if (labort_gcm) then
     320            ig=igout
     321            print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
     322            print*,'entr(ig,l)',entr(ig,l)
     323            print*,'fm(ig,l)',fm(ig,l)
     324            abort_message = 'probleme dans thermcell flux'
     325            CALL abort_gcm (modname,abort_message,1)
     326         endif
     327
     328
    298329!     enddo
    299330      endif
    … …  
    313344               detr(ig,l)=detr(ig,l)+fm(ig,l+1)
    314345               fm(ig,l+1)=0.
    315 !              print*,'fm2<0',l+1,lmax(ig)
    316346               ncorecfm2=ncorecfm2+1
    317347            endif
     348         enddo
     349
     350         labort_gcm=.false.
     351         do ig=1,ngrid
    318352            if (detr(ig,l).lt.0.) then
     353               labort_gcm=.true.
     354               igout=ig
     355            endif
     356        enddo
     357        if (labort_gcm) then
     358               ig=igout
    319359               print*,'cas 2 : ig,l,lmax(ig)',ig,l,lmax(ig)
    320360               print*,'detr(ig,l)',detr(ig,l)
    … …  
    322362               abort_message = 'probleme dans thermcell flux'
    323363               CALL abort_gcm (modname,abort_message,1)
    324             endif
    325         enddo
     364        endif
    326365!    enddo
    327366
    … …  
    388427
    389428      if (1.eq.1) then
     429      labort_gcm=.false.
    390430      do l=1,klev-1
    391431         do ig=1,ngrid
    … …  
    408448                   else
    409449                      if(l.ge.lmax(ig).and.0.eq.1) then
     450                         igout=ig
     451                         lout=l
     452                         labort_gcm=.true.
     453                      endif
     454                      entr(ig,l+1)=entr(ig,l+1)-ddd
     455                      detr(ig,l)=0.
     456                      fm(ig,l+1)=fm(ig,l)+entr(ig,l)
     457                      detr(ig,l)=0.
     458                   endif
     459                endif
     460            endif
     461         enddo
     462      enddo
     463      if (labort_gcm) then
     464                         ig=igout
     465                         l=lout
    410466                         print*,'ig,l',ig,l
    411467                         print*,'eee0',eee0
    … …  
    424480                         abort_message = 'probleme dans thermcell_flux'
    425481                         CALL abort_gcm (modname,abort_message,1)
    426                       endif
    427                       entr(ig,l+1)=entr(ig,l+1)-ddd
    428                       detr(ig,l)=0.
    429                       fm(ig,l+1)=fm(ig,l)+entr(ig,l)
    430                       detr(ig,l)=0.
    431                    endif
    432                 endif
    433             endif
    434          enddo
    435       enddo
     482      endif
    436483      endif
    437484!                 
Note: See TracChangeset for help on using the changeset viewer.