Context Navigation

integrd_p.F @ 4962

Last change on this file since 4962 was 1616, checked in by Ehouarn Millour, 13 years ago
Some cleanup around what is done during the integration step of dynamical tendencies and namely removed computation of (unused) finvmaold, thereby saving us the expense of a call to the (costly) filter at every dynamical time step. Checked (on Vargas, in seq, omp, mpi and mixed mode) that this doesn't change the GCM results, as expected. EM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 9.6 KB

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

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

Download in other formats:

Original Format