source: LMDZ6/trunk/libf/phylmd/lmdz_thermcell_dq.F90 @ 5390

Last change on this file since 5390 was 5390, checked in by yann meurdesoif, 8 days ago
  • Remove UTF8 character that inihibit fortran parsing with GPU morphosis
  • Add missing END SUBROUTINE instead of simple END, that inhibit correct parsing with regulat expression parser (quick and dirty parsing)

YM
  • 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      implicit none
163
164!=======================================================================
165!
166!   Calcul du transport verticale dans la couche limite en presence
167!   de "thermiques" explicitement representes
168!   calcul du dq/dt une fois qu'on connait les ascendances
169!
170!=======================================================================
171
172      integer ngrid,nlay,impl
173
174      real ptimestep
175      real masse(ngrid,nlay),fm(ngrid,nlay+1)
176      real entr(ngrid,nlay)
177      real q(ngrid,nlay)
178      real dq(ngrid,nlay)
179      integer lev_out                           ! niveau pour les print
180
181      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
182
183      real zzm
184
185      integer ig,k
186      real cfl
187
188      real qold(ngrid,nlay)
189      real ztimestep
190      integer niter,iter
191      CHARACTER (LEN=20) :: modname='thermcell_dq'
192      CHARACTER (LEN=80) :: abort_message
193
194
195
196! Calcul du critere CFL pour l'advection dans la subsidence
197      cfl = 0.
198      do k=1,nlay
199         do ig=1,ngrid
200            zzm=masse(ig,k)/ptimestep
201            cfl=max(cfl,fm(ig,k)/zzm)
202            if (entr(ig,k).gt.zzm) then
203               print*,'entr*dt>m,2',k,entr(ig,k)*ptimestep,masse(ig,k)
204               abort_message = 'entr dt > m, 2nd'
205               CALL abort_physic (modname,abort_message,1)
206            endif
207         enddo
208      enddo
209
210!IM 090508     print*,'CFL CFL CFL CFL ',cfl
211
212#undef CFL
213#ifdef CFL
214! On subdivise le calcul en niter pas de temps.
215      niter=int(cfl)+1
216#else
217      niter=1
218#endif
219
220      ztimestep=ptimestep/niter
221      qold=q
222
223
224do iter=1,niter
225      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
226
227!   calcul du detrainement
228      do k=1,nlay
229         do ig=1,ngrid
230            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
231!           print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
232!test
233            if (detr(ig,k).lt.0.) then
234               entr(ig,k)=entr(ig,k)-detr(ig,k)
235               detr(ig,k)=0.
236!               print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
237!     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
238            endif
239            if (fm(ig,k+1).lt.0.) then
240!               print*,'fm2<0!!!'
241            endif
242            if (entr(ig,k).lt.0.) then
243!               print*,'entr2<0!!!'
244            endif
245         enddo
246      enddo
247
248!   calcul de la valeur dans les ascendances
249      do ig=1,ngrid
250         qa(ig,1)=q(ig,1)
251      enddo
252
253      do k=2,nlay
254         do ig=1,ngrid
255            if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt.  &
256     &         1.e-5*masse(ig,k)) then
257         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
258     &         /(fm(ig,k+1)+detr(ig,k))
259            else
260               qa(ig,k)=q(ig,k)
261            endif
262            if (qa(ig,k).lt.0.) then
263!               print*,'qa<0!!!'
264            endif
265            if (q(ig,k).lt.0.) then
266!               print*,'q<0!!!'
267            endif
268         enddo
269      enddo
270
271! Calcul du flux subsident
272
273      do k=2,nlay
274         do ig=1,ngrid
275#undef centre
276#ifdef centre
277             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
278#else
279
280#define plusqueun
281#ifdef plusqueun
282! Schema avec advection sur plus qu'une maille.
283            zzm=masse(ig,k)/ztimestep
284            if (fm(ig,k)>zzm) then
285               wqd(ig,k)=zzm*q(ig,k)+(fm(ig,k)-zzm)*q(ig,k+1)
286            else
287               wqd(ig,k)=fm(ig,k)*q(ig,k)
288            endif
289#else
290            wqd(ig,k)=fm(ig,k)*q(ig,k)
291#endif
292#endif
293
294            if (wqd(ig,k).lt.0.) then
295!               print*,'wqd<0!!!'
296            endif
297         enddo
298      enddo
299      do ig=1,ngrid
300         wqd(ig,1)=0.
301         wqd(ig,nlay+1)=0.
302      enddo
303     
304
305! Calcul des tendances
306      do k=1,nlay
307         do ig=1,ngrid
308            q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
309     &               -wqd(ig,k)+wqd(ig,k+1))  &
310     &               *ztimestep/masse(ig,k)
311!            if (dq(ig,k).lt.0.) then
312!               print*,'dq<0!!!'
313!            endif
314         enddo
315      enddo
316
317
318enddo
319
320
321! Calcul des tendances
322      do k=1,nlay
323         do ig=1,ngrid
324            dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
325            q(ig,k)=qold(ig,k)
326         enddo
327      enddo
328
329      return
330    end subroutine thermcell_dq_o
331END MODULE lmdz_thermcell_dq
Note: See TracBrowser for help on using the repository browser.