source: LMDZ6/branches/Ocean_skin/libf/phylmd/calltherm.F90 @ 5008

Last change on this file since 5008 was 4368, checked in by lguez, 2 years ago

Sync latest trunk changes to Ocean_skin

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.9 KB
Line 
1!
2! $Id: calltherm.F90 4368 2022-12-05 23:01:16Z abarral $
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!!! nrlmd le 10/04/2012
12     &      ,pbl_tke,pctsrf,omega,airephy &
13     &      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
14     &      ,n2,s2,ale_bl_stat &
15     &      ,therm_tke_max,env_tke_max &
16     &      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
17     &      ,alp_bl_conv,alp_bl_stat &
18!!! fin nrlmd le 10/04/2012
19     &      ,zqla,ztva &
20#ifdef ISO         
21     &      ,xt_seri,d_xt_ajs &
22#ifdef DIAGISO         
23     &      ,q_the,xt_the &
24#endif
25#endif         
26     &   )
27
28      USE dimphy
29      USE indice_sol_mod
30      USE print_control_mod, ONLY: prt_level,lunout
31#ifdef ISO
32      use infotrac_phy, ONLY: ntiso
33#ifdef ISOVERIF
34      USE isotopes_mod, ONLY: iso_eau,iso_HDO
35      USE isotopes_verif_mod, ONLY: iso_verif_aberrant_enc_vect2D, &
36        iso_verif_egalite_vect2D
37#endif   
38#endif
39
40      implicit none
41      include "clesphys.h"
42      include "thermcell_old.h"
43
44
45!IM 140508
46      INTEGER, SAVE ::  itap
47!$OMP THREADPRIVATE(itap)
48      REAL dtime
49      LOGICAL debut
50      LOGICAL logexpr0, logexpr2(klon,klev), logexpr1(klon)
51      REAL fact(klon)
52      INTEGER nbptspb
53
54      REAL u_seri(klon,klev),v_seri(klon,klev)
55      REAL t_seri(klon,klev),q_seri(klon,klev),qmemoire(klon,klev)
56      REAL weak_inversion(klon)
57      REAL paprs(klon,klev+1)
58      REAL pplay(klon,klev)
59      REAL pphi(klon,klev)
60      real zlev(klon,klev+1)
61!test: on sort lentr et a* pour alimenter KE
62      REAL wght_th(klon,klev)
63      INTEGER lalim_conv(klon)
64      REAL zw2(klon,klev+1),fraca(klon,klev+1)
65
66!FH Update Thermiques
67      REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
68      REAL d_u_ajs(klon,klev),d_v_ajs(klon,klev)
69      real fm_therm(klon,klev+1)
70      real entr_therm(klon,klev),detr_therm(klon,klev)
71
72!********************************************************
73!     declarations
74      LOGICAL flag_bidouille_stratocu
75      real fmc_therm(klon,klev+1),zqasc(klon,klev)
76      real zqla(klon,klev)
77      real zqta(klon,klev)
78      real ztv(klon,klev),ztva(klon,klev)
79      real zpspsk(klon,klev)
80      real ztla(klon,klev)
81      real zthl(klon,klev)
82      real wmax_sec(klon)
83      real zmax_sec(klon)
84      real f_sec(klon)
85      real detrc_therm(klon,klev)
86! FH WARNING : il semble que ces save ne servent a rien
87!     save fmc_therm, detrc_therm
88      real clwcon0(klon,klev)
89      real zqsat(klon,klev)
90      real zw_sec(klon,klev+1)
91      integer lmix_sec(klon)
92      integer lmax(klon)
93      real ratqscth(klon,klev)
94      real ratqsdiff(klon,klev)
95      real zqsatth(klon,klev) 
96!nouvelles variables pour la convection
97      real ale_bl(klon)
98      real alp_bl(klon)
99      real ale(klon)
100      real alp(klon)
101!RC
102      !on garde le zmax du pas de temps precedent
103      real zmax0(klon), f0(klon)
104
105!!! nrlmd le 10/04/2012
106      real pbl_tke(klon,klev+1,nbsrf)
107      real pctsrf(klon,nbsrf)
108      real omega(klon,klev)
109      real airephy(klon)
110      real zlcl_th(klon),fraca0(klon),w0(klon),w_conv(klon)
111      real therm_tke_max0(klon),env_tke_max0(klon)
112      real n2(klon),s2(klon)
113      real ale_bl_stat(klon)
114      real therm_tke_max(klon,klev),env_tke_max(klon,klev)
115      real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
116!!! fin nrlmd le 10/04/2012
117
118!********************************************************
119
120      real, dimension(klon) :: pcon
121      real, dimension(klon,klev) :: rhobarz,wth3
122      integer,dimension(klon) :: lalim
123      real, dimension(klon,klev+1) :: fm
124      real, dimension(klon,klev) :: alim_star
125      real, dimension(klon) :: zmax
126
127
128
129
130! variables locales
131      REAL d_t_the(klon,klev), d_q_the(klon,klev)
132      REAL d_u_the(klon,klev),d_v_the(klon,klev)
133!
134      real zfm_therm(klon,klev+1),zdt
135      real zentr_therm(klon,klev),zdetr_therm(klon,klev)
136! FH A VERIFIER : SAVE INUTILES
137!      save zentr_therm,zfm_therm
138
139      character (len=20) :: modname='calltherm'
140      character (len=80) :: abort_message
141
142      integer i,k,isplit
143      logical, save :: first=.true.
144      logical :: new_thermcell
145
146#ifdef ISO
147      REAL xt_seri(ntiso,klon,klev),xtmemoire(ntiso,klon,klev)
148      REAL d_xt_ajs(ntiso,klon,klev)
149      real d_xt_the(ntiso,klon,klev)
150#ifdef DIAGISO
151      real q_the(klon,klev)
152      real xt_the(ntiso,klon,klev)
153#endif
154      real qprec(klon,klev)
155      integer ixt
156#endif
157
158
159!$OMP THREADPRIVATE(first)
160!********************************************************
161      if (first) then
162        itap=0
163        first=.false.
164      endif
165
166! Incrementer le compteur de la physique
167     itap   = itap + 1
168
169!  Modele du thermique
170!  ===================
171!         print*,'thermiques: WARNING on passe t au lieu de t_seri'
172
173
174! On prend comme valeur initiale des thermiques la valeur du pas
175! de temps precedent
176         zfm_therm(:,:)=fm_therm(:,:)
177         zdetr_therm(:,:)=detr_therm(:,:)
178         zentr_therm(:,:)=entr_therm(:,:)
179
180! On reinitialise les flux de masse a zero pour le cumul en
181! cas de splitting
182         fm_therm(:,:)=0.
183         entr_therm(:,:)=0.
184         detr_therm(:,:)=0.
185
186         ale_bl(:)=0.
187         alp_bl(:)=0.
188         if (prt_level.ge.10) then
189          print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
190         endif
191
192!   tests sur les valeurs negatives de l'eau
193         logexpr0=prt_level.ge.10
194         nbptspb=0
195         do k=1,klev
196            do i=1,klon
197! Attention teste abderr 19-03-09
198!               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
199                logexpr2(i,k)=.not.q_seri(i,k).ge.1.e-15
200               if (logexpr2(i,k)) then
201#ifdef ISO
202                qprec(i,k)=q_seri(i,k)
203#endif
204                q_seri(i,k)=1.e-15
205                nbptspb=nbptspb+1
206#ifdef ISO
207                do ixt=1,ntiso
208                  xt_seri(ixt,i,k)=1.e-15*(xt_seri(ixt,i,k)/qprec(i,k))
209                  ! xt_seri(ixt,i,k)=1.e-15*(Rdefault(index_iso(ixt)))
210                enddo
211#endif
212               endif
213!               if (logexpr0) &
214!    &             print*,'WARN eau<0 avant therm i=',i,'  k=',k  &
215!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k)
216            enddo
217         enddo
218         if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
219
220
221         new_thermcell=iflag_thermals>=15.and.iflag_thermals<=18
222#ifdef ISO
223      if (.not.new_thermcell) then
224           CALL abort_gcm('calltherm 234','isos pas prevus ici',1)
225      endif
226#ifdef ISOVERIF
227      if (iso_eau.gt.0) then
228       call iso_verif_egalite_vect2D( &
229     &           xt_seri,q_seri, &
230     &           'calltherm 174',ntiso,klon,klev)
231      endif !if (iso_eau.gt.0) then
232#endif   
233#endif
234         zdt=dtime/REAL(nsplit_thermals)
235
236
237         do isplit=1,nsplit_thermals
238
239          if (iflag_thermals>=1000) then
240            CALL thermcell_2002(klon,klev,zdt,iflag_thermals   &
241     &      ,pplay,paprs,pphi  &
242     &      ,u_seri,v_seri,t_seri,q_seri  &
243     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
244     &      ,zfm_therm,zentr_therm,fraca,zw2  &
245     &      ,r_aspect_thermals,30.,w2di_thermals  &
246     &      ,tau_thermals)
247          else if (iflag_thermals.eq.2) then
248            CALL thermcell_sec(klon,klev,zdt  &
249     &      ,pplay,paprs,pphi,zlev  &
250     &      ,u_seri,v_seri,t_seri,q_seri  &
251     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
252     &      ,zfm_therm,zentr_therm  &
253     &      ,r_aspect_thermals,30.,w2di_thermals  &
254     &      ,tau_thermals)
255          else if (iflag_thermals.eq.3) then
256            CALL thermcell(klon,klev,zdt  &
257     &      ,pplay,paprs,pphi  &
258     &      ,u_seri,v_seri,t_seri,q_seri  &
259     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
260     &      ,zfm_therm,zentr_therm  &
261     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
262     &      ,tau_thermals)
263          else if (iflag_thermals.eq.10) then
264            CALL thermcell_eau(klon,klev,zdt  &
265     &      ,pplay,paprs,pphi  &
266     &      ,u_seri,v_seri,t_seri,q_seri  &
267     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
268     &      ,zfm_therm,zentr_therm  &
269     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
270     &      ,tau_thermals)
271          else if (iflag_thermals.eq.11) then
272              abort_message = 'cas non prevu dans calltherm'
273              CALL abort_physic (modname,abort_message,1)
274          else if (iflag_thermals.eq.12) then
275            CALL calcul_sec(klon,klev,zdt  &
276     &      ,pplay,paprs,pphi,zlev  &
277     &      ,u_seri,v_seri,t_seri,q_seri  &
278     &      ,zmax_sec,wmax_sec,zw_sec,lmix_sec  &
279     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
280     &      ,tau_thermals)
281          else if (iflag_thermals==13.or.iflag_thermals==14) then
282              abort_message = 'thermcellV0_main enleve svn>2084'
283              CALL abort_physic (modname,abort_message,1)
284          else if (new_thermcell) then
285            CALL thermcell_main(itap,klon,klev,zdt  &
286     &      ,pplay,paprs,pphi,debut  &
287     &      ,u_seri,v_seri,t_seri,q_seri  &
288     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
289     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
290     &      ,ratqscth,ratqsdiff,zqsatth  &
291     &      ,zmax0,f0,zw2,fraca,ztv,zpspsk &
292     &      ,ztla,zthl,ztva &
293     &      ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
294#ifdef ISO         
295     &      ,xt_seri,d_xt_the &
296#endif         
297     &   )
298
299            CALL thermcell_alp(klon,klev,zdt  &                      ! in
300     &        ,pplay,paprs  &                                        ! in
301     &        ,zfm_therm,zentr_therm,lmax  &                         ! in
302     &        ,pbl_tke,pctsrf,omega,airephy &                        ! in
303     &        ,zw2,fraca &                                           ! in
304     &        ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &  ! in
305     &        ,ale,alp,lalim_conv,wght_th &                          ! out
306     &        ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &! out
307     &        ,n2,s2,ale_bl_stat &                                   ! out
308     &        ,therm_tke_max,env_tke_max &                           ! out
309     &        ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &          ! out
310     &        ,alp_bl_conv,alp_bl_stat &                             ! out
311     &        )
312
313           if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
314         else
315           abort_message = 'Cas des thermiques non prevu'
316           CALL abort_physic (modname,abort_message,1)
317         endif
318
319! Attention : les noms sont contre intuitif.
320! flag_bidouille_stratocu est .true. si on ne fait pas de bidouille.
321! Il aurait mieux valu avoir un nobidouille_stratocu
322! Et pour simplifier :
323! nobidouille_stratocu=.not.(iflag_thermals==13.or.iflag_thermals=15)
324! Ce serait bien de changer, mai en prenant le temps de vérifier que ca
325! fait bien ce qu'on croit.
326
327       flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18
328
329! Calcul a posteriori du niveau max des thermiques pour les schémas qui
330! ne la sortent pas.
331      if (iflag_thermals<=12.or.iflag_thermals>=1000) then
332         lmax(:)=1
333         do k=1,klev-1
334            zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1)
335         enddo
336         do k=1,klev-1
337            do i=1,klon
338               if (zfm_therm(i,k+1)>0.) lmax(i)=k
339            enddo
340         enddo
341      endif
342
343      fact(:)=0.
344      DO i=1,klon
345       logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
346       IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals)
347      ENDDO
348
349     DO k=1,klev
350!  transformation de la derivee en tendance
351            d_t_the(:,k)=d_t_the(:,k)*dtime*fact(:)
352            d_u_the(:,k)=d_u_the(:,k)*dtime*fact(:)
353            d_v_the(:,k)=d_v_the(:,k)*dtime*fact(:)
354            d_q_the(:,k)=d_q_the(:,k)*dtime*fact(:)
355            fm_therm(:,k)=fm_therm(:,k)  &
356     &      +zfm_therm(:,k)*fact(:)
357            entr_therm(:,k)=entr_therm(:,k)  &
358     &       +zentr_therm(:,k)*fact(:)
359            detr_therm(:,k)=detr_therm(:,k)  &
360     &       +zdetr_therm(:,k)*fact(:)
361#ifdef ISO
362            do ixt=1,ntiso
363              d_xt_the(ixt,:,k)=d_xt_the(ixt,:,k)*dtime*fact(:)
364            enddo
365#endif
366      ENDDO
367       fm_therm(:,klev+1)=0.
368
369
370
371!  accumulation de la tendance
372            d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
373            d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
374            d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
375            d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
376#ifdef ISO
377            d_xt_ajs(:,:,:)=d_xt_ajs(:,:,:)+d_xt_the(:,:,:)
378#endif
379
380!  incrementation des variables meteo
381            t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
382            u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
383            v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
384            qmemoire(:,:)=q_seri(:,:)
385            q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
386#ifdef ISO
387            xtmemoire(:,:,:)=xt_seri(:,:,:)
388            xt_seri(:,:,:) = xt_seri(:,:,:) + d_xt_the(:,:,:)
389#ifdef ISOVERIF
390!      write(*,*) 'calltherm 350 tmp: ajout d_xt_the'
391      if (iso_HDO.gt.0) then
392!      i=479
393!      k=4
394!      write(*,*) 'xt_seri(iso_hdo,i,k),q_seri(i,k)=', &
395!     &   xt_seri(iso_hdo,i,k),q_seri(i,k)
396!      write(*,*) 'd_xt_the(iso_hdo,i,k),d_q_the(i,k)=', &
397!     &   d_xt_the(iso_hdo,i,k),d_q_the(i,k)
398      call iso_verif_aberrant_enc_vect2D( &
399     &        xt_seri,q_seri, &
400     &        'calltherm 353, apres ajout d_xt_the',ntiso,klon,klev)
401      endif     
402#endif
403#endif
404           if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
405
406       DO i=1,klon
407            fm_therm(i,klev+1)=0.
408            ale_bl(i)=ale_bl(i)+ale(i)/REAL(nsplit_thermals)
409!            write(22,*)'ALE CALLTHERM',ale_bl(i),ale(i)
410            alp_bl(i)=alp_bl(i)+alp(i)/REAL(nsplit_thermals)
411!            write(23,*)'ALP CALLTHERM',alp_bl(i),alp(i)
412        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)
413       ENDDO
414
415!IM 060508 marche pas comme cela !!!        enddo ! isplit
416
417!   tests sur les valeurs negatives de l'eau
418         nbptspb=0
419            DO k = 1, klev
420            DO i = 1, klon
421               logexpr2(i,k)=.not.q_seri(i,k).ge.0.
422               if (logexpr2(i,k)) then
423                q_seri(i,k)=1.e-15
424                nbptspb=nbptspb+1
425#ifdef ISO
426                do ixt=1,ntiso
427                  xt_seri(ixt,i,k)=1.e-15*(xtmemoire(ixt,i,k)/qmemoire(i,k))
428                enddo
429#endif
430!                if (prt_level.ge.10) then
431!                  print*,'WARN eau<0 apres therm i=',i,'  k=',k  &
432!    &         ,' dq,q',d_q_the(i,k),q_seri(i,k),  &
433!    &         'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
434                 endif
435            ENDDO
436            ENDDO
437#ifdef ISO
438#ifdef ISOVERIF
439      if (iso_HDO.gt.0) then
440      call iso_verif_aberrant_enc_vect2D( &
441     &        xt_seri,q_seri, &
442     &        'calltherm 393, apres bidouille q<0',ntiso,klon,klev)
443      endif     
444#endif
445#endif
446
447        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb   
448! tests sur les valeurs de la temperature
449        nbptspb=0
450            DO k = 1, klev
451            DO i = 1, klon
452               logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370.
453               if (logexpr2(i,k)) nbptspb=nbptspb+1
454!              if ((t_seri(i,k).lt.50.) .or.  &
455!    &              (t_seri(i,k).gt.370.)) then
456!                 print*,'WARN temp apres therm i=',i,'  k=',k  &
457!    &         ,' t_seri',t_seri(i,k)
458!              CALL abort
459!              endif
460            ENDDO
461            ENDDO
462        IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
463         enddo ! isplit
464
465!
466!***************************************************************
467!     calcul du flux ascencant conservatif
468!            print*,'<<<<calcul flux ascendant conservatif'
469
470      fmc_therm=0.
471               do k=1,klev
472            do i=1,klon
473                  if (entr_therm(i,k).gt.0.) then
474                     fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
475                  else
476                     fmc_therm(i,k+1)=fmc_therm(i,k)
477                  endif
478                  detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1))  &
479     &                 -(fmc_therm(i,k)-fm_therm(i,k))
480               enddo
481            enddo
482     
483     
484!****************************************************************
485!     calcul de l'humidite dans l'ascendance
486!      print*,'<<<<calcul de lhumidite dans thermique'
487!CR:on ne le calcule que pour le cas sec
488      if (iflag_thermals.le.11) then     
489      do i=1,klon
490         zqasc(i,1)=q_seri(i,1)
491         do k=2,klev
492            if (fmc_therm(i,k+1).gt.1.e-6) then
493               zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1)  &
494     &              +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
495!CR:test on asseche le thermique
496!               zqasc(i,k)=zqasc(i,k)/2.
497!            else
498!               zqasc(i,k)=q_seri(i,k)
499            endif
500         enddo
501       enddo
502     
503
504!     calcul de l'eau condensee dans l'ascendance
505!             print*,'<<<<calcul de leau condensee dans thermique'
506             do i=1,klon
507                do k=1,klev
508                   clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
509                   if (clwcon0(i,k).lt.0. .or.   &
510     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
511                      clwcon0(i,k)=0.
512                   endif
513                enddo
514             enddo
515       else
516              do i=1,klon
517                do k=1,klev
518                   clwcon0(i,k)=zqla(i,k) 
519                   if (clwcon0(i,k).lt.0. .or.   &
520     &             (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
521                   clwcon0(i,k)=0.
522                   endif
523                enddo
524             enddo
525       endif
526!*******************************************************************   
527
528
529!jyg  Protection contre les temperatures nulles
530          do i=1,klon
531             do k=1,klev
532                if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0.
533             enddo
534          enddo
535
536
537      return
538
539      end
Note: See TracBrowser for help on using the repository browser.