source: LMDZ6/trunk/libf/phylmd/lmdz_thermcell_dq.f90 @ 5435

Last change on this file since 5435 was 5434, checked in by fhourdin, 12 hours ago

Superessing CPP in lmdz_*

Not possible for lmdz_thermcell_main because of isotopes
  • 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.9 KB
Line 
1MODULE lmdz_thermcell_dq
2CONTAINS
3
4      subroutine thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr,  &
5     &           masse,q,dq,qa,lev_out)
6      USE print_control_mod, ONLY: prt_level
7
8      implicit none
9
10!=======================================================================
11!
12!   Calcul du transport verticale dans la couche limite en presence
13!   de "thermiques" explicitement representes
14!   calcul du dq/dt une fois qu'on connait les ascendances
15!
16! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
17!  Introduction of an implicit computation of vertical advection in
18!  the environment of thermal plumes in thermcell_dq
19!  impl =     0 : explicit, 1 : implicit, -1 : old version
20!
21!=======================================================================
22
23! arguments
24      integer, intent(in) :: ngrid,nlay,impl
25      real, intent(in) :: ptimestep
26      real, intent(in), dimension(ngrid,nlay) :: masse
27      real, intent(inout), dimension(ngrid,nlay) :: entr,q
28      real, intent(in), dimension(ngrid,nlay+1) :: fm
29      real, intent(out), dimension(ngrid,nlay) :: dq,qa
30      integer, intent(in) :: lev_out                           ! niveau pour les print
31
32! Local
33      real, dimension(ngrid,nlay) :: detr,qold
34      real, dimension(ngrid,nlay+1) :: wqd,fqa
35      real zzm
36      integer ig,k
37      real cfl
38
39      integer niter,iter
40      CHARACTER (LEN=20) :: modname='thermcell_dq'
41      CHARACTER (LEN=80) :: abort_message
42
43
44! Old explicite scheme
45if (impl<=-1) then
46
47         call thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,  &
48     &           masse,q,dq,qa,lev_out)
49
50else
51 
52
53! Calcul du critere CFL pour l'advection dans la subsidence
54      cfl = 0.
55      do k=1,nlay
56         do ig=1,ngrid
57            zzm=masse(ig,k)/ptimestep
58            cfl=max(cfl,fm(ig,k)/zzm)
59            if (entr(ig,k).gt.zzm) then
60               print*,'entr*dt>m,1',k,entr(ig,k)*ptimestep,masse(ig,k)
61               abort_message = 'entr dt > m, 1st'
62               CALL abort_physic (modname,abort_message,1)
63            endif
64         enddo
65      enddo
66
67      qold=q
68
69
70      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
71
72!   calcul du detrainement
73      do k=1,nlay
74         do ig=1,ngrid
75            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
76!           print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
77!test
78            if (detr(ig,k).lt.0.) then
79               entr(ig,k)=entr(ig,k)-detr(ig,k)
80               detr(ig,k)=0.
81!               print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
82!     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
83            endif
84            if (fm(ig,k+1).lt.0.) then
85!               print*,'fm2<0!!!'
86            endif
87            if (entr(ig,k).lt.0.) then
88!               print*,'entr2<0!!!'
89            endif
90         enddo
91      enddo
92
93! Computation of tracer concentrations in the ascending plume
94      do ig=1,ngrid
95         qa(ig,1)=q(ig,1)
96      enddo
97
98      do k=2,nlay
99         do ig=1,ngrid
100            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
101     &         1.e-5*masse(ig,k)) then
102         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
103     &         /(fm(ig,k+1)+detr(ig,k))
104            else
105               qa(ig,k)=q(ig,k)
106            endif
107            if (qa(ig,k).lt.0.) then
108!               print*,'qa<0!!!'
109            endif
110            if (q(ig,k).lt.0.) then
111!               print*,'q<0!!!'
112            endif
113         enddo
114      enddo
115
116! Plume vertical flux
117      do k=2,nlay-1
118         fqa(:,k)=fm(:,k)*qa(:,k-1)
119      enddo
120      fqa(:,1)=0. ; fqa(:,nlay)=0.
121
122
123! Trace species evolution
124   if (impl==0) then
125      do k=1,nlay-1
126         q(:,k)=q(:,k)+(fqa(:,k)-fqa(:,k+1)-fm(:,k)*q(:,k)+fm(:,k+1)*q(:,k+1)) &
127     &               *ptimestep/masse(:,k)
128      enddo
129   else
130      do k=nlay-1,1,-1
131! FH debut de modif : le calcul ci dessous modifiait numériquement
132! la concentration quand le flux de masse etait nul car on divisait
133! puis multipliait par masse/ptimestep.
134!        q(:,k)=(masse(:,k)*q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1)) &
135!    &               /(fm(:,k)+masse(:,k)/ptimestep)
136         q(:,k)=(q(:,k)+ptimestep/masse(:,k)*(fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1))) &
137      &               /(1.+fm(:,k)*ptimestep/masse(:,k))
138! FH fin de modif.
139      enddo
140   endif
141
142! Tendencies
143      do k=1,nlay
144         do ig=1,ngrid
145            dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
146            q(ig,k)=qold(ig,k)
147         enddo
148      enddo
149
150endif ! impl=-1
151RETURN
152end subroutine thermcell_dq
153
154
155!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
156! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations
157!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
158
159      subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,  &
160     &           masse,q,dq,qa,lev_out)
161      USE print_control_mod, ONLY: prt_level
162      USE lmdz_thermcell_ini, ONLY : thermals_subsid_advect_scheme,thermals_subsid_advect_more_than_one
163      implicit none
164
165!=======================================================================
166!
167!   Calcul du transport verticale dans la couche limite en presence
168!   de "thermiques" explicitement representes
169!   calcul du dq/dt une fois qu'on connait les ascendances
170!
171!=======================================================================
172
173      integer ngrid,nlay,impl
174
175      real ptimestep
176      real masse(ngrid,nlay),fm(ngrid,nlay+1)
177      real entr(ngrid,nlay)
178      real q(ngrid,nlay)
179      real dq(ngrid,nlay)
180      integer lev_out                           ! niveau pour les print
181
182      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
183
184      real zzm
185
186      integer ig,k
187      real cfl
188
189      real qold(ngrid,nlay)
190      real ztimestep
191      integer niter,iter
192      CHARACTER (LEN=20) :: modname='thermcell_dq'
193      CHARACTER (LEN=80) :: abort_message
194
195
196
197! Calcul du critere CFL pour l'advection dans la subsidence
198      cfl = 0.
199      do k=1,nlay
200         do ig=1,ngrid
201            zzm=masse(ig,k)/ptimestep
202            cfl=max(cfl,fm(ig,k)/zzm)
203            if (entr(ig,k).gt.zzm) then
204               print*,'entr*dt>m,2',k,entr(ig,k)*ptimestep,masse(ig,k)
205               abort_message = 'entr dt > m, 2nd'
206               CALL abort_physic (modname,abort_message,1)
207            endif
208         enddo
209      enddo
210
211
212!     niter=int(cfl)+1 ! pour tourner avec un CFL différent en splitant
213      niter=1
214
215      ztimestep=ptimestep/niter
216      qold=q
217
218
219do iter=1,niter
220      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
221
222!   calcul du detrainement
223      do k=1,nlay
224         do ig=1,ngrid
225            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
226!           print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
227!test
228            if (detr(ig,k).lt.0.) then
229               entr(ig,k)=entr(ig,k)-detr(ig,k)
230               detr(ig,k)=0.
231!               print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
232!     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
233            endif
234            if (fm(ig,k+1).lt.0.) then
235!               print*,'fm2<0!!!'
236            endif
237            if (entr(ig,k).lt.0.) then
238!               print*,'entr2<0!!!'
239            endif
240         enddo
241      enddo
242
243!   calcul de la valeur dans les ascendances
244      do ig=1,ngrid
245         qa(ig,1)=q(ig,1)
246      enddo
247
248      do k=2,nlay
249         do ig=1,ngrid
250            if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt.  &
251     &         1.e-5*masse(ig,k)) then
252         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
253     &         /(fm(ig,k+1)+detr(ig,k))
254            else
255               qa(ig,k)=q(ig,k)
256            endif
257            if (qa(ig,k).lt.0.) then
258!               print*,'qa<0!!!'
259            endif
260            if (q(ig,k).lt.0.) then
261!               print*,'q<0!!!'
262            endif
263         enddo
264      enddo
265
266! Calcul du flux subsident
267
268      if ( thermals_subsid_advect_scheme == 'center' ) then
269         do k=2,nlay
270            do ig=1,ngrid
271                wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
272            enddo
273         enddo
274      else ! upstream scheme (default and recomanded)
275         do k=2,nlay
276            do ig=1,ngrid
277               zzm=masse(ig,k)/ztimestep
278               if ( fm(ig,k)<=zzm .or. thermals_subsid_advect_more_than_one == 0 ) then
279                   wqd(ig,k)=fm(ig,k)*q(ig,k)
280               else
281                  wqd(ig,k)=zzm*q(ig,k)+(fm(ig,k)-zzm)*q(ig,k+1)
282               endif
283            enddo
284         enddo
285      endif
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 subroutine thermcell_dq_o
318END MODULE lmdz_thermcell_dq
Note: See TracBrowser for help on using the repository browser.