Context Navigation

integrd_p.F @ 2499

Last change on this file since 2499 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
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 1907 2013-11-26 13:10:46Z oboucher $
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	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