Context Navigation

leapfrog_p.F @ 1541

Last change on this file since 1541 was 1520, checked in by Ehouarn Millour, 13 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 54.8 KB

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

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

Download in other formats:

Original Format