source: LMDZ5/branches/LF-private/libf/phylmd/thermcell_dq.F90 @ 5425

Last change on this file since 5425 was 1738, checked in by idelkadi, 12 years ago

Modifications for numerical stability of the boundary layer parameterizations.
Concerning yamada4.F :

  • option with asymptotic mixing length l0 imposed (iflag_pbl=8 and 9)
  • option with a new temporal scheme (iflag_pbl=10 and 11)
  • Correction for very stable PBLs (iflag_pbl=10 and 11) iflag_pbl=8 converges numerically with NPv3.1 iflag_pbl=11 -> now starts with NP from start files created by ce0l

-> the model can run with longer time-steps.

Concerning thermals :

Introduction of an implicit computation of vertical advection in
the environment of thermal plumes in thermcell_dq
impl = 0 : explicit, 1 : implicit, -1 : old version
controled by iflag_thermals =

15, 16 run with impl=-1 : numerical convergence with NPv3
17, 18 run with impl=1 : more stable

15 and 17 correspond to the activation of the stratocumulus "bidouille"

Modified routines (phylmd/):
calltherm.F90 : for managing the various options of thermcell_dq
coef_diff_turb_mod.F90 : yamada4 called for iflag_pbl<= 18 instead of 11
physiq.F : desactivation of the vertical diffusion of TKE for iflag_pbl=10
thermcellV0_main.F90 : calling thermcell_dq with implicit scheme
thermcell_dq.F90 : thermcell_dq with optional explicit/implicit scheme
thermcell_main.F90 : various option for vertical transport by thermals
yamada4.F : Yamada scheme with new options

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.2 KB
Line 
1      subroutine thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr,  &
2     &           masse,q,dq,qa,lev_out)
3      implicit none
4
5#include "iniprint.h"
6!=======================================================================
7!
8!   Calcul du transport verticale dans la couche limite en presence
9!   de "thermiques" explicitement representes
10!   calcul du dq/dt une fois qu'on connait les ascendances
11!
12! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
13!  Introduction of an implicit computation of vertical advection in
14!  the environment of thermal plumes in thermcell_dq
15!  impl =     0 : explicit, 1 : implicit, -1 : old version
16!
17!=======================================================================
18
19      integer ngrid,nlay,impl
20
21      real ptimestep
22      real masse(ngrid,nlay),fm(ngrid,nlay+1)
23      real entr(ngrid,nlay)
24      real q(ngrid,nlay)
25      real dq(ngrid,nlay)
26      integer lev_out                           ! niveau pour les print
27
28      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
29
30      real zzm
31
32      integer ig,k
33      real cfl
34
35      real qold(ngrid,nlay),fqa(ngrid,nlay+1)
36      integer niter,iter
37      CHARACTER (LEN=20) :: modname='thermcell_dq'
38      CHARACTER (LEN=80) :: abort_message
39
40
41! Old explicite scheme
42      if (impl==-1) then
43         call thermcell_dq_o(ngrid,nlay,ptimestep,fm,entr,  &
44     &           masse,q,dq,qa,lev_out)
45         return
46      endif
47
48! Calcul du critere CFL pour l'advection dans la subsidence
49      cfl = 0.
50      do k=1,nlay
51         do ig=1,ngrid
52            zzm=masse(ig,k)/ptimestep
53            cfl=max(cfl,fm(ig,k)/zzm)
54            if (entr(ig,k).gt.zzm) then
55               print*,'entr dt > m ',entr(ig,k)*ptimestep,masse(ig,k)
56               abort_message = ''
57               CALL abort_gcm (modname,abort_message,1)
58            endif
59         enddo
60      enddo
61
62      qold=q
63
64
65      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
66
67!   calcul du detrainement
68      do k=1,nlay
69         do ig=1,ngrid
70            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
71!           print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
72!test
73            if (detr(ig,k).lt.0.) then
74               entr(ig,k)=entr(ig,k)-detr(ig,k)
75               detr(ig,k)=0.
76!               print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
77!     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
78            endif
79            if (fm(ig,k+1).lt.0.) then
80!               print*,'fm2<0!!!'
81            endif
82            if (entr(ig,k).lt.0.) then
83!               print*,'entr2<0!!!'
84            endif
85         enddo
86      enddo
87
88! Computation of tracer concentrations in the ascending plume
89      do ig=1,ngrid
90         qa(ig,1)=q(ig,1)
91      enddo
92
93      do k=2,nlay
94         do ig=1,ngrid
95            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
96     &         1.e-5*masse(ig,k)) then
97         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
98     &         /(fm(ig,k+1)+detr(ig,k))
99            else
100               qa(ig,k)=q(ig,k)
101            endif
102            if (qa(ig,k).lt.0.) then
103!               print*,'qa<0!!!'
104            endif
105            if (q(ig,k).lt.0.) then
106!               print*,'q<0!!!'
107            endif
108         enddo
109      enddo
110
111! Plume vertical flux
112      do k=2,nlay-1
113         fqa(:,k)=fm(:,k)*qa(:,k-1)
114      enddo
115      fqa(:,1)=0. ; fqa(:,nlay)=0.
116
117
118! Trace species evolution
119   if (impl==0) then
120      do k=1,nlay-1
121         q(:,k)=q(:,k)+(fqa(:,k)-fqa(:,k+1)-fm(:,k)*q(:,k)+fm(:,k+1)*q(:,k+1)) &
122     &               *ptimestep/masse(:,k)
123      enddo
124   else
125      do k=nlay-1,1,-1
126         q(:,k)=(masse(:,k)*q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1)) &
127     &               /(fm(:,k)+masse(:,k)/ptimestep)
128      enddo
129   endif
130
131! Tendencies
132      do k=1,nlay
133         do ig=1,ngrid
134            dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
135            q(ig,k)=qold(ig,k)
136         enddo
137      enddo
138
139return
140end
141
142!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
143! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations
144!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
145
146      subroutine thermcell_dq_o(ngrid,nlay,ptimestep,fm,entr,  &
147     &           masse,q,dq,qa,lev_out)
148      implicit none
149
150#include "iniprint.h"
151!=======================================================================
152!
153!   Calcul du transport verticale dans la couche limite en presence
154!   de "thermiques" explicitement representes
155!   calcul du dq/dt une fois qu'on connait les ascendances
156!
157!=======================================================================
158
159      integer ngrid,nlay
160
161      real ptimestep
162      real masse(ngrid,nlay),fm(ngrid,nlay+1)
163      real entr(ngrid,nlay)
164      real q(ngrid,nlay)
165      real dq(ngrid,nlay)
166      integer lev_out                           ! niveau pour les print
167
168      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
169
170      real zzm
171
172      integer ig,k
173      real cfl
174
175      real qold(ngrid,nlay)
176      real ztimestep
177      integer niter,iter
178      CHARACTER (LEN=20) :: modname='thermcell_dq'
179      CHARACTER (LEN=80) :: abort_message
180
181
182
183! Calcul du critere CFL pour l'advection dans la subsidence
184      cfl = 0.
185      do k=1,nlay
186         do ig=1,ngrid
187            zzm=masse(ig,k)/ptimestep
188            cfl=max(cfl,fm(ig,k)/zzm)
189            if (entr(ig,k).gt.zzm) then
190               print*,'entr dt > m ',entr(ig,k)*ptimestep,masse(ig,k)
191               abort_message = ''
192               CALL abort_gcm (modname,abort_message,1)
193            endif
194         enddo
195      enddo
196
197!IM 090508     print*,'CFL CFL CFL CFL ',cfl
198
199#undef CFL
200#ifdef CFL
201! On subdivise le calcul en niter pas de temps.
202      niter=int(cfl)+1
203#else
204      niter=1
205#endif
206
207      ztimestep=ptimestep/niter
208      qold=q
209
210
211do iter=1,niter
212      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
213
214!   calcul du detrainement
215      do k=1,nlay
216         do ig=1,ngrid
217            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
218!           print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
219!test
220            if (detr(ig,k).lt.0.) then
221               entr(ig,k)=entr(ig,k)-detr(ig,k)
222               detr(ig,k)=0.
223!               print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
224!     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
225            endif
226            if (fm(ig,k+1).lt.0.) then
227!               print*,'fm2<0!!!'
228            endif
229            if (entr(ig,k).lt.0.) then
230!               print*,'entr2<0!!!'
231            endif
232         enddo
233      enddo
234
235!   calcul de la valeur dans les ascendances
236      do ig=1,ngrid
237         qa(ig,1)=q(ig,1)
238      enddo
239
240      do k=2,nlay
241         do ig=1,ngrid
242            if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt.  &
243     &         1.e-5*masse(ig,k)) then
244         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
245     &         /(fm(ig,k+1)+detr(ig,k))
246            else
247               qa(ig,k)=q(ig,k)
248            endif
249            if (qa(ig,k).lt.0.) then
250!               print*,'qa<0!!!'
251            endif
252            if (q(ig,k).lt.0.) then
253!               print*,'q<0!!!'
254            endif
255         enddo
256      enddo
257
258! Calcul du flux subsident
259
260      do k=2,nlay
261         do ig=1,ngrid
262#undef centre
263#ifdef centre
264             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
265#else
266
267#define plusqueun
268#ifdef plusqueun
269! Schema avec advection sur plus qu'une maille.
270            zzm=masse(ig,k)/ztimestep
271            if (fm(ig,k)>zzm) then
272               wqd(ig,k)=zzm*q(ig,k)+(fm(ig,k)-zzm)*q(ig,k+1)
273            else
274               wqd(ig,k)=fm(ig,k)*q(ig,k)
275            endif
276#else
277            wqd(ig,k)=fm(ig,k)*q(ig,k)
278#endif
279#endif
280
281            if (wqd(ig,k).lt.0.) then
282!               print*,'wqd<0!!!'
283            endif
284         enddo
285      enddo
286      do ig=1,ngrid
287         wqd(ig,1)=0.
288         wqd(ig,nlay+1)=0.
289      enddo
290     
291
292! Calcul des tendances
293      do k=1,nlay
294         do ig=1,ngrid
295            q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
296     &               -wqd(ig,k)+wqd(ig,k+1))  &
297     &               *ztimestep/masse(ig,k)
298!            if (dq(ig,k).lt.0.) then
299!               print*,'dq<0!!!'
300!            endif
301         enddo
302      enddo
303
304
305enddo
306
307
308! Calcul des tendances
309      do k=1,nlay
310         do ig=1,ngrid
311            dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
312            q(ig,k)=qold(ig,k)
313         enddo
314      enddo
315
316      return
317      end
Note: See TracBrowser for help on using the repository browser.