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

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.
  • 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.0 KB
Line 
1!
2! $Id: calltherm.F90 1907 2013-11-26 13:10:46Z lguez $
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.eq.1) then
177            CALL thermcell_2002(klon,klev,zdt   &
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  &
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_gcm (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_gcm (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      if (iflag_thermals<=12) 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      endif
281
282      fact(:)=0.
283      DO i=1,klon
284       logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
285       IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals)
286      ENDDO
287
288     DO k=1,klev
289!  transformation de la derivee en tendance
290            d_t_the(:,k)=d_t_the(:,k)*dtime*fact(:)
291            d_u_the(:,k)=d_u_the(:,k)*dtime*fact(:)
292            d_v_the(:,k)=d_v_the(:,k)*dtime*fact(:)
293            d_q_the(:,k)=d_q_the(:,k)*dtime*fact(:)
294            fm_therm(:,k)=fm_therm(:,k)  &
295     &      +zfm_therm(:,k)*fact(:)
296            entr_therm(:,k)=entr_therm(:,k)  &
297     &       +zentr_therm(:,k)*fact(:)
298            detr_therm(:,k)=detr_therm(:,k)  &
299     &       +zdetr_therm(:,k)*fact(:)
300      ENDDO
301       fm_therm(:,klev+1)=0.
302
303
304
305!  accumulation de la tendance
306            d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
307            d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
308            d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
309            d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
310
311!  incrementation des variables meteo
312            t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
313            u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
314            v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
315            qmemoire(:,:)=q_seri(:,:)
316            q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
317           if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
318
319       DO i=1,klon
320            fm_therm(i,klev+1)=0.
321            Ale_bl(i)=Ale_bl(i)+Ale(i)/REAL(nsplit_thermals)
322!            write(22,*)'ALE CALLTHERM',Ale_bl(i),Ale(i)
323            Alp_bl(i)=Alp_bl(i)+Alp(i)/REAL(nsplit_thermals)
324!            write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
325        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)
326       ENDDO
327
328!IM 060508 marche pas comme cela !!!        enddo ! isplit
329
330!   tests sur les valeurs negatives de l'eau
331         nbptspb=0
332            DO k = 1, klev
333            DO i = 1, klon
334               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
335               if (logexpr2(i,k)) then
336                q_seri(i,k)=1.e-15
337                nbptspb=nbptspb+1
338!                if (prt_level.ge.10) then
339!                  print*,'WARN eau<0 apres therm i=',i,'  k=',k  &
340!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k),  &
341!    &         'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
342                 endif
343            ENDDO
344            ENDDO
345        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
346! tests sur les valeurs de la temperature
347        nbptspb=0
348            DO k = 1, klev
349            DO i = 1, klon
350               logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370.
351               if (logexpr2(i,k)) nbptspb=nbptspb+1
352!              if ((t_seri(i,k).lt.50.) .or.  &
353!    &              (t_seri(i,k).gt.370.)) then
354!                 print*,'WARN temp apres therm i=',i,'  k=',k  &
355!    &         ,' t_seri',t_seri(i,k)
356!              CALL abort
357!              endif
358            ENDDO
359            ENDDO
360        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
361         enddo ! isplit
362
363!
364!***************************************************************
365!     calcul du flux ascencant conservatif
366!            print*,'<<<<calcul flux ascendant conservatif'
367
368      fmc_therm=0.
369               do k=1,klev
370            do i=1,klon
371                  if (entr_therm(i,k).gt.0.) then
372                     fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
373                  else
374                     fmc_therm(i,k+1)=fmc_therm(i,k)
375                  endif
376                  detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1))  &
377     &                 -(fmc_therm(i,k)-fm_therm(i,k))
378               enddo
379            enddo
380     
381     
382!****************************************************************
383!     calcul de l'humidite dans l'ascendance
384!      print*,'<<<<calcul de lhumidite dans thermique'
385!CR:on ne le calcule que pour le cas sec
386      if (iflag_thermals.le.11) then     
387      do i=1,klon
388         zqasc(i,1)=q_seri(i,1)
389         do k=2,klev
390            if (fmc_therm(i,k+1).gt.1.e-6) then
391               zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1)  &
392     &              +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
393!CR:test on asseche le thermique
394!               zqasc(i,k)=zqasc(i,k)/2.
395!            else
396!               zqasc(i,k)=q_seri(i,k)
397            endif
398         enddo
399       enddo
400     
401
402!     calcul de l'eau condensee dans l'ascendance
403!             print*,'<<<<calcul de leau condensee dans thermique'
404             do i=1,klon
405                do k=1,klev
406                   clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
407                   if (clwcon0(i,k).lt.0. .or.   &
408     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
409                      clwcon0(i,k)=0.
410                   endif
411                enddo
412             enddo
413       else
414              do i=1,klon
415                do k=1,klev
416                   clwcon0(i,k)=zqla(i,k) 
417                   if (clwcon0(i,k).lt.0. .or.   &
418     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
419                   clwcon0(i,k)=0.
420                   endif
421                enddo
422             enddo
423       endif
424!*******************************************************************   
425
426
427!jyg  Protection contre les temperatures nulles
428          do i=1,klon
429             do k=1,klev
430                if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0.
431             enddo
432          enddo
433
434
435      return
436
437      end
Note: See TracBrowser for help on using the repository browser.