source: LMDZ4/branches/LMDZ4-dev/libf/phylmd/calltherm.F90 @ 1211

Last change on this file since 1211 was 1142, checked in by yann meurdesoif, 16 years ago

Portage vers Vargas

Y.M + E.M

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