Context Navigation

leapfrog_p.F @ 5440

Last change on this file since 5440 was 2375, checked in by Ehouarn Millour, 9 years ago
Fix in the computation of the date; the convention is that it corresponds to the time at the end of the current dynamics or physics step (except when in backward Matsuno step where it remains unchanged as it was correctly updated during the forward part of the step). Note that the relationship between itau and date is a bit tricky as itau is incremented between Matsuno forward and backward steps (and from a leapfrog step to the next) but unchanged from Matsuno bacward step to leapfrog step. Because of change in date when calling physics, bench results will change with this revision. EM
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: 56.1 KB

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

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

Download in other formats:

Original Format