Context Navigation

convadj.F @ 1834

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

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

Download in other formats:

Original Format