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

Last change on this file since 946 was 940, checked in by Laurent Fairhead, 16 years ago

On remplace le fichier include dimphy.h par le module dimphy.F90i pour etre
coherent avec le partout
LF

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