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

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

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

  • IPSLCM5A configuration (comparison between trunk and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent
  • MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

  • IPSLCM5A configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent
  • MH07 configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files

Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

  • pour une config IPSLCM5A (comparaison tronc/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux (à part sur RN et Pb)
    • fichiers histoire égaux
  • pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux
    • fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

  • pour une config IPSLCM5A:
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc
  • pour une config MH07
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.0 KB
Line 
1!
2! $Id: calltherm.F90 1403 2010-07-01 09:02:53Z fairhead $
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
12      USE dimphy
13      implicit none
14#include "dimensions.h"
15!#include "dimphy.h"
16#include "thermcell.h"
17#include "iniprint.h"
18
19!IM 140508
20      INTEGER itap
21      REAL dtime
22      LOGICAL debut
23      LOGICAL logexpr0, logexpr2(klon,klev), logexpr1(klon)
24      REAL fact(klon)
25      INTEGER nbptspb
26
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)
34!test: on sort lentr et a* pour alimenter KE
35      REAL wght_th(klon,klev)
36      INTEGER lalim_conv(klon)
37      REAL zw2(klon,klev+1),fraca(klon,klev+1)
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)
42      real fm_therm(klon,klev+1)
43      real entr_therm(klon,klev),detr_therm(klon,klev)
44
45!********************************************************
46!     declarations
47      LOGICAL flag_bidouille_stratocu
48      real fmc_therm(klon,klev+1),zqasc(klon,klev)
49      real zqla(klon,klev)
50      real zqta(klon,klev)
51      real ztv(klon,klev)
52      real zpspsk(klon,klev)
53      real ztla(klon,klev)
54      real zthl(klon,klev)
55      real wmax_sec(klon)
56      real zmax_sec(klon)
57      real f_sec(klon)
58      real detrc_therm(klon,klev)
59! FH WARNING : il semble que ces save ne servent a rien
60!     save fmc_therm, detrc_therm
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) 
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
75      !on garde le zmax du pas de temps precedent
76      real zmax0(klon), f0(klon)
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!
84      real zfm_therm(klon,klev+1),zdt
85      real zentr_therm(klon,klev),zdetr_therm(klon,klev)
86! FH A VERIFIER : SAVE INUTILES
87!      save zentr_therm,zfm_therm
88
89      character (len=20) :: modname='calltherm'
90      character (len=80) :: abort_message
91
92      integer i,k
93      logical, save :: first=.true.
94!$OMP THREADPRIVATE(first)
95!********************************************************
96      if (first) then
97        itap=0
98        first=.false.
99      endif
100
101! Incrementer le compteur de la physique
102     itap   = itap + 1
103
104!  Modele du thermique
105!  ===================
106!         print*,'thermiques: WARNING on passe t au lieu de t_seri'
107
108
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
117         fm_therm(:,:)=0.
118         entr_therm(:,:)=0.
119         detr_therm(:,:)=0.
120
121         Ale_bl(:)=0.
122         Alp_bl(:)=0.
123         if (prt_level.ge.10) then
124          print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
125         endif
126
127!   tests sur les valeurs negatives de l'eau
128         logexpr0=prt_level.ge.10
129         nbptspb=0
130         do k=1,klev
131            do i=1,klon
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
135               if (logexpr2(i,k)) then
136                q_seri(i,k)=1.e-15
137                nbptspb=nbptspb+1
138               endif
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)
142            enddo
143         enddo
144         if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
145
146         zdt=dtime/REAL(nsplit_thermals)
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  &
156     &      ,tau_thermals,3)
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  &
164     &      ,tau_thermals,3)
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  &
172     &      ,tau_thermals,3)
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  &
180     &      ,tau_thermals,3)
181          else if (iflag_thermals.eq.11) then
182              abort_message = 'cas non prevu dans calltherm'
183              CALL abort_gcm (modname,abort_message,1)
184
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  &
191!    &      ,tau_thermals,3)
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  &
198     &      ,tau_thermals)
199          else if (iflag_thermals==13.or.iflag_thermals==14) then
200            CALL thermcellV0_main(itap,klon,klev,zdt  &
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  &
204     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
205     &      ,ratqscth,ratqsdiff,zqsatth  &
206     &      ,r_aspect_thermals,l_mix_thermals  &
207     &      ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
208     &      ,zmax0,f0,zw2,fraca)
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)
226         endif
227
228       flag_bidouille_stratocu=iflag_thermals.eq.14.or.iflag_thermals.eq.16
229
230      fact(:)=0.
231      DO i=1,klon
232       logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
233       IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals)
234      ENDDO
235
236     DO k=1,klev
237!  transformation de la derivee en tendance
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(:)
246            detr_therm(:,k)=detr_therm(:,k)  &
247     &       +zdetr_therm(:,k)*fact(:)
248      ENDDO
249       fm_therm(:,klev+1)=0.
250
251
252
253!  accumulation de la tendance
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(:,:)
258
259!  incrementation des variables meteo
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(:,:)
265           if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
266
267       DO i=1,klon
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)
269            fm_therm(i,klev+1)=0.
270            Ale_bl(i)=Ale_bl(i)+Ale(i)/REAL(nsplit_thermals)
271!            write(22,*)'ALE CALLTHERM',Ale_bl(i),Ale(i)
272            Alp_bl(i)=Alp_bl(i)+Alp(i)/REAL(nsplit_thermals)
273!            write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
274       ENDDO
275
276!IM 060508 marche pas comme cela !!!        enddo ! isplit
277
278!   tests sur les valeurs negatives de l'eau
279         nbptspb=0
280            DO k = 1, klev
281            DO i = 1, klon
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)
290                 endif
291            ENDDO
292            ENDDO
293        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
294! tests sur les valeurs de la temperature
295        nbptspb=0
296            DO k = 1, klev
297            DO i = 1, klon
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)
304!              CALL abort
305!              endif
306            ENDDO
307            ENDDO
308        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
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
375      return
376
377      end
Note: See TracBrowser for help on using the repository browser.