Context Navigation

integrd_p.F @ 1422

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

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

Download in other formats:

Original Format