Context Navigation

leapfrog_p.F @ 2404

Last change on this file since 2404 was 2375, checked in by Ehouarn Millour, 9 years ago
Fix in the computation of the date; the convention is that it corresponds to the time at the end of the current dynamics or physics step (except when in backward Matsuno step where it remains unchanged as it was correctly updated during the forward part of the step). Note that the relationship between itau and date is a bit tricky as itau is incremented between Matsuno forward and backward steps (and from a leapfrog step to the next) but unchanged from Matsuno bacward step to leapfrog step. Because of change in date when calling physics, bench results will change with this revision. EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 56.1 KB

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

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

Download in other formats:

Original Format