Context Navigation

integrd_loc.F @ 3717

Last change on this file since 3717 was 2286, checked in by Ehouarn Millour, 9 years ago

fix in dyn3d the array out of bound issue (in qminimum) that was only corrected in dyn3dmem in rev 2285.
comment out many invasive debug writes.

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

File size: 11.8 KB

Line
1	!
2	! $Id: integrd_p.F 1299 2010-01-20 14:27:21Z fairhead $
3	!
4	SUBROUTINE integrd_loc
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
9	USE mod_filtreg_p
10	USE write_field_loc
11	USE write_field
12	USE integrd_mod
13	USE infotrac, ONLY: ok_iso_verif ! ajout CRisi
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 "comconst.h"
35	#include "comgeom.h"
36	#include "comvert.h"
37	#include "logic.h"
38	#include "temps.h"
39	#include "serre.h"
40	#include "iniprint.h"
41	! include 'mpif.h'
42
43	c Arguments:
44	c ----------
45
46	INTEGER,intent(in) :: nq ! number of tracers to handle in this routine
47
48	REAL,INTENT(INOUT) :: vcov(ijb_v:ije_v,llm) ! covariant meridional wind
49	REAL,INTENT(INOUT) :: ucov(ijb_u:ije_u,llm) ! covariant zonal wind
50	REAL,INTENT(INOUT) :: teta(ijb_u:ije_u,llm) ! potential temperature
51	REAL,INTENT(INOUT) :: q(ijb_u:ije_u,llm,nq) ! advected tracers
52	REAL,INTENT(INOUT) :: ps0(ijb_u:ije_u) ! surface pressure
53	REAL,INTENT(INOUT) :: masse(ijb_u:ije_u,llm) ! atmospheric mass
54	REAL,INTENT(INOUT) :: phis(ijb_u:ije_u) ! ground geopotential !!! unused
55	! values at previous time step
56	REAL,INTENT(INOUT) :: vcovm1(ijb_v:ije_v,llm)
57	REAL,INTENT(INOUT) :: ucovm1(ijb_u:ije_u,llm)
58	REAL,INTENT(INOUT) :: tetam1(ijb_u:ije_u,llm)
59	REAL,INTENT(INOUT) :: psm1(ijb_u:ije_u)
60	REAL,INTENT(INOUT) :: massem1(ijb_u:ije_u,llm)
61	! the tendencies to add
62	REAL,INTENT(INOUT) :: dv(ijb_v:ije_v,llm)
63	REAL,INTENT(INOUT) :: du(ijb_u:ije_u,llm)
64	REAL,INTENT(INOUT) :: dteta(ijb_u:ije_u,llm)
65	REAL,INTENT(INOUT) :: dp(ijb_u:ije_u)
66	REAL,INTENT(INOUT) :: dq(ijb_u:ije_u,llm,nq) !!! unused
67	! REAL,INTENT(INOUT) ::finvmaold(ijb_u:ije_u,llm) !!! unused
68
69	c Local:
70	c ------
71
72	REAL vscr( ijb_v:ije_v ),uscr( ijb_u:ije_u )
73	REAL hscr( ijb_u:ije_u ),pscr(ijb_u:ije_u)
74	REAL massescr( ijb_u:ije_u,llm )
75	! REAL finvmasse(ijb_u:ije_u,llm)
76	REAL tpn,tps,tppn(iim),tpps(iim)
77	REAL qpn,qps,qppn(iim),qpps(iim)
78
79	INTEGER l,ij,iq,i,j
80
81	REAL SSUM
82	EXTERNAL SSUM
83	INTEGER ijb,ije,jjb,jje
84	LOGICAL :: checksum
85	LOGICAL,SAVE :: checksum_all=.TRUE.
86	INTEGER :: stop_it
87	INTEGER :: ierr
88
89	!write(,) 'integrd 88: entree, nq=',nq
90	c-----------------------------------------------------------------------
91
92	c$OMP BARRIER
93	if (pole_nord) THEN
94	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
95	DO l = 1,llm
96	DO ij = 1,iip1
97	ucov( ij , l) = 0.
98	uscr( ij ) = 0.
99	ENDDO
100	ENDDO
101	c$OMP END DO NOWAIT
102	ENDIF
103
104	if (pole_sud) THEN
105	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
106	DO l = 1,llm
107	DO ij = 1,iip1
108	ucov( ij +ip1jm, l) = 0.
109	uscr( ij +ip1jm ) = 0.
110	ENDDO
111	ENDDO
112	c$OMP END DO NOWAIT
113	ENDIF
114
115	c ............ integration de ps ..............
116
117	c CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)
118
119	ijb=ij_begin
120	ije=ij_end
121	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
122	DO l = 1,llm
123	massescr(ijb:ije,l)=masse(ijb:ije,l)
124	ENDDO
125	c$OMP END DO NOWAIT
126
127	c$OMP DO SCHEDULE(STATIC)
128	DO 2 ij = ijb,ije
129	pscr (ij) = ps0(ij)
130	ps (ij) = psm1(ij) + dt * dp(ij)
131
132	2 CONTINUE
133
134	c$OMP END DO
135	c$OMP BARRIER
136	c --> ici synchro OPENMP pour ps
137
138	checksum=.TRUE.
139	stop_it=0
140
141	c$OMP MASTER
142	!c$OMP DO SCHEDULE(STATIC)
143	DO ij = ijb,ije
144	IF( ps(ij).LT.0. ) THEN
145	IF (checksum) stop_it=ij
146	checksum=.FALSE.
147	ENDIF
148	ENDDO
149	!c$OMP END DO NOWAIT
150
151	! CALL MPI_ALLREDUCE(checksum,checksum_all,1,
152	! & MPI_LOGICAL,MPI_LOR,COMM_LMDZ,ierr)
153	IF( .NOT. checksum ) THEN
154	write(lunout,*) "integrd: ps = ", ps(stop_it)
155	write(lunout,*) " at node ij =", stop_it
156	! since ij=j+(i-1)*jjp1 , we have
157	! j=modulo(stop_it,jjp1)
158	! i=1+(stop_it-j)/jjp1
159	! write(lunout,) " lon = ",rlonv(i)180./pi, " deg",
160	! & " lat = ",rlatu(j)*180./pi, " deg"
161	call abort_gcm("integrd_loc", "negative surface pressure", 1)
162	ENDIF
163
164	c$OMP END MASTER
165	c$OMP BARRIER
166	!write(,) 'integrd 170'
167	IF (.NOT. Checksum_all) THEN
168	call WriteField_v('int_vcov',vcov)
169	call WriteField_u('int_ucov',ucov)
170	call WriteField_u('int_teta',teta)
171	call WriteField_u('int_ps0',ps0)
172	call WriteField_u('int_masse',masse)
173	call WriteField_u('int_phis',phis)
174	call WriteField_v('int_vcovm1',vcovm1)
175	call WriteField_u('int_ucovm1',ucovm1)
176	call WriteField_u('int_tetam1',tetam1)
177	call WriteField_u('int_psm1',psm1)
178	call WriteField_u('int_massem1',massem1)
179
180	call WriteField_v('int_dv',dv)
181	call WriteField_u('int_du',du)
182	call WriteField_u('int_dteta',dteta)
183	call WriteField_u('int_dp',dp)
184	! call WriteField_u('int_finvmaold',finvmaold)
185	do j=1,nq
186	call WriteField_u('int_q'//trim(int2str(j)),
187	. q(:,:,j))
188	call WriteField_u('int_dq'//trim(int2str(j)),
189	. dq(:,:,j))
190	enddo
191	call abort_gcm("integrd_loc", "", 1)
192	ENDIF
193
194
195	c
196	!write(,) 'integrd 200'
197	C$OMP MASTER
198	if (pole_nord) THEN
199
200	DO ij = 1, iim
201	tppn(ij) = aire( ij ) * ps( ij )
202	ENDDO
203	tpn = SSUM(iim,tppn,1)/apoln
204	DO ij = 1, iip1
205	ps( ij ) = tpn
206	ENDDO
207
208	ENDIF
209
210	if (pole_sud) THEN
211
212	DO ij = 1, iim
213	tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
214	ENDDO
215	tps = SSUM(iim,tpps,1)/apols
216	DO ij = 1, iip1
217	ps(ij+ip1jm) = tps
218	ENDDO
219
220	ENDIF
221	c$OMP END MASTER
222	c$OMP BARRIER
223	!write(,) 'integrd 217'
224	c
225	c ... Calcul de la nouvelle masse d'air au dernier temps integre t+1 ...
226	c
227
228	CALL pression_loc ( ip1jmp1, ap, bp, ps, p )
229
230	c$OMP BARRIER
231	CALL massdair_loc ( p , masse )
232
233	! Ehouarn : we don't use/need finvmaold and finvmasse,
234	! so might as well not compute them
235	!c CALL SCOPY( ijp1llm , masse, 1, finvmasse, 1 )
236	! ijb=ij_begin
237	! ije=ij_end
238	!
239	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
240	! DO l = 1,llm
241	! finvmasse(ijb:ije,l)=masse(ijb:ije,l)
242	! ENDDO
243	!c$OMP END DO NOWAIT
244
245	! jjb=jj_begin
246	! jje=jj_end
247	! CALL filtreg_p( finvmasse,jjb_u,jje_u,jjb,jje, jjp1, llm,
248	! & -2, 2, .TRUE., 1 )
249	c
250
251	c ............ integration de ucov, vcov, h ..............
252
253	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
254	DO 10 l = 1,llm
255
256	ijb=ij_begin
257	ije=ij_end
258	if (pole_nord) ijb=ij_begin+iip1
259	if (pole_sud) ije=ij_end-iip1
260
261	DO 4 ij = ijb,ije
262	uscr( ij ) = ucov( ij,l )
263	ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
264	4 CONTINUE
265
266	ijb=ij_begin
267	ije=ij_end
268	if (pole_sud) ije=ij_end-iip1
269
270	DO 5 ij = ijb,ije
271	vscr( ij ) = vcov( ij,l )
272	vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
273	5 CONTINUE
274
275	ijb=ij_begin
276	ije=ij_end
277
278	DO 6 ij = ijb,ije
279	hscr( ij ) = teta(ij,l)
280	teta ( ij,l ) = tetam1(ij,l) * massem1(ij,l) / masse(ij,l)
281	$ + dt * dteta(ij,l) / masse(ij,l)
282	6 CONTINUE
283
284	c .... Calcul de la valeur moyenne, unique aux poles pour teta ......
285	c
286	c
287	!write(,) 'integrd 291'
288	IF (pole_nord) THEN
289
290	DO ij = 1, iim
291	tppn(ij) = aire( ij ) * teta( ij ,l)
292	ENDDO
293	tpn = SSUM(iim,tppn,1)/apoln
294
295	DO ij = 1, iip1
296	teta( ij ,l) = tpn
297	ENDDO
298
299	ENDIF
300
301	IF (pole_sud) THEN
302
303	DO ij = 1, iim
304	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
305	ENDDO
306	tps = SSUM(iim,tpps,1)/apols
307
308	DO ij = 1, iip1
309	teta(ij+ip1jm,l) = tps
310	ENDDO
311
312	ENDIF
313	c
314
315	IF(leapf) THEN
316	c CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
317	c CALL SCOPY ( ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
318	c CALL SCOPY ( ip1jmp1, hscr(1), 1, tetam1(1, l), 1 )
319	ijb=ij_begin
320	ije=ij_end
321	ucovm1(ijb:ije,l)=uscr(ijb:ije)
322	tetam1(ijb:ije,l)=hscr(ijb:ije)
323	if (pole_sud) ije=ij_end-iip1
324	vcovm1(ijb:ije,l)=vscr(ijb:ije)
325
326	END IF
327
328	10 CONTINUE
329	c$OMP END DO NOWAIT
330
331	c
332	c ....... integration de q ......
333	c
334	ijb=ij_begin
335	ije=ij_end
336
337	if (planet_type.eq."earth") then
338	! Earth-specific treatment of first 2 tracers (water)
339	c$OMP BARRIER
340	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
341	DO l = 1, llm
342	DO ij = ijb, ije
343	deltap(ij,l) = p(ij,l) - p(ij,l+1)
344	ENDDO
345	ENDDO
346
347	c$OMP END DO NOWAIT
348	c$OMP BARRIER
349
350	if (ok_iso_verif) then
351	call check_isotopes(q,ijb,ije,'integrd 342')
352	endif !if (ok_iso_verif) then
353
354	!write(,) 'integrd 341'
355	CALL qminimum_loc( q, nq, deltap )
356	!write(,) 'integrd 343'
357
358	if (ok_iso_verif) then
359	call check_isotopes(q,ijb,ije,'integrd 346')
360	endif !if (ok_iso_verif) then
361	c
362	c ..... Calcul de la valeur moyenne, unique aux poles pour q .....
363	c
364	c$OMP BARRIER
365	IF (pole_nord) THEN
366
367	DO iq = 1, nq
368
369	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
370	DO l = 1, llm
371
372	DO ij = 1, iim
373	qppn(ij) = aire( ij ) * q( ij ,l,iq)
374	ENDDO
375	qpn = SSUM(iim,qppn,1)/apoln
376
377	DO ij = 1, iip1
378	q( ij ,l,iq) = qpn
379	ENDDO
380
381	ENDDO
382	c$OMP END DO NOWAIT
383
384	ENDDO
385
386	ENDIF
387
388	IF (pole_sud) THEN
389
390	DO iq = 1, nq
391
392	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
393	DO l = 1, llm
394
395	DO ij = 1, iim
396	qpps(ij) = aire(ij+ip1jm) * q(ij+ip1jm,l,iq)
397	ENDDO
398	qps = SSUM(iim,qpps,1)/apols
399
400	DO ij = 1, iip1
401	q(ij+ip1jm,l,iq) = qps
402	ENDDO
403
404	ENDDO
405	c$OMP END DO NOWAIT
406
407	ENDDO
408
409	ENDIF
410
411	if (ok_iso_verif) then
412	call check_isotopes(q,ijb,ije,'integrd 409')
413	endif !if (ok_iso_verif) then
414
415	! Ehouarn: forget about finvmaold
416	!c CALL SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
417
418	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
419	! DO l = 1, llm
420	! finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)
421	! ENDDO
422	!c$OMP END DO NOWAIT
423
424	endif ! of if (planet_type.eq."earth")
425
426	c
427	c
428	c ..... FIN de l'integration de q .......
429
430	15 continue
431	!write(,) 'integrd 410'
432
433	c$OMP DO SCHEDULE(STATIC)
434	DO ij=ijb,ije
435	ps0(ij)=ps(ij)
436	ENDDO
437	c$OMP END DO NOWAIT
438
439	c .................................................................
440
441
442	IF( leapf ) THEN
443	c CALL SCOPY ( ip1jmp1 , pscr , 1, psm1 , 1 )
444	c CALL SCOPY ( ip1jmp1*llm, massescr, 1, massem1, 1 )
445	c$OMP DO SCHEDULE(STATIC)
446	DO ij=ijb,ije
447	psm1(ij)=pscr(ij)
448	ENDDO
449	c$OMP END DO NOWAIT
450
451	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
452	DO l = 1, llm
453	massem1(ijb:ije,l)=massescr(ijb:ije,l)
454	ENDDO
455	c$OMP END DO NOWAIT
456	END IF
457	c$OMP BARRIER
458	RETURN
459	END

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

Download in other formats:

Original Format