source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/calltherm.F90 @ 3820

Last change on this file since 3820 was 3818, checked in by millour, 10 years ago

Some partial cleanup on uses of "dimensions.h" in physics.
At this point 3D gcm compiles and bench seems to run fine :-)
EM

File size: 15.3 KB
Line 
1!
2! $Id: calltherm.F90 1943 2014-01-22 09:51:36Z idelkadi $
3!
4      subroutine calltherm(dtime  &
5     &      ,pplay,paprs,pphi,weak_inversion  &
6     &      ,u_seri,v_seri,t_seri,q_seri,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,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
21      USE dimphy
22      USE indice_sol_mod
23      USE inifis_mod, ONLY: lunout, prt_level
24
25      implicit none
26#include "thermcell.h"
27
28
29!IM 140508
30      INTEGER, SAVE ::  itap
31!$OMP THREADPRIVATE(itap)
32      REAL dtime
33      LOGICAL debut
34      LOGICAL logexpr0, logexpr2(klon,klev), logexpr1(klon)
35      REAL fact(klon)
36      INTEGER nbptspb
37
38      REAL u_seri(klon,klev),v_seri(klon,klev)
39      REAL t_seri(klon,klev),q_seri(klon,klev),qmemoire(klon,klev)
40      REAL weak_inversion(klon)
41      REAL paprs(klon,klev+1)
42      REAL pplay(klon,klev)
43      REAL pphi(klon,klev)
44      real zlev(klon,klev+1)
45!test: on sort lentr et a* pour alimenter KE
46      REAL wght_th(klon,klev)
47      INTEGER lalim_conv(klon)
48      REAL zw2(klon,klev+1),fraca(klon,klev+1)
49
50!FH Update Thermiques
51      REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
52      REAL d_u_ajs(klon,klev),d_v_ajs(klon,klev)
53      real fm_therm(klon,klev+1)
54      real entr_therm(klon,klev),detr_therm(klon,klev)
55
56!********************************************************
57!     declarations
58      LOGICAL flag_bidouille_stratocu
59      real fmc_therm(klon,klev+1),zqasc(klon,klev)
60      real zqla(klon,klev)
61      real zqta(klon,klev)
62      real ztv(klon,klev),ztva(klon,klev)
63      real zpspsk(klon,klev)
64      real ztla(klon,klev)
65      real zthl(klon,klev)
66      real wmax_sec(klon)
67      real zmax_sec(klon)
68      real f_sec(klon)
69      real detrc_therm(klon,klev)
70! FH WARNING : il semble que ces save ne servent a rien
71!     save fmc_therm, detrc_therm
72      real clwcon0(klon,klev)
73      real zqsat(klon,klev)
74      real zw_sec(klon,klev+1)
75      integer lmix_sec(klon)
76      integer lmax(klon)
77      real ratqscth(klon,klev)
78      real ratqsdiff(klon,klev)
79      real zqsatth(klon,klev) 
80!nouvelles variables pour la convection
81      real Ale_bl(klon)
82      real Alp_bl(klon)
83      real Ale(klon)
84      real Alp(klon)
85!RC
86      !on garde le zmax du pas de temps precedent
87      real zmax0(klon), f0(klon)
88
89!!! nrlmd le 10/04/2012
90      real pbl_tke(klon,klev+1,nbsrf)
91      real pctsrf(klon,nbsrf)
92      real omega(klon,klev)
93      real airephy(klon)
94      real zlcl_th(klon),fraca0(klon),w0(klon),w_conv(klon)
95      real therm_tke_max0(klon),env_tke_max0(klon)
96      real n2(klon),s2(klon)
97      real ale_bl_stat(klon)
98      real therm_tke_max(klon,klev),env_tke_max(klon,klev)
99      real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
100!!! fin nrlmd le 10/04/2012
101
102!********************************************************
103
104
105! variables locales
106      REAL d_t_the(klon,klev), d_q_the(klon,klev)
107      REAL d_u_the(klon,klev),d_v_the(klon,klev)
108!
109      real zfm_therm(klon,klev+1),zdt
110      real zentr_therm(klon,klev),zdetr_therm(klon,klev)
111! FH A VERIFIER : SAVE INUTILES
112!      save zentr_therm,zfm_therm
113
114      character (len=20) :: modname='calltherm'
115      character (len=80) :: abort_message
116
117      integer i,k
118      logical, save :: first=.true.
119!$OMP THREADPRIVATE(first)
120!********************************************************
121      if (first) then
122        itap=0
123        first=.false.
124      endif
125
126! Incrementer le compteur de la physique
127     itap   = itap + 1
128
129!  Modele du thermique
130!  ===================
131!         print*,'thermiques: WARNING on passe t au lieu de t_seri'
132
133
134! On prend comme valeur initiale des thermiques la valeur du pas
135! de temps precedent
136         zfm_therm(:,:)=fm_therm(:,:)
137         zdetr_therm(:,:)=detr_therm(:,:)
138         zentr_therm(:,:)=entr_therm(:,:)
139
140! On reinitialise les flux de masse a zero pour le cumul en
141! cas de splitting
142         fm_therm(:,:)=0.
143         entr_therm(:,:)=0.
144         detr_therm(:,:)=0.
145
146         Ale_bl(:)=0.
147         Alp_bl(:)=0.
148         if (prt_level.ge.10) then
149          print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
150         endif
151
152!   tests sur les valeurs negatives de l'eau
153         logexpr0=prt_level.ge.10
154         nbptspb=0
155         do k=1,klev
156            do i=1,klon
157! Attention teste abderr 19-03-09
158!               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
159                logexpr2(i,k)=.not.q_seri(i,k).ge.1.e-15
160               if (logexpr2(i,k)) then
161                q_seri(i,k)=1.e-15
162                nbptspb=nbptspb+1
163               endif
164!               if (logexpr0) &
165!    &             print*,'WARN eau<0 avant therm i=',i,'  k=',k  &
166!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k)
167            enddo
168         enddo
169         if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
170
171         zdt=dtime/REAL(nsplit_thermals)
172         do isplit=1,nsplit_thermals
173
174          if (iflag_thermals>=1000) then
175            CALL thermcell_2002(klon,klev,zdt,iflag_thermals   &
176     &      ,pplay,paprs,pphi  &
177     &      ,u_seri,v_seri,t_seri,q_seri  &
178     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
179     &      ,zfm_therm,zentr_therm,fraca,zw2  &
180     &      ,r_aspect_thermals,30.,w2di_thermals  &
181     &      ,tau_thermals)
182          else if (iflag_thermals.eq.2) then
183            CALL thermcell_sec(klon,klev,zdt  &
184     &      ,pplay,paprs,pphi,zlev  &
185     &      ,u_seri,v_seri,t_seri,q_seri  &
186     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
187     &      ,zfm_therm,zentr_therm  &
188     &      ,r_aspect_thermals,30.,w2di_thermals  &
189     &      ,tau_thermals)
190          else if (iflag_thermals.eq.3) then
191            CALL thermcell(klon,klev,zdt  &
192     &      ,pplay,paprs,pphi  &
193     &      ,u_seri,v_seri,t_seri,q_seri  &
194     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
195     &      ,zfm_therm,zentr_therm  &
196     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
197     &      ,tau_thermals)
198          else if (iflag_thermals.eq.10) then
199            CALL thermcell_eau(klon,klev,zdt  &
200     &      ,pplay,paprs,pphi  &
201     &      ,u_seri,v_seri,t_seri,q_seri  &
202     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
203     &      ,zfm_therm,zentr_therm  &
204     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
205     &      ,tau_thermals)
206          else if (iflag_thermals.eq.11) then
207              abort_message = 'cas non prevu dans calltherm'
208              CALL abort_physic (modname,abort_message,1)
209
210!           CALL thermcell_pluie(klon,klev,zdt  &
211!   &      ,pplay,paprs,pphi,zlev  &
212!    &      ,u_seri,v_seri,t_seri,q_seri  &
213!    &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
214!    &      ,zfm_therm,zentr_therm,zqla  &
215!    &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
216!    &      ,tau_thermals,3)
217          else if (iflag_thermals.eq.12) then
218            CALL calcul_sec(klon,klev,zdt  &
219     &      ,pplay,paprs,pphi,zlev  &
220     &      ,u_seri,v_seri,t_seri,q_seri  &
221     &      ,zmax_sec,wmax_sec,zw_sec,lmix_sec  &
222     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
223     &      ,tau_thermals)
224          else if (iflag_thermals==13.or.iflag_thermals==14) then
225            CALL thermcellV0_main(itap,klon,klev,zdt  &
226     &      ,pplay,paprs,pphi,debut  &
227     &      ,u_seri,v_seri,t_seri,q_seri  &
228     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
229     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
230     &      ,ratqscth,ratqsdiff,zqsatth  &
231     &      ,r_aspect_thermals,l_mix_thermals  &
232     &      ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
233     &      ,zmax0,f0,zw2,fraca)
234          else if (iflag_thermals>=15.and.iflag_thermals<=18) then
235
236!            print*,'THERM iflag_thermas_ed=',iflag_thermals_ed
237            CALL thermcell_main(itap,klon,klev,zdt  &
238     &      ,pplay,paprs,pphi,debut  &
239     &      ,u_seri,v_seri,t_seri,q_seri  &
240     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
241     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
242     &      ,ratqscth,ratqsdiff,zqsatth  &
243!    &      ,r_aspect_thermals,l_mix_thermals &
244!    &      ,tau_thermals,iflag_thermals_ed,iflag_coupl &
245     &      ,Ale,Alp,lalim_conv,wght_th &
246     &      ,zmax0,f0,zw2,fraca,ztv,zpspsk &
247     &      ,ztla,zthl &
248!!! nrlmd le 10/04/2012
249     &      ,pbl_tke,pctsrf,omega,airephy &
250     &      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
251     &      ,n2,s2,ale_bl_stat &
252     &      ,therm_tke_max,env_tke_max &
253     &      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
254     &      ,alp_bl_conv,alp_bl_stat &
255!!! fin nrlmd le 10/04/2012
256     &      ,ztva )
257           if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
258         else
259           abort_message = 'Cas des thermiques non prevu'
260           CALL abort_physic (modname,abort_message,1)
261         endif
262
263! Attention : les noms sont contre intuitif.
264! flag_bidouille_stratocu est .true. si on ne fait pas de bidouille.
265! Il aurait mieux valu avoir un nobidouille_stratocu
266! Et pour simplifier :
267! nobidouille_stratocu=.not.(iflag_thermals==13.or.iflag_thermals=15)
268! Ce serait bien de changer, mai en prenant le temps de vérifier que ca
269! fait bien ce qu'on croit.
270
271       flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18
272
273! Calcul a posteriori du niveau max des thermiques pour les schémas qui
274! ne la sortent pas.
275      if (iflag_thermals<=12.or.iflag_thermals>=1000) then
276         lmax(:)=1
277         do k=1,klev-1
278            zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1)
279         enddo
280         do k=1,klev-1
281            do i=1,klon
282               if (zfm_therm(i,k+1)>0.) lmax(i)=k
283            enddo
284         enddo
285      endif
286
287      fact(:)=0.
288      DO i=1,klon
289       logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
290       IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals)
291      ENDDO
292
293     DO k=1,klev
294!  transformation de la derivee en tendance
295            d_t_the(:,k)=d_t_the(:,k)*dtime*fact(:)
296            d_u_the(:,k)=d_u_the(:,k)*dtime*fact(:)
297            d_v_the(:,k)=d_v_the(:,k)*dtime*fact(:)
298            d_q_the(:,k)=d_q_the(:,k)*dtime*fact(:)
299            fm_therm(:,k)=fm_therm(:,k)  &
300     &      +zfm_therm(:,k)*fact(:)
301            entr_therm(:,k)=entr_therm(:,k)  &
302     &       +zentr_therm(:,k)*fact(:)
303            detr_therm(:,k)=detr_therm(:,k)  &
304     &       +zdetr_therm(:,k)*fact(:)
305      ENDDO
306       fm_therm(:,klev+1)=0.
307
308
309
310!  accumulation de la tendance
311            d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
312            d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
313            d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
314            d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
315
316!  incrementation des variables meteo
317            t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
318            u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
319            v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
320            qmemoire(:,:)=q_seri(:,:)
321            q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
322           if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
323
324       DO i=1,klon
325            fm_therm(i,klev+1)=0.
326            Ale_bl(i)=Ale_bl(i)+Ale(i)/REAL(nsplit_thermals)
327!            write(22,*)'ALE CALLTHERM',Ale_bl(i),Ale(i)
328            Alp_bl(i)=Alp_bl(i)+Alp(i)/REAL(nsplit_thermals)
329!            write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
330        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)
331       ENDDO
332
333!IM 060508 marche pas comme cela !!!        enddo ! isplit
334
335!   tests sur les valeurs negatives de l'eau
336         nbptspb=0
337            DO k = 1, klev
338            DO i = 1, klon
339               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
340               if (logexpr2(i,k)) then
341                q_seri(i,k)=1.e-15
342                nbptspb=nbptspb+1
343!                if (prt_level.ge.10) then
344!                  print*,'WARN eau<0 apres therm i=',i,'  k=',k  &
345!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k),  &
346!    &         'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
347                 endif
348            ENDDO
349            ENDDO
350        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
351! tests sur les valeurs de la temperature
352        nbptspb=0
353            DO k = 1, klev
354            DO i = 1, klon
355               logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370.
356               if (logexpr2(i,k)) nbptspb=nbptspb+1
357!              if ((t_seri(i,k).lt.50.) .or.  &
358!    &              (t_seri(i,k).gt.370.)) then
359!                 print*,'WARN temp apres therm i=',i,'  k=',k  &
360!    &         ,' t_seri',t_seri(i,k)
361!              CALL abort
362!              endif
363            ENDDO
364            ENDDO
365        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
366         enddo ! isplit
367
368!
369!***************************************************************
370!     calcul du flux ascencant conservatif
371!            print*,'<<<<calcul flux ascendant conservatif'
372
373      fmc_therm=0.
374               do k=1,klev
375            do i=1,klon
376                  if (entr_therm(i,k).gt.0.) then
377                     fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
378                  else
379                     fmc_therm(i,k+1)=fmc_therm(i,k)
380                  endif
381                  detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1))  &
382     &                 -(fmc_therm(i,k)-fm_therm(i,k))
383               enddo
384            enddo
385     
386     
387!****************************************************************
388!     calcul de l'humidite dans l'ascendance
389!      print*,'<<<<calcul de lhumidite dans thermique'
390!CR:on ne le calcule que pour le cas sec
391      if (iflag_thermals.le.11) then     
392      do i=1,klon
393         zqasc(i,1)=q_seri(i,1)
394         do k=2,klev
395            if (fmc_therm(i,k+1).gt.1.e-6) then
396               zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1)  &
397     &              +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
398!CR:test on asseche le thermique
399!               zqasc(i,k)=zqasc(i,k)/2.
400!            else
401!               zqasc(i,k)=q_seri(i,k)
402            endif
403         enddo
404       enddo
405     
406
407!     calcul de l'eau condensee dans l'ascendance
408!             print*,'<<<<calcul de leau condensee dans thermique'
409             do i=1,klon
410                do k=1,klev
411                   clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
412                   if (clwcon0(i,k).lt.0. .or.   &
413     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
414                      clwcon0(i,k)=0.
415                   endif
416                enddo
417             enddo
418       else
419              do i=1,klon
420                do k=1,klev
421                   clwcon0(i,k)=zqla(i,k) 
422                   if (clwcon0(i,k).lt.0. .or.   &
423     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
424                   clwcon0(i,k)=0.
425                   endif
426                enddo
427             enddo
428       endif
429!*******************************************************************   
430
431
432!jyg  Protection contre les temperatures nulles
433          do i=1,klon
434             do k=1,klev
435                if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0.
436             enddo
437          enddo
438
439
440      return
441
442      end
Note: See TracBrowser for help on using the repository browser.