Context Navigation

leapfrog_p.F @ 1985

Last change on this file since 1985 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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: 54.5 KB

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

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

Download in other formats:

Original Format