Context Navigation

convadj.F @ 3506

Last change on this file since 3506 was 3275, checked in by afalco, 8 months ago
Pluto PCM: Changed _vap to _gas; Included surfprop.F90; callcorrk includes methane AF
File size: 9.2 KB

Line
1	subroutine convadj(ngrid,nlay,nq,ptimestep,
2	& pplay,pplev,ppopsk,
3	& pu,pv,ph,pq,
4	& pdufi,pdvfi,pdhfi,pdqfi,
5	& pduadj,pdvadj,pdhadj,pdqadj)
6
7	USE tracer_h
8	use comcstfi_mod, only: g
9	use callkeys_mod, only: tracer
10
11	implicit none
12
13	!==================================================================
14	!
15	! Purpose
16	! -------
17	! Calculates dry convective adjustment. If one tracer is N2,
18	! we take into account the molecular mass variation (e.g. when
19	! N2 condenses) to trigger convection (F. Forget 01/2005)
20	!
21	! Authors
22	! -------
23	! Original author unknown.
24	! Modif. 2005 by F. Forget.
25	!
26	!==================================================================
27
28	! ------------
29	! Declarations
30	! ------------
31
32
33	! Arguments
34	! ---------
35
36	INTEGER ngrid,nlay
37	REAL ptimestep
38	REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay)
39	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay)
40	REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay)
41	REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay)
42
43	! Tracers
44	integer nq
45	real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
46	real pdqadj(ngrid,nlay,nq)
47
48
49	! Local
50	! -----
51
52	INTEGER ig,i,l,l1,l2,jj
53	INTEGER jcnt, jadrs(ngrid)
54
55	REAL sig(nlay+1),sdsig(nlay),dsig(nlay)
56	REAL zu(ngrid,nlay),zv(ngrid,nlay)
57	REAL zh(ngrid,nlay)
58	REAL zu2(ngrid,nlay),zv2(ngrid,nlay)
59	REAL zh2(ngrid,nlay), zhc(ngrid,nlay)
60	REAL zhm,zsm,zdsm,zum,zvm,zalpha,zhmc
61
62	! Tracers
63	INTEGER iq
64	REAL zq(ngrid,nlay,nq), zq2(ngrid,nlay,nq)
65	REAL zqm(nq),zqn2m
66
67	LOGICAL vtest(ngrid),down
68
69	! for conservation test
70	real masse,cadjncons
71
72	EXTERNAL SCOPY
73
74	! --------------
75	! Initialisation
76	! --------------
77
78	DO l=1,nlay
79	DO ig=1,ngrid
80	zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep
81	zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep
82	zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep
83	ENDDO
84	ENDDO
85
86	if(tracer) then
87	DO iq =1, nq
88	DO l=1,nlay
89	DO ig=1,ngrid
90	zq(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
91	ENDDO
92	ENDDO
93	ENDDO
94	end if
95
96	CALL scopy(ngrid*nlay,zh,1,zh2,1)
97	CALL scopy(ngrid*nlay,zu,1,zu2,1)
98	CALL scopy(ngrid*nlay,zv,1,zv2,1)
99	CALL scopy(ngridnlaynq,zq,1,zq2,1)
100
101	! -----------------------------
102	! Detection of unstable columns
103	! -----------------------------
104	! If ph(above) < ph(below) we set vtest=.true.
105
106	DO ig=1,ngrid
107	vtest(ig)=.false.
108	ENDDO
109
110	CALL scopy(ngrid*nlay,zh2,1,zhc,1)
111
112	! Find out which grid points are convectively unstable
113	DO l=2,nlay
114	DO ig=1,ngrid
115	IF (zhc(ig,l).LT.zhc(ig,l-1)) THEN
116	vtest(ig)=.true.
117	ENDIF
118	ENDDO
119	ENDDO
120
121	! Make a list of them
122	jcnt=0
123	DO ig=1,ngrid
124	IF(vtest(ig)) THEN
125	jcnt=jcnt+1
126	jadrs(jcnt)=ig
127	ENDIF
128	ENDDO
129
130
131	! ---------------------------------------------------------------
132	! Adjustment of the "jcnt" unstable profiles indicated by "jadrs"
133	! ---------------------------------------------------------------
134
135	DO jj = 1, jcnt ! loop on every convective grid point
136
137	i = jadrs(jj)
138
139	! Calculate sigma in this column
140	DO l=1,nlay+1
141	sig(l)=pplev(i,l)/pplev(i,1)
142
143	ENDDO
144	DO l=1,nlay
145	dsig(l)=sig(l)-sig(l+1)
146	sdsig(l)=ppopsk(i,l)*dsig(l)
147	ENDDO
148	l2 = 1
149
150	! Test loop upwards on l2
151
152	DO
153	l2 = l2 + 1
154	IF (l2 .GT. nlay) EXIT
155	IF (zhc(i, l2).LT.zhc(i, l2-1)) THEN
156
157	! l2 is the highest level of the unstable column
158
159	l1 = l2 - 1
160	l = l1
161	zsm = sdsig(l2)
162	zdsm = dsig(l2)
163	zhm = zh2(i, l2)
164
165	! Test loop downwards
166
167	DO
168	zsm = zsm + sdsig(l)
169	zdsm = zdsm + dsig(l)
170	zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
171	zhmc = zhm
172
173	! do we have to extend the column downwards?
174
175	down = .false.
176	IF (l1 .ne. 1) then !-- and then
177	IF (zhmc.LT.zhc(i, l1-1)) then
178	down = .true.
179	END IF
180	END IF
181
182	! this could be a problem...
183
184	if (down) then
185
186	l1 = l1 - 1
187	l = l1
188
189	else
190
191	! can we extend the column upwards?
192
193	if (l2 .eq. nlay) exit
194
195	if (zhc(i, l2+1) .ge. zhmc) exit
196
197	l2 = l2 + 1
198	l = l2
199
200	end if
201
202	enddo
203
204	! New constant profile (average value)
205
206
207	zalpha=0.
208	zum=0.
209	zvm=0.
210	do iq=1,nq
211	zqm(iq) = 0.
212	end do
213	DO l = l1, l2
214	zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
215	zh2(i, l) = zhm
216	! modifs by RDW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
217	zum=zum+dsig(l)*zu2(i,l)
218	zvm=zvm+dsig(l)*zv2(i,l)
219	! zum=zum+dsig(l)*zu(i,l)
220	! zvm=zvm+dsig(l)*zv(i,l)
221	do iq=1,nq
222	zqm(iq) = zqm(iq)+dsig(l)*zq2(i,l,iq)
223	! zqm(iq) = zqm(iq)+dsig(l)*zq(i,l,iq)
224	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
225
226	! to conserve tracers/ KE, we must calculate zum, zvm and zqm using
227	! the up-to-date column values. If we do not do this, there are cases
228	! where convection stops at one level and starts at the next where we
229	! can break conservation of stuff (particularly tracers) significantly.
230
231	end do
232	ENDDO
233	zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
234	zum=zum/(sig(l1)-sig(l2+1))
235	zvm=zvm/(sig(l1)-sig(l2+1))
236	do iq=1,nq
237	zqm(iq) = zqm(iq)/(sig(l1)-sig(l2+1))
238	end do
239
240	IF(zalpha.GT.1.) THEN
241	zalpha=1.
242	ELSE
243	! IF(zalpha.LT.0.) STOP
244	IF(zalpha.LT.1.e-4) zalpha=1.e-4
245	ENDIF
246
247	DO l=l1,l2
248	zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
249	zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
250	do iq=1,nq
251	! zq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq))
252	zq2(i,l,iq)=zqm(iq)
253	end do
254	ENDDO
255
256
257	l2 = l2 + 1
258
259	END IF ! End of l1 to l2 instability treatment
260	! We now continue to test from l2 upwards
261
262	ENDDO ! End of upwards loop on l2
263
264
265	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
266	! check conservation
267	cadjncons=0.0
268
269	if(cadjncons.lt.-1.e-6)then
270	print*,'convadj has just crashed...'
271	print*,'i = ',i
272	print*,'l1 = ',l1
273	print*,'l2 = ',l2
274	print*,'cadjncons = ',cadjncons
275	do l = 1, nlay
276	print*,'dsig = ',dsig(l)
277	end do
278	do l = 1, nlay
279	print*,'dsig = ',dsig(l)
280	end do
281	do l = 1, nlay
282	print*,'sig = ',sig(l)
283	end do
284	do l = 1, nlay
285	print*,'pplay(ig,:) = ',pplay(i,l)
286	end do
287	do l = 1, nlay+1
288	print*,'pplev(ig,:) = ',pplev(i,l)
289	end do
290	do l = 1, nlay
291	print*,'ph(ig,:) = ',ph(i,l)
292	end do
293	do l = 1, nlay
294	print*,'ph(ig,:) = ',ph(i,l)
295	end do
296	do l = 1, nlay
297	print*,'ph(ig,:) = ',ph(i,l)
298	end do
299	do l = 1, nlay
300	print*,'zh(ig,:) = ',zh(i,l)
301	end do
302	do l = 1, nlay
303	print*,'zh2(ig,:) = ',zh2(i,l)
304	end do
305	! do l = 1, nlay
306	! print*,'zq(ig,:,vap) = ',zq(i,l,igcm_h2o_gas)
307	! end do
308	! do l = 1, nlay
309	! print*,'zq2(ig,:,vap) = ',zq2(i,l,igcm_h2o_gas)
310	! end do
311	! print*,'zqm(vap) = ',zqm(igcm_h2o_gas)
312	print*,'jadrs=',jadrs
313
314	call abort
315	endif
316	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
317
318
319	ENDDO
320
321	DO l=1,nlay
322	DO ig=1,ngrid
323	pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep
324	pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep
325	pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep
326	ENDDO
327	ENDDO
328
329	if(tracer) then
330	do iq=1, nq
331	do l=1,nlay
332	DO ig=1,ngrid
333	pdqadj(ig,l,iq)=(zq2(ig,l,iq)-zq(ig,l,iq))/ptimestep
334	end do
335	end do
336	end do
337	end if
338
339
340	! output
341	! if (ngrid.eq.1) then
342	! ig=1
343	! iq =1
344	! write(,)'**** l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq)'
345	! do l=nlay,1,-1
346	! write(,) l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq)
347	! end do
348	! end if
349
350
351	return
352	end

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format