source: LMDZ6/trunk/libf/phylmd/calltherm_mod.F90 @ 5773

Last change on this file since 5773 was 5765, checked in by rkazeroni, 13 days ago

For GPU porting:

  • Add one Fortran comment for each ported param to specify possible names for horizontal variables used in the param
  • Add Fortran comment to exclude a file and its routines from GPU porting

Note:

  • The "gpum horizontal" comment is placed in the param main. It also applies to dependencies of the main.
  • The "gpum nogpu" excludes all modules and routines of the file from GPU porting. This is used for routines not targeted for GPU computing.
  • Property svn:eol-style set to native
File size: 18.5 KB
Line 
1!
2! $Id: calltherm.F90 5390 2024-12-05 16:09:25Z ymeurdesoif $
3!
4!$gpum horizontal klon ngrid
5MODULE calltherm_mod
6
7CONTAINS
8
9      subroutine calltherm(itap, dtime  &
10     &      ,pplay,paprs,pphi,weak_inversion  &
11     &      ,u_seri_,v_seri_,t_seri_,q_seri_,t_env,q_env,zqsat,debut  &
12     &      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs  &
13     &    ,fm_therm,entr_therm,detr_therm,zqasc,clwcon0,lmax,ratqscth,&
14     &       ratqsdiff,zqsatth,ale_bl,alp_bl,lalim_conv,wght_th, &
15     &       zmax0,f0,zw2,fraca,ztv,zpspsk,ztla,zthl &
16!!! nrlmd le 10/04/2012
17     &      ,pbl_tke,pctsrf,omega,airephy &
18     &      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
19     &      ,n2,s2,strig,zcong,ale_bl_stat &
20     &      ,therm_tke_max,env_tke_max &
21     &      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
22     &      ,alp_bl_conv,alp_bl_stat &
23!!! fin nrlmd le 10/04/2012
24     &      ,zqla,ztva &
25#ifdef ISO         
26     &      ,xt_seri,d_xt_ajs &
27#ifdef DIAGISO         
28     &      ,q_the,xt_the &
29#endif
30#endif         
31     &   )
32
33      USE dimphy
34      USE indice_sol_mod
35      USE print_control_mod, ONLY: prt_level,lunout
36      USE lmdz_thermcell_alp, ONLY: thermcell_alp
37      USE lmdz_thermcell_main, ONLY: thermcell_main
38      USE lmdz_thermcell_old, ONLY: thermcell_old=>thermcell, thermcell_2002, thermcell_eau, calcul_sec, thermcell_sec
39#ifdef ISO
40      use infotrac_phy, ONLY: ntiso
41#ifdef ISOVERIF
42      USE isotopes_mod, ONLY: iso_eau,iso_HDO
43      USE isotopes_verif_mod, ONLY: iso_verif_aberrant_enc_vect2D, &
44        iso_verif_egalite_vect2D
45#endif   
46#endif
47      USE clesphys_mod_h
48      USE thermcell_old_mod_h, ONLY: r_aspect_thermals, l_mix_thermals, w2di_thermals
49      implicit none
50
51      INTEGER, INTENT(IN) :: itap
52!IM 140508
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),PARAMETER :: modname='calltherm'
149      character (len=80) :: abort_message
150
151      integer i,k,isplit
152      logical :: new_thermcell
153
154#ifdef ISO
155      REAL xt_seri(ntiso,klon,klev),xtmemoire(ntiso,klon,klev)
156      REAL d_xt_ajs(ntiso,klon,klev)
157      real d_xt_the(ntiso,klon,klev)
158#ifdef DIAGISO
159      real q_the(klon,klev)
160      real xt_the(ntiso,klon,klev)
161#endif
162      real qprec(klon,klev)
163      integer ixt
164#endif
165
166
167      u_seri(:,:)=u_seri_(:,:)
168      v_seri(:,:)=v_seri_(:,:)
169      t_seri(:,:)=t_seri_(:,:)
170      q_seri(:,:)=q_seri_(:,:)
171
172
173!  Modele du thermique
174!  ===================
175!         print*,'thermiques: WARNING on passe t au lieu de t_seri'
176
177
178! On prend comme valeur initiale des thermiques la valeur du pas
179! de temps precedent
180         zfm_therm(:,:)=fm_therm(:,:)
181         zdetr_therm(:,:)=detr_therm(:,:)
182         zentr_therm(:,:)=entr_therm(:,:)
183
184! On reinitialise les flux de masse a zero pour le cumul en
185! cas de splitting
186         fm_therm(:,:)=0.
187         entr_therm(:,:)=0.
188         detr_therm(:,:)=0.
189
190         ale_bl(:)=0.
191         alp_bl(:)=0.
192         if (prt_level.ge.10) then
193          print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
194         endif
195
196!   tests sur les valeurs negatives de l'eau
197         logexpr0=prt_level.ge.10
198         nbptspb=0
199         do k=1,klev
200            do i=1,klon
201! Attention teste abderr 19-03-09
202!               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
203                logexpr2(i,k)=.not.q_seri(i,k).ge.1.e-15
204               if (logexpr2(i,k)) then
205#ifdef ISO
206                qprec(i,k)=q_seri(i,k)
207#endif
208                q_seri(i,k)=1.e-15
209                nbptspb=nbptspb+1
210#ifdef ISO
211                do ixt=1,ntiso
212                  xt_seri(ixt,i,k)=1.e-15*(xt_seri(ixt,i,k)/qprec(i,k))
213                  ! xt_seri(ixt,i,k)=1.e-15*(Rdefault(index_iso(ixt)))
214                enddo
215#endif
216               endif
217!               if (logexpr0) &
218!    &             print*,'WARN eau<0 avant therm i=',i,'  k=',k  &
219!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k)
220            enddo
221         enddo
222         if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
223
224
225         new_thermcell=iflag_thermals>=15.and.iflag_thermals<=18
226#ifdef ISO
227      if (.not.new_thermcell) then
228           CALL abort_physic('calltherm 234','isos pas prevus ici',1)
229      endif
230#ifdef ISOVERIF
231      if (iso_eau.gt.0) then
232       call iso_verif_egalite_vect2D( &
233     &           xt_seri,q_seri, &
234     &           'calltherm 174',ntiso,klon,klev)
235      endif !if (iso_eau.gt.0) then
236#endif   
237#endif
238         zdt=dtime/REAL(nsplit_thermals)
239
240
241         do isplit=1,nsplit_thermals
242
243          if (iflag_thermals>=1000) then
244            CALL thermcell_2002(klon,klev,zdt,iflag_thermals   &
245     &      ,pplay,paprs,pphi  &
246     &      ,u_seri,v_seri,t_seri,q_seri  &
247     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
248     &      ,zfm_therm,zentr_therm,fraca,zw2  &
249     &      ,r_aspect_thermals,30.,w2di_thermals  &
250     &      ,tau_thermals)
251          else if (iflag_thermals.eq.2) then
252            CALL thermcell_sec(klon,klev,zdt  &
253     &      ,pplay,paprs,pphi,zlev  &
254     &      ,u_seri,v_seri,t_seri,q_seri  &
255     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
256     &      ,zfm_therm,zentr_therm  &
257     &      ,r_aspect_thermals,30.,w2di_thermals  &
258     &      ,tau_thermals)
259          else if (iflag_thermals.eq.3) then
260            CALL thermcell_old(klon,klev,zdt  &
261     &      ,pplay,paprs,pphi  &
262     &      ,u_seri,v_seri,t_seri,q_seri  &
263     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
264     &      ,zfm_therm,zentr_therm  &
265     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
266     &      ,tau_thermals)
267          else if (iflag_thermals.eq.10) then
268            CALL thermcell_eau(klon,klev,zdt  &
269     &      ,pplay,paprs,pphi  &
270     &      ,u_seri,v_seri,t_seri,q_seri  &
271     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
272     &      ,zfm_therm,zentr_therm  &
273     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
274     &      ,tau_thermals)
275          else if (iflag_thermals.eq.11) then
276              abort_message = 'cas non prevu dans calltherm'
277              CALL abort_physic (modname,abort_message,1)
278          else if (iflag_thermals.eq.12) then
279            CALL calcul_sec(klon,klev,zdt  &
280     &      ,pplay,paprs,pphi,zlev  &
281     &      ,u_seri,v_seri,t_seri,q_seri  &
282     &      ,zmax_sec,wmax_sec,zw_sec,lmix_sec  &
283     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
284     &      ,tau_thermals)
285          else if (iflag_thermals==13.or.iflag_thermals==14) then
286              abort_message = 'thermcellV0_main enleve svn>2084'
287              CALL abort_physic (modname,abort_message,1)
288          else if (new_thermcell) then
289            CALL thermcell_main(itap,klon,klev,zdt  &
290     &      ,pplay,paprs,pphi,debut  &
291     &      ,u_seri,v_seri,t_seri,q_seri,t_env,q_env  &
292     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
293     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
294     &      ,ratqscth,ratqsdiff,zqsatth  &
295     &      ,zmax0,f0,zw2,fraca,ztv,zpspsk &
296     &      ,ztla,zthl,ztva &
297     &      ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax,zcong &
298#ifdef ISO         
299     &      ,xt_seri,d_xt_the &
300#endif         
301     &   )
302
303            CALL thermcell_alp(klon,klev,zdt  &                      ! in
304     &        ,pplay,paprs  &                                        ! in
305     &        ,zfm_therm,zentr_therm,lmax  &                         ! in
306     &        ,pbl_tke,pctsrf,omega,airephy &                        ! in
307     &        ,zw2,fraca &                                           ! in
308     &        ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &  ! in
309     &        ,zcong,ale,alp,lalim_conv,wght_th &                          ! out
310     &        ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &! out
311     &        ,n2,s2,strig,ale_bl_stat &                                   ! out
312     &        ,therm_tke_max,env_tke_max &                           ! out
313     &        ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &          ! out
314     &        ,alp_bl_conv,alp_bl_stat &                             ! out
315     &        )
316
317           if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
318         else
319           abort_message = 'Cas des thermiques non prevu'
320           CALL abort_physic (modname,abort_message,1)
321         endif
322
323! Attention : les noms sont contre intuitif.
324! flag_bidouille_stratocu est .true. si on ne fait pas de bidouille.
325! Il aurait mieux valu avoir un nobidouille_stratocu
326! Et pour simplifier :
327! nobidouille_stratocu=.not.(iflag_thermals==13.or.iflag_thermals=15)
328! Ce serait bien de changer, mai en prenant le temps de vérifier que ca
329! fait bien ce qu'on croit.
330
331       flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18
332
333! Calcul a posteriori du niveau max des thermiques pour les schémas qui
334! ne la sortent pas.
335      if (iflag_thermals<=12.or.iflag_thermals>=1000) then
336         lmax(:)=1
337         do k=1,klev-1
338            zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1)
339         enddo
340         do k=1,klev-1
341            do i=1,klon
342               if (zfm_therm(i,k+1)>0.) lmax(i)=k
343            enddo
344         enddo
345      endif
346
347      fact(:)=0.
348      DO i=1,klon
349       logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
350       IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals)
351      ENDDO
352
353     DO k=1,klev
354!  transformation de la derivee en tendance
355            d_t_the(:,k)=d_t_the(:,k)*dtime*fact(:)
356            d_u_the(:,k)=d_u_the(:,k)*dtime*fact(:)
357            d_v_the(:,k)=d_v_the(:,k)*dtime*fact(:)
358            d_q_the(:,k)=d_q_the(:,k)*dtime*fact(:)
359            fm_therm(:,k)=fm_therm(:,k)  &
360     &      +zfm_therm(:,k)*fact(:)
361            entr_therm(:,k)=entr_therm(:,k)  &
362     &       +zentr_therm(:,k)*fact(:)
363            detr_therm(:,k)=detr_therm(:,k)  &
364     &       +zdetr_therm(:,k)*fact(:)
365#ifdef ISO
366            do ixt=1,ntiso
367              d_xt_the(ixt,:,k)=d_xt_the(ixt,:,k)*dtime*fact(:)
368            enddo
369#endif
370      ENDDO
371       fm_therm(:,klev+1)=0.
372
373
374
375!  accumulation de la tendance
376            d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
377            d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
378            d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
379            d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
380#ifdef ISO
381            d_xt_ajs(:,:,:)=d_xt_ajs(:,:,:)+d_xt_the(:,:,:)
382#endif
383
384!  incrementation des variables meteo
385            t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
386            u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
387            v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
388            qmemoire(:,:)=q_seri(:,:)
389            q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
390#ifdef ISO
391            xtmemoire(:,:,:)=xt_seri(:,:,:)
392            xt_seri(:,:,:) = xt_seri(:,:,:) + d_xt_the(:,:,:)
393#ifdef ISOVERIF
394!      write(*,*) 'calltherm 350 tmp: ajout d_xt_the'
395      if (iso_HDO.gt.0) then
396!      i=479
397!      k=4
398!      write(*,*) 'xt_seri(iso_hdo,i,k),q_seri(i,k)=', &
399!     &   xt_seri(iso_hdo,i,k),q_seri(i,k)
400!      write(*,*) 'd_xt_the(iso_hdo,i,k),d_q_the(i,k)=', &
401!     &   d_xt_the(iso_hdo,i,k),d_q_the(i,k)
402      call iso_verif_aberrant_enc_vect2D( &
403     &        xt_seri,q_seri, &
404     &        'calltherm 353, apres ajout d_xt_the',ntiso,klon,klev)
405      endif     
406#endif
407#endif
408           if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
409
410       DO i=1,klon
411            fm_therm(i,klev+1)=0.
412            ale_bl(i)=ale_bl(i)+ale(i)/REAL(nsplit_thermals)
413!            write(22,*)'ALE CALLTHERM',ale_bl(i),ale(i)
414            alp_bl(i)=alp_bl(i)+alp(i)/REAL(nsplit_thermals)
415!            write(23,*)'ALP CALLTHERM',alp_bl(i),alp(i)
416        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)
417       ENDDO
418
419!IM 060508 marche pas comme cela !!!        enddo ! isplit
420
421!   tests sur les valeurs negatives de l'eau
422         nbptspb=0
423            DO k = 1, klev
424            DO i = 1, klon
425               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
426               if (logexpr2(i,k)) then
427                q_seri(i,k)=1.e-15
428                nbptspb=nbptspb+1
429#ifdef ISO
430                do ixt=1,ntiso
431                  xt_seri(ixt,i,k)=1.e-15*(xtmemoire(ixt,i,k)/qmemoire(i,k))
432                enddo
433#endif
434!                if (prt_level.ge.10) then
435!                  print*,'WARN eau<0 apres therm i=',i,'  k=',k  &
436!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k),  &
437!    &         'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
438                 endif
439            ENDDO
440            ENDDO
441#ifdef ISO
442#ifdef ISOVERIF
443      if (iso_HDO.gt.0) then
444      call iso_verif_aberrant_enc_vect2D( &
445     &        xt_seri,q_seri, &
446     &        'calltherm 393, apres bidouille q<0',ntiso,klon,klev)
447      endif     
448#endif
449#endif
450
451        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
452! tests sur les valeurs de la temperature
453        nbptspb=0
454            DO k = 1, klev
455            DO i = 1, klon
456               logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370.
457               if (logexpr2(i,k)) nbptspb=nbptspb+1
458!              if ((t_seri(i,k).lt.50.) .or.  &
459!    &              (t_seri(i,k).gt.370.)) then
460!                 print*,'WARN temp apres therm i=',i,'  k=',k  &
461!    &         ,' t_seri',t_seri(i,k)
462!              CALL abort
463!              endif
464            ENDDO
465            ENDDO
466        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
467         enddo ! isplit
468
469!
470!***************************************************************
471!     calcul du flux ascencant conservatif
472!            print*,'<<<<calcul flux ascendant conservatif'
473
474      fmc_therm=0.
475               do k=1,klev
476            do i=1,klon
477                  if (entr_therm(i,k).gt.0.) then
478                     fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
479                  else
480                     fmc_therm(i,k+1)=fmc_therm(i,k)
481                  endif
482                  detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1))  &
483     &                 -(fmc_therm(i,k)-fm_therm(i,k))
484               enddo
485            enddo
486     
487     
488!****************************************************************
489!     calcul de l'humidite dans l'ascendance
490!      print*,'<<<<calcul de lhumidite dans thermique'
491!CR:on ne le calcule que pour le cas sec
492      if (iflag_thermals.le.11) then     
493      do i=1,klon
494         zqasc(i,1)=q_seri(i,1)
495         do k=2,klev
496            if (fmc_therm(i,k+1).gt.1.e-6) then
497               zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1)  &
498     &              +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
499!CR:test on asseche le thermique
500!               zqasc(i,k)=zqasc(i,k)/2.
501!            else
502!               zqasc(i,k)=q_seri(i,k)
503            endif
504         enddo
505       enddo
506     
507
508!     calcul de l'eau condensee dans l'ascendance
509!             print*,'<<<<calcul de leau condensee dans thermique'
510             do i=1,klon
511                do k=1,klev
512                   clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
513                   if (clwcon0(i,k).lt.0. .or.   &
514     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
515                      clwcon0(i,k)=0.
516                   endif
517                enddo
518             enddo
519       else
520              do i=1,klon
521                do k=1,klev
522                   clwcon0(i,k)=zqla(i,k) 
523                   if (clwcon0(i,k).lt.0. .or.   &
524     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
525                   clwcon0(i,k)=0.
526                   endif
527                enddo
528             enddo
529       endif
530!*******************************************************************   
531
532
533!jyg  Protection contre les temperatures nulles
534          do i=1,klon
535             do k=1,klev
536                if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0.
537             enddo
538          enddo
539
540
541      return
542
543      end subroutine calltherm
544
545END MODULE calltherm_mod
Note: See TracBrowser for help on using the repository browser.