Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

convadj.F @ 38

Last change on this file since 38 was 38, checked in by emillour, 14 years ago
Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs. EM
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	! Modif. 2010 by R. Wordsworth
23	!
24	!==================================================================
25
26	! ------------
27	! Declarations
28	! ------------
29
30	#include "dimensions.h"
31	#include "dimphys.h"
32	#include "comcstfi.h"
33	#include "callkeys.h"
34	#include "tracer.h"
35
36
37	! Arguments
38	! ---------
39
40	INTEGER,intent(in) :: ngrid,nlay
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(ngridmx)
63
64	REAL sig(nlayermx+1),sdsig(nlayermx),dsig(nlayermx)
65	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
66	REAL zh(ngridmx,nlayermx)
67	REAL zu2(ngridmx,nlayermx),zv2(ngridmx,nlayermx)
68	REAL zh2(ngridmx,nlayermx), zhc(ngridmx,nlayermx)
69	REAL zhm,zsm,zdsm,zum,zvm,zalpha,zhmc
70
71	! Tracers
72	INTEGER iq,ico2
73	save ico2
74	REAL zq(ngridmx,nlayermx,nqmx), zq2(ngridmx,nlayermx,nqmx)
75	REAL zqm(nqmx),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(ngridmx),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	IF (firstcall) THEN
93	IF(ngrid.NE.ngridmx) THEN
94	PRINT*
95	PRINT*,'STOP in convadj'
96	PRINT*,'ngrid =',ngrid
97	PRINT*,'ngridmx =',ngridmx
98	ENDIF
99	ico2=0
100	if (tracer) then
101	! Prepare Special treatment if one of the tracers is CO2 gas
102	do iq=1,nqmx
103	if (noms(iq).eq."co2") then
104	ico2=iq
105	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
106	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
107	! Compute A and B coefficient use to compute
108	! mean molecular mass Mair defined by
109	! 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
110	! 1/Mair = A*q(ico2) + B
111	A =(1/m_co2 - 1/m_noco2)
112	B=1/m_noco2
113	end if
114	enddo
115	endif
116	firstcall=.false.
117	ENDIF ! of IF (firstcall)
118
119	DO l=1,nlay
120	DO ig=1,ngrid
121	zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep
122	zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep
123	zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep
124	ENDDO
125	ENDDO
126
127	if(tracer) then
128	DO iq =1, nq
129	DO l=1,nlay
130	DO ig=1,ngrid
131	zq(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
132	ENDDO
133	ENDDO
134	ENDDO
135	end if
136
137	zh2(:,:)=zh(:,:)
138	zu2(:,:)=zu(:,:)
139	zv2(:,:)=zv(:,:)
140	zq2(:,:,:)=zq(:,:,:)
141
142	! -----------------------------
143	! Detection of unstable columns
144	! -----------------------------
145	! If ph(above) < ph(below) we set vtest=.true.
146
147	DO ig=1,ngrid
148	vtest(ig)=.false.
149	ENDDO
150
151	if (ico2.ne.0) then
152	! Special case if one of the tracers is CO2 gas
153	DO l=1,nlay
154	DO ig=1,ngrid
155	zhc(ig,l) = zh2(ig,l)(Azq2(ig,l,ico2)+B)
156	ENDDO
157	ENDDO
158	else
159	zhc(:,:)=zh2(:,:)
160	end if
161
162	! Find out which grid points are convectively unstable
163	DO l=2,nlay
164	DO ig=1,ngrid
165	IF(zhc(ig,l).LT.zhc(ig,l-1)) vtest(ig)=.true.
166	ENDDO
167	ENDDO
168
169	! Make a list of them
170	jcnt=0
171	DO ig=1,ngrid
172	IF(vtest(ig)) THEN
173	jcnt=jcnt+1
174	jadrs(jcnt)=ig
175	ENDIF
176	ENDDO
177
178
179	! ---------------------------------------------------------------
180	! Adjustment of the "jcnt" unstable profiles indicated by "jadrs"
181	! ---------------------------------------------------------------
182
183	DO jj = 1, jcnt ! loop on every convective grid point
184
185	i = jadrs(jj)
186
187	! Calculate sigma in this column
188	DO l=1,nlay+1
189	sig(l)=pplev(i,l)/pplev(i,1)
190
191	ENDDO
192	DO l=1,nlay
193	dsig(l)=sig(l)-sig(l+1)
194	sdsig(l)=ppopsk(i,l)*dsig(l)
195	ENDDO
196	l2 = 1
197
198	! Test loop upwards on l2
199
200	DO
201	l2 = l2 + 1
202	IF (l2 .GT. nlay) EXIT
203	IF (zhc(i, l2) .LT. zhc(i, l2-1)) THEN
204
205	! l2 is the highest level of the unstable column
206
207	l1 = l2 - 1
208	l = l1
209	zsm = sdsig(l2)
210	zdsm = dsig(l2)
211	zhm = zh2(i, l2)
212	if(ico2.ne.0) zqco2m = zq2(i,l2,ico2)
213
214	! Test loop downwards
215
216	DO
217	zsm = zsm + sdsig(l)
218	zdsm = zdsm + dsig(l)
219	zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
220	if(ico2.ne.0) then
221	zqco2m =
222	& zqco2m + dsig(l) * (zq2(i,l,ico2) - zqco2m) / zdsm
223	zhmc = zhm(Azqco2m+B)
224	else
225	zhmc = zhm
226	end if
227
228	! do we have to extend the column downwards?
229
230	down = .false.
231	IF (l1 .ne. 1) then !-- and then
232	IF (zhmc .lt. zhc(i, l1-1)) then
233	down = .true.
234	END IF
235	END IF
236
237	! this could be a problem...
238
239	if (down) then
240
241	l1 = l1 - 1
242	l = l1
243
244	else
245
246	! can we extend the column upwards?
247
248	if (l2 .eq. nlay) exit
249
250	if (zhc(i, l2+1) .ge. zhmc) exit
251
252	l2 = l2 + 1
253	l = l2
254
255	end if ! of if (down)
256
257	enddo ! of do ! Test loop forward
258
259	! New constant profile (average value)
260
261
262	zalpha=0.
263	zum=0.
264	zvm=0.
265	do iq=1,nq
266	zqm(iq) = 0.
267	end do
268	DO l = l1, l2
269	if(ico2.ne.0) then
270	zalpha=zalpha+
271	& ABS(zhc(i,l)/(A+Bzqco2m) -zhm)dsig(l)
272	else
273	zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
274	endif
275	zh2(i, l) = zhm
276	! modifs by RDW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
277	zum=zum+dsig(l)*zu2(i,l)
278	zvm=zvm+dsig(l)*zv2(i,l)
279	! zum=zum+dsig(l)*zu(i,l)
280	! zvm=zvm+dsig(l)*zv(i,l)
281	do iq=1,nq
282	zqm(iq) = zqm(iq)+dsig(l)*zq2(i,l,iq)
283	! zqm(iq) = zqm(iq)+dsig(l)*zq(i,l,iq)
284	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
285
286	! to conserve tracers/ KE, we must calculate zum, zvm and zqm using
287	! the up-to-date column values. If we do not do this, there are cases
288	! where convection stops at one level and starts at the next where we
289	! can break conservation of stuff (particularly tracers) significantly.
290
291	end do
292	ENDDO ! of DO l = l1, l2
293	zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
294	zum=zum/(sig(l1)-sig(l2+1))
295	zvm=zvm/(sig(l1)-sig(l2+1))
296	do iq=1,nq
297	zqm(iq) = zqm(iq)/(sig(l1)-sig(l2+1))
298	end do
299
300	IF(zalpha.GT.1.) THEN
301	zalpha=1.
302	ELSE
303	! IF(zalpha.LT.0.) STOP
304	IF(zalpha.LT.1.e-4) zalpha=1.e-4
305	ENDIF
306
307	DO l=l1,l2
308	zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
309	zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
310	do iq=1,nq
311	! zq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq))
312	zq2(i,l,iq)=zqm(iq)
313	end do
314	ENDDO
315	if (ico2.ne.0) then
316	DO l=l1, l2
317	zhc(i,l) = zh2(i,l)(Azq2(i,l,ico2)+B)
318	ENDDO
319	end if
320
321
322	l2 = l2 + 1
323
324	END IF ! End of l1 to l2 instability treatment
325	! We now continue to test from l2 upwards
326
327	ENDDO ! End of upwards loop on l2
328
329
330	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
331	! check conservation
332	cadjncons=0.0
333	if(water)then
334	iq = igcm_h2o_vap
335	do l = 1, nlay
336	masse = (pplev(i,l) - pplev(i,l+1))/g
337	cadjncons = cadjncons +
338	& masse*(zq2(i,l,iq)-zq(i,l,iq))/ptimestep
339	end do
340	endif
341
342	if(cadjncons.lt.-1.e-6)then
343	print*,'convadj has just crashed...'
344	print*,'i = ',i
345	print*,'l1 = ',l1
346	print*,'l2 = ',l2
347	print*,'cadjncons = ',cadjncons
348	do l = 1, nlay
349	print*,'dsig = ',dsig(l)
350	end do
351	do l = 1, nlay
352	print*,'dsig = ',dsig(l)
353	end do
354	do l = 1, nlay
355	print*,'sig = ',sig(l)
356	end do
357	do l = 1, nlay
358	print*,'pplay(ig,:) = ',pplay(i,l)
359	end do
360	do l = 1, nlay+1
361	print*,'pplev(ig,:) = ',pplev(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	ENDDO ! of DO jj = 1, jcnt ! loop on every convective grid point
386
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(,)'**convadj: 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
420

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

Download in other formats:

Original Format