Context Navigation

leapfrog_p.F @ 2014

Last change on this file since 2014 was 1959, checked in by emillour, 7 years ago
Common (parallel) dynamics: Add some initializations of global arrays. This is not a bug fix, strictly speaking, as results are unchanged. But in some cases computations extend into halos and these cases are identified as errors with recent versions of ifort (2018) with the -init=arrays,snan debugging option. EM
File size: 65.8 KB

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

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

Download in other formats:

Original Format