source: lmdz_wrf/trunk/WRFV3/lmdz/thermcell_dq.F90 @ 1544

Last change on this file since 1544 was 142, checked in by lfita, 10 years ago

Incorporing all the work done in LMDZ_WRFmeas about:

1.- Chekcing for NaNs?
2.- Instabilities in thermcell (negative SQRT & robustness IF)

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