My Project
 All Classes Files Functions Variables Macros
calltherm.F90
Go to the documentation of this file.
1 !
2 ! $Id: calltherm.F90 1776 2013-06-19 15:47:17Z 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 !!! 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  & )
20 
21  USE dimphy
22  implicit none
23 #include "dimensions.h"
24 !#include "dimphy.h"
25 #include "thermcell.h"
26 #include "iniprint.h"
27 
28 !!! nrlmd le 10/04/2012
29 #include "indicesol.h"
30 !!! fin nrlmd le 10/04/2012
31 
32 !IM 140508
33  INTEGER, SAVE :: itap
34 !$OMP THREADPRIVATE(itap)
35  REAL dtime
36  LOGICAL debut
37  LOGICAL logexpr0, logexpr2(klon,klev), logexpr1(klon)
38  REAL fact(klon)
39  INTEGER nbptspb
40 
41  REAL u_seri(klon,klev),v_seri(klon,klev)
42  REAL t_seri(klon,klev),q_seri(klon,klev),qmemoire(klon,klev)
43  REAL weak_inversion(klon)
44  REAL paprs(klon,klev+1)
45  REAL pplay(klon,klev)
46  REAL pphi(klon,klev)
47  real zlev(klon,klev+1)
48 !test: on sort lentr et a* pour alimenter KE
49  REAL wght_th(klon,klev)
50  INTEGER lalim_conv(klon)
51  REAL zw2(klon,klev+1),fraca(klon,klev+1)
52 
53 !FH Update Thermiques
54  REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
55  REAL d_u_ajs(klon,klev),d_v_ajs(klon,klev)
56  real fm_therm(klon,klev+1)
57  real entr_therm(klon,klev),detr_therm(klon,klev)
58 
59 !********************************************************
60 ! declarations
61  LOGICAL flag_bidouille_stratocu
62  real fmc_therm(klon,klev+1),zqasc(klon,klev)
63  real zqla(klon,klev)
64  real zqta(klon,klev)
65  real ztv(klon,klev)
66  real zpspsk(klon,klev)
67  real ztla(klon,klev)
68  real zthl(klon,klev)
69  real wmax_sec(klon)
70  real zmax_sec(klon)
71  real f_sec(klon)
72  real detrc_therm(klon,klev)
73 ! FH WARNING : il semble que ces save ne servent a rien
74 ! save fmc_therm, detrc_therm
75  real clwcon0(klon,klev)
76  real zqsat(klon,klev)
77  real zw_sec(klon,klev+1)
78  integer lmix_sec(klon)
79  integer lmax(klon)
80  real ratqscth(klon,klev)
81  real ratqsdiff(klon,klev)
82  real zqsatth(klon,klev)
83 !nouvelles variables pour la convection
84  real ale_bl(klon)
85  real alp_bl(klon)
86  real ale(klon)
87  real alp(klon)
88 !RC
89  !on garde le zmax du pas de temps precedent
90  real zmax0(klon), f0(klon)
91 
92 !!! nrlmd le 10/04/2012
93  real pbl_tke(klon,klev+1,nbsrf)
94  real pctsrf(klon,nbsrf)
95  real omega(klon,klev)
96  real airephy(klon)
97  real zlcl_th(klon),fraca0(klon),w0(klon),w_conv(klon)
98  real therm_tke_max0(klon),env_tke_max0(klon)
99  real n2(klon),s2(klon)
100  real ale_bl_stat(klon)
101  real therm_tke_max(klon,klev),env_tke_max(klon,klev)
102  real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
103 !!! fin nrlmd le 10/04/2012
104 
105 !********************************************************
106 
107 
108 ! variables locales
109  REAL d_t_the(klon,klev), d_q_the(klon,klev)
110  REAL d_u_the(klon,klev),d_v_the(klon,klev)
111 !
112  real zfm_therm(klon,klev+1),zdt
113  real zentr_therm(klon,klev),zdetr_therm(klon,klev)
114 ! FH A VERIFIER : SAVE INUTILES
115 ! save zentr_therm,zfm_therm
116 
117  character (len=20) :: modname='calltherm'
118  character (len=80) :: abort_message
119 
120  integer i,k
121  logical, save :: first=.true.
122 !$OMP THREADPRIVATE(first)
123 !********************************************************
124  if (first) then
125  itap=0
126  first=.false.
127  endif
128 
129 ! Incrementer le compteur de la physique
130  itap = itap + 1
131 
132 ! Modele du thermique
133 ! ===================
134 ! print*,'thermiques: WARNING on passe t au lieu de t_seri'
135 
136 
137 ! On prend comme valeur initiale des thermiques la valeur du pas
138 ! de temps precedent
139  zfm_therm(:,:)=fm_therm(:,:)
140  zdetr_therm(:,:)=detr_therm(:,:)
141  zentr_therm(:,:)=entr_therm(:,:)
142 
143 ! On reinitialise les flux de masse a zero pour le cumul en
144 ! cas de splitting
145  fm_therm(:,:)=0.
146  entr_therm(:,:)=0.
147  detr_therm(:,:)=0.
148 
149  ale_bl(:)=0.
150  alp_bl(:)=0.
151  if (prt_level.ge.10) then
152  print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
153  endif
154 
155 ! tests sur les valeurs negatives de l'eau
156  logexpr0=prt_level.ge.10
157  nbptspb=0
158  do k=1,klev
159  do i=1,klon
160 ! Attention teste abderr 19-03-09
161 ! logexpr2(i,k)=.not.q_seri(i,k).ge.0.
162  logexpr2(i,k)=.not.q_seri(i,k).ge.1.e-15
163  if (logexpr2(i,k)) then
164  q_seri(i,k)=1.e-15
165  nbptspb=nbptspb+1
166  endif
167 ! if (logexpr0) &
168 ! & print*,'WARN eau<0 avant therm i=',i,' k=',k &
169 ! & ,' dq,q',d_q_the(i,k),q_seri(i,k)
170  enddo
171  enddo
172  if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
173 
174  zdt=dtime/REAL(nsplit_thermals)
175  do isplit=1,nsplit_thermals
176 
177  if (iflag_thermals.eq.1) then
178  CALL thermcell_2002(klon,klev,zdt &
179  & ,pplay,paprs,pphi &
180  & ,u_seri,v_seri,t_seri,q_seri &
181  & ,d_u_the,d_v_the,d_t_the,d_q_the &
182  & ,zfm_therm,zentr_therm &
183  & ,r_aspect_thermals,30.,w2di_thermals &
184  & ,tau_thermals)
185  else if (iflag_thermals.eq.2) then
186  CALL thermcell_sec(klon,klev,zdt &
187  & ,pplay,paprs,pphi,zlev &
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 &
191  & ,r_aspect_thermals,30.,w2di_thermals &
192  & ,tau_thermals)
193  else if (iflag_thermals.eq.3) then
194  CALL thermcell(klon,klev,zdt &
195  & ,pplay,paprs,pphi &
196  & ,u_seri,v_seri,t_seri,q_seri &
197  & ,d_u_the,d_v_the,d_t_the,d_q_the &
198  & ,zfm_therm,zentr_therm &
199  & ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
200  & ,tau_thermals)
201  else if (iflag_thermals.eq.10) then
202  CALL thermcell_eau(klon,klev,zdt &
203  & ,pplay,paprs,pphi &
204  & ,u_seri,v_seri,t_seri,q_seri &
205  & ,d_u_the,d_v_the,d_t_the,d_q_the &
206  & ,zfm_therm,zentr_therm &
207  & ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
208  & ,tau_thermals)
209  else if (iflag_thermals.eq.11) then
210  abort_message = 'cas non prevu dans calltherm'
211  CALL abort_gcm(modname,abort_message,1)
212 
213 ! CALL thermcell_pluie(klon,klev,zdt &
214 ! & ,pplay,paprs,pphi,zlev &
215 ! & ,u_seri,v_seri,t_seri,q_seri &
216 ! & ,d_u_the,d_v_the,d_t_the,d_q_the &
217 ! & ,zfm_therm,zentr_therm,zqla &
218 ! & ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
219 ! & ,tau_thermals,3)
220  else if (iflag_thermals.eq.12) then
221  CALL calcul_sec(klon,klev,zdt &
222  & ,pplay,paprs,pphi,zlev &
223  & ,u_seri,v_seri,t_seri,q_seri &
224  & ,zmax_sec,wmax_sec,zw_sec,lmix_sec &
225  & ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
226  & ,tau_thermals)
227  else if (iflag_thermals==13.or.iflag_thermals==14) then
228  CALL thermcellv0_main(itap,klon,klev,zdt &
229  & ,pplay,paprs,pphi,debut &
230  & ,u_seri,v_seri,t_seri,q_seri &
231  & ,d_u_the,d_v_the,d_t_the,d_q_the &
232  & ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax &
233  & ,ratqscth,ratqsdiff,zqsatth &
234  & ,r_aspect_thermals,l_mix_thermals &
235  & ,tau_thermals,ale,alp,lalim_conv,wght_th &
236  & ,zmax0,f0,zw2,fraca)
237  else if (iflag_thermals>=15.and.iflag_thermals<=18) then
238 
239 ! print*,'THERM iflag_thermas_ed=',iflag_thermals_ed
240  CALL thermcell_main(itap,klon,klev,zdt &
241  & ,pplay,paprs,pphi,debut &
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,zdetr_therm,zqasc,zqla,lmax &
245  & ,ratqscth,ratqsdiff,zqsatth &
246 ! & ,r_aspect_thermals,l_mix_thermals &
247 ! & ,tau_thermals,iflag_thermals_ed,iflag_coupl &
248  & ,ale,alp,lalim_conv,wght_th &
249  & ,zmax0,f0,zw2,fraca,ztv,zpspsk &
250  & ,ztla,zthl &
251 !!! nrlmd le 10/04/2012
252  & ,pbl_tke,pctsrf,omega,airephy &
253  & ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
254  & ,n2,s2,ale_bl_stat &
255  & ,therm_tke_max,env_tke_max &
256  & ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
257  & ,alp_bl_conv,alp_bl_stat &
258 !!! fin nrlmd le 10/04/2012
259  & )
260  if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
261  else
262  abort_message = 'Cas des thermiques non prevu'
263  CALL abort_gcm(modname,abort_message,1)
264  endif
265 
266 ! Attention : les noms sont contre intuitif.
267 ! flag_bidouille_stratocu est .true. si on ne fait pas de bidouille.
268 ! Il aurait mieux valu avoir un nobidouille_stratocu
269 ! Et pour simplifier :
270 ! nobidouille_stratocu=.not.(iflag_thermals==13.or.iflag_thermals=15)
271 ! Ce serait bien de changer, mai en prenant le temps de vérifier que ca
272 ! fait bien ce qu'on croit.
273 
274  flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18
275 
276  if (iflag_thermals<=12) then
277  lmax=1
278  do k=1,klev-1
279  zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1)
280  enddo
281  endif
282 
283  fact(:)=0.
284  DO i=1,klon
285  logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
286  IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals)
287  ENDDO
288 
289  DO k=1,klev
290 ! transformation de la derivee en tendance
291  d_t_the(:,k)=d_t_the(:,k)*dtime*fact(:)
292  d_u_the(:,k)=d_u_the(:,k)*dtime*fact(:)
293  d_v_the(:,k)=d_v_the(:,k)*dtime*fact(:)
294  d_q_the(:,k)=d_q_the(:,k)*dtime*fact(:)
295  fm_therm(:,k)=fm_therm(:,k) &
296  & +zfm_therm(:,k)*fact(:)
297  entr_therm(:,k)=entr_therm(:,k) &
298  & +zentr_therm(:,k)*fact(:)
299  detr_therm(:,k)=detr_therm(:,k) &
300  & +zdetr_therm(:,k)*fact(:)
301  ENDDO
302  fm_therm(:,klev+1)=0.
303 
304 
305 
306 ! accumulation de la tendance
307  d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
308  d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
309  d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
310  d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
311 
312 ! incrementation des variables meteo
313  t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
314  u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
315  v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
316  qmemoire(:,:)=q_seri(:,:)
317  q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
318  if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
319 
320  DO i=1,klon
321  fm_therm(i,klev+1)=0.
322  ale_bl(i)=ale_bl(i)+ale(i)/REAL(nsplit_thermals)
323 ! write(22,*)'ALE CALLTHERM',Ale_bl(i),Ale(i)
324  alp_bl(i)=alp_bl(i)+alp(i)/REAL(nsplit_thermals)
325 ! write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
326  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)
327  ENDDO
328 
329 !IM 060508 marche pas comme cela !!! enddo ! isplit
330 
331 ! tests sur les valeurs negatives de l'eau
332  nbptspb=0
333  DO k = 1, klev
334  DO i = 1, klon
335  logexpr2(i,k)=.not.q_seri(i,k).ge.0.
336  if (logexpr2(i,k)) then
337  q_seri(i,k)=1.e-15
338  nbptspb=nbptspb+1
339 ! if (prt_level.ge.10) then
340 ! print*,'WARN eau<0 apres therm i=',i,' k=',k &
341 ! & ,' dq,q',d_q_the(i,k),q_seri(i,k), &
342 ! & 'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
343  endif
344  ENDDO
345  ENDDO
346  IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
347 ! tests sur les valeurs de la temperature
348  nbptspb=0
349  DO k = 1, klev
350  DO i = 1, klon
351  logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370.
352  if (logexpr2(i,k)) nbptspb=nbptspb+1
353 ! if ((t_seri(i,k).lt.50.) .or. &
354 ! & (t_seri(i,k).gt.370.)) then
355 ! print*,'WARN temp apres therm i=',i,' k=',k &
356 ! & ,' t_seri',t_seri(i,k)
357 ! CALL abort
358 ! endif
359  ENDDO
360  ENDDO
361  IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
362  enddo ! isplit
363 
364 !
365 !***************************************************************
366 ! calcul du flux ascencant conservatif
367 ! print*,'<<<<calcul flux ascendant conservatif'
368 
369  fmc_therm=0.
370  do k=1,klev
371  do i=1,klon
372  if (entr_therm(i,k).gt.0.) then
373  fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
374  else
375  fmc_therm(i,k+1)=fmc_therm(i,k)
376  endif
377  detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1)) &
378  & -(fmc_therm(i,k)-fm_therm(i,k))
379  enddo
380  enddo
381 
382 
383 !****************************************************************
384 ! calcul de l'humidite dans l'ascendance
385 ! print*,'<<<<calcul de lhumidite dans thermique'
386 !CR:on ne le calcule que pour le cas sec
387  if (iflag_thermals.le.11) then
388  do i=1,klon
389  zqasc(i,1)=q_seri(i,1)
390  do k=2,klev
391  if (fmc_therm(i,k+1).gt.1.e-6) then
392  zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1) &
393  & +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
394 !CR:test on asseche le thermique
395 ! zqasc(i,k)=zqasc(i,k)/2.
396 ! else
397 ! zqasc(i,k)=q_seri(i,k)
398  endif
399  enddo
400  enddo
401 
402 
403 ! calcul de l'eau condensee dans l'ascendance
404 ! print*,'<<<<calcul de leau condensee dans thermique'
405  do i=1,klon
406  do k=1,klev
407  clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
408  if (clwcon0(i,k).lt.0. .or. &
409  & (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
410  clwcon0(i,k)=0.
411  endif
412  enddo
413  enddo
414  else
415  do i=1,klon
416  do k=1,klev
417  clwcon0(i,k)=zqla(i,k)
418  if (clwcon0(i,k).lt.0. .or. &
419  & (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
420  clwcon0(i,k)=0.
421  endif
422  enddo
423  enddo
424  endif
425 !*******************************************************************
426 
427 
428 !jyg Protection contre les temperatures nulles
429  do i=1,klon
430  do k=1,klev
431  if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0.
432  enddo
433  enddo
434 
435 
436  return
437 
438  end