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

Last change on this file since 3814 was 3814, checked in by ymipsl, 10 years ago

remove all dynamic dependency in LMDZ physics except for the include "dimensions.h"

YM

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