LMDZ
thermcell_dq.F90
Go to the documentation of this file.
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 
145 return
146 end
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 
217 do 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 
311 enddo
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
dq
do llm!au dessus on relaxe vers profil init!on fait l hypothese que dans ce il n y a plus d eau liq au dessus!donc la relaxation en thetal et qt devient relaxation en tempe et qv l dq1 relax dq(l, 1)
prt_level
!FH On elimine toutes les clefs physiques dans la dynamique prt_level
Definition: write_paramLMDZ_dyn.h:9
thermcell_dq
subroutine thermcell_dq(ngrid, nlay, impl, ptimestep, fm, entr, masse, q, dq, qa, lev_out)
Definition: thermcell_dq.F90:3
print_control_mod
Definition: print_control_mod.F90:2
abort_physic
subroutine abort_physic(modname, message, ierr)
Definition: abort_physic.F90:3
thermcell_dq_o
subroutine thermcell_dq_o(ngrid, nlay, impl, ptimestep, fm, entr, masse, q, dq, qa, lev_out)
Definition: thermcell_dq.F90:154