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

Last change on this file since 3979 was 1428, checked in by Laurent Fairhead, 14 years ago

Problems with some initialisations in the new physics package. These are
temporary fixes. Further investigations needed. Trac ticket issued. JYG


Problèmes d'initialisations dans la nouvelle physique. Ces corrections
sont temporaires, il faut trouver la vraie raison des plantages. Un ticket
trac est créé. JYG

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