source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/thermcell_dq.F90 @ 2404

Last change on this file since 2404 was 1985, checked in by fhourdin, 11 years ago

Modification de du schéma implicite des thermique pour que la concentration

soit inchangée en cas de flux de masse uniformément nul.

Implicit scheme for thermal plume modified to keep unchanged the concentration

when the mass flux is zero.

  • 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: 8.6 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,impl,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,1',k,entr(ig,k)*ptimestep,masse(ig,k)
56               abort_message = 'entr dt > m, 1st'
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! FH debut de modif : le calcul ci dessous modifiait numériquement
127! la concentration quand le flux de masse etait nul car on divisait
128! puis multipliait par masse/ptimestep.
129!        q(:,k)=(masse(:,k)*q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1)) &
130!    &               /(fm(:,k)+masse(:,k)/ptimestep)
131         q(:,k)=(q(:,k)+ptimestep/masse(:,k)*(fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1))) &
132      &               /(1.+fm(:,k)*ptimestep/masse(:,k))
133! FH fin de modif.
134      enddo
135   endif
136
137! Tendencies
138      do k=1,nlay
139         do ig=1,ngrid
140            dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
141            q(ig,k)=qold(ig,k)
142         enddo
143      enddo
144
145return
146end
147
148!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
149! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations
150!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
151
152      subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,  &
153     &           masse,q,dq,qa,lev_out)
154      implicit none
155
156#include "iniprint.h"
157!=======================================================================
158!
159!   Calcul du transport verticale dans la couche limite en presence
160!   de "thermiques" explicitement representes
161!   calcul du dq/dt une fois qu'on connait les ascendances
162!
163!=======================================================================
164
165      integer ngrid,nlay,impl
166
167      real ptimestep
168      real masse(ngrid,nlay),fm(ngrid,nlay+1)
169      real entr(ngrid,nlay)
170      real q(ngrid,nlay)
171      real dq(ngrid,nlay)
172      integer lev_out                           ! niveau pour les print
173
174      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
175
176      real zzm
177
178      integer ig,k
179      real cfl
180
181      real qold(ngrid,nlay)
182      real ztimestep
183      integer niter,iter
184      CHARACTER (LEN=20) :: modname='thermcell_dq'
185      CHARACTER (LEN=80) :: abort_message
186
187
188
189! Calcul du critere CFL pour l'advection dans la subsidence
190      cfl = 0.
191      do k=1,nlay
192         do ig=1,ngrid
193            zzm=masse(ig,k)/ptimestep
194            cfl=max(cfl,fm(ig,k)/zzm)
195            if (entr(ig,k).gt.zzm) then
196               print*,'entr*dt>m,2',k,entr(ig,k)*ptimestep,masse(ig,k)
197               abort_message = 'entr dt > m, 2nd'
198               CALL abort_gcm (modname,abort_message,1)
199            endif
200         enddo
201      enddo
202
203!IM 090508     print*,'CFL CFL CFL CFL ',cfl
204
205#undef CFL
206#ifdef CFL
207! On subdivise le calcul en niter pas de temps.
208      niter=int(cfl)+1
209#else
210      niter=1
211#endif
212
213      ztimestep=ptimestep/niter
214      qold=q
215
216
217do iter=1,niter
218      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
219
220!   calcul du detrainement
221      do k=1,nlay
222         do ig=1,ngrid
223            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
224!           print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
225!test
226            if (detr(ig,k).lt.0.) then
227               entr(ig,k)=entr(ig,k)-detr(ig,k)
228               detr(ig,k)=0.
229!               print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
230!     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
231            endif
232            if (fm(ig,k+1).lt.0.) then
233!               print*,'fm2<0!!!'
234            endif
235            if (entr(ig,k).lt.0.) then
236!               print*,'entr2<0!!!'
237            endif
238         enddo
239      enddo
240
241!   calcul de la valeur dans les ascendances
242      do ig=1,ngrid
243         qa(ig,1)=q(ig,1)
244      enddo
245
246      do k=2,nlay
247         do ig=1,ngrid
248            if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt.  &
249     &         1.e-5*masse(ig,k)) then
250         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
251     &         /(fm(ig,k+1)+detr(ig,k))
252            else
253               qa(ig,k)=q(ig,k)
254            endif
255            if (qa(ig,k).lt.0.) then
256!               print*,'qa<0!!!'
257            endif
258            if (q(ig,k).lt.0.) then
259!               print*,'q<0!!!'
260            endif
261         enddo
262      enddo
263
264! Calcul du flux subsident
265
266      do k=2,nlay
267         do ig=1,ngrid
268#undef centre
269#ifdef centre
270             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
271#else
272
273#define plusqueun
274#ifdef plusqueun
275! Schema avec advection sur plus qu'une maille.
276            zzm=masse(ig,k)/ztimestep
277            if (fm(ig,k)>zzm) then
278               wqd(ig,k)=zzm*q(ig,k)+(fm(ig,k)-zzm)*q(ig,k+1)
279            else
280               wqd(ig,k)=fm(ig,k)*q(ig,k)
281            endif
282#else
283            wqd(ig,k)=fm(ig,k)*q(ig,k)
284#endif
285#endif
286
287            if (wqd(ig,k).lt.0.) then
288!               print*,'wqd<0!!!'
289            endif
290         enddo
291      enddo
292      do ig=1,ngrid
293         wqd(ig,1)=0.
294         wqd(ig,nlay+1)=0.
295      enddo
296     
297
298! Calcul des tendances
299      do k=1,nlay
300         do ig=1,ngrid
301            q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
302     &               -wqd(ig,k)+wqd(ig,k+1))  &
303     &               *ztimestep/masse(ig,k)
304!            if (dq(ig,k).lt.0.) then
305!               print*,'dq<0!!!'
306!            endif
307         enddo
308      enddo
309
310
311enddo
312
313
314! Calcul des tendances
315      do k=1,nlay
316         do ig=1,ngrid
317            dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
318            q(ig,k)=qold(ig,k)
319         enddo
320      enddo
321
322      return
323      end
Note: See TracBrowser for help on using the repository browser.