Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

convadj.F90 @ 3662

Last change on this file since 3662 was 3662, checked in by gmilcareck, 4 months ago
Generic PCM: Minor changing for the virtual potential temperature correction. And convadj.F becomes convadj.F90 . GM
File size: 10.9 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_theta_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),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_theta_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(:,:)(1.e0+(1.e0/epsi_generic-1.e0)zq(:,:,igcm_generic_vap))
117	endif
118	ENDDO
119	zhc(:,:)=zvh(:,:)
120	else
121	zhc(:,:)=zh2(:,:)
122	endif
123
124	! Find out which grid points are convectively unstable
125	DO l=2,nlay
126	DO ig=1,ngrid
127	IF (zhc(ig,l).LT.zhc(ig,l-1)) THEN
128	vtest(ig)=.true.
129	ENDIF
130	ENDDO
131	ENDDO
132
133	! Make a list of them
134	jcnt=0
135	DO ig=1,ngrid
136	IF(vtest(ig)) THEN
137	jcnt=jcnt+1
138	jadrs(jcnt)=ig
139	ENDIF
140	ENDDO
141
142
143	! ---------------------------------------------------------------
144	! Adjustment of the "jcnt" unstable profiles indicated by "jadrs"
145	! ---------------------------------------------------------------
146
147	DO jj = 1, jcnt ! loop on every convective grid point
148
149	i = jadrs(jj)
150
151	! Calculate sigma in this column
152	DO l=1,nlay+1
153	sig(l)=pplev(i,l)/pplev(i,1)
154
155	ENDDO
156	DO l=1,nlay
157	dsig(l)=sig(l)-sig(l+1)
158	sdsig(l)=ppopsk(i,l)*dsig(l)
159	ENDDO
160	l2 = 1
161
162	! Test loop upwards on l2
163
164	DO
165	l2 = l2 + 1
166	IF (l2 .GT. nlay) EXIT
167	IF (zhc(i, l2).LT.zhc(i, l2-1)) THEN
168
169	! l2 is the highest level of the unstable column
170
171	l1 = l2 - 1
172	l = l1
173	zsm = sdsig(l2)
174	zdsm = dsig(l2)
175	zhm = zh2(i, l2)
176	if (generic_condensation .and. virtual_theta_correction) then
177	zhmc = zh2(i, l2)(1.e0+(1.e0/epsi_generic-1.e0)zq(i,l2,igcm_generic_vap))
178	endif
179
180	! Test loop downwards
181
182	DO
183	zsm = zsm + sdsig(l)
184	zdsm = zdsm + dsig(l)
185
186	if (generic_condensation .and. virtual_theta_correction) then
187	zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
188	zhmc = zhmc + sdsig(l) * (zh2(i, l)(1.e0+(1.e0/epsi_generic-1.e0)zq(i,l,igcm_generic_vap)) - zhmc) / zsm
189	else
190	zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
191	zhmc = zhm
192	endif
193
194	! do we have to extend the column downwards?
195
196	down = .false.
197	IF (l1 .ne. 1) then !-- and then
198	IF (zhmc.LT.zhc(i, l1-1)) then
199	down = .true.
200	END IF
201	END IF
202
203	! this could be a problem...
204
205	if (down) then
206
207	l1 = l1 - 1
208	l = l1
209
210	else
211
212	! can we extend the column upwards?
213
214	if (l2 .eq. nlay) exit
215
216	if (zhc(i, l2+1) .ge. zhmc) exit
217
218	l2 = l2 + 1
219	l = l2
220
221	end if
222
223	enddo
224
225	! New constant profile (average value)
226
227
228	zalpha=0.
229	zum=0.
230	zvm=0.
231	do iq=1,nq
232	zqm(iq) = 0.
233	end do
234	DO l = l1, l2
235	zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
236	zh2(i, l) = zhm
237	! modifs by RDW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
238	zum=zum+dsig(l)*zu2(i,l)
239	zvm=zvm+dsig(l)*zv2(i,l)
240	! zum=zum+dsig(l)*zu(i,l)
241	! zvm=zvm+dsig(l)*zv(i,l)
242	do iq=1,nq
243	zqm(iq) = zqm(iq)+dsig(l)*zq2(i,l,iq)
244	! zqm(iq) = zqm(iq)+dsig(l)*zq(i,l,iq)
245	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
246
247	! to conserve tracers/ KE, we must calculate zum, zvm and zqm using
248	! the up-to-date column values. If we do not do this, there are cases
249	! where convection stops at one level and starts at the next where we
250	! can break conservation of stuff (particularly tracers) significantly.
251
252	end do
253	ENDDO
254	zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
255	zum=zum/(sig(l1)-sig(l2+1))
256	zvm=zvm/(sig(l1)-sig(l2+1))
257	do iq=1,nq
258	zqm(iq) = zqm(iq)/(sig(l1)-sig(l2+1))
259	end do
260
261	IF(zalpha.GT.1.) THEN
262	zalpha=1.
263	ELSE
264	! IF(zalpha.LT.0.) STOP
265	IF(zalpha.LT.1.e-4) zalpha=1.e-4
266	ENDIF
267
268	DO l=l1,l2
269	zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
270	zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
271	do iq=1,nq
272	! zq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq))
273	zq2(i,l,iq)=zqm(iq)
274	end do
275	ENDDO
276
277
278	l2 = l2 + 1
279
280	END IF ! End of l1 to l2 instability treatment
281	! We now continue to test from l2 upwards
282
283	ENDDO ! End of upwards loop on l2
284
285
286	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
287	! check conservation
288	cadjncons=0.0
289	if(water)then
290	do l = 1, nlay
291	masse = (pplev(i,l) - pplev(i,l+1))/g
292	iq = igcm_h2o_vap
293	cadjncons = cadjncons + masse*(zq2(i,l,iq)-zq(i,l,iq))/ptimestep
294	end do
295	endif
296
297	if(cadjncons.lt.-1.e-6)then
298	print*,'convadj has just crashed...'
299	print*,'i = ',i
300	print*,'l1 = ',l1
301	print*,'l2 = ',l2
302	print*,'cadjncons = ',cadjncons
303	do l = 1, nlay
304	print*,'dsig = ',dsig(l)
305	end do
306	do l = 1, nlay
307	print*,'dsig = ',dsig(l)
308	end do
309	do l = 1, nlay
310	print*,'sig = ',sig(l)
311	end do
312	do l = 1, nlay
313	print*,'pplay(ig,:) = ',pplay(i,l)
314	end do
315	do l = 1, nlay+1
316	print*,'pplev(ig,:) = ',pplev(i,l)
317	end do
318	do l = 1, nlay
319	print*,'ph(ig,:) = ',ph(i,l)
320	end do
321	do l = 1, nlay
322	print*,'ph(ig,:) = ',ph(i,l)
323	end do
324	do l = 1, nlay
325	print*,'ph(ig,:) = ',ph(i,l)
326	end do
327	do l = 1, nlay
328	print*,'zh(ig,:) = ',zh(i,l)
329	end do
330	do l = 1, nlay
331	print*,'zh2(ig,:) = ',zh2(i,l)
332	end do
333	do l = 1, nlay
334	print*,'zq(ig,:,vap) = ',zq(i,l,igcm_h2o_vap)
335	end do
336	do l = 1, nlay
337	print*,'zq2(ig,:,vap) = ',zq2(i,l,igcm_h2o_vap)
338	end do
339	print*,'zqm(vap) = ',zqm(igcm_h2o_vap)
340	print*,'jadrs=',jadrs
341
342	call abort
343	endif
344	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
345
346
347	ENDDO
348
349	pdhadj(:,:)=(zh2(:,:)-zh(:,:))/ptimestep
350	pduadj(:,:)=(zu2(:,:)-zu(:,:))/ptimestep
351	pdvadj(:,:)=(zv2(:,:)-zv(:,:))/ptimestep
352
353	if(tracer) then
354	pdqadj(:,:,:)=(zq2(:,:,:)-zq(:,:,:))/ptimestep
355	end if
356
357	end subroutine convadj
358
359	end module convadj_mod

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

Download in other formats:

Original Format