Context Navigation

leapfrog_p.F @ 2280

Last change on this file since 2280 was 2039, checked in by fhourdin, 11 years ago

Inclusion d'un appel supplementaire a geopot dans leapfrog dans

les versions dyn3d et dyn3dpar pour garantir la convergence
des 3 versions dyn3d/dyn3dpar/dyn3dmem.
La convergence fonctionne avec un compilateur gfortran 4.6.3 / openmpi
distribué sous ubuntu avec la nouvelle physique (NPv3.2) et le guidage.
options : gfortran -fdefault-real-8 -DNC_DOUBLE -O3
Une seule exception : si on guide avec dyn3dpar et openMP.

Correction du guidage dans dyn3dmem.

un print supprime dans cv3_toutines.F90

Rm : On ne devrait perdre la convergence numérique avec les précédentes
versions que pour la nouvelle physique puisque le géopotentiel
n'est utilisé dans la physique que par les theermiques et sisvat.

=====================================================================

Added call to geopot in leapfrog in dyn3d and dyn3dpar in order to

insure numerical convergence with dyn3dmem.
The convergence dyn3d/dyn3dpar/dyn3dmem is satistied with
gfortran 4.6.3 / openmpi, new physics (NPv3.2) and nudging.
options : gfortran -fdefault-real-8 -DNC_DOUBLE -O3
Only one exception : with nudging & dynd3dpar & openMP.

Bug fixing for nudging in dyn3dmem.

print supressed in cv3_routines.F90

Rm : Numerical convergence with previous releases should be lost for new
physics only since the geopotential is used only by thermals and
sisvat in the physics.

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: 55.8 KB

Line
1	!
2	! $Id: leapfrog_p.F 2039 2014-05-08 00:11:19Z musat $
3	!
4	c
5	c
6
7	SUBROUTINE leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
8	& time_0)
9
10	use exner_hyb_m, only: exner_hyb
11	use exner_milieu_m, only: exner_milieu
12	use exner_hyb_p_m, only: exner_hyb_p
13	use exner_milieu_p_m, only: exner_milieu_p
14	USE misc_mod
15	USE parallel_lmdz
16	USE times
17	USE mod_hallo
18	USE Bands
19	USE Write_Field
20	USE Write_Field_p
21	USE vampir
22	USE timer_filtre, ONLY : print_filtre_timer
23	USE infotrac, ONLY: nqtot
24	USE guide_p_mod, ONLY : guide_main
25	USE getparam
26	USE control_mod, ONLY: nday, day_step, planet_type, offline,
27	& iconser, iphysiq, iperiod, dissip_period,
28	& iecri, ip_ebil_dyn, ok_dynzon, ok_dyn_ins,
29	& periodav, ok_dyn_ave, output_grads_dyn,
30	& iapp_tracvl
31	IMPLICIT NONE
32
33	c ...... Version du 10/01/98 ..........
34
35	c avec coordonnees verticales hybrides
36	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
37
38	c=======================================================================
39	c
40	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
41	c -------
42	c
43	c Objet:
44	c ------
45	c
46	c GCM LMD nouvelle grille
47	c
48	c=======================================================================
49	c
50	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
51	c et possibilite d'appeler une fonction f(y) a derivee tangente
52	c hyperbolique a la place de la fonction a derivee sinusoidale.
53
54	c ... Possibilite de choisir le shema pour l'advection de
55	c q , en modifiant iadv dans traceur.def (10/02) .
56	c
57	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
58	c Pour Van-Leer iadv=10
59	c
60	c-----------------------------------------------------------------------
61	c Declarations:
62	c -------------
63
64	#include "dimensions.h"
65	#include "paramet.h"
66	#include "comconst.h"
67	#include "comdissnew.h"
68	#include "comvert.h"
69	#include "comgeom.h"
70	#include "logic.h"
71	#include "temps.h"
72	#include "ener.h"
73	#include "description.h"
74	#include "serre.h"
75	!#include "com_io_dyn.h"
76	#include "iniprint.h"
77	#include "academic.h"
78
79	INTEGER,PARAMETER :: longcles = 20
80	REAL,INTENT(IN) :: clesphy0( longcles ) ! not used
81	REAL,INTENT(IN) :: time_0 ! not used
82
83	c dynamical variables:
84	REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
85	REAL,INTENT(INOUT) :: vcov(ip1jm,llm) ! meridional covariant wind
86	REAL,INTENT(INOUT) :: teta(ip1jmp1,llm) ! potential temperature
87	REAL,INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
88	REAL,INTENT(INOUT) :: masse(ip1jmp1,llm) ! air mass
89	REAL,INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
90	REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
91
92	REAL,SAVE :: p (ip1jmp1,llmp1 ) ! interlayer pressure
93	REAL,SAVE :: pks(ip1jmp1) ! exner at the surface
94	REAL,SAVE :: pk(ip1jmp1,llm) ! exner at mid-layer
95	REAL,SAVE :: pkf(ip1jmp1,llm) ! filtered exner at mid-layer
96	REAL,SAVE :: phi(ip1jmp1,llm) ! geopotential
97	REAL,SAVE :: w(ip1jmp1,llm) ! vertical velocity
98
99	real zqmin,zqmax
100
101	c variables dynamiques intermediaire pour le transport
102	REAL,SAVE :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
103
104	c variables dynamiques au pas -1
105	REAL,SAVE :: vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
106	REAL,SAVE :: tetam1(ip1jmp1,llm),psm1(ip1jmp1)
107	REAL,SAVE :: massem1(ip1jmp1,llm)
108
109	c tendances dynamiques
110	REAL,SAVE :: dv(ip1jm,llm),du(ip1jmp1,llm)
111	REAL,SAVE :: dteta(ip1jmp1,llm),dp(ip1jmp1)
112	REAL,DIMENSION(:,:,:), ALLOCATABLE, SAVE :: dq
113
114	c tendances de la dissipation
115	REAL,SAVE :: dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
116	REAL,SAVE :: dtetadis(ip1jmp1,llm)
117
118	c tendances physiques
119	REAL,SAVE :: dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
120	REAL,SAVE :: dtetafi(ip1jmp1,llm)
121	REAL,SAVE :: dpfi(ip1jmp1)
122	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi
123
124	c variables pour le fichier histoire
125	REAL dtav ! intervalle de temps elementaire
126
127	REAL tppn(iim),tpps(iim),tpn,tps
128	c
129	INTEGER itau,itaufinp1,iav
130	! INTEGER iday ! jour julien
131	REAL time
132
133	REAL SSUM
134	! REAL,SAVE :: finvmaold(ip1jmp1,llm)
135
136	cym LOGICAL lafin
137	LOGICAL :: lafin
138	INTEGER ij,iq,l
139	INTEGER ik
140
141	real time_step, t_wrt, t_ops
142
143	! jD_cur: jour julien courant
144	! jH_cur: heure julienne courante
145	REAL :: jD_cur, jH_cur
146	INTEGER :: an, mois, jour
147	REAL :: secondes
148
149	logical :: physic
150	LOGICAL first,callinigrads
151
152	data callinigrads/.true./
153	character*10 string10
154
155	REAL,SAVE :: flxw(ip1jmp1,llm) ! flux de masse verticale
156
157	c+jld variables test conservation energie
158	REAL,SAVE :: ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
159	C Tendance de la temp. potentiel d (theta)/ d t due a la
160	C tansformation d'energie cinetique en energie thermique
161	C cree par la dissipation
162	REAL,SAVE :: dtetaecdt(ip1jmp1,llm)
163	REAL,SAVE :: vcont(ip1jm,llm),ucont(ip1jmp1,llm)
164	REAL,SAVE :: vnat(ip1jm,llm),unat(ip1jmp1,llm)
165	REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
166	CHARACTER*15 ztit
167	! INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
168	! SAVE ip_ebil_dyn
169	! DATA ip_ebil_dyn/0/
170	c-jld
171
172	character*80 dynhist_file, dynhistave_file
173	character(len=*),parameter :: modname="leapfrog"
174	character*80 abort_message
175
176
177	logical,PARAMETER :: dissip_conservative=.TRUE.
178
179	INTEGER testita
180	PARAMETER (testita = 9)
181
182	logical , parameter :: flag_verif = .false.
183
184	c declaration liees au parallelisme
185	INTEGER :: ierr
186	LOGICAL :: FirstCaldyn
187	LOGICAL :: FirstPhysic
188	INTEGER :: ijb,ije,j,i
189	type(Request) :: TestRequest
190	type(Request) :: Request_Dissip
191	type(Request) :: Request_physic
192	REAL,SAVE :: dvfi_tmp(iip1,llm),dufi_tmp(iip1,llm)
193	REAL,SAVE :: dtetafi_tmp(iip1,llm)
194	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi_tmp
195	REAL,SAVE :: dpfi_tmp(iip1)
196
197	INTEGER :: true_itau
198	INTEGER :: iapptrac
199	INTEGER :: AdjustCount
200	! INTEGER :: var_time
201	LOGICAL :: ok_start_timer=.FALSE.
202	LOGICAL, SAVE :: firstcall=.TRUE.
203
204	c$OMP MASTER
205	ItCount=0
206	c$OMP END MASTER
207	true_itau=0
208	FirstCaldyn=.TRUE.
209	FirstPhysic=.TRUE.
210	iapptrac=0
211	AdjustCount = 0
212	lafin=.false.
213
214	if (nday>=0) then
215	itaufin = nday*day_step
216	else
217	itaufin = -nday
218	endif
219
220	itaufinp1 = itaufin +1
221
222	itau = 0
223	physic=.true.
224	if (iflag_phys==0.or.iflag_phys==2) physic=.false.
225	! iday = day_ini+itau/day_step
226	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
227	! IF(time.GT.1.) THEN
228	! time = time-1.
229	! iday = iday+1
230	! ENDIF
231
232	c Allocate variables depending on dynamic variable nqtot
233	c$OMP MASTER
234	IF (firstcall) THEN
235	firstcall=.FALSE.
236	ALLOCATE(dq(ip1jmp1,llm,nqtot))
237	ALLOCATE(dqfi(ip1jmp1,llm,nqtot))
238	ALLOCATE(dqfi_tmp(iip1,llm,nqtot))
239	END IF
240	c$OMP END MASTER
241	c$OMP BARRIER
242
243	c-----------------------------------------------------------------------
244	c On initialise la pression et la fonction d'Exner :
245	c --------------------------------------------------
246
247	c$OMP MASTER
248	dq(:,:,:)=0.
249	CALL pression ( ip1jmp1, ap, bp, ps, p )
250	if (pressure_exner) then
251	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
252	else
253	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
254	endif
255	c$OMP END MASTER
256	c-----------------------------------------------------------------------
257	c Debut de l'integration temporelle:
258	c ----------------------------------
259	c et du parallelisme !!
260
261	1 CONTINUE ! Matsuno Forward step begins here
262
263	jD_cur = jD_ref + day_ini - day_ref + &
264	& itau/day_step
265	jH_cur = jH_ref + start_time + &
266	& mod(itau,day_step)/float(day_step)
267	if (jH_cur > 1.0 ) then
268	jD_cur = jD_cur +1.
269	jH_cur = jH_cur -1.
270	endif
271
272
273	#ifdef CPP_IOIPSL
274	if (ok_guide) then
275	!$OMP MASTER
276	call guide_main(itau,ucov,vcov,teta,q,masse,ps)
277	!$OMP END MASTER
278	!$OMP BARRIER
279	endif
280	#endif
281
282	c
283	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
284	c CALL test_period ( ucov,vcov,teta,q,p,phis )
285	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
286	c ENDIF
287	c
288	cym CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
289	cym CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
290	cym CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
291	cym CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
292	cym CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
293
294	if (FirstCaldyn) then
295	c$OMP MASTER
296	ucovm1=ucov
297	vcovm1=vcov
298	tetam1= teta
299	massem1= masse
300	psm1= ps
301
302	! Ehouarn: finvmaold is actually not used
303	! finvmaold = masse
304	! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
305	c$OMP END MASTER
306	c$OMP BARRIER
307	else
308	! Save fields obtained at previous time step as '...m1'
309	ijb=ij_begin
310	ije=ij_end
311
312	c$OMP MASTER
313	psm1 (ijb:ije) = ps (ijb:ije)
314	c$OMP END MASTER
315
316	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
317	DO l=1,llm
318	ije=ij_end
319	ucovm1 (ijb:ije,l) = ucov (ijb:ije,l)
320	tetam1 (ijb:ije,l) = teta (ijb:ije,l)
321	massem1 (ijb:ije,l) = masse (ijb:ije,l)
322	! finvmaold(ijb:ije,l)=masse(ijb:ije,l)
323
324	if (pole_sud) ije=ij_end-iip1
325	vcovm1(ijb:ije,l) = vcov (ijb:ije,l)
326
327
328	ENDDO
329	c$OMP ENDDO
330
331	! Ehouarn: finvmaold not used
332	! CALL filtreg_p ( finvmaold ,jj_begin,jj_end,jjp1,
333	! . llm, -2,2, .TRUE., 1 )
334
335	endif ! of if (FirstCaldyn)
336
337	forward = .TRUE.
338	leapf = .FALSE.
339	dt = dtvr
340
341	c ... P.Le Van .26/04/94 ....
342
343	cym CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
344	cym CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
345
346	cym ne sert a rien
347	cym call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
348
349	2 CONTINUE ! Matsuno backward or leapfrog step begins here
350
351	c$OMP MASTER
352	ItCount=ItCount+1
353	if (MOD(ItCount,1)==1) then
354	debug=.true.
355	else
356	debug=.false.
357	endif
358	c$OMP END MASTER
359	c-----------------------------------------------------------------------
360
361	c date:
362	c -----
363
364
365	c gestion des appels de la physique et des dissipations:
366	c ------------------------------------------------------
367	c
368	c ... P.Le Van ( 6/02/95 ) ....
369
370	apphys = .FALSE.
371	statcl = .FALSE.
372	conser = .FALSE.
373	apdiss = .FALSE.
374
375	IF( purmats ) THEN
376	! Purely Matsuno time stepping
377	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
378	IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward )
379	s apdiss = .TRUE.
380	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
381	s .and. physic ) apphys = .TRUE.
382	ELSE
383	! Leapfrog/Matsuno time stepping
384	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
385	IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
386	s apdiss = .TRUE.
387	IF( MOD(itau+1,iphysiq).EQ.0.AND.physic) apphys=.TRUE.
388	END IF
389
390	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
391	! supress dissipation step
392	if (llm.eq.1) then
393	apdiss=.false.
394	endif
395
396	cym ---> Pour le moment
397	cym apphys = .FALSE.
398	statcl = .FALSE.
399	conser = .FALSE. ! ie: no output of control variables to stdout in //
400
401	if (firstCaldyn) then
402	c$OMP MASTER
403	call SetDistrib(jj_Nb_Caldyn)
404	c$OMP END MASTER
405	c$OMP BARRIER
406	firstCaldyn=.FALSE.
407	cym call InitTime
408	c$OMP MASTER
409	call Init_timer
410	c$OMP END MASTER
411	endif
412
413	c$OMP MASTER
414	IF (ok_start_timer) THEN
415	CALL InitTime
416	ok_start_timer=.FALSE.
417	ENDIF
418	c$OMP END MASTER
419
420	if (Adjust) then
421	c$OMP MASTER
422	AdjustCount=AdjustCount+1
423	if (iapptrac==iapp_tracvl .and. (forward. OR . leapf)
424	& .and. itau/iphysiq>2 .and. Adjustcount>30) then
425	AdjustCount=0
426	call allgather_timer_average
427
428	if (prt_level > 9) then
429
430	print ,'********************************'
431	print ,'*** TIMER CALDYN ****'
432	do i=0,mpi_size-1
433	print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
434	& ' : temps moyen :',
435	& timer_average(jj_nb_caldyn(i),timer_caldyn,i),
436	& '+-',timer_delta(jj_nb_caldyn(i),timer_caldyn,i)
437	enddo
438
439	print ,'********************************'
440	print ,'*** TIMER VANLEER ****'
441	do i=0,mpi_size-1
442	print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
443	& ' : temps moyen :',
444	& timer_average(jj_nb_vanleer(i),timer_vanleer,i),
445	& '+-',timer_delta(jj_nb_vanleer(i),timer_vanleer,i)
446	enddo
447
448	print ,'********************************'
449	print ,'*** TIMER DISSIP ****'
450	do i=0,mpi_size-1
451	print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
452	& ' : temps moyen :',
453	& timer_average(jj_nb_dissip(i),timer_dissip,i),
454	& '+-',timer_delta(jj_nb_dissip(i),timer_dissip,i)
455	enddo
456
457	if (mpi_rank==0) call WriteBands
458
459	endif
460
461	call AdjustBands_caldyn
462	if (mpi_rank==0) call WriteBands
463
464	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
465	& jj_Nb_caldyn,0,0,TestRequest)
466	call Register_SwapFieldHallo(ucovm1,ucovm1,ip1jmp1,llm,
467	& jj_Nb_caldyn,0,0,TestRequest)
468	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
469	& jj_Nb_caldyn,0,0,TestRequest)
470	call Register_SwapFieldHallo(vcovm1,vcovm1,ip1jm,llm,
471	& jj_Nb_caldyn,0,0,TestRequest)
472	call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
473	& jj_Nb_caldyn,0,0,TestRequest)
474	call Register_SwapFieldHallo(tetam1,tetam1,ip1jmp1,llm,
475	& jj_Nb_caldyn,0,0,TestRequest)
476	call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
477	& jj_Nb_caldyn,0,0,TestRequest)
478	call Register_SwapFieldHallo(massem1,massem1,ip1jmp1,llm,
479	& jj_Nb_caldyn,0,0,TestRequest)
480	call Register_SwapFieldHallo(ps,ps,ip1jmp1,1,
481	& jj_Nb_caldyn,0,0,TestRequest)
482	call Register_SwapFieldHallo(psm1,psm1,ip1jmp1,1,
483	& jj_Nb_caldyn,0,0,TestRequest)
484	call Register_SwapFieldHallo(pkf,pkf,ip1jmp1,llm,
485	& jj_Nb_caldyn,0,0,TestRequest)
486	call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
487	& jj_Nb_caldyn,0,0,TestRequest)
488	call Register_SwapFieldHallo(pks,pks,ip1jmp1,1,
489	& jj_Nb_caldyn,0,0,TestRequest)
490	call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
491	& jj_Nb_caldyn,0,0,TestRequest)
492	call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
493	& jj_Nb_caldyn,0,0,TestRequest)
494	! call Register_SwapFieldHallo(finvmaold,finvmaold,ip1jmp1,llm,
495	! & jj_Nb_caldyn,0,0,TestRequest)
496
497	do j=1,nqtot
498	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
499	& jj_nb_caldyn,0,0,TestRequest)
500	enddo
501
502	call SetDistrib(jj_nb_caldyn)
503	call SendRequest(TestRequest)
504	call WaitRequest(TestRequest)
505
506	call AdjustBands_dissip
507	call AdjustBands_physic
508
509	endif
510	c$OMP END MASTER
511	endif
512
513
514
515	c-----------------------------------------------------------------------
516	c calcul des tendances dynamiques:
517	c --------------------------------
518	c$OMP BARRIER
519	c$OMP MASTER
520	call VTb(VThallo)
521	c$OMP END MASTER
522
523	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,TestRequest)
524	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,TestRequest)
525	call Register_Hallo(teta,ip1jmp1,llm,1,1,1,1,TestRequest)
526	call Register_Hallo(ps,ip1jmp1,1,1,2,2,1,TestRequest)
527	call Register_Hallo(pkf,ip1jmp1,llm,1,1,1,1,TestRequest)
528	call Register_Hallo(pk,ip1jmp1,llm,1,1,1,1,TestRequest)
529	call Register_Hallo(pks,ip1jmp1,1,1,1,1,1,TestRequest)
530	call Register_Hallo(p,ip1jmp1,llmp1,1,1,1,1,TestRequest)
531
532	c do j=1,nqtot
533	c call Register_Hallo(q(1,1,j),ip1jmp1,llm,1,1,1,1,
534	c * TestRequest)
535	c enddo
536
537	call SendRequest(TestRequest)
538	c$OMP BARRIER
539	call WaitRequest(TestRequest)
540
541	c$OMP MASTER
542	call VTe(VThallo)
543	c$OMP END MASTER
544	c$OMP BARRIER
545
546	if (debug) then
547	!$OMP BARRIER
548	!$OMP MASTER
549	call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
550	call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
551	call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
552	call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
553	call WriteField_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
554	call WriteField_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
555	call WriteField_p('pks',reshape(pks,(/iip1,jmp1/)))
556	call WriteField_p('pkf',reshape(pkf,(/iip1,jmp1,llm/)))
557	call WriteField_p('phis',reshape(phis,(/iip1,jmp1/)))
558	do j=1,nqtot
559	call WriteField_p('q'//trim(int2str(j)),
560	. reshape(q(:,:,j),(/iip1,jmp1,llm/)))
561	enddo
562	!$OMP END MASTER
563	c$OMP BARRIER
564	endif
565
566
567	True_itau=True_itau+1
568
569	c$OMP MASTER
570	IF (prt_level>9) THEN
571	WRITE(lunout,*)"leapfrog_p: Iteration No",True_itau
572	ENDIF
573
574
575	call start_timer(timer_caldyn)
576
577	! compute geopotential phi()
578	CALL geopot_p ( ip1jmp1, teta , pk , pks, phis , phi )
579
580
581	call VTb(VTcaldyn)
582	c$OMP END MASTER
583	! var_time=time+iday-day_ini
584
585	c$OMP BARRIER
586	! CALL FTRACE_REGION_BEGIN("caldyn")
587	time = jD_cur + jH_cur
588	CALL caldyn_p
589	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
590	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
591
592	! CALL FTRACE_REGION_END("caldyn")
593
594	c$OMP MASTER
595	call VTe(VTcaldyn)
596	c$OMP END MASTER
597
598	cc$OMP BARRIER
599	cc$OMP MASTER
600	! call WriteField_p('du',reshape(du,(/iip1,jmp1,llm/)))
601	! call WriteField_p('dv',reshape(dv,(/iip1,jjm,llm/)))
602	! call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
603	! call WriteField_p('dp',reshape(dp,(/iip1,jmp1/)))
604	! call WriteField_p('w',reshape(w,(/iip1,jmp1,llm/)))
605	! call WriteField_p('pbaru',reshape(pbaru,(/iip1,jmp1,llm/)))
606	! call WriteField_p('pbarv',reshape(pbarv,(/iip1,jjm,llm/)))
607	! call WriteField_p('p',reshape(p,(/iip1,jmp1,llmp1/)))
608	! call WriteField_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
609	! call WriteField_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
610	cc$OMP END MASTER
611
612	c-----------------------------------------------------------------------
613	c calcul des tendances advection des traceurs (dont l'humidite)
614	c -------------------------------------------------------------
615
616	IF( forward. OR . leapf ) THEN
617	! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
618	CALL caladvtrac_p(q,pbaru,pbarv,
619	* p, masse, dq, teta,
620	. flxw,pk, iapptrac)
621
622	C Stokage du flux de masse pour traceurs OFF-LINE
623	IF (offline .AND. .NOT. adjust) THEN
624	CALL fluxstokenc_p(pbaru,pbarv,masse,teta,phi,phis,
625	. dtvr, itau)
626	ENDIF
627
628	ENDIF ! of IF( forward. OR . leapf )
629
630	c-----------------------------------------------------------------------
631	c integrations dynamique et traceurs:
632	c ----------------------------------
633
634	c$OMP MASTER
635	call VTb(VTintegre)
636	c$OMP END MASTER
637	c call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
638	c call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
639	c call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
640	c call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
641	c call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
642	c call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
643	c call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
644	c call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
645	cc$OMP PARALLEL DEFAULT(SHARED)
646	c$OMP BARRIER
647	! CALL FTRACE_REGION_BEGIN("integrd")
648
649	CALL integrd_p ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
650	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
651	! $ finvmaold )
652
653	! CALL FTRACE_REGION_END("integrd")
654	c$OMP BARRIER
655	cc$OMP MASTER
656	c call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
657	c call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
658	c call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
659	c call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
660	c call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
661	c call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
662	c call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
663	c call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
664	c
665	c call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
666	c do j=1,nqtot
667	c call WriteField_p('q'//trim(int2str(j)),
668	c . reshape(q(:,:,j),(/iip1,jmp1,llm/)))
669	c call WriteField_p('dq'//trim(int2str(j)),
670	c . reshape(dq(:,:,j),(/iip1,jmp1,llm/)))
671	c enddo
672	cc$OMP END MASTER
673
674
675	c$OMP MASTER
676	call VTe(VTintegre)
677	c$OMP END MASTER
678	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
679	c
680	c-----------------------------------------------------------------------
681	c calcul des tendances physiques:
682	c -------------------------------
683	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
684	c
685	IF( purmats ) THEN
686	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
687	ELSE
688	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
689	ENDIF
690
691	cc$OMP END PARALLEL
692
693	c
694	c
695	IF( apphys ) THEN
696	c
697	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
698	c
699	cc$OMP PARALLEL DEFAULT(SHARED)
700	cc$OMP+ PRIVATE(rdaym_ini,rdayvrai,ijb,ije)
701
702	c$OMP MASTER
703	call suspend_timer(timer_caldyn)
704
705	if (prt_level >= 10) then
706	write(lunout,*)
707	& 'leapfrog_p: Entree dans la physique : Iteration No ',true_itau
708	endif
709	c$OMP END MASTER
710
711	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
712
713	c$OMP BARRIER
714	if (pressure_exner) then
715	CALL exner_hyb_p( ip1jmp1, ps, p,pks, pk, pkf )
716	else
717	CALL exner_milieu_p( ip1jmp1, ps, p, pks, pk, pkf )
718	endif
719	! Appel a geopot ajoute le 2014/05/08 pour garantir la convergence numerique
720	! avec dyn3dmem
721	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
722	c$OMP BARRIER
723	jD_cur = jD_ref + day_ini - day_ref
724	$ + itau/day_step
725
726	IF (planet_type .eq."generic") THEN
727	! AS: we make jD_cur to be pday
728	jD_cur = int(day_ini + itau/day_step)
729	ENDIF
730
731	jH_cur = jH_ref + start_time + &
732	& mod(itau,day_step)/float(day_step)
733	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
734	if (jH_cur > 1.0 ) then
735	jD_cur = jD_cur +1.
736	jH_cur = jH_cur -1.
737	endif
738
739	c rajout debug
740	c lafin = .true.
741
742
743	c Inbterface avec les routines de phylmd (phymars ... )
744	c -----------------------------------------------------
745
746	c+jld
747
748	c Diagnostique de conservation de l'energie : initialisation
749	IF (ip_ebil_dyn.ge.1 ) THEN
750	ztit='bil dyn'
751	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
752	IF (planet_type.eq."earth") THEN
753	#ifdef CPP_EARTH
754	CALL diagedyn(ztit,2,1,1,dtphys
755	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
756	#endif
757	ENDIF
758	ENDIF
759	c-jld
760	c$OMP BARRIER
761	c$OMP MASTER
762	call VTb(VThallo)
763	c$OMP END MASTER
764
765	call SetTag(Request_physic,800)
766
767	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
768	* jj_Nb_physic,2,2,Request_physic)
769
770	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
771	* jj_Nb_physic,2,2,Request_physic)
772
773	call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
774	* jj_Nb_physic,2,2,Request_physic)
775
776	call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
777	* jj_Nb_physic,1,2,Request_physic)
778
779	call Register_SwapFieldHallo(ps,ps,ip1jmp1,1,
780	* jj_Nb_physic,2,2,Request_physic)
781
782	call Register_SwapFieldHallo(p,p,ip1jmp1,llmp1,
783	* jj_Nb_physic,2,2,Request_physic)
784
785	call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
786	* jj_Nb_physic,2,2,Request_physic)
787
788	call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
789	* jj_Nb_physic,2,2,Request_physic)
790
791	call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
792	* jj_Nb_physic,2,2,Request_physic)
793
794	call Register_SwapFieldHallo(w,w,ip1jmp1,llm,
795	* jj_Nb_physic,2,2,Request_physic)
796
797	c call SetDistrib(jj_nb_vanleer)
798	do j=1,nqtot
799
800	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
801	* jj_Nb_physic,2,2,Request_physic)
802	enddo
803
804	call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
805	* jj_Nb_physic,2,2,Request_physic)
806
807	call SendRequest(Request_Physic)
808	c$OMP BARRIER
809	call WaitRequest(Request_Physic)
810
811	c$OMP BARRIER
812	c$OMP MASTER
813	call SetDistrib(jj_nb_Physic)
814	call VTe(VThallo)
815
816	call VTb(VTphysiq)
817	c$OMP END MASTER
818	c$OMP BARRIER
819
820	cc$OMP MASTER
821	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
822	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
823	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
824	c call WriteField_p('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
825	c call WriteField_p('pkfi',reshape(pk,(/iip1,jmp1,llm/)))
826	cc$OMP END MASTER
827	cc$OMP BARRIER
828	! CALL FTRACE_REGION_BEGIN("calfis")
829	CALL calfis_p(lafin ,jD_cur, jH_cur,
830	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
831	$ du,dv,dteta,dq,
832	$ flxw,
833	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
834	! CALL FTRACE_REGION_END("calfis")
835	ijb=ij_begin
836	ije=ij_end
837	if ( .not. pole_nord) then
838	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
839	DO l=1,llm
840	dufi_tmp(1:iip1,l) = dufi(ijb:ijb+iim,l)
841	dvfi_tmp(1:iip1,l) = dvfi(ijb:ijb+iim,l)
842	dtetafi_tmp(1:iip1,l)= dtetafi(ijb:ijb+iim,l)
843	dqfi_tmp(1:iip1,l,:) = dqfi(ijb:ijb+iim,l,:)
844	ENDDO
845	c$OMP END DO NOWAIT
846
847	c$OMP MASTER
848	dpfi_tmp(1:iip1) = dpfi(ijb:ijb+iim)
849	c$OMP END MASTER
850	endif ! of if ( .not. pole_nord)
851
852	c$OMP BARRIER
853	c$OMP MASTER
854	call SetDistrib(jj_nb_Physic_bis)
855
856	call VTb(VThallo)
857	c$OMP END MASTER
858	c$OMP BARRIER
859
860	call Register_Hallo(dufi,ip1jmp1,llm,
861	* 1,0,0,1,Request_physic)
862
863	call Register_Hallo(dvfi,ip1jm,llm,
864	* 1,0,0,1,Request_physic)
865
866	call Register_Hallo(dtetafi,ip1jmp1,llm,
867	* 1,0,0,1,Request_physic)
868
869	call Register_Hallo(dpfi,ip1jmp1,1,
870	* 1,0,0,1,Request_physic)
871
872	do j=1,nqtot
873	call Register_Hallo(dqfi(1,1,j),ip1jmp1,llm,
874	* 1,0,0,1,Request_physic)
875	enddo
876
877	call SendRequest(Request_Physic)
878	c$OMP BARRIER
879	call WaitRequest(Request_Physic)
880
881	c$OMP BARRIER
882	c$OMP MASTER
883	call VTe(VThallo)
884
885	call SetDistrib(jj_nb_Physic)
886	c$OMP END MASTER
887	c$OMP BARRIER
888	ijb=ij_begin
889	if (.not. pole_nord) then
890
891	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
892	DO l=1,llm
893	dufi(ijb:ijb+iim,l) = dufi(ijb:ijb+iim,l)+dufi_tmp(1:iip1,l)
894	dvfi(ijb:ijb+iim,l) = dvfi(ijb:ijb+iim,l)+dvfi_tmp(1:iip1,l)
895	dtetafi(ijb:ijb+iim,l) = dtetafi(ijb:ijb+iim,l)
896	& +dtetafi_tmp(1:iip1,l)
897	dqfi(ijb:ijb+iim,l,:) = dqfi(ijb:ijb+iim,l,:)
898	& + dqfi_tmp(1:iip1,l,:)
899	ENDDO
900	c$OMP END DO NOWAIT
901
902	c$OMP MASTER
903	dpfi(ijb:ijb+iim) = dpfi(ijb:ijb+iim)+ dpfi_tmp(1:iip1)
904	c$OMP END MASTER
905
906	endif ! of if (.not. pole_nord)
907	c$OMP BARRIER
908	cc$OMP MASTER
909	c call WriteField_p('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
910	c call WriteField_p('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
911	c call WriteField_p('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
912	c call WriteField_p('dpfi',reshape(dpfi,(/iip1,jmp1/)))
913	cc$OMP END MASTER
914	c
915	c do j=1,nqtot
916	c call WriteField_p('dqfi'//trim(int2str(j)),
917	c . reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
918	c enddo
919
920	c ajout des tendances physiques:
921	c ------------------------------
922	CALL addfi_p( dtphys, leapf, forward ,
923	$ ucov, vcov, teta , q ,ps ,
924	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
925	! since addfi updates ps(), also update p(), masse() and pk()
926	CALL pression_p(ip1jmp1,ap,bp,ps,p)
927	c$OMP BARRIER
928	CALL massdair_p(p,masse)
929	c$OMP BARRIER
930	if (pressure_exner) then
931	CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk,pkf)
932	else
933	CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk,pkf)
934	endif
935	c$OMP BARRIER
936
937	IF (ok_strato) THEN
938	CALL top_bound_p(vcov,ucov,teta,masse,dtphys)
939	ENDIF
940
941	c$OMP BARRIER
942	c$OMP MASTER
943	call VTe(VTphysiq)
944
945	call VTb(VThallo)
946	c$OMP END MASTER
947
948	call SetTag(Request_physic,800)
949	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
950	* jj_Nb_caldyn,Request_physic)
951
952	call Register_SwapField(vcov,vcov,ip1jm,llm,
953	* jj_Nb_caldyn,Request_physic)
954
955	call Register_SwapField(teta,teta,ip1jmp1,llm,
956	* jj_Nb_caldyn,Request_physic)
957
958	call Register_SwapField(masse,masse,ip1jmp1,llm,
959	* jj_Nb_caldyn,Request_physic)
960
961	call Register_SwapField(ps,ps,ip1jmp1,1,
962	* jj_Nb_caldyn,Request_physic)
963
964	call Register_SwapField(p,p,ip1jmp1,llmp1,
965	* jj_Nb_caldyn,Request_physic)
966
967	call Register_SwapField(pk,pk,ip1jmp1,llm,
968	* jj_Nb_caldyn,Request_physic)
969
970	call Register_SwapField(phis,phis,ip1jmp1,1,
971	* jj_Nb_caldyn,Request_physic)
972
973	call Register_SwapField(phi,phi,ip1jmp1,llm,
974	* jj_Nb_caldyn,Request_physic)
975
976	call Register_SwapField(w,w,ip1jmp1,llm,
977	* jj_Nb_caldyn,Request_physic)
978
979	do j=1,nqtot
980
981	call Register_SwapField(q(1,1,j),q(1,1,j),ip1jmp1,llm,
982	* jj_Nb_caldyn,Request_physic)
983
984	enddo
985
986	call SendRequest(Request_Physic)
987	c$OMP BARRIER
988	call WaitRequest(Request_Physic)
989
990	c$OMP BARRIER
991	c$OMP MASTER
992	call VTe(VThallo)
993	call SetDistrib(jj_Nb_caldyn)
994	c$OMP END MASTER
995	c$OMP BARRIER
996	c
997	c Diagnostique de conservation de l'energie : difference
998	IF (ip_ebil_dyn.ge.1 ) THEN
999	ztit='bil phys'
1000	CALL diagedyn(ztit,2,1,1,dtphys
1001	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
1002	ENDIF
1003
1004	cc$OMP MASTER
1005	c if (debug) then
1006	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
1007	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
1008	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
1009	c endif
1010	cc$OMP END MASTER
1011
1012
1013	c-jld
1014	c$OMP MASTER
1015	call resume_timer(timer_caldyn)
1016	if (FirstPhysic) then
1017	ok_start_timer=.TRUE.
1018	FirstPhysic=.false.
1019	endif
1020	c$OMP END MASTER
1021	ENDIF ! of IF( apphys )
1022
1023	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
1024	! Academic case : Simple friction and Newtonan relaxation
1025	! -------------------------------------------------------
1026	c$OMP MASTER
1027	if (FirstPhysic) then
1028	ok_start_timer=.TRUE.
1029	FirstPhysic=.false.
1030	endif
1031	c$OMP END MASTER
1032
1033	ijb=ij_begin
1034	ije=ij_end
1035	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1036	do l=1,llm
1037	teta(ijb:ije,l)=teta(ijb:ije,l)-dtvr*
1038	& (teta(ijb:ije,l)-tetarappel(ijb:ije,l))*
1039	& (knewt_g+knewt_t(l)*clat4(ijb:ije))
1040	enddo ! of do l=1,llm
1041	!$OMP END DO
1042
1043	!$OMP MASTER
1044	if (planet_type.eq."giant") then
1045	! add an intrinsic heat flux at the base of the atmosphere
1046	teta(ijb:ije,1) = teta(ijb:ije,1)
1047	& + dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1)
1048	endif
1049	!$OMP END MASTER
1050	!$OMP BARRIER
1051
1052	call Register_Hallo(ucov,ip1jmp1,llm,0,1,1,0,Request_Physic)
1053	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Physic)
1054	call SendRequest(Request_Physic)
1055	c$OMP BARRIER
1056	call WaitRequest(Request_Physic)
1057	c$OMP BARRIER
1058	call friction_p(ucov,vcov,dtvr)
1059	!$OMP BARRIER
1060
1061	! Sponge layer (if any)
1062	IF (ok_strato) THEN
1063	CALL top_bound_p(vcov,ucov,teta,masse,dtvr)
1064	!$OMP BARRIER
1065	ENDIF ! of IF (ok_strato)
1066	ENDIF ! of IF(iflag_phys.EQ.2)
1067
1068
1069	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
1070	c$OMP BARRIER
1071	if (pressure_exner) then
1072	CALL exner_hyb_p( ip1jmp1, ps, p, pks, pk, pkf )
1073	else
1074	CALL exner_milieu_p( ip1jmp1, ps, p, pks, pk, pkf )
1075	endif
1076	c$OMP BARRIER
1077	CALL massdair_p(p,masse)
1078	c$OMP BARRIER
1079
1080	cc$OMP END PARALLEL
1081
1082	c-----------------------------------------------------------------------
1083	c dissipation horizontale et verticale des petites echelles:
1084	c ----------------------------------------------------------
1085
1086	IF(apdiss) THEN
1087	cc$OMP PARALLEL DEFAULT(SHARED)
1088	cc$OMP+ PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
1089	c$OMP MASTER
1090	call suspend_timer(timer_caldyn)
1091
1092	c print*,'Entree dans la dissipation : Iteration No ',true_itau
1093	c calcul de l'energie cinetique avant dissipation
1094	c print *,'Passage dans la dissipation'
1095
1096	call VTb(VThallo)
1097	c$OMP END MASTER
1098
1099	c$OMP BARRIER
1100
1101	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
1102	* jj_Nb_dissip,1,1,Request_dissip)
1103
1104	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
1105	* jj_Nb_dissip,1,1,Request_dissip)
1106
1107	call Register_SwapField(teta,teta,ip1jmp1,llm,
1108	* jj_Nb_dissip,Request_dissip)
1109
1110	call Register_SwapField(p,p,ip1jmp1,llmp1,
1111	* jj_Nb_dissip,Request_dissip)
1112
1113	call Register_SwapField(pk,pk,ip1jmp1,llm,
1114	* jj_Nb_dissip,Request_dissip)
1115
1116	call SendRequest(Request_dissip)
1117	c$OMP BARRIER
1118	call WaitRequest(Request_dissip)
1119
1120	c$OMP BARRIER
1121	c$OMP MASTER
1122	call SetDistrib(jj_Nb_dissip)
1123	call VTe(VThallo)
1124	call VTb(VTdissipation)
1125	call start_timer(timer_dissip)
1126	c$OMP END MASTER
1127	c$OMP BARRIER
1128
1129	call covcont_p(llm,ucov,vcov,ucont,vcont)
1130	call enercin_p(vcov,ucov,vcont,ucont,ecin0)
1131
1132	c dissipation
1133
1134	! CALL FTRACE_REGION_BEGIN("dissip")
1135	CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
1136	! CALL FTRACE_REGION_END("dissip")
1137
1138	ijb=ij_begin
1139	ije=ij_end
1140	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1141	DO l=1,llm
1142	ucov(ijb:ije,l)=ucov(ijb:ije,l)+dudis(ijb:ije,l)
1143	ENDDO
1144	c$OMP END DO NOWAIT
1145	if (pole_sud) ije=ije-iip1
1146	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1147	DO l=1,llm
1148	vcov(ijb:ije,l)=vcov(ijb:ije,l)+dvdis(ijb:ije,l)
1149	ENDDO
1150	c$OMP END DO NOWAIT
1151
1152	c teta=teta+dtetadis
1153
1154
1155	c------------------------------------------------------------------------
1156	if (dissip_conservative) then
1157	C On rajoute la tendance due a la transform. Ec -> E therm. cree
1158	C lors de la dissipation
1159	c$OMP BARRIER
1160	c$OMP MASTER
1161	call suspend_timer(timer_dissip)
1162	call VTb(VThallo)
1163	c$OMP END MASTER
1164	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Request_Dissip)
1165	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Dissip)
1166	call SendRequest(Request_Dissip)
1167	c$OMP BARRIER
1168	call WaitRequest(Request_Dissip)
1169	c$OMP MASTER
1170	call VTe(VThallo)
1171	call resume_timer(timer_dissip)
1172	c$OMP END MASTER
1173	c$OMP BARRIER
1174	call covcont_p(llm,ucov,vcov,ucont,vcont)
1175	call enercin_p(vcov,ucov,vcont,ucont,ecin)
1176
1177	ijb=ij_begin
1178	ije=ij_end
1179	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1180	do l=1,llm
1181	do ij=ijb,ije
1182	dtetaecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
1183	dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
1184	enddo
1185	enddo
1186	c$OMP END DO NOWAIT
1187	endif ! of if (dissip_conservative)
1188
1189	ijb=ij_begin
1190	ije=ij_end
1191	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1192	do l=1,llm
1193	do ij=ijb,ije
1194	teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
1195	enddo
1196	enddo
1197	c$OMP END DO NOWAIT
1198	c------------------------------------------------------------------------
1199
1200
1201	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
1202	c ... Calcul de la valeur moyenne, unique de h aux poles .....
1203	c
1204
1205	ijb=ij_begin
1206	ije=ij_end
1207
1208	if (pole_nord) then
1209	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1210	DO l = 1, llm
1211	DO ij = 1,iim
1212	tppn(ij) = aire( ij ) * teta( ij ,l)
1213	ENDDO
1214	tpn = SSUM(iim,tppn,1)/apoln
1215
1216	DO ij = 1, iip1
1217	teta( ij ,l) = tpn
1218	ENDDO
1219	ENDDO
1220	c$OMP END DO NOWAIT
1221
1222	if (1 == 0) then
1223	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
1224	!!! 2) should probably not be here anyway
1225	!!! but are kept for those who would want to revert to previous behaviour
1226	c$OMP MASTER
1227	DO ij = 1,iim
1228	tppn(ij) = aire( ij ) * ps ( ij )
1229	ENDDO
1230	tpn = SSUM(iim,tppn,1)/apoln
1231
1232	DO ij = 1, iip1
1233	ps( ij ) = tpn
1234	ENDDO
1235	c$OMP END MASTER
1236	endif ! of if (1 == 0)
1237	endif ! of of (pole_nord)
1238
1239	if (pole_sud) then
1240	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1241	DO l = 1, llm
1242	DO ij = 1,iim
1243	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
1244	ENDDO
1245	tps = SSUM(iim,tpps,1)/apols
1246
1247	DO ij = 1, iip1
1248	teta(ij+ip1jm,l) = tps
1249	ENDDO
1250	ENDDO
1251	c$OMP END DO NOWAIT
1252
1253	if (1 == 0) then
1254	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
1255	!!! 2) should probably not be here anyway
1256	!!! but are kept for those who would want to revert to previous behaviour
1257	c$OMP MASTER
1258	DO ij = 1,iim
1259	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
1260	ENDDO
1261	tps = SSUM(iim,tpps,1)/apols
1262
1263	DO ij = 1, iip1
1264	ps(ij+ip1jm) = tps
1265	ENDDO
1266	c$OMP END MASTER
1267	endif ! of if (1 == 0)
1268	endif ! of if (pole_sud)
1269
1270
1271	c$OMP BARRIER
1272	c$OMP MASTER
1273	call VTe(VTdissipation)
1274
1275	call stop_timer(timer_dissip)
1276
1277	call VTb(VThallo)
1278	c$OMP END MASTER
1279	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
1280	* jj_Nb_caldyn,Request_dissip)
1281
1282	call Register_SwapField(vcov,vcov,ip1jm,llm,
1283	* jj_Nb_caldyn,Request_dissip)
1284
1285	call Register_SwapField(teta,teta,ip1jmp1,llm,
1286	* jj_Nb_caldyn,Request_dissip)
1287
1288	call Register_SwapField(p,p,ip1jmp1,llmp1,
1289	* jj_Nb_caldyn,Request_dissip)
1290
1291	call Register_SwapField(pk,pk,ip1jmp1,llm,
1292	* jj_Nb_caldyn,Request_dissip)
1293
1294	call SendRequest(Request_dissip)
1295	c$OMP BARRIER
1296	call WaitRequest(Request_dissip)
1297
1298	c$OMP BARRIER
1299	c$OMP MASTER
1300	call SetDistrib(jj_Nb_caldyn)
1301	call VTe(VThallo)
1302	call resume_timer(timer_caldyn)
1303	c print *,'fin dissipation'
1304	c$OMP END MASTER
1305	c$OMP BARRIER
1306	END IF ! of IF(apdiss)
1307
1308	cc$OMP END PARALLEL
1309
1310	c ajout debug
1311	c IF( lafin ) then
1312	c abort_message = 'Simulation finished'
1313	c call abort_gcm(modname,abort_message,0)
1314	c ENDIF
1315
1316	c ********************************************************************
1317	c ********************************************************************
1318	c .... fin de l'integration dynamique et physique pour le pas itau ..
1319	c ********************************************************************
1320	c ********************************************************************
1321
1322	c preparation du pas d'integration suivant ......
1323	cym call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1324	cym call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1325	c$OMP MASTER
1326	call stop_timer(timer_caldyn)
1327	c$OMP END MASTER
1328	IF (itau==itaumax) then
1329	c$OMP MASTER
1330	call allgather_timer_average
1331
1332	if (mpi_rank==0) then
1333
1334	print ,'********************************'
1335	print ,'*** TIMER CALDYN ****'
1336	do i=0,mpi_size-1
1337	print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
1338	& ' : temps moyen :',
1339	& timer_average(jj_nb_caldyn(i),timer_caldyn,i)
1340	enddo
1341
1342	print ,'********************************'
1343	print ,'*** TIMER VANLEER ****'
1344	do i=0,mpi_size-1
1345	print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
1346	& ' : temps moyen :',
1347	& timer_average(jj_nb_vanleer(i),timer_vanleer,i)
1348	enddo
1349
1350	print ,'********************************'
1351	print ,'*** TIMER DISSIP ****'
1352	do i=0,mpi_size-1
1353	print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
1354	& ' : temps moyen :',
1355	& timer_average(jj_nb_dissip(i),timer_dissip,i)
1356	enddo
1357
1358	print ,'********************************'
1359	print ,'*** TIMER PHYSIC ****'
1360	do i=0,mpi_size-1
1361	print *,'proc',i,' : Nb Bandes :',jj_nb_physic(i),
1362	& ' : temps moyen :',
1363	& timer_average(jj_nb_physic(i),timer_physic,i)
1364	enddo
1365
1366	endif
1367
1368	print ,'Taille du Buffer MPI (REAL8)',MaxBufferSize
1369	print ,'Taille du Buffer MPI utilise (REAL8)',MaxBufferSize_Used
1370	print *, 'Temps total ecoule sur la parallelisation :',DiffTime()
1371	print *, 'Temps CPU ecoule sur la parallelisation :',DiffCpuTime()
1372	CALL print_filtre_timer
1373	call fin_getparam
1374	call finalize_parallel
1375	c$OMP END MASTER
1376	c$OMP BARRIER
1377	RETURN
1378	ENDIF ! of IF (itau==itaumax)
1379
1380	IF ( .NOT.purmats ) THEN
1381	c ........................................................
1382	c .............. schema matsuno + leapfrog ..............
1383	c ........................................................
1384
1385	IF(forward. OR. leapf) THEN
1386	itau= itau + 1
1387	! iday= day_ini+itau/day_step
1388	! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1389	! IF(time.GT.1.) THEN
1390	! time = time-1.
1391	! iday = iday+1
1392	! ENDIF
1393	ENDIF
1394
1395
1396	IF( itau. EQ. itaufinp1 ) then
1397
1398	if (flag_verif) then
1399	write(79,*) 'ucov',ucov
1400	write(80,*) 'vcov',vcov
1401	write(81,*) 'teta',teta
1402	write(82,*) 'ps',ps
1403	write(83,*) 'q',q
1404	WRITE(85,*) 'q1 = ',q(:,:,1)
1405	WRITE(86,*) 'q3 = ',q(:,:,3)
1406	endif
1407
1408
1409	c$OMP MASTER
1410	call fin_getparam
1411	call finalize_parallel
1412	c$OMP END MASTER
1413	abort_message = 'Simulation finished'
1414	call abort_gcm(modname,abort_message,0)
1415	RETURN
1416	ENDIF
1417	c-----------------------------------------------------------------------
1418	c ecriture du fichier histoire moyenne:
1419	c -------------------------------------
1420
1421	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1422	c$OMP BARRIER
1423	IF(itau.EQ.itaufin) THEN
1424	iav=1
1425	ELSE
1426	iav=0
1427	ENDIF
1428	#ifdef CPP_IOIPSL
1429	IF (ok_dynzon) THEN
1430	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
1431	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1432	ENDIF !ok_dynzon
1433	#endif
1434	IF (ok_dyn_ave) THEN
1435	!$OMP MASTER
1436	#ifdef CPP_IOIPSL
1437	! Ehouarn: Gather fields and make master send to output
1438	call Gather_Field(vcov,ip1jm,llm,0)
1439	call Gather_Field(ucov,ip1jmp1,llm,0)
1440	call Gather_Field(teta,ip1jmp1,llm,0)
1441	call Gather_Field(pk,ip1jmp1,llm,0)
1442	call Gather_Field(phi,ip1jmp1,llm,0)
1443	do iq=1,nqtot
1444	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1445	enddo
1446	call Gather_Field(masse,ip1jmp1,llm,0)
1447	call Gather_Field(ps,ip1jmp1,1,0)
1448	call Gather_Field(phis,ip1jmp1,1,0)
1449	if (mpi_rank==0) then
1450	CALL writedynav(itau,vcov,
1451	& ucov,teta,pk,phi,q,masse,ps,phis)
1452	endif
1453	#endif
1454	!$OMP END MASTER
1455	ENDIF ! of IF (ok_dyn_ave)
1456	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
1457
1458	c-----------------------------------------------------------------------
1459	c ecriture de la bande histoire:
1460	c ------------------------------
1461
1462	IF( MOD(itau,iecri).EQ.0) THEN
1463	! Ehouarn: output only during LF or Backward Matsuno
1464	if (leapf.or.(.not.leapf.and.(.not.forward))) then
1465	c$OMP BARRIER
1466	c$OMP MASTER
1467	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1468
1469	cym unat=0.
1470
1471	ijb=ij_begin
1472	ije=ij_end
1473
1474	if (pole_nord) then
1475	ijb=ij_begin+iip1
1476	unat(1:iip1,:)=0.
1477	endif
1478
1479	if (pole_sud) then
1480	ije=ij_end-iip1
1481	unat(ij_end-iip1+1:ij_end,:)=0.
1482	endif
1483
1484	do l=1,llm
1485	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1486	enddo
1487
1488	ijb=ij_begin
1489	ije=ij_end
1490	if (pole_sud) ije=ij_end-iip1
1491
1492	do l=1,llm
1493	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1494	enddo
1495
1496	#ifdef CPP_IOIPSL
1497	if (ok_dyn_ins) then
1498	! Ehouarn: Gather fields and make master write to output
1499	call Gather_Field(vcov,ip1jm,llm,0)
1500	call Gather_Field(ucov,ip1jmp1,llm,0)
1501	call Gather_Field(teta,ip1jmp1,llm,0)
1502	call Gather_Field(phi,ip1jmp1,llm,0)
1503	do iq=1,nqtot
1504	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1505	enddo
1506	call Gather_Field(masse,ip1jmp1,llm,0)
1507	call Gather_Field(ps,ip1jmp1,1,0)
1508	call Gather_Field(phis,ip1jmp1,1,0)
1509	if (mpi_rank==0) then
1510	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1511	endif
1512	! CALL writehist_p(histid,histvid, itau,vcov,
1513	! & ucov,teta,phi,q,masse,ps,phis)
1514	! or use writefield_p
1515	! call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1516	! call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1517	! call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
1518	! call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
1519	endif ! of if (ok_dyn_ins)
1520	#endif
1521	! For some Grads outputs of fields
1522	if (output_grads_dyn) then
1523	! Ehouarn: hope this works the way I think it does:
1524	call Gather_Field(unat,ip1jmp1,llm,0)
1525	call Gather_Field(vnat,ip1jm,llm,0)
1526	call Gather_Field(teta,ip1jmp1,llm,0)
1527	call Gather_Field(ps,ip1jmp1,1,0)
1528	do iq=1,nqtot
1529	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1530	enddo
1531	if (mpi_rank==0) then
1532	#include "write_grads_dyn.h"
1533	endif
1534	endif ! of if (output_grads_dyn)
1535	c$OMP END MASTER
1536	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
1537	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1538
1539	IF(itau.EQ.itaufin) THEN
1540
1541	c$OMP BARRIER
1542	c$OMP MASTER
1543
1544	! if (planet_type.eq."earth") then
1545	! Write an Earth-format restart file
1546	CALL dynredem1_p("restart.nc",0.0,
1547	& vcov,ucov,teta,q,masse,ps)
1548	! endif ! of if (planet_type.eq."earth")
1549
1550	! CLOSE(99)
1551	c$OMP END MASTER
1552	ENDIF ! of IF (itau.EQ.itaufin)
1553
1554	c-----------------------------------------------------------------------
1555	c gestion de l'integration temporelle:
1556	c ------------------------------------
1557
1558	IF( MOD(itau,iperiod).EQ.0 ) THEN
1559	GO TO 1
1560	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
1561
1562	IF( forward ) THEN
1563	c fin du pas forward et debut du pas backward
1564
1565	forward = .FALSE.
1566	leapf = .FALSE.
1567	GO TO 2
1568
1569	ELSE
1570	c fin du pas backward et debut du premier pas leapfrog
1571
1572	leapf = .TRUE.
1573	dt = 2.*dtvr
1574	GO TO 2
1575	END IF
1576	ELSE
1577
1578	c ...... pas leapfrog .....
1579
1580	leapf = .TRUE.
1581	dt = 2.*dtvr
1582	GO TO 2
1583	END IF ! of IF (MOD(itau,iperiod).EQ.0)
1584	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
1585
1586
1587	ELSE ! of IF (.not.purmats)
1588
1589	c ........................................................
1590	c .............. schema matsuno ...............
1591	c ........................................................
1592	IF( forward ) THEN
1593
1594	itau = itau + 1
1595	! iday = day_ini+itau/day_step
1596	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1597	!
1598	! IF(time.GT.1.) THEN
1599	! time = time-1.
1600	! iday = iday+1
1601	! ENDIF
1602
1603	forward = .FALSE.
1604	IF( itau. EQ. itaufinp1 ) then
1605	c$OMP MASTER
1606	call fin_getparam
1607	call finalize_parallel
1608	c$OMP END MASTER
1609	abort_message = 'Simulation finished'
1610	call abort_gcm(modname,abort_message,0)
1611	RETURN
1612	ENDIF
1613	GO TO 2
1614
1615	ELSE ! of IF(forward) i.e. backward step
1616
1617	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1618	IF(itau.EQ.itaufin) THEN
1619	iav=1
1620	ELSE
1621	iav=0
1622	ENDIF
1623	#ifdef CPP_IOIPSL
1624	IF (ok_dynzon) THEN
1625	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
1626	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1627	END IF !ok_dynzon
1628	#endif
1629	IF (ok_dyn_ave) THEN
1630	!$OMP MASTER
1631	#ifdef CPP_IOIPSL
1632	! Ehouarn: Gather fields and make master send to output
1633	call Gather_Field(vcov,ip1jm,llm,0)
1634	call Gather_Field(ucov,ip1jmp1,llm,0)
1635	call Gather_Field(teta,ip1jmp1,llm,0)
1636	call Gather_Field(pk,ip1jmp1,llm,0)
1637	call Gather_Field(phi,ip1jmp1,llm,0)
1638	do iq=1,nqtot
1639	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1640	enddo
1641	call Gather_Field(masse,ip1jmp1,llm,0)
1642	call Gather_Field(ps,ip1jmp1,1,0)
1643	call Gather_Field(phis,ip1jmp1,1,0)
1644	if (mpi_rank==0) then
1645	CALL writedynav(itau,vcov,
1646	& ucov,teta,pk,phi,q,masse,ps,phis)
1647	endif
1648	#endif
1649	!$OMP END MASTER
1650	ENDIF ! of IF (ok_dyn_ave)
1651
1652	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
1653
1654
1655	IF(MOD(itau,iecri ).EQ.0) THEN
1656	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
1657	c$OMP BARRIER
1658	c$OMP MASTER
1659	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1660
1661	cym unat=0.
1662	ijb=ij_begin
1663	ije=ij_end
1664
1665	if (pole_nord) then
1666	ijb=ij_begin+iip1
1667	unat(1:iip1,:)=0.
1668	endif
1669
1670	if (pole_sud) then
1671	ije=ij_end-iip1
1672	unat(ij_end-iip1+1:ij_end,:)=0.
1673	endif
1674
1675	do l=1,llm
1676	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1677	enddo
1678
1679	ijb=ij_begin
1680	ije=ij_end
1681	if (pole_sud) ije=ij_end-iip1
1682
1683	do l=1,llm
1684	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1685	enddo
1686
1687	#ifdef CPP_IOIPSL
1688	if (ok_dyn_ins) then
1689	! Ehouarn: Gather fields and make master send to output
1690	call Gather_Field(vcov,ip1jm,llm,0)
1691	call Gather_Field(ucov,ip1jmp1,llm,0)
1692	call Gather_Field(teta,ip1jmp1,llm,0)
1693	call Gather_Field(phi,ip1jmp1,llm,0)
1694	do iq=1,nqtot
1695	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1696	enddo
1697	call Gather_Field(masse,ip1jmp1,llm,0)
1698	call Gather_Field(ps,ip1jmp1,1,0)
1699	call Gather_Field(phis,ip1jmp1,1,0)
1700	if (mpi_rank==0) then
1701	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1702	endif
1703	! CALL writehist_p(histid, histvid, itau,vcov ,
1704	! & ucov,teta,phi,q,masse,ps,phis)
1705	endif ! of if (ok_dyn_ins)
1706	#endif
1707	! For some Grads output (but does it work?)
1708	if (output_grads_dyn) then
1709	call Gather_Field(unat,ip1jmp1,llm,0)
1710	call Gather_Field(vnat,ip1jm,llm,0)
1711	call Gather_Field(teta,ip1jmp1,llm,0)
1712	call Gather_Field(ps,ip1jmp1,1,0)
1713	do iq=1,nqtot
1714	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1715	enddo
1716	c
1717	if (mpi_rank==0) then
1718	#include "write_grads_dyn.h"
1719	endif
1720	endif ! of if (output_grads_dyn)
1721
1722	c$OMP END MASTER
1723	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1724
1725	IF(itau.EQ.itaufin) THEN
1726	! if (planet_type.eq."earth") then
1727	c$OMP MASTER
1728	CALL dynredem1_p("restart.nc",0.0,
1729	. vcov,ucov,teta,q,masse,ps)
1730	c$OMP END MASTER
1731	! endif ! of if (planet_type.eq."earth")
1732	ENDIF ! of IF(itau.EQ.itaufin)
1733
1734	forward = .TRUE.
1735	GO TO 1
1736
1737	ENDIF ! of IF (forward)
1738
1739	END IF ! of IF(.not.purmats)
1740	c$OMP MASTER
1741	call fin_getparam
1742	call finalize_parallel
1743	c$OMP END MASTER
1744	RETURN
1745	END

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

Download in other formats:

Original Format