Context Navigation

convadj.F @ 3529

Last change on this file since 3529 was 1668, checked in by jvatant, 8 years ago
Another round of cleaning of dust and dummy tracers JVO
File size: 9.1 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	USE tracer_h
9	use comcstfi_mod, only: g
10	use callkeys_mod, only: tracer
11
12	implicit none
13
14	!==================================================================
15	!
16	! Purpose
17	! -------
18	! Calculates dry convective adjustment.
19	!
20	! Authors
21	! -------
22	! Original author unknown.
23	! Modif. 2005 by F. Forget.
24	!
25	!==================================================================
26
27	! ------------
28	! Declarations
29	! ------------
30
31
32	! Arguments
33	! ---------
34
35	INTEGER ngrid,nlay
36	REAL ptimestep
37	REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay)
38	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay)
39	REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay)
40	REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay)
41
42	! Tracers
43	integer nq
44	real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
45	real pdqadj(ngrid,nlay,nq)
46
47
48	! Local
49	! -----
50
51	INTEGER ig,i,l,l1,l2,jj
52	INTEGER jcnt, jadrs(ngrid)
53
54	REAL sig(nlay+1),sdsig(nlay),dsig(nlay)
55	REAL zu(ngrid,nlay),zv(ngrid,nlay)
56	REAL zh(ngrid,nlay)
57	REAL zu2(ngrid,nlay),zv2(ngrid,nlay)
58	REAL zh2(ngrid,nlay), zhc(ngrid,nlay)
59	REAL zhm,zsm,zdsm,zum,zvm,zalpha,zhmc
60
61	! Tracers
62	INTEGER iq
63	REAL zq(ngrid,nlay,nq), zq2(ngrid,nlay,nq)
64	REAL zqm(nq)
65	real A , B
66	save A, B
67	!$OMP THREADPRIVATE(A,B)
68
69	real mtot1, mtot2 , mm1, mm2
70	integer l1ref, l2ref
71	LOGICAL vtest(ngrid),down, firstcall
72	save firstcall
73	data firstcall/.true./
74	!$OMP THREADPRIVATE(firstcall)
75
76	! for conservation test
77	real masse,cadjncons
78
79	EXTERNAL SCOPY
80
81	! --------------
82	! Initialisation
83	! --------------
84
85	IF (firstcall) THEN
86	firstcall=.false.
87	ENDIF ! of IF (firstcall)
88
89	DO l=1,nlay
90	DO ig=1,ngrid
91	zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep
92	zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep
93	zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep
94	ENDDO
95	ENDDO
96
97	if(tracer) then
98	DO iq =1, nq
99	DO l=1,nlay
100	DO ig=1,ngrid
101	zq(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
102	ENDDO
103	ENDDO
104	ENDDO
105	end if
106
107	CALL scopy(ngrid*nlay,zh,1,zh2,1)
108	CALL scopy(ngrid*nlay,zu,1,zu2,1)
109	CALL scopy(ngrid*nlay,zv,1,zv2,1)
110	CALL scopy(ngridnlaynq,zq,1,zq2,1)
111
112	! -----------------------------
113	! Detection of unstable columns
114	! -----------------------------
115	! If ph(above) < ph(below) we set vtest=.true.
116
117	DO ig=1,ngrid
118	vtest(ig)=.false.
119	ENDDO
120
121	CALL scopy(ngrid*nlay,zh2,1,zhc,1)
122
123	! Find out which grid points are convectively unstable
124	DO l=2,nlay
125	DO ig=1,ngrid
126	IF(zhc(ig,l).LT.zhc(ig,l-1)) vtest(ig)=.true.
127	ENDDO
128	ENDDO
129
130	! Make a list of them
131	jcnt=0
132	DO ig=1,ngrid
133	IF(vtest(ig)) THEN
134	jcnt=jcnt+1
135	jadrs(jcnt)=ig
136	ENDIF
137	ENDDO
138
139
140	! ---------------------------------------------------------------
141	! Adjustment of the "jcnt" unstable profiles indicated by "jadrs"
142	! ---------------------------------------------------------------
143
144	DO jj = 1, jcnt ! loop on every convective grid point
145
146	i = jadrs(jj)
147
148	! Calculate sigma in this column
149	DO l=1,nlay+1
150	sig(l)=pplev(i,l)/pplev(i,1)
151
152	ENDDO
153	DO l=1,nlay
154	dsig(l)=sig(l)-sig(l+1)
155	sdsig(l)=ppopsk(i,l)*dsig(l)
156	ENDDO
157	l2 = 1
158
159	! Test loop upwards on l2
160
161	DO
162	l2 = l2 + 1
163	IF (l2 .GT. nlay) EXIT
164	IF (zhc(i, l2) .LT. zhc(i, l2-1)) THEN
165
166	! l2 is the highest level of the unstable column
167
168	l1 = l2 - 1
169	l = l1
170	zsm = sdsig(l2)
171	zdsm = dsig(l2)
172	zhm = zh2(i, l2)
173
174	! Test loop downwards
175
176	DO
177	zsm = zsm + sdsig(l)
178	zdsm = zdsm + dsig(l)
179	zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
180	zhmc = zhm
181
182	! do we have to extend the column downwards?
183
184	down = .false.
185	IF (l1 .ne. 1) then !-- and then
186	IF (zhmc .lt. zhc(i, l1-1)) then
187	down = .true.
188	END IF
189	END IF
190
191	! this could be a problem...
192
193	if (down) then
194
195	l1 = l1 - 1
196	l = l1
197
198	else
199
200	! can we extend the column upwards?
201
202	if (l2 .eq. nlay) exit
203
204	if (zhc(i, l2+1) .ge. zhmc) exit
205
206	l2 = l2 + 1
207	l = l2
208
209	end if
210
211	enddo
212
213	! New constant profile (average value)
214
215
216	zalpha=0.
217	zum=0.
218	zvm=0.
219	do iq=1,nq
220	zqm(iq) = 0.
221	end do
222	DO l = l1, l2
223	zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
224	zh2(i, l) = zhm
225	! modifs by RDW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
226	zum=zum+dsig(l)*zu2(i,l)
227	zvm=zvm+dsig(l)*zv2(i,l)
228	! zum=zum+dsig(l)*zu(i,l)
229	! zvm=zvm+dsig(l)*zv(i,l)
230	do iq=1,nq
231	zqm(iq) = zqm(iq)+dsig(l)*zq2(i,l,iq)
232	! zqm(iq) = zqm(iq)+dsig(l)*zq(i,l,iq)
233	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
234
235	! to conserve tracers/ KE, we must calculate zum, zvm and zqm using
236	! the up-to-date column values. If we do not do this, there are cases
237	! where convection stops at one level and starts at the next where we
238	! can break conservation of stuff (particularly tracers) significantly.
239
240	end do
241	ENDDO
242	zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
243	zum=zum/(sig(l1)-sig(l2+1))
244	zvm=zvm/(sig(l1)-sig(l2+1))
245	do iq=1,nq
246	zqm(iq) = zqm(iq)/(sig(l1)-sig(l2+1))
247	end do
248
249	IF(zalpha.GT.1.) THEN
250	zalpha=1.
251	ELSE
252	! IF(zalpha.LT.0.) STOP
253	IF(zalpha.LT.1.e-4) zalpha=1.e-4
254	ENDIF
255
256	DO l=l1,l2
257	zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
258	zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
259	do iq=1,nq
260	! zq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq))
261	zq2(i,l,iq)=zqm(iq)
262	end do
263	ENDDO
264
265
266	l2 = l2 + 1
267
268	END IF ! End of l1 to l2 instability treatment
269	! We now continue to test from l2 upwards
270
271	ENDDO ! End of upwards loop on l2
272
273
274	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
275	! check conservation
276	cadjncons=0.0
277
278	if(cadjncons.lt.-1.e-6)then
279	print*,'convadj has just crashed...'
280	print*,'i = ',i
281	print*,'l1 = ',l1
282	print*,'l2 = ',l2
283	print*,'cadjncons = ',cadjncons
284	do l = 1, nlay
285	print*,'dsig = ',dsig(l)
286	end do
287	do l = 1, nlay
288	print*,'dsig = ',dsig(l)
289	end do
290	do l = 1, nlay
291	print*,'sig = ',sig(l)
292	end do
293	do l = 1, nlay
294	print*,'pplay(ig,:) = ',pplay(i,l)
295	end do
296	do l = 1, nlay+1
297	print*,'pplev(ig,:) = ',pplev(i,l)
298	end do
299	do l = 1, nlay
300	print*,'ph(ig,:) = ',ph(i,l)
301	end do
302	do l = 1, nlay
303	print*,'ph(ig,:) = ',ph(i,l)
304	end do
305	do l = 1, nlay
306	print*,'ph(ig,:) = ',ph(i,l)
307	end do
308	do l = 1, nlay
309	print*,'zh(ig,:) = ',zh(i,l)
310	end do
311	do l = 1, nlay
312	print*,'zh2(ig,:) = ',zh2(i,l)
313	end do
314	print*,'jadrs=',jadrs
315
316	call abort
317	endif
318	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
319
320
321	ENDDO
322
323	DO l=1,nlay
324	DO ig=1,ngrid
325	pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep
326	pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep
327	pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep
328	ENDDO
329	ENDDO
330
331	if(tracer) then
332	do iq=1, nq
333	do l=1,nlay
334	DO ig=1,ngrid
335	pdqadj(ig,l,iq)=(zq2(ig,l,iq)-zq(ig,l,iq))/ptimestep
336	end do
337	end do
338	end do
339	end if
340
341
342	! output
343	! if (ngrid.eq.1) then
344	! ig=1
345	! iq =1
346	! write(,)'**** l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq)'
347	! do l=nlay,1,-1
348	! write(,) l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq)
349	! end do
350	! end if
351
352
353	return
354	end

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

Download in other formats:

Original Format