Context Navigation

leapfrog_p.F @ 2124

Last change on this file since 2124 was 2083, checked in by Ehouarn Millour, 10 years ago

Minor fix in dyn3dpar/leapfrog_p.F , should call geopot_p and not geopot
Added a sanity check in iniacademic
Added flag "resetvarc" to trigger a reset of initial values in sortvarc
Removed "sortvarc0" since the job can now be done with "resetvarc" and having set flag resertvarc to true.

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

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

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

Download in other formats:

Original Format