Context Navigation

convadj.F @ 2276

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

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

Download in other formats:

Original Format