source: LMDZ6/branches/Amaury_dev/libf/phylmd/calltherm.F90 @ 5441

Last change on this file since 5441 was 5231, checked in by abarral, 3 months ago

Merge r5217

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