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

Location:
LMDZ4/trunk/libf/phylmd
Files:
1 added
11 edited

Legend:

Unmodified
Added
Removed

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

    r941 r972  
    1818use phys_state_var_mod
    1919IMPLICIT none
     20#include "iniprint.h"
    2021
    2122! Arguments :
    … …  
    3233INTEGER jadrs(klon*klev), jbad
    3334INTEGER jqadrs(klon*klev), jqbad
     35INTEGER kadrs(klon*klev)
     36INTEGER kqadrs(klon*klev)
    3437
    3538integer debug_level
    36 
     39logical, save :: first=.true.
     40INTEGER, SAVE :: itap
    3741!======================================================================
    3842! Initialisations
    3943
    4044debug_level=10
    41 
     45     if (first) then
     46        itap=0
     47        first=.false.
     48     endif
     49! Incrementer le compteur de la physique
     50     itap   = itap + 1
    4251!======================================================================
    4352! Ajout des tendances sur le vent et l'eau liquide
    … …  
    6271            jbad = jbad + 1
    6372            jadrs(jbad) = i
     73            kadrs(jbad) = k
    6474            ENDIF
    6575            IF ( zq<0. .or. zq>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
    6676            jqbad = jqbad + 1
    6777            jqadrs(jqbad) = i
     78            kqadrs(jqbad) = k
    6879            ENDIF
    6980            t_seri(i,k)=zt
    … …  
    97108         print*,'l    T     dT       Q     dQ    '
    98109         DO k = 1, klev
    99             write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
     110           zq=q_seri(i,k)+zdq(i,k)
     111           if (zq.lt.1.e-15) then
     112              if (q_seri(i,k).lt.1.e-15) then
     113!              print*,' cas q_seri<1.e-15 i k q_seri zq zdq :',i,k,q_seri(i,k),zq,zdq(i,k)
     114               q_seri(i,k)=1.e-15
     115               zdq(i,k)=(1.e-15-q_seri(i,k))
     116              endif
     117           endif
    100118!           zq=q_seri(i,k)+zdq(i,k)
    101119!           if (zq.lt.1.e-15) then
    … …  
    103121!           endif
    104122         ENDDO
    105          call print_debug_phys(i,debug_level,text)
    106123      ENDDO
    107124ENDIF
    108125!
    109126
     127!IM ajout memes tests pour reverifier les jbad, jqbad beg
     128      jbad=0
     129      jqbad=0
     130      DO k = 1, klev
     131         DO i = 1, klon
     132            IF ( t_seri(i,k)>370. .or. t_seri(i,k)<130. .or. abs(zdt(i,k))>50. ) then
     133            jbad = jbad + 1
     134            jadrs(jbad) = i
     135!           if(prt_level.ge.10) THEN
     136!             print*,'cas2 i k t_seri zdt',i,k,t_seri(i,k),zdt(i,k)
     137!           endif
     138            ENDIF
     139            IF ( q_seri(i,k)<0. .or. q_seri(i,k)>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
     140            jqbad = jqbad + 1
     141            jqadrs(jqbad) = i
     142            kqadrs(jqbad) = k
     143!           if(prt_level.ge.10) THEN
     144!             print*,'cas2 i k q_seri zdq',i,k,q_seri(i,k),zdq(i,k)
     145!           endif
     146            ENDIF
     147         ENDDO
     148      ENDDO
     149IF (jbad .GT. 0) THEN
     150      DO j = 1, jbad
     151         i=jadrs(j)
     152         k=kadrs(j)
     153         print*,'PLANTAGE2 POUR LE POINT i itap rlon rlat txt jbad zdt t',i,itap,rlon(i),rlat(i),text,jbad, &
     154       &        zdt(i,k),t_seri(i,k)-zdt(i,k)
     155!        if(prt_level.ge.10) THEN
     156         if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN
     157          print*,'l    T     dT       Q     dQ    '
     158          DO k = 1, klev
     159             write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
     160          ENDDO
     161          call print_debug_phys(i,debug_level,text)
     162         endif
     163      ENDDO
     164ENDIF
     165!
     166IF (jqbad .GT. 0) THEN
     167      DO j = 1, jqbad
     168         i=jqadrs(j)
     169         k=kqadrs(j)
     170         print*,'WARNING  : EAU2 POUR LE POINT i itap rlon rlat txt jqbad zdq q zdql ql',i,itap,rlon(i),rlat(i),text,jqbad,&
     171       &        zdq(i,k), q_seri(i,k)-zdq(i,k), zdql(i,k), ql_seri(i,k)-zdql(i,k)
     172!        if(prt_level.ge.10) THEN
     173         if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN
     174          print*,'l    T     dT       Q     dQ    '
     175          DO k = 1, klev
     176            write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
     177          ENDDO
     178          call print_debug_phys(i,debug_level,text)
     179         endif
     180      ENDDO
     181ENDIF
     182
     183      CALL hgardfou(t_seri,ftsol,text)
    110184      RETURN
    111185      END

  • LMDZ4/trunk/libf/phylmd/cv3a_compress.F

    r879 r972  
    123123
    124124      if (nn.ne.ncum) then
    125          print*,'strange! nn not equal to ncum: ',nn,ncum
     125         print*,'WARNING nn not equal to ncum: ',nn,ncum
    126126         stop
    127127      endif
    … …  
    150150 150  continue
    151151
     152      if (nn.ne.ncum) then
     153         print*,'WARNING nn not equal to ncum: ',nn,ncum
     154         stop
     155      endif
     156
    152157      RETURN
    153158      END

  • LMDZ4/trunk/libf/phylmd/hgardfou.F

    r941 r972  
    5151           ok = .FALSE.
    5252           DO i = 1, jbad
    53              PRINT *,'i,k,temperature =',jadrs(i),k,zt(jadrs(i)),
    54      $       rlon(jadrs(i)),rlat(jadrs(i))
     53            PRINT *,'i,k,temperature rlon rlat=',jadrs(i),k,zt(jadrs(i))
     54     $      ,rlon(jadrs(i)),rlat(jadrs(i))
    5555           ENDDO
    5656         ENDIF
    … …  
    7070           ok = .FALSE.
    7171           DO i = 1, jbad
    72              PRINT *,'i,k,temperature =',jadrs(i),k,zt(jadrs(i)),
    73      $       rlon(jadrs(i)),rlat(jadrs(i))
     72            PRINT *,'i,k,temperature rlon rlat=',jadrs(i),k,zt(jadrs(i))
     73     $      ,rlon(jadrs(i)),rlat(jadrs(i))
    7474           ENDDO
    7575         ENDIF

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

    r889 r972  
    880880!****************************************************************************************
    881881! H and Q
    882        print *,'pbl_surface: ok_flux_surf=',ok_flux_surf
    883        print *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
     882!      print *,'pbl_surface: ok_flux_surf=',ok_flux_surf
     883!      print *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
    884884       if (ok_flux_surf) then
    885885          y_flux_t1(:) =  fsens

  • LMDZ4/trunk/libf/phylmd/thermcell.h

    r879 r972  
    11      integer iflag_thermals,nsplit_thermals
    2       real r_aspect_thermals,l_mix_thermals,tho_thermals
     2      real r_aspect_thermals,l_mix_thermals,tau_thermals
    33      integer w2di_thermals,isplit
    44      integer iflag_coupl,iflag_clos,iflag_wake
    55
    66      common/ctherm1/iflag_thermals,nsplit_thermals
    7       common/ctherm2/r_aspect_thermals,l_mix_thermals,tho_thermals
     7      common/ctherm2/r_aspect_thermals,l_mix_thermals,tau_thermals
    88      common/ctherm3/w2di_thermals
    99      common/ctherm4/iflag_coupl,iflag_clos,iflag_wake

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

    r938 r972  
    11      SUBROUTINE thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
    2      &   zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,f0,lev_out)
     2     &   zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
    33
    44!-------------------------------------------------------------------------
    … …  
    2424
    2525      REAL f(ngrid)
    26       REAL f0(ngrid)
    2726
    2827      do ig=1,ngrid
    … …  
    5251             f(ig)=wmax_sec(ig)/(max(500.,zmax_sec(ig))*r_aspect  &
    5352     &             *alim_star2(ig))
    54              f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
    55      &                     zmax_sec(ig))*wmax_sec(ig))
     53!            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
     54!    &                     zmax_sec(ig))*wmax_sec(ig))
    5655             else
    5756             f(ig)=wmax(ig)/zdenom
    58             f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
    59      &                     zmax(ig))*wmax(ig))
     57!            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
     58!     &                     zmax(ig))*wmax(ig))
    6059             endif
    6160         endif
    62          f0(ig)=f(ig)
     61!         f0(ig)=f(ig)
    6362      enddo
    6463      if (prt_level.ge.1) print*,'apres fermeture'

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

    r938 r972  
    2323      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
    2424
     25      real zzm
     26
    2527      integer ig,k
     28      real cfl
     29
     30      real qold(ngrid,nlay)
     31      real ztimestep
     32      integer niter,iter
     33
     34
     35
     36! Calcul du critere CFL pour l'advection dans la subsidence
     37      do k=1,nlay
     38         do ig=1,ngrid
     39            zzm=masse(ig,k)/ptimestep
     40            cfl=max(cfl,fm(ig,k)/zzm)
     41            if (entr(ig,k).gt.zzm) then
     42               print*,'entr dt > m ',entr(ig,k)*ptimestep,masse(ig,k)
     43               stop
     44            endif
     45         enddo
     46      enddo
     47
     48!IM 090508     print*,'CFL CFL CFL CFL ',cfl
     49
     50#undef CFL
     51#ifdef CFL
     52! On subdivise le calcul en niter pas de temps.
     53      niter=int(cfl)+1
     54#else
     55      niter=1
     56#endif
     57
     58      ztimestep=ptimestep/niter
     59      qold=q
     60
     61
     62do iter=1,niter
     63      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
    2664
    2765!   calcul du detrainement
    28 
    29       if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
    30 
    3166      do k=1,nlay
    3267         do ig=1,ngrid
    … …  
    5691      do k=2,nlay
    5792         do ig=1,ngrid
    58             if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
     93            if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt.  &
    5994     &         1.e-5*masse(ig,k)) then
    6095         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
    … …  
    72107      enddo
    73108
     109! Calcul du flux subsident
     110
    74111      do k=2,nlay
    75112         do ig=1,ngrid
    76 !             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
     113#undef centre
     114#ifdef centre
     115             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
     116#else
     117
     118#define plusqueun
     119#ifdef plusqueun
     120! Schema avec advection sur plus qu'une maille.
     121            zzm=masse(ig,k)/ztimestep
     122            if (fm(ig,k)>zzm) then
     123               wqd(ig,k)=zzm*q(ig,k)+(fm(ig,k)-zzm)*q(ig,k+1)
     124            else
     125               wqd(ig,k)=fm(ig,k)*q(ig,k)
     126            endif
     127#else
    77128            wqd(ig,k)=fm(ig,k)*q(ig,k)
     129#endif
     130#endif
     131
    78132            if (wqd(ig,k).lt.0.) then
    79133!               print*,'wqd<0!!!'
    … …  
    86140      enddo
    87141     
     142
     143! Calcul des tendances
    88144      do k=1,nlay
    89145         do ig=1,ngrid
    90             dq(ig,k)=(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
     146            q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
    91147     &               -wqd(ig,k)+wqd(ig,k+1))  &
    92      &               /masse(ig,k)
     148     &               *ztimestep/masse(ig,k)
    93149!            if (dq(ig,k).lt.0.) then
    94150!               print*,'dq<0!!!'
    … …  
    97153      enddo
    98154
     155
     156enddo
     157
     158
     159! Calcul des tendances
     160      do k=1,nlay
     161         do ig=1,ngrid
     162            dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
     163            q(ig,k)=qold(ig,k)
     164         enddo
     165      enddo
     166
    99167      return
    100168      end

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

    r938 r972  
    5454      enddo
    5555
     56      print*,'WARNING on initialise gamma(1:ngrid,1)=0.'
     57      gamma(1:ngrid,1)=0.
    5658      do k=2,nlay
    5759         do ig=1,ngrid
    … …  
    6062!   On itère sur la valeur du coeff de freinage.
    6163!              gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
     64!IM 060508 beg
     65!             if(0.5*(fraca(ig,k+1)+fraca(ig,k)).LT.0.) THEN
     66!              print*,'th_dv2 ig k fraca(:,k) fraca(:k+1)', &
     67!    &         ig,k,fraca(ig,k),fraca(ig,k+1)
     68!             endif
     69!             if(larga(ig).EQ.0.) THEN
     70!              print*,'th_dv2 ig larga=0.',ig
     71!             endif
     72              if(larga(ig).GT.0.) THEN
     73!IM 060508 end
    6274               gamma0=masse(ig,k)  &
    6375     &         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )  &
    6476     &         *0.5/larga(ig)  &
    6577     &         *1.
     78!IM 060508 beg
     79              else
     80               if(prt_level.GE.10) print*,'WARNING cas ELSE on initialise gamma0=0.'
     81               gamma0=0.
     82              endif !(larga(ig).GT.0.) THEN
     83!IM 060508 end
    6684!    s         *0.5
    6785!              gamma0=0.
    … …  
    117135      do k=1,nlay
    118136         do ig=1,ngrid
     137!IM
     138         if(prt_level.GE.10) print*,'th_dv2 ig k gamma entr detr ua ue va ve wud wvd masse',ig,k,gamma(ig,k), &
     139     &   entr(ig,k),detr(ig,k),ua(ig,k),ue(ig,k),va(ig,k),ve(ig,k),wud(ig,k),wvd(ig,k),wud(ig,k+1),wvd(ig,k+1), &
     140     &   masse(ig,k)
     141!
    119142            du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)  &
    120143     &               -(entr(ig,k)+gamma(ig,k))*ue(ig,k)  &

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

    r938 r972  
    1       SUBROUTINE thermcell_init(ngrid,nlay,ztv,zlev,  &
     1      SUBROUTINE thermcell_init(ngrid,nlay,ztv,zlay,zlev,  &
    22     &                  lalim,lmin,alim_star,alim_star_tot,lev_out)
    33
    … …  
    1212      INTEGER ngrid,nlay
    1313      REAL ztv(ngrid,nlay)
    14       REAL zlev(ngrid,nlay)
     14      REAL zlay(ngrid,nlay)
     15      REAL zlev(ngrid,nlay+1)
    1516!arguments de sortie
    1617      INTEGER lalim(ngrid)
    … …  
    2021      integer lev_out                           ! niveau pour les print
    2122     
     23      REAL zzalim(ngrid)
    2224!CR: ponderation entrainement des couches instables
    2325!def des alim_star tels que alim=f*alim_star     
    … …  
    5860      enddo
    5961!
     62      zzalim(:)=0.
     63      do l=1,nlay-1
     64         do ig=1,ngrid
     65             if (l<lalim(ig)) then
     66                zzalim(ig)=zzalim(ig)+zlay(ig,l)*(ztv(ig,l)-ztv(ig,l+1))
     67             endif
     68          enddo
     69      enddo
     70      do ig=1,ngrid
     71          if (lalim(ig)>1) then
     72             zzalim(ig)=zlay(ig,1)+zzalim(ig)/(ztv(ig,1)-ztv(ig,lalim(ig)))
     73          else
     74             zzalim(ig)=zlay(ig,1)
     75          endif
     76      enddo
     77
     78      if(prt_level.GE.10) print*,'ZZALIM LALIM ',zzalim,lalim,zlay(1,lalim(1))
     79
    6080! definition de l'entrainement des couches
     81      if (1.eq.1) then
    6182      do l=1,nlay-1
    6283         do ig=1,ngrid
    … …  
    6990         enddo
    7091      enddo
     92      else
     93      do l=1,nlay-1
     94         do ig=1,ngrid
     95            if (ztv(ig,l).gt.ztv(ig,l+1).and.  &
     96     &          l.ge.lmin(ig).and.l.lt.lalim(ig)) then
     97             alim_star(ig,l)=max(3.*zzalim(ig)-zlay(ig,l),0.) &
     98     &        *(zlev(ig,l+1)-zlev(ig,l))
     99            endif
     100         enddo
     101      enddo
     102      endif
    71103     
    72104! pas de thermique si couche 1 stable

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

    r940 r972  
    22! $Header$
    33!
    4       SUBROUTINE thermcell_main(ngrid,nlay,ptimestep  &
     4      SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep  &
    55     &                  ,pplay,pplev,pphi,debut  &
    66     &                  ,pu,pv,pt,po  &
    77     &                  ,pduadj,pdvadj,pdtadj,pdoadj  &
    8      &                  ,fm0,entr0,zqla,lmax  &
     8     &                  ,fm0,entr0,detr0,zqla,lmax  &
    99     &                  ,ratqscth,ratqsdiff,zqsatth  &
    10      &                  ,r_aspect,l_mix,w2di,tho &
     10     &                  ,r_aspect,l_mix,tau_thermals &
    1111     &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
    1212     &                  ,zmax0, f0)
    1313
    14       use dimphy
     14      USE dimphy
    1515      IMPLICIT NONE
    1616
    … …  
    5050!   ----------
    5151
    52       INTEGER ngrid,nlay,w2di,tho
     52!IM 140508
     53      INTEGER itap
     54
     55      INTEGER ngrid,nlay,w2di
     56      real tau_thermals
    5357      real ptimestep,l_mix,r_aspect
    5458      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
    … …  
    6165!   local:
    6266!   ------
     67
     68      integer icount
     69      data icount/0/
     70      save icount
    6371
    6472      integer,save :: igout=1
    … …  
    7987!FH/IM     save zmax0
    8088
     89      real lambda
     90
    8191      real zlev(klon,klev+1),zlay(klon,klev)
    8292      real deltaz(klon,klev)
    83       REAL zh(klon,klev),zdhadj(klon,klev)
     93      REAL zh(klon,klev)
    8494      real zthl(klon,klev),zdthladj(klon,klev)
    8595      REAL ztv(klon,klev)
    … …  
    97107      real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
    98108      real q2(klon,klev)
    99 !      common/comtherm/thetath2,wth2
     109! FH probleme de dimensionnement avec l'allocation dynamique
     110!     common/comtherm/thetath2,wth2
    100111   
    101112      real ratqscth(klon,klev)
    … …  
    109120
    110121      logical sorties
    111       real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
     122      real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
    112123      real zpspsk(klon,klev)
    113124
    114125      real wmax(klon)
    115126      real wmax_sec(klon)
    116       real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
    117       real detr0(klon,klev)
    118       real fm(klon,klev+1),entr(klon,klev)
     127      real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
     128      real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
    119129
    120130      real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
    … …  
    164174      seuil=0.25
    165175
     176      if (debut)  then
     177         fm0=0.
     178         entr0=0.
     179         detr0=0.
     180      endif
     181
     182      fm=0. ; entr=0. ; detr=0.
     183
     184      icount=icount+1
     185
     186!IM 090508 beg
     187!print*,'====================================================================='
     188!print*,'====================================================================='
     189!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
     190!print*,'====================================================================='
     191!print*,'====================================================================='
     192!IM 090508 end
     193
    166194      if (prt_level.ge.1) print*,'thermcell_main V4'
    167195
    … …  
    176204!Initialisation
    177205!
    178       do ig=1,klon     
     206!    IF (1.eq.0) THEN
     207!     do ig=1,klon     
    179208!FH/IM 130308     if ((debut).or.((.not.debut).and.(f0(ig).lt.1.e-10))) then
    180       if ((.not.debut).and.(f0(ig).lt.1.e-10)) then
    181             f0(ig)=1.e-5
    182             zmax0(ig)=40.
     209!     if ((.not.debut).and.(f0(ig).lt.1.e-10)) then
     210!           f0(ig)=1.e-5
     211!           zmax0(ig)=40.
    183212!v1d        therm=.false.
    184       endif
    185       enddo
    186 
     213!     endif
     214!     enddo
     215!    ENDIF !(1.eq.0) THEN
     216     print*,'WARNING thermcell_main f0=max(f0,1.e-2)'
     217     do ig=1,klon
     218      if (prt_level.ge.20) then
     219       print*,'th_main ig f0',ig,f0(ig)
     220      endif
     221         f0(ig)=max(f0(ig),1.e-2)
     222!IMmarche pas ?!       if (f0(ig)<1.e-2) f0(ig)=1.e-2
     223     enddo
    187224
    188225!-----------------------------------------------------------------------
    … …  
    238275         rho(:,l)=pplay(:,l)/(zpspsk(:,l)*RD*ztv(:,l))
    239276      enddo
     277
     278!IM
     279      print*,'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
     280      rhobarz(:,1)=rho(:,1)
    240281
    241282      do l=2,nlay
    … …  
    309350!
    310351      entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
    311       CALL thermcell_init(ngrid,nlay,ztv,zlev,  &
     352      CALL thermcell_init(ngrid,nlay,ztv,zlay,zlev,  &
    312353     &                    lalim,lmin,alim_star,alim_star_tot,lev_out)
    313354
    … …  
    318359      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_init'
    319360      if (prt_level.ge.10) then
    320          write(lunout,*) 'Dans thermcell_main 1'
    321          write(lunout,*) 'lmin ',lmin(igout)
    322          write(lunout,*) 'lalim ',lalim(igout)
    323          write(lunout,*) ' ig l alim_star thetav'
    324          write(lunout,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) &
     361         write(lunout1,*) 'Dans thermcell_main 1'
     362         write(lunout1,*) 'lmin ',lmin(igout)
     363         write(lunout1,*) 'lalim ',lalim(igout)
     364         write(lunout1,*) ' ig l alim_star thetav'
     365         write(lunout1,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) &
    325366     &   ,ztv(igout,l),l=1,lalim(igout)+4)
    326367      endif
    … …  
    346387      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
    347388      if (prt_level.ge.10) then
    348          write(lunout,*) 'Dans thermcell_main 1b'
    349          write(lunout,*) 'lmin ',lmin(igout)
    350          write(lunout,*) 'lalim ',lalim(igout)
    351          write(lunout,*) ' ig l alim_star entr_star detr_star f_star '
    352          write(lunout,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
     389         write(lunout1,*) 'Dans thermcell_main 1b'
     390         write(lunout1,*) 'lmin ',lmin(igout)
     391         write(lunout1,*) 'lalim ',lalim(igout)
     392         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
     393         write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
    353394     &    ,l=1,lalim(igout)+4)
    354395      endif
    … …  
    360401!--------------------------------------------------------------------------------
    361402!
    362       CALL thermcell_plume(ngrid,nlay,ztv,zthl,po,zl,rhobarz,  &
     403      if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
     404!IM 140508   CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
     405      CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
    363406     &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,  &
    364407     &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
    365408     &           ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out)
     409      if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
     410
    366411      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
    367412      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
    … …  
    369414      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
    370415      if (prt_level.ge.10) then
    371          write(lunout,*) 'Dans thermcell_main 2'
    372          write(lunout,*) 'lmin ',lmin(igout)
    373          write(lunout,*) 'lalim ',lalim(igout)
    374          write(lunout,*) ' ig l alim_star entr_star detr_star f_star '
    375          write(lunout,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
     416         write(lunout1,*) 'Dans thermcell_main 2'
     417         write(lunout1,*) 'lmin ',lmin(igout)
     418         write(lunout1,*) 'lalim ',lalim(igout)
     419         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
     420         write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
    376421     &    ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
    377422      endif
    … …  
    397442
    398443      CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
    399      &   zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,f0,lev_out)
     444     &   zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
    400445
    401446      if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
     447
     448      if (tau_thermals>1.) then
     449         lambda=exp(-ptimestep/tau_thermals)
     450         f0=(1.-lambda)*f+lambda*f0
     451      else
     452         f0=f
     453      endif
     454
     455! Test valable seulement en 1D mais pas genant
     456      if (.not. (f0(1).ge.0.) ) then
     457           stop'Dans thermcell_main'
     458      endif
    402459
    403460!-------------------------------------------------------------------------------
    … …  
    405462!-------------------------------------------------------------------------------
    406463
    407       CALL thermcell_flux(ngrid,nlay,ptimestep,masse, &
     464      CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
    408465     &       lalim,lmax,alim_star,  &
    409466     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
    410      &       detr,zqla,zmax,lev_out,lunout,igout)
     467     &       detr,zqla,lev_out,lunout1,igout)
     468!IM 060508    &       detr,zqla,zmax,lev_out,lunout,igout)
    411469
    412470      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
    … …  
    414472      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
    415473
    416 !c------------------------------------------------------------------
    417 !   calcul du transport vertical
    418 !------------------------------------------------------------------
    419 
    420       if (w2di.eq.1) then
    421          fm0=fm0+ptimestep*(fm-fm0)/float(tho)
    422          entr0=entr0+ptimestep*(entr-entr0)/float(tho)
     474!------------------------------------------------------------------
     475!   On ne prend pas directement les profils issus des calculs precedents
     476!   mais on s'autorise genereusement une relaxation vers ceci avec
     477!   une constante de temps tau_thermals (typiquement 1800s).
     478!------------------------------------------------------------------
     479
     480      if (tau_thermals>1.) then
     481         lambda=exp(-ptimestep/tau_thermals)
     482         fm0=(1.-lambda)*fm+lambda*fm0
     483         entr0=(1.-lambda)*entr+lambda*entr0
     484!        detr0=(1.-lambda)*detr+lambda*detr0
    423485      else
    424486         fm0=fm
    … …  
    426488         detr0=detr
    427489      endif
     490
     491!c------------------------------------------------------------------
     492!   calcul du transport vertical
     493!------------------------------------------------------------------
    428494
    429495      call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
    … …  
    453519!------------------------------------------------------------------
    454520
     521!IM 090508 
    455522      if (1.eq.1) then
     523!IM 070508 vers. _dq       
     524!     if (1.eq.0) then
    456525
    457526
    … …  
    461530         call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse  &
    462531     &    ,fraca,zmax  &
    463      &    ,zu,zv,pduadj,pdvadj,zua,zva,igout,lev_out)
     532     &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
     533!IM 050508    &    ,zu,zv,pduadj,pdvadj,zua,zva,igout,lev_out)
    464534      else
    465535
    … …  
    478548      enddo
    479549
    480 !    print*,'14 OK convect8'
     550      if (prt_level.ge.1) print*,'14 OK convect8'
    481551!------------------------------------------------------------------
    482552!   Calculs de diagnostiques pour les sorties
    … …  
    485555     
    486556      if (sorties) then
    487 !      print*,'14a OK convect8'
     557      if (prt_level.ge.1) print*,'14a OK convect8'
    488558! calcul du niveau de condensation
    489559! initialisation
    … …  
    505575         enddo
    506576      enddo
    507 !    print*,'14b OK convect8'
     577      if (prt_level.ge.1) print*,'14b OK convect8'
    508578      do k=nlay,1,-1
    509579         do ig=1,ngrid
    … …  
    514584         enddo
    515585      enddo
    516 !    print*,'14c OK convect8'
     586      if (prt_level.ge.1) print*,'14c OK convect8'
    517587!calcul des moments
    518588!initialisation
    … …  
    526596         enddo
    527597      enddo     
    528 !     print*,'14d OK convect8'
     598      if (prt_level.ge.1) print*,'14d OK convect8'
     599      print*,'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
    529600      do l=1,nlay
    530601         do ig=1,ngrid
    531602            zf=fraca(ig,l)
    532603            zf2=zf/(1.-zf)
     604!
     605      if (prt_level.ge.10) print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
     606!
     607      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)
    533608            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))**2
    534             wth2(ig,l)=zf2*(zw2(ig,l))**2
     609            if(zw2(ig,l).gt.1.e-10) then
     610             wth2(ig,l)=zf2*(zw2(ig,l))**2
     611            else
     612             wth2(ig,l)=0.
     613            endif
    535614!           print*,'wth2=',wth2(ig,l)
    536615            wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))  &
    537616     &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
     617      if (prt_level.ge.10) print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
    538618            q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
    539619!test: on calcul q2/po=ratqsc
    … …  
    548628      n_int(ig)=0
    549629      enddo
    550       do ig=1,ngrid
    551 !      do l=nivcon(ig),lmax(ig)
    552       do l=1,lmax(ig)
    553       alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)
    554       ale_int(ig)=ale_int(ig)+0.5*zw2(ig,l)**2
    555       n_int(ig)=n_int(ig)+1
     630!
     631      do l=1,nlay
     632      do ig=1,ngrid
     633       if(l.LE.lmax(ig)) THEN
     634        alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)
     635        ale_int(ig)=ale_int(ig)+0.5*zw2(ig,l)**2
     636        n_int(ig)=n_int(ig)+1
     637       endif
    556638      enddo
    557639      enddo
    … …  
    638720
    639721!calcul du ratqscdiff
    640 !    print*,'14e OK convect8'
     722      if (prt_level.ge.1) print*,'14e OK convect8'
    641723      var=0.
    642724      vardiff=0.
    … …  
    647729         enddo
    648730      enddo
    649 !    print*,'14f OK convect8'
     731      if (prt_level.ge.1) print*,'14f OK convect8'
    650732      do ig=1,ngrid
    651733          do l=1,lalim(ig)
    … …  
    658740          enddo
    659741      enddo
    660 !    print*,'14g OK convect8'
     742      if (prt_level.ge.1) print*,'14g OK convect8'
    661743      do l=1,nlay
    662744         do ig=1,ngrid
    … …  
    679761
    680762
    681 ! #define troisD
     763#define troisD
     764#undef troisD
    682765      if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
    683766#ifdef troisD
    … …  
    689772      if (prt_level.ge.1) print*,'thermcell_main FIN  OK'
    690773
     774!     if(icount.eq.501) stop'au pas 301 dans thermcell_main'
    691775      return
    692776      end
    … …  
    717801!  test sur la hauteur des thermiques ...
    718802         do i=1,klon
    719             if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
     803!IMtemp           if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
     804           if (prt_level.ge.10) then
    720805               print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
    721806               print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
    … …  
    724809               enddo
    725810!              stop
    726             endif
     811           endif
    727812         enddo
    728813

  • 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.