source: lmdz_wrf/trunk/WRFV3/lmdz/calltherm.F90

Last change on this file was 179, checked in by lfita, 10 years ago

No Ale, Alp, ... values because pphi is at the intermediate level!!!
Removing raw prints of checkig for Ale, Alp

File size: 15.3 KB
Line 
1!
2! $Id: calltherm.F90 1795 2013-07-18 08:20:28Z emillour $
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! L. Fita, LMD July 2014. Initializing tendencies and variables (should
129!  be done?)
130!      d_t_the = 0.
131!      d_q_the = 0.
132!      d_u_the = 0.
133!      d_v_the = 0.
134!      Alp = 0.
135!      Ale = 0.
136
137! Incrementer le compteur de la physique
138     itap   = itap + 1
139
140!  Modele du thermique
141!  ===================
142!         print*,'thermiques: WARNING on passe t au lieu de t_seri'
143
144! On prend comme valeur initiale des thermiques la valeur du pas
145! de temps precedent
146         zfm_therm(:,:)=fm_therm(:,:)
147         zdetr_therm(:,:)=detr_therm(:,:)
148         zentr_therm(:,:)=entr_therm(:,:)
149
150! On reinitialise les flux de masse a zero pour le cumul en
151! cas de splitting
152         fm_therm(:,:)=0.
153         entr_therm(:,:)=0.
154         detr_therm(:,:)=0.
155
156         Ale_bl(:)=0.
157         Alp_bl(:)=0.
158         if (prt_level.ge.10) then
159          print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
160         endif
161
162!   tests sur les valeurs negatives de l'eau
163         logexpr0=prt_level.ge.10
164         nbptspb=0
165         do k=1,klev
166            do i=1,klon
167! Attention teste abderr 19-03-09
168!               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
169                logexpr2(i,k)=.not.q_seri(i,k).ge.1.e-15
170               if (logexpr2(i,k)) then
171                q_seri(i,k)=1.e-15
172                nbptspb=nbptspb+1
173               endif
174!               if (logexpr0) &
175!    &             print*,'WARN eau<0 avant therm i=',i,'  k=',k  &
176!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k)
177            enddo
178         enddo
179         if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
180
181         zdt=dtime/REAL(nsplit_thermals)
182         do isplit=1,nsplit_thermals
183
184          if (iflag_thermals.eq.1) then
185            CALL thermcell_2002(klon,klev,zdt   &
186     &      ,pplay,paprs,pphi  &
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.2) then
193            CALL thermcell_sec(klon,klev,zdt  &
194     &      ,pplay,paprs,pphi,zlev  &
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,30.,w2di_thermals  &
199     &      ,tau_thermals)
200          else if (iflag_thermals.eq.3) then
201            CALL thermcell(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.10) then
209            CALL thermcell_eau(klon,klev,zdt  &
210     &      ,pplay,paprs,pphi  &
211     &      ,u_seri,v_seri,t_seri,q_seri  &
212     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
213     &      ,zfm_therm,zentr_therm  &
214     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
215     &      ,tau_thermals)
216          else if (iflag_thermals.eq.11) then
217              abort_message = 'cas non prevu dans calltherm'
218              CALL abort_gcm (modname,abort_message,1)
219
220!           CALL thermcell_pluie(klon,klev,zdt  &
221!   &      ,pplay,paprs,pphi,zlev  &
222!    &      ,u_seri,v_seri,t_seri,q_seri  &
223!    &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
224!    &      ,zfm_therm,zentr_therm,zqla  &
225!    &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
226!    &      ,tau_thermals,3)
227          else if (iflag_thermals.eq.12) then
228            CALL calcul_sec(klon,klev,zdt  &
229     &      ,pplay,paprs,pphi,zlev  &
230     &      ,u_seri,v_seri,t_seri,q_seri  &
231     &      ,zmax_sec,wmax_sec,zw_sec,lmix_sec  &
232     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
233     &      ,tau_thermals)
234          else if (iflag_thermals==13.or.iflag_thermals==14) then
235            CALL thermcellV0_main(itap,klon,klev,zdt  &
236     &      ,pplay,paprs,pphi,debut  &
237     &      ,u_seri,v_seri,t_seri,q_seri  &
238     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
239     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
240     &      ,ratqscth,ratqsdiff,zqsatth  &
241     &      ,r_aspect_thermals,l_mix_thermals  &
242     &      ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
243     &      ,zmax0,f0,zw2,fraca)
244          else if (iflag_thermals>=15.and.iflag_thermals<=18) then
245
246!            print*,'THERM iflag_thermas_ed=',iflag_thermals_ed
247
248            CALL thermcell_main(itap,klon,klev,zdt  &
249     &      ,pplay,paprs,pphi,debut  &
250     &      ,u_seri,v_seri,t_seri,q_seri  &
251     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
252     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
253     &      ,ratqscth,ratqsdiff,zqsatth  &
254!    &      ,r_aspect_thermals,l_mix_thermals &
255!    &      ,tau_thermals,iflag_thermals_ed,iflag_coupl &
256     &      ,Ale,Alp,lalim_conv,wght_th &
257     &      ,zmax0,f0,zw2,fraca,ztv,zpspsk &
258     &      ,ztla,zthl &
259!!! nrlmd le 10/04/2012
260     &      ,pbl_tke,pctsrf,omega,airephy &
261     &      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
262     &      ,n2,s2,ale_bl_stat &
263     &      ,therm_tke_max,env_tke_max &
264     &      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
265     &      ,alp_bl_conv,alp_bl_stat &
266!!! fin nrlmd le 10/04/2012
267     &      ,ztva )
268           if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
269         else
270           abort_message = 'Cas des thermiques non prevu'
271           CALL abort_gcm (modname,abort_message,1)
272         endif
273
274! Attention : les noms sont contre intuitif.
275! flag_bidouille_stratocu est .true. si on ne fait pas de bidouille.
276! Il aurait mieux valu avoir un nobidouille_stratocu
277! Et pour simplifier :
278! nobidouille_stratocu=.not.(iflag_thermals==13.or.iflag_thermals=15)
279! Ce serait bien de changer, mai en prenant le temps de vérifier que ca
280! fait bien ce qu'on croit.
281       flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18
282
283      if (iflag_thermals<=12) then
284         lmax=1
285         do k=1,klev-1
286            zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1)
287         enddo
288      endif
289
290      fact(:)=0.
291      DO i=1,klon
292       logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
293       IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals)
294      ENDDO
295
296     DO k=1,klev
297!  transformation de la derivee en tendance
298            d_t_the(:,k)=d_t_the(:,k)*dtime*fact(:)
299            d_u_the(:,k)=d_u_the(:,k)*dtime*fact(:)
300            d_v_the(:,k)=d_v_the(:,k)*dtime*fact(:)
301            d_q_the(:,k)=d_q_the(:,k)*dtime*fact(:)
302            fm_therm(:,k)=fm_therm(:,k)  &
303     &      +zfm_therm(:,k)*fact(:)
304            entr_therm(:,k)=entr_therm(:,k)  &
305     &       +zentr_therm(:,k)*fact(:)
306            detr_therm(:,k)=detr_therm(:,k)  &
307     &       +zdetr_therm(:,k)*fact(:)
308      ENDDO
309       fm_therm(:,klev+1)=0.
310
311!  accumulation de la tendance
312            d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
313            d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
314            d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
315            d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
316
317!  incrementation des variables meteo
318            t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
319            u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
320            v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
321            qmemoire(:,:)=q_seri(:,:)
322            q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
323           if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
324
325       DO i=1,klon
326            fm_therm(i,klev+1)=0.
327            Ale_bl(i)=Ale_bl(i)+Ale(i)/REAL(nsplit_thermals)
328!            write(22,*)'ALE CALLTHERM',Ale_bl(i),Ale(i)
329            Alp_bl(i)=Alp_bl(i)+Alp(i)/REAL(nsplit_thermals)
330!            write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
331        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)
332       ENDDO
333
334!IM 060508 marche pas comme cela !!!        enddo ! isplit
335
336!   tests sur les valeurs negatives de l'eau
337         nbptspb=0
338            DO k = 1, klev
339            DO i = 1, klon
340               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
341               if (logexpr2(i,k)) then
342                q_seri(i,k)=1.e-15
343                nbptspb=nbptspb+1
344!                if (prt_level.ge.10) then
345!                  print*,'WARN eau<0 apres therm i=',i,'  k=',k  &
346!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k),  &
347!    &         'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
348                 endif
349            ENDDO
350            ENDDO
351        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
352! tests sur les valeurs de la temperature
353        nbptspb=0
354            DO k = 1, klev
355            DO i = 1, klon
356               logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370.
357               if (logexpr2(i,k)) nbptspb=nbptspb+1
358!              if ((t_seri(i,k).lt.50.) .or.  &
359!    &              (t_seri(i,k).gt.370.)) then
360!                 print*,'WARN temp apres therm i=',i,'  k=',k  &
361!    &         ,' t_seri',t_seri(i,k)
362!              CALL abort
363!              endif
364            ENDDO
365            ENDDO
366        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
367         enddo ! isplit
368
369!
370!***************************************************************
371!     calcul du flux ascencant conservatif
372!            print*,'<<<<calcul flux ascendant conservatif'
373
374      fmc_therm=0.
375               do k=1,klev
376            do i=1,klon
377                  if (entr_therm(i,k).gt.0.) then
378                     fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
379                  else
380                     fmc_therm(i,k+1)=fmc_therm(i,k)
381                  endif
382                  detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1))  &
383     &                 -(fmc_therm(i,k)-fm_therm(i,k))
384               enddo
385            enddo
386     
387     
388!****************************************************************
389!     calcul de l'humidite dans l'ascendance
390!      print*,'<<<<calcul de lhumidite dans thermique'
391!CR:on ne le calcule que pour le cas sec
392      if (iflag_thermals.le.11) then     
393      do i=1,klon
394         zqasc(i,1)=q_seri(i,1)
395         do k=2,klev
396            if (fmc_therm(i,k+1).gt.1.e-6) then
397               zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1)  &
398     &              +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
399!CR:test on asseche le thermique
400!               zqasc(i,k)=zqasc(i,k)/2.
401!            else
402!               zqasc(i,k)=q_seri(i,k)
403            endif
404         enddo
405       enddo
406     
407
408!     calcul de l'eau condensee dans l'ascendance
409!             print*,'<<<<calcul de leau condensee dans thermique'
410             do i=1,klon
411                do k=1,klev
412                   clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
413                   if (clwcon0(i,k).lt.0. .or.   &
414     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
415                      clwcon0(i,k)=0.
416                   endif
417                enddo
418             enddo
419       else
420              do i=1,klon
421                do k=1,klev
422                   clwcon0(i,k)=zqla(i,k) 
423                   if (clwcon0(i,k).lt.0. .or.   &
424     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
425                   clwcon0(i,k)=0.
426                   endif
427                enddo
428             enddo
429       endif
430!*******************************************************************   
431
432
433!jyg  Protection contre les temperatures nulles
434          do i=1,klon
435             do k=1,klev
436                if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0.
437             enddo
438          enddo
439
440
441      return
442
443      end
Note: See TracBrowser for help on using the repository browser.