source: LMDZ6/branches/LMDZ-ECRAD/libf/phylmd/thermcell_dq.F90 @ 4109

Last change on this file since 4109 was 2311, checked in by Ehouarn Millour, 9 years ago

Further modifications to enforce physics/dynamics separation:

  • moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
  • created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
  • created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
  • reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

EM

  • 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      USE print_control_mod, ONLY: prt_level
4      implicit none
5
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_physic (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      USE print_control_mod, ONLY: prt_level
155      implicit none
156
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_physic (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.