My Project
 All Classes Files Functions Variables Macros
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  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 
41 ! Old explicite scheme
42  if (impl==-1) then
43  call thermcell_dq_o(ngrid,nlay,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 ',entr(ig,k)*ptimestep,masse(ig,k)
56  abort_message = ''
57  CALL abort_gcm(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  q(:,k)=(masse(:,k)*q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1)) &
127  & /(fm(:,k)+masse(:,k)/ptimestep)
128  enddo
129  endif
130 
131 ! Tendencies
132  do k=1,nlay
133  do ig=1,ngrid
134  dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
135  q(ig,k)=qold(ig,k)
136  enddo
137  enddo
138 
139 return
140 end
141 
142 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
143 ! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations
144 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
145 
146  subroutine thermcell_dq_o(ngrid,nlay,ptimestep,fm,entr, &
147  & masse,q,dq,qa,lev_out)
148  implicit none
149 
150 #include "iniprint.h"
151 !=======================================================================
152 !
153 ! Calcul du transport verticale dans la couche limite en presence
154 ! de "thermiques" explicitement representes
155 ! calcul du dq/dt une fois qu'on connait les ascendances
156 !
157 !=======================================================================
158 
159  integer ngrid,nlay
160 
161  real ptimestep
162  real masse(ngrid,nlay),fm(ngrid,nlay+1)
163  real entr(ngrid,nlay)
164  real q(ngrid,nlay)
165  real dq(ngrid,nlay)
166  integer lev_out ! niveau pour les print
167 
168  real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
169 
170  real zzm
171 
172  integer ig,k
173  real cfl
174 
175  real qold(ngrid,nlay)
176  real ztimestep
177  integer niter,iter
178  CHARACTER (LEN=20) :: modname='thermcell_dq'
179  CHARACTER (LEN=80) :: abort_message
180 
181 
182 
183 ! Calcul du critere CFL pour l'advection dans la subsidence
184  cfl = 0.
185  do k=1,nlay
186  do ig=1,ngrid
187  zzm=masse(ig,k)/ptimestep
188  cfl=max(cfl,fm(ig,k)/zzm)
189  if (entr(ig,k).gt.zzm) then
190  print*,'entr dt > m ',entr(ig,k)*ptimestep,masse(ig,k)
191  abort_message = ''
192  CALL abort_gcm(modname,abort_message,1)
193  endif
194  enddo
195  enddo
196 
197 !IM 090508 print*,'CFL CFL CFL CFL ',cfl
198 
199 #undef CFL
200 #ifdef CFL
201 ! On subdivise le calcul en niter pas de temps.
202  niter=int(cfl)+1
203 #else
204  niter=1
205 #endif
206 
207  ztimestep=ptimestep/niter
208  qold=q
209 
210 
211 do iter=1,niter
212  if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
213 
214 ! calcul du detrainement
215  do k=1,nlay
216  do ig=1,ngrid
217  detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
218 ! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
219 !test
220  if (detr(ig,k).lt.0.) then
221  entr(ig,k)=entr(ig,k)-detr(ig,k)
222  detr(ig,k)=0.
223 ! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
224 ! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
225  endif
226  if (fm(ig,k+1).lt.0.) then
227 ! print*,'fm2<0!!!'
228  endif
229  if (entr(ig,k).lt.0.) then
230 ! print*,'entr2<0!!!'
231  endif
232  enddo
233  enddo
234 
235 ! calcul de la valeur dans les ascendances
236  do ig=1,ngrid
237  qa(ig,1)=q(ig,1)
238  enddo
239 
240  do k=2,nlay
241  do ig=1,ngrid
242  if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt. &
243  & 1.e-5*masse(ig,k)) then
244  qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k)) &
245  & /(fm(ig,k+1)+detr(ig,k))
246  else
247  qa(ig,k)=q(ig,k)
248  endif
249  if (qa(ig,k).lt.0.) then
250 ! print*,'qa<0!!!'
251  endif
252  if (q(ig,k).lt.0.) then
253 ! print*,'q<0!!!'
254  endif
255  enddo
256  enddo
257 
258 ! Calcul du flux subsident
259 
260  do k=2,nlay
261  do ig=1,ngrid
262 #undef centre
263 #ifdef centre
264  wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
265 #else
266 
267 #define plusqueun
268 #ifdef plusqueun
269 ! Schema avec advection sur plus qu'une maille.
270  zzm=masse(ig,k)/ztimestep
271  if (fm(ig,k)>zzm) then
272  wqd(ig,k)=zzm*q(ig,k)+(fm(ig,k)-zzm)*q(ig,k+1)
273  else
274  wqd(ig,k)=fm(ig,k)*q(ig,k)
275  endif
276 #else
277  wqd(ig,k)=fm(ig,k)*q(ig,k)
278 #endif
279 #endif
280 
281  if (wqd(ig,k).lt.0.) then
282 ! print*,'wqd<0!!!'
283  endif
284  enddo
285  enddo
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 
305 enddo
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