Context Navigation

convadj.F @ 1036

Last change on this file since 1036 was 1036, checked in by emillour, 11 years ago

Mars GCM: (a first step towards using parallel dynamics)

IMPORTANT CHANGE: Implemented dynamic tracers. It is no longer necessary to compile the model with the '-t #' option, number of tracers is simply read from tracer.def file (as before). Adapted makegcm_* scripts (and co.) accordingly. Technical aspects of the switch to dynamic tracers are:
- advtrac.h (in dyn3d) removed and replaced by module infotrac.F
- tracer.h (in phymars) removed and replaced by module tracer_mod.F90 (which contains nqmx, the number of tracers, etc. and can be used anywhere in the physics).
Included some side cleanups: removed unused files (in dyn3d) anldoppler2.F, anl_mcdstats.F and anl_stats-diag.F, and all the unecessary dimensions.* files in grid/dimension.
Checked that changes are clean and that GCM yields identical results (in debug mode) to previous svn version.

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

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

Download in other formats:

Original Format