Context Navigation

convadj.F @ 3300

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

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

Download in other formats:

Original Format