Context Navigation

convadj.F @ 3523

Last change on this file since 3523 was 3322, checked in by emillour, 7 months ago
Generic PCM: Cleanup and cosmetics on convective adjustment routine. Make it a module in the process. EM
File size: 10.7 KB

Line
1	module convadj_mod
2
3	implicit none
4
5	contains
6
7	subroutine convadj(ngrid,nlay,nq,ptimestep,
8	& pplay,pplev,ppopsk,
9	& pu,pv,ph,pq,
10	& pdufi,pdvfi,pdhfi,pdqfi,
11	& pduadj,pdvadj,pdhadj,pdqadj)
12
13	use tracer_h, only: igcm_h2o_vap
14	use comcstfi_mod, only: g
15	use generic_cloud_common_h, only: epsi_generic
16	use generic_tracer_index_mod, only: generic_tracer_index
17	use callkeys_mod, only: tracer,water,generic_condensation,
18	& virtual_correction
19
20	implicit none
21
22	!==================================================================
23	!
24	! Purpose
25	! -------
26	! Compute dry convective adjustment.
27	! See old reference paper: Hourdin et al. JAS 1993
28	! "Meteorological Variability and the Annual Surface
29	! Pressure Cycle on Mars"
30	! https://doi.org/10.1175/1520-0469(1993)050%3C3625:MVATAS%3E2.0.CO;2
31	!==================================================================
32
33	! Arguments
34	! ---------
35
36	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
37	INTEGER,INTENT(IN) :: nlay ! number of atmospheric layers
38	INTEGER,INTENT(IN) :: nq ! number of tracers
39	REAL,INTENT(IN) :: ptimestep ! physics time step (s)
40	REAL,INTENT(IN) :: pplay(ngrid,nlay) ! mi-layer pressure (Pa)
41	REAL,INTENT(IN) :: pplev(ngrid,nlay+1) ! inter-layer pressure (Pa)
42	REAL,INTENT(IN) :: ppopsk(ngrid,nlay) ! Exner (wrt surface pressure)
43	! fields from dynamics
44	REAL,INTENT(IN) :: ph(ngrid,nlay) ! potential temperature (K)
45	REAL,INTENT(IN) :: pu(ngrid,nlay) ! zonal wind (m/s)
46	REAL,INTENT(IN) :: pv(ngrid,nlay) ! meridioanl wind (m/s)
47	REAL,INTENT(IN) :: pq(ngrid,nlay,nq) ! tracers (../kg_air)
48	! tendencies due to previous physical processes
49	REAL,INTENT(IN) :: pdufi(ngrid,nlay) ! on zonal wind (m/s/s)
50	REAL,INTENT(IN) :: pdvfi(ngrid,nlay) ! on meridional wind (m/s/s)
51	REAL,INTENT(IN) :: pdhfi(ngrid,nlay)! on potential temperature (/K/s)
52	REAL,INTENT(IN) :: pdqfi(ngrid,nlay,nq) ! on tracers (../kg_air/s)
53	! tendencies due to convetive adjustement
54	REAL,INTENT(OUT) :: pduadj(ngrid,nlay) ! on zonal wind (m/s/s)
55	REAL,INTENT(OUT) :: pdvadj(ngrid,nlay) ! on meridinal wind (m/s/s)
56	REAL,INTENT(OUT) :: pdhadj(ngrid,nlay) ! on potential temperature (/K/s)
57	REAL,INTENT(OUT) :: pdqadj(ngrid,nlay,nq) ! on traceurs (../kg_air/s)
58
59
60	! Local
61	! -----
62
63	INTEGER ig,i,l,l1,l2,jj
64	INTEGER jcnt, jadrs(ngrid)
65
66	REAL sig(nlay+1),sdsig(nlay),dsig(nlay)
67	REAL zu(ngrid,nlay),zv(ngrid,nlay)
68	REAL zh(ngrid,nlay), zvh(ngrid,nlay)
69	REAL zu2(ngrid,nlay),zv2(ngrid,nlay)
70	REAL zh2(ngrid,nlay),zvh2(ngrid,nlay),zhc(ngrid,nlay)
71	REAL zhm,zsm,zdsm,zum,zvm,zalpha,zhmc
72
73	! Tracers
74	INTEGER iq
75	REAL zq(ngrid,nlay,nq), zq2(ngrid,nlay,nq)
76	REAL zqm(nq)
77
78	integer igcm_generic_vap, igcm_generic_ice
79	logical call_ice_vap_generic
80
81	LOGICAL vtest(ngrid),down
82
83	! for conservation test
84	real masse,cadjncons
85
86	! --------------
87	! Initialisation
88	! --------------
89	zh(:,:)=ph(:,:)+pdhfi(:,:)*ptimestep
90	zu(:,:)=pu(:,:)+pdufi(:,:)*ptimestep
91	zv(:,:)=pv(:,:)+pdvfi(:,:)*ptimestep
92
93	if(tracer) then
94	zq(:,:,:)=pq(:,:,:)+pdqfi(:,:,:)*ptimestep
95	end if
96
97	zh2(:,:)=zh(:,:)
98	zu2(:,:)=zu(:,:)
99	zv2(:,:)=zv(:,:)
100	zq2(:,:,:)=zq(:,:,:)
101
102	! -----------------------------
103	! Detection of unstable columns
104	! -----------------------------
105	! If ph(above) < ph(below) we set vtest=.true.
106
107	DO ig=1,ngrid
108	vtest(ig)=.false.
109	ENDDO
110
111	if((generic_condensation) .and. (virtual_correction)) THEN
112	DO iq=1,nq
113	call generic_tracer_index(nq,iq,igcm_generic_vap,
114	& igcm_generic_ice,call_ice_vap_generic)
115	if(call_ice_vap_generic) then
116	zvh(:,:)=zh(:,:)*
117	& (1.+zq(:,:,igcm_generic_vap)/epsi_generic)/
118	& (1.+zq(:,:,igcm_generic_vap))
119	endif
120	ENDDO
121	zvh2(:,:)=zvh(:,:)
122	zhc(:,:)=zvh2(:,:)
123	else
124	zhc(:,:)=zh2(:,:)
125	endif
126
127	! Find out which grid points are convectively unstable
128	DO l=2,nlay
129	DO ig=1,ngrid
130	IF (zhc(ig,l).LT.zhc(ig,l-1)) THEN
131	vtest(ig)=.true.
132	ENDIF
133	ENDDO
134	ENDDO
135
136	! Make a list of them
137	jcnt=0
138	DO ig=1,ngrid
139	IF(vtest(ig)) THEN
140	jcnt=jcnt+1
141	jadrs(jcnt)=ig
142	ENDIF
143	ENDDO
144
145
146	! ---------------------------------------------------------------
147	! Adjustment of the "jcnt" unstable profiles indicated by "jadrs"
148	! ---------------------------------------------------------------
149
150	DO jj = 1, jcnt ! loop on every convective grid point
151
152	i = jadrs(jj)
153
154	! Calculate sigma in this column
155	DO l=1,nlay+1
156	sig(l)=pplev(i,l)/pplev(i,1)
157
158	ENDDO
159	DO l=1,nlay
160	dsig(l)=sig(l)-sig(l+1)
161	sdsig(l)=ppopsk(i,l)*dsig(l)
162	ENDDO
163	l2 = 1
164
165	! Test loop upwards on l2
166
167	DO
168	l2 = l2 + 1
169	IF (l2 .GT. nlay) EXIT
170	IF (zhc(i, l2).LT.zhc(i, l2-1)) THEN
171
172	! l2 is the highest level of the unstable column
173
174	l1 = l2 - 1
175	l = l1
176	zsm = sdsig(l2)
177	zdsm = dsig(l2)
178	zhm = zh2(i, l2)
179
180	! Test loop downwards
181
182	DO
183	zsm = zsm + sdsig(l)
184	zdsm = zdsm + dsig(l)
185	zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
186
187	if (generic_condensation .and. virtual_correction) then
188	zhmc = zhm*
189	& (1.+zq(i,l,igcm_generic_vap)/epsi_generic)/
190	& (1.+zq(i,l,igcm_generic_vap))
191	else
192	zhmc = zhm
193	endif
194
195	! do we have to extend the column downwards?
196
197	down = .false.
198	IF (l1 .ne. 1) then !-- and then
199	IF (zhmc.LT.zhc(i, l1-1)) then
200	down = .true.
201	END IF
202	END IF
203
204	! this could be a problem...
205
206	if (down) then
207
208	l1 = l1 - 1
209	l = l1
210
211	else
212
213	! can we extend the column upwards?
214
215	if (l2 .eq. nlay) exit
216
217	if (zhc(i, l2+1) .ge. zhmc) exit
218
219	l2 = l2 + 1
220	l = l2
221
222	end if
223
224	enddo
225
226	! New constant profile (average value)
227
228
229	zalpha=0.
230	zum=0.
231	zvm=0.
232	do iq=1,nq
233	zqm(iq) = 0.
234	end do
235	DO l = l1, l2
236	zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
237	zh2(i, l) = zhm
238	! modifs by RDW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
239	zum=zum+dsig(l)*zu2(i,l)
240	zvm=zvm+dsig(l)*zv2(i,l)
241	! zum=zum+dsig(l)*zu(i,l)
242	! zvm=zvm+dsig(l)*zv(i,l)
243	do iq=1,nq
244	zqm(iq) = zqm(iq)+dsig(l)*zq2(i,l,iq)
245	! zqm(iq) = zqm(iq)+dsig(l)*zq(i,l,iq)
246	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
247
248	! to conserve tracers/ KE, we must calculate zum, zvm and zqm using
249	! the up-to-date column values. If we do not do this, there are cases
250	! where convection stops at one level and starts at the next where we
251	! can break conservation of stuff (particularly tracers) significantly.
252
253	end do
254	ENDDO
255	zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
256	zum=zum/(sig(l1)-sig(l2+1))
257	zvm=zvm/(sig(l1)-sig(l2+1))
258	do iq=1,nq
259	zqm(iq) = zqm(iq)/(sig(l1)-sig(l2+1))
260	end do
261
262	IF(zalpha.GT.1.) THEN
263	zalpha=1.
264	ELSE
265	! IF(zalpha.LT.0.) STOP
266	IF(zalpha.LT.1.e-4) zalpha=1.e-4
267	ENDIF
268
269	DO l=l1,l2
270	zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
271	zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
272	do iq=1,nq
273	! zq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq))
274	zq2(i,l,iq)=zqm(iq)
275	end do
276	ENDDO
277
278
279	l2 = l2 + 1
280
281	END IF ! End of l1 to l2 instability treatment
282	! We now continue to test from l2 upwards
283
284	ENDDO ! End of upwards loop on l2
285
286
287	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
288	! check conservation
289	cadjncons=0.0
290	if(water)then
291	do l = 1, nlay
292	masse = (pplev(i,l) - pplev(i,l+1))/g
293	iq = igcm_h2o_vap
294	cadjncons = cadjncons +
295	& masse*(zq2(i,l,iq)-zq(i,l,iq))/ptimestep
296	end do
297	endif
298
299	if(cadjncons.lt.-1.e-6)then
300	print*,'convadj has just crashed...'
301	print*,'i = ',i
302	print*,'l1 = ',l1
303	print*,'l2 = ',l2
304	print*,'cadjncons = ',cadjncons
305	do l = 1, nlay
306	print*,'dsig = ',dsig(l)
307	end do
308	do l = 1, nlay
309	print*,'dsig = ',dsig(l)
310	end do
311	do l = 1, nlay
312	print*,'sig = ',sig(l)
313	end do
314	do l = 1, nlay
315	print*,'pplay(ig,:) = ',pplay(i,l)
316	end do
317	do l = 1, nlay+1
318	print*,'pplev(ig,:) = ',pplev(i,l)
319	end do
320	do l = 1, nlay
321	print*,'ph(ig,:) = ',ph(i,l)
322	end do
323	do l = 1, nlay
324	print*,'ph(ig,:) = ',ph(i,l)
325	end do
326	do l = 1, nlay
327	print*,'ph(ig,:) = ',ph(i,l)
328	end do
329	do l = 1, nlay
330	print*,'zh(ig,:) = ',zh(i,l)
331	end do
332	do l = 1, nlay
333	print*,'zh2(ig,:) = ',zh2(i,l)
334	end do
335	do l = 1, nlay
336	print*,'zq(ig,:,vap) = ',zq(i,l,igcm_h2o_vap)
337	end do
338	do l = 1, nlay
339	print*,'zq2(ig,:,vap) = ',zq2(i,l,igcm_h2o_vap)
340	end do
341	print*,'zqm(vap) = ',zqm(igcm_h2o_vap)
342	print*,'jadrs=',jadrs
343
344	call abort
345	endif
346	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
347
348
349	ENDDO
350
351	pdhadj(:,:)=(zh2(:,:)-zh(:,:))/ptimestep
352	pduadj(:,:)=(zu2(:,:)-zu(:,:))/ptimestep
353	pdvadj(:,:)=(zv2(:,:)-zv(:,:))/ptimestep
354
355	if(tracer) then
356	pdqadj(:,:,:)=(zq2(:,:,:)-zq(:,:,:))/ptimestep
357	end if
358
359	end subroutine convadj
360
361	end module convadj_mod

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

Download in other formats:

Original Format