Context Navigation

convadj.F @ 2156

Last change on this file since 2156 was 1779, checked in by aslmd, 7 years ago
LMDZ.MARS (purely comments) marked the absolute firstcalls not supposed to change with runtime (e.g. not domain-related). this is most of them. those firstcall can stay local and do not need to be linked with the caller's general firstcall.
File size: 11.6 KB

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

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

Download in other formats:

Original Format