Context Navigation

convadj.F @ 2176

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

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

Download in other formats:

Original Format