Context Navigation

convadj.F @ 374

Last change on this file since 374 was 253, checked in by emillour, 13 years ago

Generic GCM

Massive update to version 0.7

EM+RW

File size: 11.6 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,
6	& pdqadj)
7
8	implicit none
9
10	!==================================================================
11	!
12	! Purpose
13	! -------
14	! Calculates dry convective adjustment. If one tracer is CO2,
15	! we take into account the molecular mass variation (e.g. when
16	! CO2 condenses) to trigger convection (F. Forget 01/2005)
17	!
18	! Authors
19	! -------
20	! Original author unknown.
21	! Modif. 2005 by F. Forget.
22	!
23	!==================================================================
24
25	! ------------
26	! Declarations
27	! ------------
28
29	#include "dimensions.h"
30	#include "dimphys.h"
31	#include "comcstfi.h"
32	#include "callkeys.h"
33	#include "tracer.h"
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(ngridmx)
57
58	REAL sig(nlayermx+1),sdsig(nlayermx),dsig(nlayermx)
59	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
60	REAL zh(ngridmx,nlayermx)
61	REAL zu2(ngridmx,nlayermx),zv2(ngridmx,nlayermx)
62	REAL zh2(ngridmx,nlayermx), zhc(ngridmx,nlayermx)
63	REAL zhm,zsm,zdsm,zum,zvm,zalpha,zhmc
64
65	! Tracers
66	INTEGER iq,ico2
67	save ico2
68	REAL zq(ngridmx,nlayermx,nqmx), zq2(ngridmx,nlayermx,nqmx)
69	REAL zqm(nqmx),zqco2m
70	real m_co2, m_noco2, A , B
71	save A, B
72
73	real mtot1, mtot2 , mm1, mm2
74	integer l1ref, l2ref
75	LOGICAL vtest(ngridmx),down,firstcall
76	save firstcall
77	data firstcall/.true./
78
79	! for conservation test
80	real masse,cadjncons
81
82	EXTERNAL SCOPY
83
84	! --------------
85	! Initialisation
86	! --------------
87
88	IF (firstcall) THEN
89	IF(ngrid.NE.ngridmx) THEN
90	PRINT*
91	PRINT*,'STOP in convadj'
92	PRINT*,'ngrid =',ngrid
93	PRINT*,'ngridmx =',ngridmx
94	ENDIF
95	ico2=0
96	if (tracer) then
97	! Prepare Special treatment if one of the tracers is CO2 gas
98	do iq=1,nqmx
99	if (noms(iq).eq."co2") then
100	print*,'dont go there'
101	stop
102	ico2=iq
103	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
104	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
105	! Compute A and B coefficient use to compute
106	! mean molecular mass Mair defined by
107	! 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
108	! 1/Mair = A*q(ico2) + B
109	A =(1/m_co2 - 1/m_noco2)
110	B=1/m_noco2
111	end if
112	enddo
113	endif
114	firstcall=.false.
115	ENDIF ! of IF (firstcall)
116
117	DO l=1,nlay
118	DO ig=1,ngrid
119	zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep
120	zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep
121	zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep
122	ENDDO
123	ENDDO
124
125	if(tracer) then
126	DO iq =1, nq
127	DO l=1,nlay
128	DO ig=1,ngrid
129	zq(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
130	ENDDO
131	ENDDO
132	ENDDO
133	end if
134
135	CALL scopy(ngrid*nlay,zh,1,zh2,1)
136	CALL scopy(ngrid*nlay,zu,1,zu2,1)
137	CALL scopy(ngrid*nlay,zv,1,zv2,1)
138	CALL scopy(ngridnlaynq,zq,1,zq2,1)
139
140	! -----------------------------
141	! Detection of unstable columns
142	! -----------------------------
143	! If ph(above) < ph(below) we set vtest=.true.
144
145	DO ig=1,ngrid
146	vtest(ig)=.false.
147	ENDDO
148
149	if (ico2.ne.0) then
150	! Special case if one of the tracers is CO2 gas
151	DO l=1,nlay
152	DO ig=1,ngrid
153	zhc(ig,l) = zh2(ig,l)(Azq2(ig,l,ico2)+B)
154	ENDDO
155	ENDDO
156	else
157	CALL scopy(ngrid*nlay,zh2,1,zhc,1)
158	end if
159
160	! Find out which grid points are convectively unstable
161	DO l=2,nlay
162	DO ig=1,ngrid
163	IF(zhc(ig,l).LT.zhc(ig,l-1)) vtest(ig)=.true.
164	ENDDO
165	ENDDO
166
167	! Make a list of them
168	jcnt=0
169	DO ig=1,ngrid
170	IF(vtest(ig)) THEN
171	jcnt=jcnt+1
172	jadrs(jcnt)=ig
173	ENDIF
174	ENDDO
175
176
177	! ---------------------------------------------------------------
178	! Adjustment of the "jcnt" unstable profiles indicated by "jadrs"
179	! ---------------------------------------------------------------
180
181	DO jj = 1, jcnt ! loop on every convective grid point
182
183	i = jadrs(jj)
184
185	! Calculate sigma in this column
186	DO l=1,nlay+1
187	sig(l)=pplev(i,l)/pplev(i,1)
188
189	ENDDO
190	DO l=1,nlay
191	dsig(l)=sig(l)-sig(l+1)
192	sdsig(l)=ppopsk(i,l)*dsig(l)
193	ENDDO
194	l2 = 1
195
196	! Test loop upwards on l2
197
198	DO
199	l2 = l2 + 1
200	IF (l2 .GT. nlay) EXIT
201	IF (zhc(i, l2) .LT. zhc(i, l2-1)) THEN
202
203	! l2 is the highest level of the unstable column
204
205	l1 = l2 - 1
206	l = l1
207	zsm = sdsig(l2)
208	zdsm = dsig(l2)
209	zhm = zh2(i, l2)
210	if(ico2.ne.0) zqco2m = zq2(i,l2,ico2)
211
212	! Test loop downwards
213
214	DO
215	zsm = zsm + sdsig(l)
216	zdsm = zdsm + dsig(l)
217	zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
218	if(ico2.ne.0) then
219	zqco2m =
220	& zqco2m + dsig(l) * (zq2(i,l,ico2) - zqco2m) / zdsm
221	zhmc = zhm(Azqco2m+B)
222	else
223	zhmc = zhm
224	end if
225
226	! do we have to extend the column downwards?
227
228	down = .false.
229	IF (l1 .ne. 1) then !-- and then
230	IF (zhmc .lt. zhc(i, l1-1)) then
231	down = .true.
232	END IF
233	END IF
234
235	! this could be a problem...
236
237	if (down) then
238
239	l1 = l1 - 1
240	l = l1
241
242	else
243
244	! can we extend the column upwards?
245
246	if (l2 .eq. nlay) exit
247
248	if (zhc(i, l2+1) .ge. zhmc) exit
249
250	l2 = l2 + 1
251	l = l2
252
253	end if
254
255	enddo
256
257	! New constant profile (average value)
258
259
260	zalpha=0.
261	zum=0.
262	zvm=0.
263	do iq=1,nq
264	zqm(iq) = 0.
265	end do
266	DO l = l1, l2
267	if(ico2.ne.0) then
268	zalpha=zalpha+
269	& ABS(zhc(i,l)/(A+Bzqco2m) -zhm)dsig(l)
270	else
271	zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
272	endif
273	zh2(i, l) = zhm
274	! modifs by RDW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
275	zum=zum+dsig(l)*zu2(i,l)
276	zvm=zvm+dsig(l)*zv2(i,l)
277	! zum=zum+dsig(l)*zu(i,l)
278	! zvm=zvm+dsig(l)*zv(i,l)
279	do iq=1,nq
280	zqm(iq) = zqm(iq)+dsig(l)*zq2(i,l,iq)
281	! zqm(iq) = zqm(iq)+dsig(l)*zq(i,l,iq)
282	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
283
284	! to conserve tracers/ KE, we must calculate zum, zvm and zqm using
285	! the up-to-date column values. If we do not do this, there are cases
286	! where convection stops at one level and starts at the next where we
287	! can break conservation of stuff (particularly tracers) significantly.
288
289	end do
290	ENDDO
291	zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
292	zum=zum/(sig(l1)-sig(l2+1))
293	zvm=zvm/(sig(l1)-sig(l2+1))
294	do iq=1,nq
295	zqm(iq) = zqm(iq)/(sig(l1)-sig(l2+1))
296	end do
297
298	IF(zalpha.GT.1.) THEN
299	zalpha=1.
300	ELSE
301	! IF(zalpha.LT.0.) STOP
302	IF(zalpha.LT.1.e-4) zalpha=1.e-4
303	ENDIF
304
305	DO l=l1,l2
306	zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
307	zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
308	do iq=1,nq
309	! zq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq))
310	zq2(i,l,iq)=zqm(iq)
311	end do
312	ENDDO
313	if (ico2.ne.0) then
314	DO l=l1, l2
315	zhc(i,l) = zh2(i,l)(Azq2(i,l,ico2)+B)
316	ENDDO
317	end if
318
319
320	l2 = l2 + 1
321
322	END IF ! End of l1 to l2 instability treatment
323	! We now continue to test from l2 upwards
324
325	ENDDO ! End of upwards loop on l2
326
327
328	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
329	! check conservation
330	cadjncons=0.0
331	if(water)then
332	do l = 1, nlay
333	masse = (pplev(i,l) - pplev(i,l+1))/g
334	iq = igcm_h2o_vap
335	cadjncons = cadjncons +
336	& masse*(zq2(i,l,iq)-zq(i,l,iq))/ptimestep
337	end do
338	endif
339
340	if(cadjncons.lt.-1.e-6)then
341	print*,'convadj has just crashed...'
342	print*,'i = ',i
343	print*,'l1 = ',l1
344	print*,'l2 = ',l2
345	print*,'cadjncons = ',cadjncons
346	do l = 1, nlay
347	print*,'dsig = ',dsig(l)
348	end do
349	do l = 1, nlay
350	print*,'dsig = ',dsig(l)
351	end do
352	do l = 1, nlay
353	print*,'sig = ',sig(l)
354	end do
355	do l = 1, nlay
356	print*,'pplay(ig,:) = ',pplay(i,l)
357	end do
358	do l = 1, nlay+1
359	print*,'pplev(ig,:) = ',pplev(i,l)
360	end do
361	do l = 1, nlay
362	print*,'ph(ig,:) = ',ph(i,l)
363	end do
364	do l = 1, nlay
365	print*,'ph(ig,:) = ',ph(i,l)
366	end do
367	do l = 1, nlay
368	print*,'ph(ig,:) = ',ph(i,l)
369	end do
370	do l = 1, nlay
371	print*,'zh(ig,:) = ',zh(i,l)
372	end do
373	do l = 1, nlay
374	print*,'zh2(ig,:) = ',zh2(i,l)
375	end do
376	do l = 1, nlay
377	print*,'zq(ig,:,vap) = ',zq(i,l,igcm_h2o_vap)
378	end do
379	do l = 1, nlay
380	print*,'zq2(ig,:,vap) = ',zq2(i,l,igcm_h2o_vap)
381	end do
382	print*,'zqm(vap) = ',zqm(igcm_h2o_vap)
383	print*,'jadrs=',jadrs
384
385	call abort
386	endif
387	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
388
389
390	ENDDO
391
392	DO l=1,nlay
393	DO ig=1,ngrid
394	pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep
395	pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep
396	pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep
397	ENDDO
398	ENDDO
399
400	if(tracer) then
401	do iq=1, nq
402	do l=1,nlay
403	do ig=1, ngrid
404	pdqadj(ig,l,iq)=(zq2(ig,l,iq)-zq(ig,l,iq))/ptimestep
405	end do
406	end do
407	end do
408	end if
409
410
411	! output
412	! if (ngrid.eq.1) then
413	! ig=1
414	! iq =1
415	! write(,)'**** l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq)'
416	! do l=nlay,1,-1
417	! write(,) l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq)
418	! end do
419	! end if
420
421
422	return
423	end

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

Download in other formats:

Original Format