Context Navigation

integrd_p.F @ 3532

Last change on this file since 3532 was 2126, checked in by emillour, 6 years ago

Common dynamics:
Some work to enforce total tracer mass conservation in the dynamics.
Still to be further studied and validated.
For now these changes are triggered by setting a "force_conserv_tracer"
flag to ".true." in run.def (default is ".false." to not change anything
with respect to previous versions).
When force_conserv_tracer=.true. then:

Rescale tracer mass in caladvtrac after tracer advection computations
Recompute q ratios once atmospheric mass has been updated in integrd

These steps technically ensure total tracer mass conservation but it
might be the tracer advection scheme and/or time-stepping updating
sequence of fields that should be rethought or fixed.
EM

File size: 10.4 KB

Line
1	!
2	! $Id: integrd_p.F 1616 2012-02-17 11:59:00Z emillour $
3	!
4	SUBROUTINE integrd_p
5	$ ( nq,vcovm1,ucovm1,tetam1,psm1,massem1,
6	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps0,masse,phis) !,finvmaold)
7	USE parallel_lmdz, ONLY: ij_begin, ij_end, pole_nord, pole_sud,
8	& omp_chunk
9	USE control_mod, only : planet_type,force_conserv_tracer
10	USE comvert_mod, ONLY: ap,bp
11	USE comconst_mod, ONLY: pi
12	USE logic_mod, ONLY: leapf
13	USE temps_mod, ONLY: dt
14	IMPLICIT NONE
15
16
17	c=======================================================================
18	c
19	c Auteur: P. Le Van
20	c -------
21	c
22	c objet:
23	c ------
24	c
25	c Incrementation des tendances dynamiques
26	c
27	c=======================================================================
28	c-----------------------------------------------------------------------
29	c Declarations:
30	c -------------
31
32	#include "dimensions.h"
33	#include "paramet.h"
34	#include "comgeom.h"
35	#include "iniprint.h"
36
37	c Arguments:
38	c ----------
39
40	integer,intent(in) :: nq ! number of tracers to handle in this routine
41	real,intent(inout) :: vcov(ip1jm,llm) ! covariant meridional wind
42	real,intent(inout) :: ucov(ip1jmp1,llm) ! covariant zonal wind
43	real,intent(inout) :: teta(ip1jmp1,llm) ! potential temperature
44	real,intent(inout) :: q(ip1jmp1,llm,nq) ! advected tracers
45	real,intent(inout) :: ps0(ip1jmp1) ! surface pressure
46	real,intent(inout) :: masse(ip1jmp1,llm) ! atmospheric mass
47	real,intent(in) :: phis(ip1jmp1) ! ground geopotential !!! unused
48	! values at previous time step
49	real,intent(inout) :: vcovm1(ip1jm,llm)
50	real,intent(inout) :: ucovm1(ip1jmp1,llm)
51	real,intent(inout) :: tetam1(ip1jmp1,llm)
52	real,intent(inout) :: psm1(ip1jmp1)
53	real,intent(inout) :: massem1(ip1jmp1,llm)
54	! the tendencies to add
55	real,intent(in) :: dv(ip1jm,llm)
56	real,intent(in) :: du(ip1jmp1,llm)
57	real,intent(in) :: dteta(ip1jmp1,llm)
58	real,intent(in) :: dp(ip1jmp1)
59	real,intent(in) :: dq(ip1jmp1,llm,nq) !!! unused
60	! real,intent(out) :: finvmaold(ip1jmp1,llm) !!! unused
61
62	c Local:
63	c ------
64
65	REAL vscr( ip1jm ),uscr( ip1jmp1 ),hscr( ip1jmp1 ),pscr(ip1jmp1)
66	REAL massescr( ip1jmp1,llm )
67	REAL :: massratio(ip1jmp1,llm)
68	! REAL finvmasse(ip1jmp1,llm)
69	REAL,SAVE :: p(ip1jmp1,llmp1)
70	REAL tpn,tps,tppn(iim),tpps(iim)
71	REAL qpn,qps,qppn(iim),qpps(iim)
72	REAL,SAVE :: deltap( ip1jmp1,llm )
73
74	INTEGER l,ij,iq,i,j
75
76	REAL SSUM
77	EXTERNAL SSUM
78	INTEGER ijb,ije,jjb,jje
79	REAL,SAVE :: ps(ip1jmp1)=0
80	LOGICAL :: checksum
81	INTEGER :: stop_it
82	c-----------------------------------------------------------------------
83	c$OMP BARRIER
84	if (pole_nord) THEN
85	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
86	DO l = 1,llm
87	DO ij = 1,iip1
88	ucov( ij , l) = 0.
89	uscr( ij ) = 0.
90	ENDDO
91	ENDDO
92	c$OMP END DO NOWAIT
93	ENDIF
94
95	if (pole_sud) THEN
96	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
97	DO l = 1,llm
98	DO ij = 1,iip1
99	ucov( ij +ip1jm, l) = 0.
100	uscr( ij +ip1jm ) = 0.
101	ENDDO
102	ENDDO
103	c$OMP END DO NOWAIT
104	ENDIF
105
106	c ............ integration de ps ..............
107
108	c CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)
109
110	ijb=ij_begin
111	ije=ij_end
112	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
113	DO l = 1,llm
114	massescr(ijb:ije,l)=masse(ijb:ije,l)
115	ENDDO
116	c$OMP END DO NOWAIT
117
118	c$OMP DO SCHEDULE(STATIC)
119	DO 2 ij = ijb,ije
120	pscr (ij) = ps0(ij)
121	ps (ij) = psm1(ij) + dt * dp(ij)
122	2 CONTINUE
123	c$OMP END DO
124	c$OMP BARRIER
125	c --> ici synchro OPENMP pour ps
126
127	checksum=.TRUE.
128	stop_it=0
129
130	c$OMP DO SCHEDULE(STATIC)
131	DO ij = ijb,ije
132	IF( ps(ij).LT.0. ) THEN
133	IF (checksum) stop_it=ij
134	checksum=.FALSE.
135	ENDIF
136	ENDDO
137	c$OMP END DO NOWAIT
138
139	IF( .NOT. checksum ) THEN
140	write(lunout,*) "integrd: negative surface pressure ",
141	& ps(stop_it)
142	write(lunout,*) " at node ij =", stop_it
143	! since ij=j+(i-1)*jjp1 , we have
144	j=modulo(stop_it,jjp1)
145	i=1+(stop_it-j)/jjp1
146	write(lunout,) " lon = ",rlonv(i)180./pi, " deg",
147	& " lat = ",rlatu(j)*180./pi, " deg"
148	write(lunout,*) " psm1(ij)=",psm1(stop_it)," dt=",dt,
149	& " dp(ij)=",dp(stop_it)
150	call abort_gcm("integrd_p", "negative surface pressure", 1)
151	ENDIF
152
153	c
154	C$OMP MASTER
155	if (pole_nord) THEN
156
157	DO ij = 1, iim
158	tppn(ij) = aire( ij ) * ps( ij )
159	ENDDO
160	tpn = SSUM(iim,tppn,1)/apoln
161	DO ij = 1, iip1
162	ps( ij ) = tpn
163	ENDDO
164
165	ENDIF
166
167	if (pole_sud) THEN
168
169	DO ij = 1, iim
170	tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
171	ENDDO
172	tps = SSUM(iim,tpps,1)/apols
173	DO ij = 1, iip1
174	ps(ij+ip1jm) = tps
175	ENDDO
176
177	ENDIF
178	c$OMP END MASTER
179	c$OMP BARRIER
180	c
181	c ... Calcul de la nouvelle masse d'air au dernier temps integre t+1 ...
182	c
183
184	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
185	c$OMP BARRIER
186	CALL massdair_p ( p , masse )
187
188	! Ehouarn : we don't use/need finvmaold and finvmasse,
189	! so might as well not compute them
190	!c CALL SCOPY( ijp1llm , masse, 1, finvmasse, 1 )
191	! ijb=ij_begin
192	! ije=ij_end
193	!
194	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
195	! DO l = 1,llm
196	! finvmasse(ijb:ije,l)=masse(ijb:ije,l)
197	! ENDDO
198	!c$OMP END DO NOWAIT
199	!
200	! jjb=jj_begin
201	! jje=jj_end
202	! CALL filtreg_p( finvmasse,jjb,jje, jjp1, llm, -2, 2, .TRUE., 1 )
203	c
204
205	c ............ integration de ucov, vcov, h ..............
206
207	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
208	DO 10 l = 1,llm
209
210	ijb=ij_begin
211	ije=ij_end
212	if (pole_nord) ijb=ij_begin+iip1
213	if (pole_sud) ije=ij_end-iip1
214
215	DO 4 ij = ijb,ije
216	uscr( ij ) = ucov( ij,l )
217	ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
218	4 CONTINUE
219
220	ijb=ij_begin
221	ije=ij_end
222	if (pole_sud) ije=ij_end-iip1
223
224	DO 5 ij = ijb,ije
225	vscr( ij ) = vcov( ij,l )
226	vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
227	5 CONTINUE
228
229	ijb=ij_begin
230	ije=ij_end
231
232	DO 6 ij = ijb,ije
233	hscr( ij ) = teta(ij,l)
234	teta ( ij,l ) = tetam1(ij,l) * massem1(ij,l) / masse(ij,l)
235	$ + dt * dteta(ij,l) / masse(ij,l)
236	6 CONTINUE
237
238	c .... Calcul de la valeur moyenne, unique aux poles pour teta ......
239	c
240	c
241	IF (pole_nord) THEN
242
243	DO ij = 1, iim
244	tppn(ij) = aire( ij ) * teta( ij ,l)
245	ENDDO
246	tpn = SSUM(iim,tppn,1)/apoln
247
248	DO ij = 1, iip1
249	teta( ij ,l) = tpn
250	ENDDO
251
252	ENDIF
253
254	IF (pole_sud) THEN
255
256	DO ij = 1, iim
257	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
258	ENDDO
259	tps = SSUM(iim,tpps,1)/apols
260
261	DO ij = 1, iip1
262	teta(ij+ip1jm,l) = tps
263	ENDDO
264
265	ENDIF
266	c
267
268	IF(leapf) THEN
269	c CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
270	c CALL SCOPY ( ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
271	c CALL SCOPY ( ip1jmp1, hscr(1), 1, tetam1(1, l), 1 )
272	ijb=ij_begin
273	ije=ij_end
274	ucovm1(ijb:ije,l)=uscr(ijb:ije)
275	tetam1(ijb:ije,l)=hscr(ijb:ije)
276	if (pole_sud) ije=ij_end-iip1
277	vcovm1(ijb:ije,l)=vscr(ijb:ije)
278
279	END IF
280
281	10 CONTINUE
282	c$OMP END DO NOWAIT
283
284	c
285	c ....... integration de q ......
286	c
287	ijb=ij_begin
288	ije=ij_end
289
290	if (planet_type.eq."earth") then
291	! Earth-specific treatment of first 2 tracers (water)
292	c$OMP BARRIER
293	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
294	DO l = 1, llm
295	DO ij = ijb, ije
296	deltap(ij,l) = p(ij,l) - p(ij,l+1)
297	ENDDO
298	ENDDO
299	c$OMP END DO NOWAIT
300	c$OMP BARRIER
301
302	CALL qminimum_p( q, nq, deltap )
303	c
304	c ..... Calcul de la valeur moyenne, unique aux poles pour q .....
305	c
306	c$OMP BARRIER
307	IF (pole_nord) THEN
308
309	DO iq = 1, nq
310
311	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
312	DO l = 1, llm
313
314	DO ij = 1, iim
315	qppn(ij) = aire( ij ) * q( ij ,l,iq)
316	ENDDO
317	qpn = SSUM(iim,qppn,1)/apoln
318
319	DO ij = 1, iip1
320	q( ij ,l,iq) = qpn
321	ENDDO
322
323	ENDDO
324	c$OMP END DO NOWAIT
325
326	ENDDO
327
328	ENDIF
329
330	IF (pole_sud) THEN
331
332	DO iq = 1, nq
333
334	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
335	DO l = 1, llm
336
337	DO ij = 1, iim
338	qpps(ij) = aire(ij+ip1jm) * q(ij+ip1jm,l,iq)
339	ENDDO
340	qps = SSUM(iim,qpps,1)/apols
341
342	DO ij = 1, iip1
343	q(ij+ip1jm,l,iq) = qps
344	ENDDO
345
346	ENDDO
347	c$OMP END DO NOWAIT
348
349	ENDDO
350
351	ENDIF
352
353	! Ehouarn: forget about finvmaold
354	!c CALL SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
355	!
356	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
357	! DO l = 1, llm
358	! finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)
359	! ENDDO
360	!c$OMP END DO NOWAIT
361
362	endif ! of if (planet_type.eq."earth")
363
364
365	if (force_conserv_tracer) then
366	! Ehouarn: try to keep total amont of tracers fixed
367	! by acounting for mass change in each cell
368	massratio(ijb:ije,1:llm)=massescr(ijb:ije,1:llm)
369	& /masse(ijb:ije,1:llm)
370	do iq=1,nq
371	q(ijb:ije,1:llm,iq)=q(ijb:ije,1:llm,iq)
372	& *massratio(ijb:ije,1:llm)
373	enddo
374	endif ! of if (force_conserv_tracer)
375	c
376	c
377	c ..... FIN de l'integration de q .......
378
379	15 continue
380
381	c$OMP DO SCHEDULE(STATIC)
382	DO ij=ijb,ije
383	ps0(ij)=ps(ij)
384	ENDDO
385	c$OMP END DO NOWAIT
386
387	c .................................................................
388
389
390	IF( leapf ) THEN
391	c CALL SCOPY ( ip1jmp1 , pscr , 1, psm1 , 1 )
392	c CALL SCOPY ( ip1jmp1*llm, massescr, 1, massem1, 1 )
393	c$OMP DO SCHEDULE(STATIC)
394	DO ij=ijb,ije
395	psm1(ij)=pscr(ij)
396	ENDDO
397	c$OMP END DO NOWAIT
398
399	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
400	DO l = 1, llm
401	massem1(ijb:ije,l)=massescr(ijb:ije,l)
402	ENDDO
403	c$OMP END DO NOWAIT
404	END IF
405	c$OMP BARRIER
406	RETURN
407	END

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

Download in other formats:

Original Format