Directory: | ./ |
---|---|
File: | phys/thermcell_dq.f90 |
Date: | 2022-01-11 19:19:34 |
Exec | Total | Coverage | |
---|---|---|---|
Lines: | 35 | 90 | 38.9% |
Branches: | 42 | 94 | 44.7% |
Line | Branch | Exec | Source |
---|---|---|---|
1 | 2880 | 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 | 5760 | 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 | 5760 | 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 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2880 times.
|
2880 | 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 |
2/2✓ Branch 0 taken 112320 times.
✓ Branch 1 taken 2880 times.
|
115200 | do k=1,nlay |
51 |
2/2✓ Branch 0 taken 111646080 times.
✓ Branch 1 taken 112320 times.
|
111761280 | do ig=1,ngrid |
52 | 111646080 | zzm=masse(ig,k)/ptimestep | |
53 | cfl=max(cfl,fm(ig,k)/zzm) | ||
54 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 111646080 times.
|
111758400 | 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 |
4/4✓ Branch 0 taken 112320 times.
✓ Branch 1 taken 2880 times.
✓ Branch 2 taken 111646080 times.
✓ Branch 3 taken 112320 times.
|
111761280 | qold=q |
63 | |||
64 | |||
65 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2880 times.
|
2880 | if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0' |
66 | |||
67 | ! calcul du detrainement | ||
68 |
2/2✓ Branch 0 taken 112320 times.
✓ Branch 1 taken 2880 times.
|
115200 | do k=1,nlay |
69 |
2/2✓ Branch 0 taken 111646080 times.
✓ Branch 1 taken 112320 times.
|
111761280 | do ig=1,ngrid |
70 | 111646080 | 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 |
2/2✓ Branch 0 taken 48966 times.
✓ Branch 1 taken 111597114 times.
|
111646080 | if (detr(ig,k).lt.0.) then |
74 | 48966 | entr(ig,k)=entr(ig,k)-detr(ig,k) | |
75 | 48966 | 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 | 112320 | 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 |
2/2✓ Branch 0 taken 2862720 times.
✓ Branch 1 taken 2880 times.
|
2865600 | do ig=1,ngrid |
90 | 2865600 | qa(ig,1)=q(ig,1) | |
91 | enddo | ||
92 | |||
93 |
2/2✓ Branch 0 taken 109440 times.
✓ Branch 1 taken 2880 times.
|
112320 | do k=2,nlay |
94 |
2/2✓ Branch 0 taken 108783360 times.
✓ Branch 1 taken 109440 times.
|
108895680 | do ig=1,ngrid |
95 |
2/2✓ Branch 0 taken 7820046 times.
✓ Branch 1 taken 100963314 times.
|
108783360 | 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 | 7820046 | & /(fm(ig,k+1)+detr(ig,k)) | |
99 | else | ||
100 | 100963314 | qa(ig,k)=q(ig,k) | |
101 | endif | ||
102 | if (qa(ig,k).lt.0.) then | ||
103 | ! print*,'qa<0!!!' | ||
104 | endif | ||
105 | 109440 | if (q(ig,k).lt.0.) then | |
106 | ! print*,'q<0!!!' | ||
107 | endif | ||
108 | enddo | ||
109 | enddo | ||
110 | |||
111 | ! Plume vertical flux | ||
112 |
2/2✓ Branch 0 taken 106560 times.
✓ Branch 1 taken 2880 times.
|
109440 | do k=2,nlay-1 |
113 |
2/2✓ Branch 0 taken 105920640 times.
✓ Branch 1 taken 106560 times.
|
106030080 | fqa(:,k)=fm(:,k)*qa(:,k-1) |
114 | enddo | ||
115 |
4/4✓ Branch 0 taken 2862720 times.
✓ Branch 1 taken 2880 times.
✓ Branch 2 taken 2862720 times.
✓ Branch 3 taken 2880 times.
|
5728320 | fqa(:,1)=0. ; fqa(:,nlay)=0. |
116 | |||
117 | |||
118 | ! Trace species evolution | ||
119 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2880 times.
|
2880 | 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 |
2/2✓ Branch 0 taken 2880 times.
✓ Branch 1 taken 109440 times.
|
112320 | 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 |
2/2✓ Branch 0 taken 108783360 times.
✓ Branch 1 taken 109440 times.
|
108895680 | & /(1.+fm(:,k)*ptimestep/masse(:,k)) |
133 | ! FH fin de modif. | ||
134 | enddo | ||
135 | endif | ||
136 | |||
137 | ! Tendencies | ||
138 |
2/2✓ Branch 0 taken 2880 times.
✓ Branch 1 taken 112320 times.
|
115200 | do k=1,nlay |
139 |
2/2✓ Branch 0 taken 111646080 times.
✓ Branch 1 taken 112320 times.
|
111761280 | do ig=1,ngrid |
140 | 111646080 | dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep | |
141 | 111758400 | 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 | niter=1 | ||
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 | |||
263 | ! Schema avec advection sur plus qu'une maille. | ||
264 | ✗ | zzm=masse(ig,k)/ztimestep | |
265 | ✗ | if (fm(ig,k)>zzm) then | |
266 | ✗ | wqd(ig,k)=zzm*q(ig,k)+(fm(ig,k)-zzm)*q(ig,k+1) | |
267 | else | ||
268 | ✗ | wqd(ig,k)=fm(ig,k)*q(ig,k) | |
269 | endif | ||
270 | |||
271 | ✗ | if (wqd(ig,k).lt.0.) then | |
272 | ! print*,'wqd<0!!!' | ||
273 | endif | ||
274 | enddo | ||
275 | enddo | ||
276 | ✗ | do ig=1,ngrid | |
277 | ✗ | wqd(ig,1)=0. | |
278 | ✗ | wqd(ig,nlay+1)=0. | |
279 | enddo | ||
280 | |||
281 | |||
282 | ! Calcul des tendances | ||
283 | ✗ | do k=1,nlay | |
284 | ✗ | do ig=1,ngrid | |
285 | q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k) & | ||
286 | & -wqd(ig,k)+wqd(ig,k+1)) & | ||
287 | ✗ | & *ztimestep/masse(ig,k) | |
288 | ! if (dq(ig,k).lt.0.) then | ||
289 | ! print*,'dq<0!!!' | ||
290 | ! endif | ||
291 | enddo | ||
292 | enddo | ||
293 | |||
294 | |||
295 | enddo | ||
296 | |||
297 | |||
298 | ! Calcul des tendances | ||
299 | ✗ | do k=1,nlay | |
300 | ✗ | do ig=1,ngrid | |
301 | ✗ | dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep | |
302 | ✗ | q(ig,k)=qold(ig,k) | |
303 | enddo | ||
304 | enddo | ||
305 | |||
306 | ✗ | return | |
307 | end | ||
308 |