source: LMDZ6/trunk/libf/phylmd/calltherm.F90 @ 5322

Last change on this file since 5322 was 5291, checked in by abarral, 3 weeks ago

Move thermcell_old.h iotd.h to module

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.6 KB
RevLine 
[878]1!
[1403]2! $Id: calltherm.F90 5291 2024-10-28 15:44:39Z abarral $
[878]3!
4      subroutine calltherm(dtime  &
5     &      ,pplay,paprs,pphi,weak_inversion  &
[4690]6     &      ,u_seri_,v_seri_,t_seri_,q_seri_,t_env,q_env,zqsat,debut  &
[878]7     &      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs  &
[4837]8     &    ,fm_therm,entr_therm,detr_therm,zqasc,clwcon0,lmax,ratqscth,&
[4089]9     &       ratqsdiff,zqsatth,ale_bl,alp_bl,lalim_conv,wght_th, &
[1638]10     &       zmax0,f0,zw2,fraca,ztv,zpspsk,ztla,zthl &
11!!! nrlmd le 10/04/2012
12     &      ,pbl_tke,pctsrf,omega,airephy &
13     &      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
[4843]14     &      ,n2,s2,strig,zcong,ale_bl_stat &
[1638]15     &      ,therm_tke_max,env_tke_max &
16     &      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
17     &      ,alp_bl_conv,alp_bl_stat &
18!!! fin nrlmd le 10/04/2012
[4089]19     &      ,zqla,ztva &
20#ifdef ISO         
21     &      ,xt_seri,d_xt_ajs &
22#ifdef DIAGISO         
23     &      ,q_the,xt_the &
24#endif
25#endif         
26     &   )
[878]27
[940]28      USE dimphy
[1785]29      USE indice_sol_mod
[2311]30      USE print_control_mod, ONLY: prt_level,lunout
[4590]31      USE lmdz_thermcell_alp, ONLY: thermcell_alp
32      USE lmdz_thermcell_main, ONLY: thermcell_main
33      USE lmdz_thermcell_old, ONLY: thermcell, thermcell_2002, thermcell_eau, calcul_sec, thermcell_sec
[4089]34#ifdef ISO
[4143]35      use infotrac_phy, ONLY: ntiso
[4089]36#ifdef ISOVERIF
37      USE isotopes_mod, ONLY: iso_eau,iso_HDO
38      USE isotopes_verif_mod, ONLY: iso_verif_aberrant_enc_vect2D, &
39        iso_verif_egalite_vect2D
40#endif   
41#endif
[5282]42      USE clesphys_mod_h
[5291]43      USE thermcell_old_mod_h, ONLY: r_aspect_thermals, l_mix_thermals, w2di_thermals
[878]44      implicit none
45
[1638]46
[973]47!IM 140508
[1776]48      INTEGER, SAVE ::  itap
49!$OMP THREADPRIVATE(itap)
[878]50      REAL dtime
51      LOGICAL debut
[973]52      LOGICAL logexpr0, logexpr2(klon,klev), logexpr1(klon)
53      REAL fact(klon)
54      INTEGER nbptspb
55
[4678]56      REAL, DIMENSION(klon,klev), INTENT(IN) :: u_seri_,v_seri_
[4837]57      REAL, DIMENSION(klon,klev), INTENT(IN) :: t_seri_,q_seri_
58      REAL, DIMENSION(klon,klev), INTENT(IN) :: t_env,q_env
[4678]59      REAL, DIMENSION(klon,klev) :: u_seri,v_seri
60      REAL, DIMENSION(klon,klev) :: t_seri,q_seri
61      REAL, DIMENSION(klon,klev) :: qmemoire
[878]62      REAL weak_inversion(klon)
63      REAL paprs(klon,klev+1)
64      REAL pplay(klon,klev)
65      REAL pphi(klon,klev)
66      real zlev(klon,klev+1)
[879]67!test: on sort lentr et a* pour alimenter KE
68      REAL wght_th(klon,klev)
69      INTEGER lalim_conv(klon)
[1026]70      REAL zw2(klon,klev+1),fraca(klon,klev+1)
[878]71
72!FH Update Thermiques
73      REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
74      REAL d_u_ajs(klon,klev),d_v_ajs(klon,klev)
[973]75      real fm_therm(klon,klev+1)
76      real entr_therm(klon,klev),detr_therm(klon,klev)
[878]77
78!********************************************************
79!     declarations
[1403]80      LOGICAL flag_bidouille_stratocu
[973]81      real fmc_therm(klon,klev+1),zqasc(klon,klev)
[878]82      real zqla(klon,klev)
[1790]83      real ztv(klon,klev),ztva(klon,klev)
[1403]84      real zpspsk(klon,klev)
85      real ztla(klon,klev)
86      real zthl(klon,klev)
[878]87      real wmax_sec(klon)
[4843]88      real zcong(klon)
[878]89      real zmax_sec(klon)
90      real f_sec(klon)
[973]91      real detrc_therm(klon,klev)
92! FH WARNING : il semble que ces save ne servent a rien
93!     save fmc_therm, detrc_therm
[878]94      real clwcon0(klon,klev)
95      real zqsat(klon,klev)
96      real zw_sec(klon,klev+1)
97      integer lmix_sec(klon)
98      integer lmax(klon)
99      real ratqscth(klon,klev)
100      real ratqsdiff(klon,klev)
101      real zqsatth(klon,klev) 
[879]102!nouvelles variables pour la convection
[4089]103      real ale_bl(klon)
104      real alp_bl(klon)
105      real ale(klon)
106      real alp(klon)
[879]107!RC
[927]108      !on garde le zmax du pas de temps precedent
109      real zmax0(klon), f0(klon)
[1638]110
111!!! nrlmd le 10/04/2012
112      real pbl_tke(klon,klev+1,nbsrf)
113      real pctsrf(klon,nbsrf)
114      real omega(klon,klev)
115      real airephy(klon)
116      real zlcl_th(klon),fraca0(klon),w0(klon),w_conv(klon)
117      real therm_tke_max0(klon),env_tke_max0(klon)
[4843]118      real n2(klon),s2(klon),strig(klon)
[1638]119      real ale_bl_stat(klon)
120      real therm_tke_max(klon,klev),env_tke_max(klon,klev)
121      real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
122!!! fin nrlmd le 10/04/2012
123
[878]124!********************************************************
125
[4089]126      real, dimension(klon) :: pcon
127      real, dimension(klon,klev) :: rhobarz,wth3
128      integer,dimension(klon) :: lalim
129      real, dimension(klon,klev+1) :: fm
130      real, dimension(klon,klev) :: alim_star
131      real, dimension(klon) :: zmax
[878]132
[4089]133
134
135
[878]136! variables locales
137      REAL d_t_the(klon,klev), d_q_the(klon,klev)
138      REAL d_u_the(klon,klev),d_v_the(klon,klev)
139!
[973]140      real zfm_therm(klon,klev+1),zdt
141      real zentr_therm(klon,klev),zdetr_therm(klon,klev)
142! FH A VERIFIER : SAVE INUTILES
[940]143!      save zentr_therm,zfm_therm
[973]144
[1403]145      character (len=20) :: modname='calltherm'
146      character (len=80) :: abort_message
147
[4089]148      integer i,k,isplit
[973]149      logical, save :: first=.true.
[4089]150      logical :: new_thermcell
151
152#ifdef ISO
[4143]153      REAL xt_seri(ntiso,klon,klev),xtmemoire(ntiso,klon,klev)
154      REAL d_xt_ajs(ntiso,klon,klev)
155      real d_xt_the(ntiso,klon,klev)
[4089]156#ifdef DIAGISO
157      real q_the(klon,klev)
[4143]158      real xt_the(ntiso,klon,klev)
[4089]159#endif
160      real qprec(klon,klev)
161      integer ixt
162#endif
163
164
[987]165!$OMP THREADPRIVATE(first)
[878]166!********************************************************
[973]167      if (first) then
168        itap=0
169        first=.false.
170      endif
[878]171
[4678]172      u_seri(:,:)=u_seri_(:,:)
173      v_seri(:,:)=v_seri_(:,:)
174      t_seri(:,:)=t_seri_(:,:)
175      q_seri(:,:)=q_seri_(:,:)
176
[973]177! Incrementer le compteur de la physique
178     itap   = itap + 1
179
[878]180!  Modele du thermique
181!  ===================
182!         print*,'thermiques: WARNING on passe t au lieu de t_seri'
183
184
[973]185! On prend comme valeur initiale des thermiques la valeur du pas
186! de temps precedent
187         zfm_therm(:,:)=fm_therm(:,:)
188         zdetr_therm(:,:)=detr_therm(:,:)
189         zentr_therm(:,:)=entr_therm(:,:)
190
191! On reinitialise les flux de masse a zero pour le cumul en
192! cas de splitting
[878]193         fm_therm(:,:)=0.
194         entr_therm(:,:)=0.
[973]195         detr_therm(:,:)=0.
196
[4089]197         ale_bl(:)=0.
198         alp_bl(:)=0.
[938]199         if (prt_level.ge.10) then
200          print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
201         endif
[878]202
203!   tests sur les valeurs negatives de l'eau
[973]204         logexpr0=prt_level.ge.10
205         nbptspb=0
[878]206         do k=1,klev
207            do i=1,klon
[1146]208! Attention teste abderr 19-03-09
209!               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
210                logexpr2(i,k)=.not.q_seri(i,k).ge.1.e-15
[973]211               if (logexpr2(i,k)) then
[4089]212#ifdef ISO
213                qprec(i,k)=q_seri(i,k)
214#endif
[938]215                q_seri(i,k)=1.e-15
[973]216                nbptspb=nbptspb+1
[4089]217#ifdef ISO
[4143]218                do ixt=1,ntiso
[4089]219                  xt_seri(ixt,i,k)=1.e-15*(xt_seri(ixt,i,k)/qprec(i,k))
220                  ! xt_seri(ixt,i,k)=1.e-15*(Rdefault(index_iso(ixt)))
221                enddo
222#endif
[878]223               endif
[973]224!               if (logexpr0) &
225!    &             print*,'WARN eau<0 avant therm i=',i,'  k=',k  &
226!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k)
[878]227            enddo
228         enddo
[973]229         if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
[878]230
[4089]231
232         new_thermcell=iflag_thermals>=15.and.iflag_thermals<=18
233#ifdef ISO
234      if (.not.new_thermcell) then
[5217]235           CALL abort_physic('calltherm 234','isos pas prevus ici',1)
[4089]236      endif
237#ifdef ISOVERIF
238      if (iso_eau.gt.0) then
239       call iso_verif_egalite_vect2D( &
240     &           xt_seri,q_seri, &
[4143]241     &           'calltherm 174',ntiso,klon,klev)
[4089]242      endif !if (iso_eau.gt.0) then
243#endif   
244#endif
[1403]245         zdt=dtime/REAL(nsplit_thermals)
[4089]246
247
[878]248         do isplit=1,nsplit_thermals
249
[1943]250          if (iflag_thermals>=1000) then
251            CALL thermcell_2002(klon,klev,zdt,iflag_thermals   &
[878]252     &      ,pplay,paprs,pphi  &
253     &      ,u_seri,v_seri,t_seri,q_seri  &
254     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
[1943]255     &      ,zfm_therm,zentr_therm,fraca,zw2  &
[878]256     &      ,r_aspect_thermals,30.,w2di_thermals  &
[1517]257     &      ,tau_thermals)
[878]258          else if (iflag_thermals.eq.2) then
259            CALL thermcell_sec(klon,klev,zdt  &
260     &      ,pplay,paprs,pphi,zlev  &
261     &      ,u_seri,v_seri,t_seri,q_seri  &
262     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
263     &      ,zfm_therm,zentr_therm  &
264     &      ,r_aspect_thermals,30.,w2di_thermals  &
[1517]265     &      ,tau_thermals)
[878]266          else if (iflag_thermals.eq.3) then
267            CALL thermcell(klon,klev,zdt  &
268     &      ,pplay,paprs,pphi  &
269     &      ,u_seri,v_seri,t_seri,q_seri  &
270     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
271     &      ,zfm_therm,zentr_therm  &
272     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
[1517]273     &      ,tau_thermals)
[878]274          else if (iflag_thermals.eq.10) then
275            CALL thermcell_eau(klon,klev,zdt  &
276     &      ,pplay,paprs,pphi  &
277     &      ,u_seri,v_seri,t_seri,q_seri  &
278     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
279     &      ,zfm_therm,zentr_therm  &
280     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
[1517]281     &      ,tau_thermals)
[878]282          else if (iflag_thermals.eq.11) then
[1403]283              abort_message = 'cas non prevu dans calltherm'
[2311]284              CALL abort_physic (modname,abort_message,1)
[878]285          else if (iflag_thermals.eq.12) then
286            CALL calcul_sec(klon,klev,zdt  &
287     &      ,pplay,paprs,pphi,zlev  &
288     &      ,u_seri,v_seri,t_seri,q_seri  &
289     &      ,zmax_sec,wmax_sec,zw_sec,lmix_sec  &
290     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
[973]291     &      ,tau_thermals)
[1403]292          else if (iflag_thermals==13.or.iflag_thermals==14) then
[4089]293              abort_message = 'thermcellV0_main enleve svn>2084'
294              CALL abort_physic (modname,abort_message,1)
295          else if (new_thermcell) then
[1403]296            CALL thermcell_main(itap,klon,klev,zdt  &
297     &      ,pplay,paprs,pphi,debut  &
[4690]298     &      ,u_seri,v_seri,t_seri,q_seri,t_env,q_env  &
[1403]299     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
300     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
301     &      ,ratqscth,ratqsdiff,zqsatth  &
302     &      ,zmax0,f0,zw2,fraca,ztv,zpspsk &
[4089]303     &      ,ztla,zthl,ztva &
[4843]304     &      ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax,zcong &
[4089]305#ifdef ISO         
306     &      ,xt_seri,d_xt_the &
307#endif         
308     &   )
309
310            CALL thermcell_alp(klon,klev,zdt  &                      ! in
311     &        ,pplay,paprs  &                                        ! in
312     &        ,zfm_therm,zentr_therm,lmax  &                         ! in
313     &        ,pbl_tke,pctsrf,omega,airephy &                        ! in
314     &        ,zw2,fraca &                                           ! in
315     &        ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &  ! in
[4843]316     &        ,zcong,ale,alp,lalim_conv,wght_th &                          ! out
[4089]317     &        ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &! out
[4843]318     &        ,n2,s2,strig,ale_bl_stat &                                   ! out
[4089]319     &        ,therm_tke_max,env_tke_max &                           ! out
320     &        ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &          ! out
321     &        ,alp_bl_conv,alp_bl_stat &                             ! out
322     &        )
323
[1403]324           if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
325         else
326           abort_message = 'Cas des thermiques non prevu'
[2311]327           CALL abort_physic (modname,abort_message,1)
[878]328         endif
329
[1638]330! Attention : les noms sont contre intuitif.
331! flag_bidouille_stratocu est .true. si on ne fait pas de bidouille.
332! Il aurait mieux valu avoir un nobidouille_stratocu
333! Et pour simplifier :
334! nobidouille_stratocu=.not.(iflag_thermals==13.or.iflag_thermals=15)
335! Ce serait bien de changer, mai en prenant le temps de vérifier que ca
336! fait bien ce qu'on croit.
[878]337
[1738]338       flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18
[1638]339
[1943]340! Calcul a posteriori du niveau max des thermiques pour les schémas qui
341! ne la sortent pas.
342      if (iflag_thermals<=12.or.iflag_thermals>=1000) then
343         lmax(:)=1
[1638]344         do k=1,klev-1
345            zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1)
346         enddo
[1943]347         do k=1,klev-1
348            do i=1,klon
349               if (zfm_therm(i,k+1)>0.) lmax(i)=k
350            enddo
351         enddo
[1638]352      endif
353
[973]354      fact(:)=0.
[878]355      DO i=1,klon
[1403]356       logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
357       IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals)
[973]358      ENDDO
[878]359
[973]360     DO k=1,klev
[878]361!  transformation de la derivee en tendance
[973]362            d_t_the(:,k)=d_t_the(:,k)*dtime*fact(:)
363            d_u_the(:,k)=d_u_the(:,k)*dtime*fact(:)
364            d_v_the(:,k)=d_v_the(:,k)*dtime*fact(:)
365            d_q_the(:,k)=d_q_the(:,k)*dtime*fact(:)
366            fm_therm(:,k)=fm_therm(:,k)  &
367     &      +zfm_therm(:,k)*fact(:)
368            entr_therm(:,k)=entr_therm(:,k)  &
369     &       +zentr_therm(:,k)*fact(:)
[1026]370            detr_therm(:,k)=detr_therm(:,k)  &
371     &       +zdetr_therm(:,k)*fact(:)
[4089]372#ifdef ISO
[4143]373            do ixt=1,ntiso
[4089]374              d_xt_the(ixt,:,k)=d_xt_the(ixt,:,k)*dtime*fact(:)
375            enddo
376#endif
[973]377      ENDDO
378       fm_therm(:,klev+1)=0.
[878]379
380
381
382!  accumulation de la tendance
[973]383            d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
384            d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
385            d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
386            d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
[4089]387#ifdef ISO
388            d_xt_ajs(:,:,:)=d_xt_ajs(:,:,:)+d_xt_the(:,:,:)
389#endif
[878]390
391!  incrementation des variables meteo
[973]392            t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
393            u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
394            v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
395            qmemoire(:,:)=q_seri(:,:)
396            q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
[4089]397#ifdef ISO
398            xtmemoire(:,:,:)=xt_seri(:,:,:)
399            xt_seri(:,:,:) = xt_seri(:,:,:) + d_xt_the(:,:,:)
400#ifdef ISOVERIF
401!      write(*,*) 'calltherm 350 tmp: ajout d_xt_the'
402      if (iso_HDO.gt.0) then
403!      i=479
404!      k=4
405!      write(*,*) 'xt_seri(iso_hdo,i,k),q_seri(i,k)=', &
406!     &   xt_seri(iso_hdo,i,k),q_seri(i,k)
407!      write(*,*) 'd_xt_the(iso_hdo,i,k),d_q_the(i,k)=', &
408!     &   d_xt_the(iso_hdo,i,k),d_q_the(i,k)
409      call iso_verif_aberrant_enc_vect2D( &
410     &        xt_seri,q_seri, &
[4143]411     &        'calltherm 353, apres ajout d_xt_the',ntiso,klon,klev)
[4089]412      endif     
413#endif
414#endif
[1403]415           if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
[878]416
[879]417       DO i=1,klon
418            fm_therm(i,klev+1)=0.
[4089]419            ale_bl(i)=ale_bl(i)+ale(i)/REAL(nsplit_thermals)
420!            write(22,*)'ALE CALLTHERM',ale_bl(i),ale(i)
421            alp_bl(i)=alp_bl(i)+alp(i)/REAL(nsplit_thermals)
422!            write(23,*)'ALP CALLTHERM',alp_bl(i),alp(i)
423        if(prt_level.GE.10) print*,'calltherm i alp_bl alp ale_bl ale',i,alp_bl(i),alp(i),ale_bl(i),ale(i)
[879]424       ENDDO
425
[973]426!IM 060508 marche pas comme cela !!!        enddo ! isplit
427
[878]428!   tests sur les valeurs negatives de l'eau
[973]429         nbptspb=0
[878]430            DO k = 1, klev
431            DO i = 1, klon
[973]432               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
433               if (logexpr2(i,k)) then
434                q_seri(i,k)=1.e-15
435                nbptspb=nbptspb+1
[4089]436#ifdef ISO
[4143]437                do ixt=1,ntiso
[4089]438                  xt_seri(ixt,i,k)=1.e-15*(xtmemoire(ixt,i,k)/qmemoire(i,k))
439                enddo
440#endif
[973]441!                if (prt_level.ge.10) then
442!                  print*,'WARN eau<0 apres therm i=',i,'  k=',k  &
443!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k),  &
444!    &         'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
[938]445                 endif
[878]446            ENDDO
447            ENDDO
[4089]448#ifdef ISO
449#ifdef ISOVERIF
450      if (iso_HDO.gt.0) then
451      call iso_verif_aberrant_enc_vect2D( &
452     &        xt_seri,q_seri, &
[4143]453     &        'calltherm 393, apres bidouille q<0',ntiso,klon,klev)
[4089]454      endif     
455#endif
456#endif
457
[973]458        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
[878]459! tests sur les valeurs de la temperature
[973]460        nbptspb=0
[878]461            DO k = 1, klev
462            DO i = 1, klon
[973]463               logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370.
464               if (logexpr2(i,k)) nbptspb=nbptspb+1
465!              if ((t_seri(i,k).lt.50.) .or.  &
466!    &              (t_seri(i,k).gt.370.)) then
467!                 print*,'WARN temp apres therm i=',i,'  k=',k  &
468!    &         ,' t_seri',t_seri(i,k)
[878]469!              CALL abort
[973]470!              endif
[878]471            ENDDO
472            ENDDO
[973]473        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
[878]474         enddo ! isplit
475
476!
477!***************************************************************
478!     calcul du flux ascencant conservatif
479!            print*,'<<<<calcul flux ascendant conservatif'
480
481      fmc_therm=0.
482               do k=1,klev
483            do i=1,klon
484                  if (entr_therm(i,k).gt.0.) then
485                     fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
486                  else
487                     fmc_therm(i,k+1)=fmc_therm(i,k)
488                  endif
489                  detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1))  &
490     &                 -(fmc_therm(i,k)-fm_therm(i,k))
491               enddo
492            enddo
493     
494     
495!****************************************************************
496!     calcul de l'humidite dans l'ascendance
497!      print*,'<<<<calcul de lhumidite dans thermique'
498!CR:on ne le calcule que pour le cas sec
499      if (iflag_thermals.le.11) then     
500      do i=1,klon
501         zqasc(i,1)=q_seri(i,1)
502         do k=2,klev
503            if (fmc_therm(i,k+1).gt.1.e-6) then
504               zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1)  &
505     &              +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
506!CR:test on asseche le thermique
507!               zqasc(i,k)=zqasc(i,k)/2.
508!            else
509!               zqasc(i,k)=q_seri(i,k)
510            endif
511         enddo
512       enddo
513     
514
515!     calcul de l'eau condensee dans l'ascendance
516!             print*,'<<<<calcul de leau condensee dans thermique'
517             do i=1,klon
518                do k=1,klev
519                   clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
520                   if (clwcon0(i,k).lt.0. .or.   &
521     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
522                      clwcon0(i,k)=0.
523                   endif
524                enddo
525             enddo
526       else
527              do i=1,klon
528                do k=1,klev
529                   clwcon0(i,k)=zqla(i,k) 
530                   if (clwcon0(i,k).lt.0. .or.   &
531     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
532                   clwcon0(i,k)=0.
533                   endif
534                enddo
535             enddo
536       endif
537!*******************************************************************   
538
539
[1428]540!jyg  Protection contre les temperatures nulles
541          do i=1,klon
542             do k=1,klev
543                if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0.
544             enddo
545          enddo
546
547
[878]548      return
549
550      end
Note: See TracBrowser for help on using the repository browser.