source: LMDZ4/trunk/libf/phylmd/calltherm.F90 @ 1085

Last change on this file since 1085 was 1026, checked in by lmdzadmin, 16 years ago

Modifs thermiques
FH/IM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.8 KB
Line 
1!
2! $Header$
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)
11
12      USE dimphy
13      implicit none
14#include "dimensions.h"
15!#include "dimphy.h"
16#include "thermcell.h"
17#include "iniprint.h"
18
19!  A inclure eventuellement dans les fichiers de configuration
20      data r_aspect_thermals,l_mix_thermals/2.,30./
21      data w2di_thermals/1/
22
23!IM 140508
24      INTEGER itap
25      REAL dtime
26      LOGICAL debut
27      LOGICAL logexpr0, logexpr2(klon,klev), logexpr1(klon)
28      REAL fact(klon)
29      INTEGER nbptspb
30
31      REAL u_seri(klon,klev),v_seri(klon,klev)
32      REAL t_seri(klon,klev),q_seri(klon,klev),qmemoire(klon,klev)
33      REAL weak_inversion(klon)
34      REAL paprs(klon,klev+1)
35      REAL pplay(klon,klev)
36      REAL pphi(klon,klev)
37      real zlev(klon,klev+1)
38!test: on sort lentr et a* pour alimenter KE
39      REAL wght_th(klon,klev)
40      INTEGER lalim_conv(klon)
41      REAL zw2(klon,klev+1),fraca(klon,klev+1)
42
43!FH Update Thermiques
44      REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
45      REAL d_u_ajs(klon,klev),d_v_ajs(klon,klev)
46      real fm_therm(klon,klev+1)
47      real entr_therm(klon,klev),detr_therm(klon,klev)
48
49!********************************************************
50!     declarations
51      real fmc_therm(klon,klev+1),zqasc(klon,klev)
52      real zqla(klon,klev)
53      real zqta(klon,klev)
54      real wmax_sec(klon)
55      real zmax_sec(klon)
56      real f_sec(klon)
57      real detrc_therm(klon,klev)
58! FH WARNING : il semble que ces save ne servent a rien
59!     save fmc_therm, detrc_therm
60      real clwcon0(klon,klev)
61      real zqsat(klon,klev)
62      real zw_sec(klon,klev+1)
63      integer lmix_sec(klon)
64      integer lmax(klon)
65      real ratqscth(klon,klev)
66      real ratqsdiff(klon,klev)
67      real zqsatth(klon,klev) 
68!nouvelles variables pour la convection
69      real Ale_bl(klon)
70      real Alp_bl(klon)
71      real Ale(klon)
72      real Alp(klon)
73!RC
74      !on garde le zmax du pas de temps precedent
75      real zmax0(klon), f0(klon)
76!********************************************************
77
78
79! variables locales
80      REAL d_t_the(klon,klev), d_q_the(klon,klev)
81      REAL d_u_the(klon,klev),d_v_the(klon,klev)
82!
83      real zfm_therm(klon,klev+1),zdt
84      real zentr_therm(klon,klev),zdetr_therm(klon,klev)
85! FH A VERIFIER : SAVE INUTILES
86!      save zentr_therm,zfm_therm
87
88      integer i,k
89      logical, save :: first=.true.
90!$OMP THREADPRIVATE(first)
91!********************************************************
92      if (first) then
93        itap=0
94        first=.false.
95      endif
96
97! Incrementer le compteur de la physique
98     itap   = itap + 1
99
100!  Modele du thermique
101!  ===================
102!         print*,'thermiques: WARNING on passe t au lieu de t_seri'
103
104
105! On prend comme valeur initiale des thermiques la valeur du pas
106! de temps precedent
107         zfm_therm(:,:)=fm_therm(:,:)
108         zdetr_therm(:,:)=detr_therm(:,:)
109         zentr_therm(:,:)=entr_therm(:,:)
110
111! On reinitialise les flux de masse a zero pour le cumul en
112! cas de splitting
113         fm_therm(:,:)=0.
114         entr_therm(:,:)=0.
115         detr_therm(:,:)=0.
116
117         Ale_bl(:)=0.
118         Alp_bl(:)=0.
119         if (prt_level.ge.10) then
120          print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
121         endif
122
123!   tests sur les valeurs negatives de l'eau
124         logexpr0=prt_level.ge.10
125         nbptspb=0
126         do k=1,klev
127            do i=1,klon
128               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
129               if (logexpr2(i,k)) then
130                q_seri(i,k)=1.e-15
131                nbptspb=nbptspb+1
132               endif
133!               if (logexpr0) &
134!    &             print*,'WARN eau<0 avant therm i=',i,'  k=',k  &
135!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k)
136            enddo
137         enddo
138         if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
139
140         zdt=dtime/float(nsplit_thermals)
141         do isplit=1,nsplit_thermals
142
143          if (iflag_thermals.eq.1) then
144            CALL thermcell_2002(klon,klev,zdt   &
145     &      ,pplay,paprs,pphi  &
146     &      ,u_seri,v_seri,t_seri,q_seri  &
147     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
148     &      ,zfm_therm,zentr_therm  &
149     &      ,r_aspect_thermals,30.,w2di_thermals  &
150     &      ,tau_thermals,3)
151          else if (iflag_thermals.eq.2) then
152            CALL thermcell_sec(klon,klev,zdt  &
153     &      ,pplay,paprs,pphi,zlev  &
154     &      ,u_seri,v_seri,t_seri,q_seri  &
155     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
156     &      ,zfm_therm,zentr_therm  &
157     &      ,r_aspect_thermals,30.,w2di_thermals  &
158     &      ,tau_thermals,3)
159          else if (iflag_thermals.eq.3) then
160            CALL thermcell(klon,klev,zdt  &
161     &      ,pplay,paprs,pphi  &
162     &      ,u_seri,v_seri,t_seri,q_seri  &
163     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
164     &      ,zfm_therm,zentr_therm  &
165     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
166     &      ,tau_thermals,3)
167          else if (iflag_thermals.eq.10) then
168            CALL thermcell_eau(klon,klev,zdt  &
169     &      ,pplay,paprs,pphi  &
170     &      ,u_seri,v_seri,t_seri,q_seri  &
171     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
172     &      ,zfm_therm,zentr_therm  &
173     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
174     &      ,tau_thermals,3)
175          else if (iflag_thermals.eq.11) then
176            stop'cas non prevu dans calltherm'
177!           CALL thermcell_pluie(klon,klev,zdt  &
178!   &      ,pplay,paprs,pphi,zlev  &
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,zqla  &
182!    &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
183!    &      ,tau_thermals,3)
184          else if (iflag_thermals.eq.12) then
185            CALL calcul_sec(klon,klev,zdt  &
186     &      ,pplay,paprs,pphi,zlev  &
187     &      ,u_seri,v_seri,t_seri,q_seri  &
188     &      ,zmax_sec,wmax_sec,zw_sec,lmix_sec  &
189     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
190     &      ,tau_thermals)
191          else if (iflag_thermals.ge.13) then
192            CALL thermcell_main(itap,klon,klev,zdt  &
193     &      ,pplay,paprs,pphi,debut  &
194     &      ,u_seri,v_seri,t_seri,q_seri  &
195     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
196     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
197     &      ,ratqscth,ratqsdiff,zqsatth  &
198     &      ,r_aspect_thermals,l_mix_thermals  &
199     &      ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
200     &      ,zmax0,f0,zw2,fraca)
201         endif
202
203
204      fact(:)=0.
205      DO i=1,klon
206       logexpr1(i)=iflag_thermals.lt.14.or.weak_inversion(i).gt.0.5
207       IF(logexpr1(i)) fact(i)=1./float(nsplit_thermals)
208      ENDDO
209
210     DO k=1,klev
211!  transformation de la derivee en tendance
212            d_t_the(:,k)=d_t_the(:,k)*dtime*fact(:)
213            d_u_the(:,k)=d_u_the(:,k)*dtime*fact(:)
214            d_v_the(:,k)=d_v_the(:,k)*dtime*fact(:)
215            d_q_the(:,k)=d_q_the(:,k)*dtime*fact(:)
216            fm_therm(:,k)=fm_therm(:,k)  &
217     &      +zfm_therm(:,k)*fact(:)
218            entr_therm(:,k)=entr_therm(:,k)  &
219     &       +zentr_therm(:,k)*fact(:)
220            detr_therm(:,k)=detr_therm(:,k)  &
221     &       +zdetr_therm(:,k)*fact(:)
222      ENDDO
223       fm_therm(:,klev+1)=0.
224
225
226
227!  accumulation de la tendance
228            d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
229            d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
230            d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
231            d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
232
233!  incrementation des variables meteo
234            t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
235            u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
236            v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
237            qmemoire(:,:)=q_seri(:,:)
238            q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
239
240       DO i=1,klon
241        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)
242            fm_therm(i,klev+1)=0.
243            Ale_bl(i)=Ale_bl(i)+Ale(i)/float(nsplit_thermals)
244!            write(22,*)'ALE CALLTHERM',Ale_bl(i),Ale(i)
245            Alp_bl(i)=Alp_bl(i)+Alp(i)/float(nsplit_thermals)
246!            write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
247       ENDDO
248
249!IM 060508 marche pas comme cela !!!        enddo ! isplit
250
251!   tests sur les valeurs negatives de l'eau
252         nbptspb=0
253            DO k = 1, klev
254            DO i = 1, klon
255               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
256               if (logexpr2(i,k)) then
257                q_seri(i,k)=1.e-15
258                nbptspb=nbptspb+1
259!                if (prt_level.ge.10) then
260!                  print*,'WARN eau<0 apres therm i=',i,'  k=',k  &
261!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k),  &
262!    &         'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
263                 endif
264!       stop
265            ENDDO
266            ENDDO
267        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
268! tests sur les valeurs de la temperature
269        nbptspb=0
270            DO k = 1, klev
271            DO i = 1, klon
272               logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370.
273               if (logexpr2(i,k)) nbptspb=nbptspb+1
274!              if ((t_seri(i,k).lt.50.) .or.  &
275!    &              (t_seri(i,k).gt.370.)) then
276!                 print*,'WARN temp apres therm i=',i,'  k=',k  &
277!    &         ,' t_seri',t_seri(i,k)
278!              CALL abort
279!              endif
280            ENDDO
281            ENDDO
282        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
283         enddo ! isplit
284
285!
286!***************************************************************
287!     calcul du flux ascencant conservatif
288!            print*,'<<<<calcul flux ascendant conservatif'
289
290      fmc_therm=0.
291               do k=1,klev
292            do i=1,klon
293                  if (entr_therm(i,k).gt.0.) then
294                     fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
295                  else
296                     fmc_therm(i,k+1)=fmc_therm(i,k)
297                  endif
298                  detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1))  &
299     &                 -(fmc_therm(i,k)-fm_therm(i,k))
300               enddo
301            enddo
302     
303     
304!****************************************************************
305!     calcul de l'humidite dans l'ascendance
306!      print*,'<<<<calcul de lhumidite dans thermique'
307!CR:on ne le calcule que pour le cas sec
308      if (iflag_thermals.le.11) then     
309      do i=1,klon
310         zqasc(i,1)=q_seri(i,1)
311         do k=2,klev
312            if (fmc_therm(i,k+1).gt.1.e-6) then
313               zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1)  &
314     &              +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
315!CR:test on asseche le thermique
316!               zqasc(i,k)=zqasc(i,k)/2.
317!            else
318!               zqasc(i,k)=q_seri(i,k)
319            endif
320         enddo
321       enddo
322     
323
324!     calcul de l'eau condensee dans l'ascendance
325!             print*,'<<<<calcul de leau condensee dans thermique'
326             do i=1,klon
327                do k=1,klev
328                   clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
329                   if (clwcon0(i,k).lt.0. .or.   &
330     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
331                      clwcon0(i,k)=0.
332                   endif
333                enddo
334             enddo
335       else
336              do i=1,klon
337                do k=1,klev
338                   clwcon0(i,k)=zqla(i,k) 
339                   if (clwcon0(i,k).lt.0. .or.   &
340     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
341                   clwcon0(i,k)=0.
342                   endif
343                enddo
344             enddo
345       endif
346!*******************************************************************   
347
348
349      return
350
351      end
Note: See TracBrowser for help on using the repository browser.