Ignore:
Timestamp:
Oct 11, 2007, 3:43:42 PM (17 years ago)
Author:
Laurent Fairhead
Message:

Mise a jour de la physique avec thermiques avec la version de FH d'aout 2007
LF

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/trunk/libf/phytherm/thermcell_main.F90

    r814 r852  
     1!
     2! $Header$
     3!
    14      SUBROUTINE thermcell_main(ngrid,nlay,ptimestep  &
    25     &                  ,pplay,pplev,pphi,debut  &
    … …  
    5558!   ------
    5659
    57       integer,save :: igout=871
     60!      integer,save :: igout=4259
     61      integer,save :: igout=2928
    5862      integer,save :: lunout=6
    59       integer,save :: lev_out=0
     63      integer,save :: lev_out=10
    6064
    6165      INTEGER ig,k,l,ll
    … …  
    129133      real zlevinter(klon)
    130134      logical debut
     135       real seuil
    131136
    132137!
    … …  
    142147!   ---------------
    143148!
     149
     150      seuil=0.25
     151
    144152      if (lev_out.ge.1) print*,'thermcell_main V4'
    145153
    … …  
    228236
    229237!------------------------------------------------------------------
    230 !   Calcul de w2, carre de w a partir de la cape
    231 !
    232 !   Indicages:
    233 !   l'ascendance provenant du niveau k traverse l'interface l avec
    234 !   une vitesse wa(k,l).
    235 !
    236 !                       --------------------
    237 !
    238 !                       + + + + + + + + + +
    239 !
    240 !  wa(k,l)   ----       --------------------    l
    241 !             /\
    242 !            /||\       + + + + + + + + + +
    243 !             ||
    244 !             ||        --------------------
    245 !             ||
    246 !             ||        + + + + + + + + + +
    247 !             ||
    248 !             ||        --------------------
    249 !             ||__
    250 !             |___      + + + + + + + + + +     k
    251 !
    252 !                       --------------------
    253 !
    254 !
    255 !
     238!
     239!             /|\
     240!    --------  |  F_k+1 -------   
     241!                              ----> D_k
     242!             /|\              <---- E_k , A_k
     243!    --------  |  F_k ---------
     244!                              ----> D_k-1
     245!                              <---- E_k-1 , A_k-1
     246!
     247!
     248!
     249!
     250!
     251!    ---------------------------
     252!
     253!    ----- F_lmax+1=0 ----------         \
     254!            lmax     (zmax)              |
     255!    ---------------------------          |
     256!                                         |
     257!    ---------------------------          |
     258!                                         |
     259!    ---------------------------          |
     260!                                         |
     261!    ---------------------------          |
     262!                                         |
     263!    ---------------------------          |
     264!                                         |  E
     265!    ---------------------------          |  D
     266!                                         |
     267!    ---------------------------          |
     268!                                         |
     269!    ---------------------------  \       |
     270!            lalim                 |      |
     271!    ---------------------------   |      |
     272!                                  |      |
     273!    ---------------------------   |      |
     274!                                  | A    |
     275!    ---------------------------   |      |
     276!                                  |      |
     277!    ---------------------------   |      |
     278!    lmin  (=1 pour le moment)     |      |
     279!    ----- F_lmin=0 ------------  /      /
     280!
     281!    ---------------------------
     282!    //////////////////////////
     283!
     284!
     285!=============================================================================
     286!  Calculs initiaux ne faisant pas intervenir les changements de phase
     287!=============================================================================
     288
    256289!------------------------------------------------------------------
    257 !definition du profil d alimentation a partir de la flottabilite:
    258 !alim_star, alim_star_tot, lalim et lmin
     290!  1. alim_star est le profil vertical de l'alimentation à la base du
     291!     panache thermique, calculé à partir de la flotabilité de l'air sec
     292!  2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
    259293!------------------------------------------------------------------
    260294!
    … …  
    262296      CALL thermcell_init(ngrid,nlay,ztv,zlev,  &
    263297     &                    lalim,lmin,alim_star,alim_star_tot,lev_out)
     298
     299call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin  ')
     300call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ')
     301
    264302
    265303      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_init'
    … …  
    268306         write(lunout,*) 'lmin ',lmin(igout)
    269307         write(lunout,*) 'lalim ',lalim(igout)
     308         write(lunout,*) ' ig l alim_star thetav'
     309         write(lunout,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) &
     310     &   ,ztv(igout,l),l=1,lalim(igout)+4)
    270311      endif
    271312
    272 !-----------------------------------------------------------------------------------
    273 !calcul des caracteristiques du thermique sec pour le calcul de detr et la fermeture
    274 !-----------------------------------------------------------------------------------
     313!-----------------------------------------------------------------------------
     314!  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
     315!     panache sec conservatif (e=d=0) alimente selon alim_star
     316!     Il s'agit d'un calcul de type CAPE
     317!     zmax_sec est utilisé pour déterminer la géométrie du thermique.
     318!------------------------------------------------------------------------------
    275319!
    276320      CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
    277321     &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)
    278322
     323call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
     324call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
     325
    279326      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_dry'
     327      if (lev_out.ge.10) then
     328         write(lunout,*) 'Dans thermcell_main 1b'
     329         write(lunout,*) 'lmin ',lmin(igout)
     330         write(lunout,*) 'lalim ',lalim(igout)
     331         write(lunout,*) ' ig l alim_star entr_star detr_star f_star '
     332         write(lunout,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
     333     &    ,l=1,lalim(igout)+4)
     334      endif
    280335
    281336
    … …  
    289344     &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
    290345     &           ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out)
     346      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
     347      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
    291348
    292349      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_plume'
    … …  
    297354         write(lunout,*) ' ig l alim_star entr_star detr_star f_star '
    298355         write(lunout,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
    299      &    ,f_star(igout,l+1),l=1,lmax(igout))
     356     &    ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
    300357      endif
    301358
    … …  
    307364     &           zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
    308365
     366
     367      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
     368      call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
     369      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
     370      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
     371
    309372      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_height'
    310373
    … …  
    322385!-------------------------------------------------------------------------------
    323386
    324       CALL thermcell_flux(ngrid,nlay,ptimestep,masse, &
     387      CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
    325388     &       lalim,lmax,alim_star,  &
    326389     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
    … …  
    328391
    329392      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_flux'
     393      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
     394      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
    330395
    331396!c------------------------------------------------------------------
    … …  
    391456      enddo
    392457     
     458      print*,'14a OK convect8'
    393459! calcul du niveau de condensation
    394460! initialisation
    … …  
    410476         enddo
    411477      enddo
     478      print*,'14b OK convect8'
    412479      do k=nlay,1,-1
    413480         do ig=1,ngrid
    … …  
    418485         enddo
    419486      enddo
     487      print*,'14c OK convect8'
    420488!calcul des moments
    421489!initialisation
    … …  
    429497         enddo
    430498      enddo     
     499      print*,'14d OK convect8'
    431500      do l=1,nlay
    432501         do ig=1,ngrid
    … …  
    440509            q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
    441510!test: on calcul q2/po=ratqsc
    442             ratqscth(ig,l)=sqrt(q2(ig,l))/(po(ig,l)*1000.)
     511            ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
    443512         enddo
    444513      enddo
    445514!calcul du ratqscdiff
     515      print*,'14e OK convect8'
    446516      var=0.
    447517      vardiff=0.
    … …  
    452522         enddo
    453523      enddo
     524      print*,'14f OK convect8'
    454525      do ig=1,ngrid
    455526          do l=1,lalim(ig)
    … …  
    462533          enddo
    463534      enddo
     535      print*,'14g OK convect8'
    464536      do l=1,nlay
    465537         do ig=1,ngrid
    … …  
    496568
    497569!-----------------------------------------------------------------------------
     570
     571      subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
     572
     573      integer klon,klev
     574      real pplev(klon,klev+1),pplay(klon,klev)
     575      real ztv(klon,klev)
     576      real po(klon,klev)
     577      real ztva(klon,klev)
     578      real zqla(klon,klev)
     579      real f_star(klon,klev)
     580      real zw2(klon,klev)
     581      integer long(klon)
     582      real seuil
     583      character*21 comment
     584
     585      print*,'TEST ',comment
     586!  test sur la hauteur des thermiques ...
     587         do i=1,klon
     588            if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
     589               print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
     590               print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
     591               do k=1,klev
     592                  write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)
     593               enddo
     594!              stop
     595            endif
     596         enddo
     597
     598
     599      return
     600      end
     601
Note: See TracChangeset for help on using the changeset viewer.