Context Navigation

convadj.F @ 3737

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

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

Download in other formats:

Original Format