source: LMDZ5/trunk/libf/phylmd/calltherm.F90 @ 2479

Last change on this file since 2479 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

EM

  • 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: 15.3 KB
Line 
1!
2! $Id: calltherm.F90 2346 2015-08-21 15:13:46Z 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 print_control_mod, ONLY: prt_level,lunout
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.