Context Navigation

leapfrog_p.F @ 1678

Last change on this file since 1678 was 1676, checked in by aslmd, 12 years ago

Repare le probleme "un jour sans fin" des versions precedentes pour le modele generique. Pas d'impact sur modele terrestre.

NB: changements egalement faits sur dyn3dmem en prevision de son utilisation

Fix the bug "Groundhog Day" which was a problem for previous versions using generic physics. No impact whatsoever on terrestrial model.

NB: changes also done on dyn3dmem which is planned to be used soon

Author : AS equipe planeto

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

File size: 55.7 KB

Rev	Line
[764]	1	!
[1279]	2	! $Id: leapfrog_p.F 1676 2012-11-08 20:03:30Z emillour $
[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
	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 ------------------------------
[1000]	906	IF (ok_strato) THEN
[1279]	907	CALL top_bound_p( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
[1000]	908	ENDIF
	909
[1146]	910	CALL addfi_p( dtphys, leapf, forward ,
[630]	911	$ ucov, vcov, teta , q ,ps ,
	912	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	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	! -------------------------------------------------------
[630]	996	ijb=ij_begin
	997	ije=ij_end
[1403]	998	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	999	do l=1,llm
[1454]	1000	teta(ijb:ije,l)=teta(ijb:ije,l)-dtvr*
	1001	& (teta(ijb:ije,l)-tetarappel(ijb:ije,l))*
	1002	& (knewt_g+knewt_t(l)*clat4(ijb:ije))
	1003	enddo ! of do l=1,llm
[1403]	1004	!$OMP END DO
[630]	1005
[1505]	1006	!$OMP MASTER
	1007	if (planet_type.eq."giant") then
	1008	! add an intrinsic heat flux at the base of the atmosphere
	1009	teta(ijb:ije,1) = teta(ijb:ije,1)
	1010	& + dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1)
	1011	endif
	1012	!$OMP END MASTER
	1013	!$OMP BARRIER
	1014
[630]	1015	call Register_Hallo(ucov,ip1jmp1,llm,0,1,1,0,Request_Physic)
	1016	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Physic)
	1017	call SendRequest(Request_Physic)
[1279]	1018	c$OMP BARRIER
[630]	1019	call WaitRequest(Request_Physic)
[1403]	1020	c$OMP BARRIER
[1454]	1021	call friction_p(ucov,vcov,dtvr)
[1403]	1022	!$OMP BARRIER
[1454]	1023
	1024	! Sponge layer (if any)
	1025	IF (ok_strato) THEN
	1026	! set dufi,dvfi,... to zero
	1027	ijb=ij_begin
	1028	ije=ij_end
	1029	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1030	do l=1,llm
	1031	dufi(ijb:ije,l)=0
	1032	dtetafi(ijb:ije,l)=0
	1033	dqfi(ijb:ije,l,1:nqtot)=0
	1034	enddo
	1035	!$OMP END DO
[1492]	1036	!$OMP MASTER
[1454]	1037	dpfi(ijb:ije)=0
[1492]	1038	!$OMP END MASTER
[1454]	1039	ijb=ij_begin
	1040	ije=ij_end
	1041	if (pole_sud) ije=ije-iip1
	1042	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1043	do l=1,llm
	1044	dvfi(ijb:ije,l)=0
	1045	enddo
	1046	!$OMP END DO
	1047
	1048	CALL top_bound_p(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
	1049	CALL addfi_p( dtvr, leapf, forward ,
	1050	$ ucov, vcov, teta , q ,ps ,
	1051	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	1052	!$OMP BARRIER
	1053	ENDIF ! of IF (ok_strato)
[1146]	1054	ENDIF ! of IF(iflag_phys.EQ.2)
[630]	1055
	1056
[985]	1057	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
[774]	1058	c$OMP BARRIER
[1625]	1059	if (pressure_exner) then
[1520]	1060	CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
[1625]	1061	else
[1520]	1062	CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
	1063	endif
[774]	1064	c$OMP BARRIER
[630]	1065
[774]	1066	cc$OMP END PARALLEL
[630]	1067
	1068	c-----------------------------------------------------------------------
	1069	c dissipation horizontale et verticale des petites echelles:
	1070	c ----------------------------------------------------------
	1071
	1072	IF(apdiss) THEN
[774]	1073	cc$OMP PARALLEL DEFAULT(SHARED)
	1074	cc$OMP+ PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
[764]	1075	c$OMP MASTER
[630]	1076	call suspend_timer(timer_caldyn)
[949]	1077
	1078	c print*,'Entree dans la dissipation : Iteration No ',true_itau
[630]	1079	c calcul de l'energie cinetique avant dissipation
[949]	1080	c print *,'Passage dans la dissipation'
[630]	1081
	1082	call VTb(VThallo)
[764]	1083	c$OMP END MASTER
[630]	1084
[764]	1085	c$OMP BARRIER
[985]	1086
[630]	1087	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
	1088	* jj_Nb_dissip,1,1,Request_dissip)
	1089
	1090	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
	1091	* jj_Nb_dissip,1,1,Request_dissip)
	1092
	1093	call Register_SwapField(teta,teta,ip1jmp1,llm,
	1094	* jj_Nb_dissip,Request_dissip)
	1095
	1096	call Register_SwapField(p,p,ip1jmp1,llmp1,
	1097	* jj_Nb_dissip,Request_dissip)
	1098
	1099	call Register_SwapField(pk,pk,ip1jmp1,llm,
	1100	* jj_Nb_dissip,Request_dissip)
	1101
	1102	call SendRequest(Request_dissip)
[985]	1103	c$OMP BARRIER
[630]	1104	call WaitRequest(Request_dissip)
[985]	1105
	1106	c$OMP BARRIER
	1107	c$OMP MASTER
[630]	1108	call SetDistrib(jj_Nb_dissip)
[949]	1109	call VTe(VThallo)
[630]	1110	call VTb(VTdissipation)
[949]	1111	call start_timer(timer_dissip)
[764]	1112	c$OMP END MASTER
	1113	c$OMP BARRIER
	1114
[630]	1115	call covcont_p(llm,ucov,vcov,ucont,vcont)
	1116	call enercin_p(vcov,ucov,vcont,ucont,ecin0)
	1117
	1118	c dissipation
	1119
[985]	1120	! CALL FTRACE_REGION_BEGIN("dissip")
[630]	1121	CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
[985]	1122	! CALL FTRACE_REGION_END("dissip")
[764]	1123
[949]	1124	ijb=ij_begin
	1125	ije=ij_end
	1126	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1127	DO l=1,llm
	1128	ucov(ijb:ije,l)=ucov(ijb:ije,l)+dudis(ijb:ije,l)
[764]	1129	ENDDO
[949]	1130	c$OMP END DO NOWAIT
	1131	if (pole_sud) ije=ije-iip1
	1132	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1133	DO l=1,llm
	1134	vcov(ijb:ije,l)=vcov(ijb:ije,l)+dvdis(ijb:ije,l)
[764]	1135	ENDDO
[949]	1136	c$OMP END DO NOWAIT
[764]	1137
[630]	1138	c teta=teta+dtetadis
	1139
	1140
	1141	c------------------------------------------------------------------------
	1142	if (dissip_conservative) then
	1143	C On rajoute la tendance due a la transform. Ec -> E therm. cree
	1144	C lors de la dissipation
[764]	1145	c$OMP BARRIER
	1146	c$OMP MASTER
[630]	1147	call suspend_timer(timer_dissip)
[949]	1148	call VTb(VThallo)
[985]	1149	c$OMP END MASTER
[630]	1150	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Request_Dissip)
	1151	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Dissip)
	1152	call SendRequest(Request_Dissip)
[985]	1153	c$OMP BARRIER
[630]	1154	call WaitRequest(Request_Dissip)
[985]	1155	c$OMP MASTER
[630]	1156	call VTe(VThallo)
	1157	call resume_timer(timer_dissip)
[764]	1158	c$OMP END MASTER
[949]	1159	c$OMP BARRIER
[630]	1160	call covcont_p(llm,ucov,vcov,ucont,vcont)
	1161	call enercin_p(vcov,ucov,vcont,ucont,ecin)
	1162
[949]	1163	ijb=ij_begin
	1164	ije=ij_end
	1165	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1166	do l=1,llm
	1167	do ij=ijb,ije
	1168	dtetaecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
	1169	dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
[630]	1170	enddo
[949]	1171	enddo
	1172	c$OMP END DO NOWAIT
[1403]	1173	endif ! of if (dissip_conservative)
[630]	1174
	1175	ijb=ij_begin
	1176	ije=ij_end
[949]	1177	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1178	do l=1,llm
	1179	do ij=ijb,ije
[630]	1180	teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
[949]	1181	enddo
	1182	enddo
	1183	c$OMP END DO NOWAIT
[630]	1184	c------------------------------------------------------------------------
	1185
	1186
	1187	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
	1188	c ... Calcul de la valeur moyenne, unique de h aux poles .....
	1189	c
	1190
	1191	ijb=ij_begin
	1192	ije=ij_end
	1193
	1194	if (pole_nord) then
[764]	1195	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	1196	DO l = 1, llm
	1197	DO ij = 1,iim
	1198	tppn(ij) = aire( ij ) * teta( ij ,l)
	1199	ENDDO
	1200	tpn = SSUM(iim,tppn,1)/apoln
	1201
	1202	DO ij = 1, iip1
	1203	teta( ij ,l) = tpn
	1204	ENDDO
	1205	ENDDO
[764]	1206	c$OMP END DO NOWAIT
	1207
[1614]	1208	if (1 == 0) then
	1209	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
	1210	!!! 2) should probably not be here anyway
	1211	!!! but are kept for those who would want to revert to previous behaviour
[764]	1212	c$OMP MASTER
[630]	1213	DO ij = 1,iim
	1214	tppn(ij) = aire( ij ) * ps ( ij )
	1215	ENDDO
	1216	tpn = SSUM(iim,tppn,1)/apoln
	1217
	1218	DO ij = 1, iip1
	1219	ps( ij ) = tpn
	1220	ENDDO
[764]	1221	c$OMP END MASTER
[1614]	1222	endif ! of if (1 == 0)
	1223	endif ! of of (pole_nord)
[630]	1224
	1225	if (pole_sud) then
[764]	1226	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	1227	DO l = 1, llm
	1228	DO ij = 1,iim
	1229	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	1230	ENDDO
	1231	tps = SSUM(iim,tpps,1)/apols
	1232
	1233	DO ij = 1, iip1
	1234	teta(ij+ip1jm,l) = tps
	1235	ENDDO
	1236	ENDDO
[764]	1237	c$OMP END DO NOWAIT
	1238
[1614]	1239	if (1 == 0) then
	1240	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
	1241	!!! 2) should probably not be here anyway
	1242	!!! but are kept for those who would want to revert to previous behaviour
[764]	1243	c$OMP MASTER
[630]	1244	DO ij = 1,iim
	1245	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
	1246	ENDDO
	1247	tps = SSUM(iim,tpps,1)/apols
	1248
	1249	DO ij = 1, iip1
	1250	ps(ij+ip1jm) = tps
	1251	ENDDO
[764]	1252	c$OMP END MASTER
[1614]	1253	endif ! of if (1 == 0)
	1254	endif ! of if (pole_sud)
[630]	1255
[764]	1256
	1257	c$OMP BARRIER
	1258	c$OMP MASTER
[630]	1259	call VTe(VTdissipation)
	1260
	1261	call stop_timer(timer_dissip)
[949]	1262
[630]	1263	call VTb(VThallo)
[985]	1264	c$OMP END MASTER
[630]	1265	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
	1266	* jj_Nb_caldyn,Request_dissip)
	1267
	1268	call Register_SwapField(vcov,vcov,ip1jm,llm,
	1269	* jj_Nb_caldyn,Request_dissip)
	1270
	1271	call Register_SwapField(teta,teta,ip1jmp1,llm,
	1272	* jj_Nb_caldyn,Request_dissip)
	1273
	1274	call Register_SwapField(p,p,ip1jmp1,llmp1,
	1275	* jj_Nb_caldyn,Request_dissip)
	1276
	1277	call Register_SwapField(pk,pk,ip1jmp1,llm,
	1278	* jj_Nb_caldyn,Request_dissip)
	1279
	1280	call SendRequest(Request_dissip)
[985]	1281	c$OMP BARRIER
[630]	1282	call WaitRequest(Request_dissip)
[985]	1283
	1284	c$OMP BARRIER
	1285	c$OMP MASTER
[630]	1286	call SetDistrib(jj_Nb_caldyn)
	1287	call VTe(VThallo)
	1288	call resume_timer(timer_caldyn)
[961]	1289	c print *,'fin dissipation'
[764]	1290	c$OMP END MASTER
[985]	1291	c$OMP BARRIER
[1403]	1292	END IF ! of IF(apdiss)
[630]	1293
[774]	1294	cc$OMP END PARALLEL
	1295
[630]	1296	c ajout debug
	1297	c IF( lafin ) then
	1298	c abort_message = 'Simulation finished'
	1299	c call abort_gcm(modname,abort_message,0)
	1300	c ENDIF
	1301
	1302	c ********************************************************************
	1303	c ********************************************************************
	1304	c .... fin de l'integration dynamique et physique pour le pas itau ..
	1305	c ********************************************************************
	1306	c ********************************************************************
	1307
	1308	c preparation du pas d'integration suivant ......
	1309	cym call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	1310	cym call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
[774]	1311	c$OMP MASTER
[630]	1312	call stop_timer(timer_caldyn)
[774]	1313	c$OMP END MASTER
[630]	1314	IF (itau==itaumax) then
[774]	1315	c$OMP MASTER
[630]	1316	call allgather_timer_average
	1317
	1318	if (mpi_rank==0) then
	1319
	1320	print ,'********************************'
[949]	1321	print ,'*** TIMER CALDYN ****'
	1322	do i=0,mpi_size-1
	1323	print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
[630]	1324	& ' : temps moyen :',
	1325	& timer_average(jj_nb_caldyn(i),timer_caldyn,i)
	1326	enddo
	1327
	1328	print ,'********************************'
[949]	1329	print ,'*** TIMER VANLEER ****'
	1330	do i=0,mpi_size-1
	1331	print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
[630]	1332	& ' : temps moyen :',
	1333	& timer_average(jj_nb_vanleer(i),timer_vanleer,i)
	1334	enddo
	1335
	1336	print ,'********************************'
[949]	1337	print ,'*** TIMER DISSIP ****'
	1338	do i=0,mpi_size-1
	1339	print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
[630]	1340	& ' : temps moyen :',
	1341	& timer_average(jj_nb_dissip(i),timer_dissip,i)
	1342	enddo
[949]	1343
	1344	print ,'********************************'
	1345	print ,'*** TIMER PHYSIC ****'
	1346	do i=0,mpi_size-1
	1347	print *,'proc',i,' : Nb Bandes :',jj_nb_physic(i),
[630]	1348	& ' : temps moyen :',
	1349	& timer_average(jj_nb_physic(i),timer_physic,i)
	1350	enddo
[949]	1351
[630]	1352	endif
	1353
	1354	print ,'Taille du Buffer MPI (REAL8)',MaxBufferSize
	1355	print ,'Taille du Buffer MPI utilise (REAL8)',MaxBufferSize_Used
	1356	print *, 'Temps total ecoule sur la parallelisation :',DiffTime()
	1357	print *, 'Temps CPU ecoule sur la parallelisation :',DiffCpuTime()
[985]	1358	CALL print_filtre_timer
[1279]	1359	call fin_getparam
[630]	1360	call finalize_parallel
[774]	1361	c$OMP END MASTER
[985]	1362	c$OMP BARRIER
[949]	1363	RETURN
[630]	1364	ENDIF
	1365
	1366	IF ( .NOT.purmats ) THEN
	1367	c ........................................................
	1368	c .............. schema matsuno + leapfrog ..............
	1369	c ........................................................
	1370
	1371	IF(forward. OR. leapf) THEN
	1372	itau= itau + 1
[1279]	1373	! iday= day_ini+itau/day_step
[1403]	1374	! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
[1279]	1375	! IF(time.GT.1.) THEN
	1376	! time = time-1.
	1377	! iday = iday+1
	1378	! ENDIF
[630]	1379	ENDIF
	1380
	1381
[1286]	1382	IF( itau. EQ. itaufinp1 ) then
[630]	1383
[1286]	1384	if (flag_verif) then
	1385	write(79,*) 'ucov',ucov
	1386	write(80,*) 'vcov',vcov
	1387	write(81,*) 'teta',teta
	1388	write(82,*) 'ps',ps
	1389	write(83,*) 'q',q
	1390	WRITE(85,*) 'q1 = ',q(:,:,1)
	1391	WRITE(86,*) 'q3 = ',q(:,:,3)
	1392	endif
	1393
	1394
[774]	1395	c$OMP MASTER
[1279]	1396	call fin_getparam
[764]	1397	call finalize_parallel
[774]	1398	c$OMP END MASTER
[630]	1399	abort_message = 'Simulation finished'
	1400	call abort_gcm(modname,abort_message,0)
[985]	1401	RETURN
[630]	1402	ENDIF
	1403	c-----------------------------------------------------------------------
	1404	c ecriture du fichier histoire moyenne:
	1405	c -------------------------------------
	1406
	1407	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
[774]	1408	c$OMP BARRIER
[630]	1409	IF(itau.EQ.itaufin) THEN
	1410	iav=1
	1411	ELSE
	1412	iav=0
	1413	ENDIF
	1414	#ifdef CPP_IOIPSL
[1146]	1415	IF (ok_dynzon) THEN
[774]	1416	call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
[985]	1417	call SendRequest(TestRequest)
	1418	c$OMP BARRIER
[630]	1419	call WaitRequest(TestRequest)
[1000]	1420	c$OMP BARRIER
[985]	1421	c$OMP MASTER
[1279]	1422	! CALL writedynav_p(histaveid, itau,vcov ,
	1423	! , ucov,teta,pk,phi,q,masse,ps,phis)
[1146]	1424
	1425	c ATTENTION!!! bilan_dyn_p ne marche probablement pas avec OpenMP
	1426	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
	1427	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[985]	1428	c$OMP END MASTER
[1146]	1429	ENDIF !ok_dynzon
[630]	1430	#endif
[1403]	1431	IF (ok_dyn_ave) THEN
	1432	!$OMP MASTER
	1433	#ifdef CPP_IOIPSL
	1434	! Ehouarn: Gather fields and make master send to output
	1435	call Gather_Field(vcov,ip1jm,llm,0)
	1436	call Gather_Field(ucov,ip1jmp1,llm,0)
	1437	call Gather_Field(teta,ip1jmp1,llm,0)
	1438	call Gather_Field(pk,ip1jmp1,llm,0)
	1439	call Gather_Field(phi,ip1jmp1,llm,0)
	1440	do iq=1,nqtot
	1441	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
	1442	enddo
	1443	call Gather_Field(masse,ip1jmp1,llm,0)
	1444	call Gather_Field(ps,ip1jmp1,1,0)
	1445	call Gather_Field(phis,ip1jmp1,1,0)
	1446	if (mpi_rank==0) then
	1447	CALL writedynav(itau,vcov,
	1448	& ucov,teta,pk,phi,q,masse,ps,phis)
	1449	endif
	1450	#endif
	1451	!$OMP END MASTER
	1452	ENDIF ! of IF (ok_dyn_ave)
	1453	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
[630]	1454
	1455	c-----------------------------------------------------------------------
	1456	c ecriture de la bande histoire:
	1457	c ------------------------------
	1458
[1403]	1459	IF( MOD(itau,iecri).EQ.0) THEN
	1460	! Ehouarn: output only during LF or Backward Matsuno
	1461	if (leapf.or.(.not.leapf.and.(.not.forward))) then
[774]	1462	c$OMP BARRIER
	1463	c$OMP MASTER
[1146]	1464	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
[774]	1465
[949]	1466	cym unat=0.
	1467
[1146]	1468	ijb=ij_begin
	1469	ije=ij_end
[949]	1470
[1146]	1471	if (pole_nord) then
	1472	ijb=ij_begin+iip1
	1473	unat(1:iip1,:)=0.
	1474	endif
[949]	1475
[1146]	1476	if (pole_sud) then
	1477	ije=ij_end-iip1
	1478	unat(ij_end-iip1+1:ij_end,:)=0.
	1479	endif
[949]	1480
[1146]	1481	do l=1,llm
	1482	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
	1483	enddo
[630]	1484
[1146]	1485	ijb=ij_begin
	1486	ije=ij_end
	1487	if (pole_sud) ije=ij_end-iip1
[949]	1488
[1146]	1489	do l=1,llm
	1490	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
	1491	enddo
[949]	1492
[630]	1493	#ifdef CPP_IOIPSL
[1403]	1494	if (ok_dyn_ins) then
	1495	! Ehouarn: Gather fields and make master write to output
	1496	call Gather_Field(vcov,ip1jm,llm,0)
	1497	call Gather_Field(ucov,ip1jmp1,llm,0)
	1498	call Gather_Field(teta,ip1jmp1,llm,0)
	1499	call Gather_Field(phi,ip1jmp1,llm,0)
	1500	do iq=1,nqtot
	1501	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
	1502	enddo
	1503	call Gather_Field(masse,ip1jmp1,llm,0)
	1504	call Gather_Field(ps,ip1jmp1,1,0)
	1505	call Gather_Field(phis,ip1jmp1,1,0)
	1506	if (mpi_rank==0) then
	1507	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	1508	endif
[1279]	1509	! CALL writehist_p(histid,histvid, itau,vcov,
	1510	! & ucov,teta,phi,q,masse,ps,phis)
[1403]	1511	! or use writefield_p
	1512	! call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	1513	! call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
	1514	! call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
	1515	! call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
	1516	endif ! of if (ok_dyn_ins)
[630]	1517	#endif
[1146]	1518	! For some Grads outputs of fields
	1519	if (output_grads_dyn) then
	1520	! Ehouarn: hope this works the way I think it does:
	1521	call Gather_Field(unat,ip1jmp1,llm,0)
	1522	call Gather_Field(vnat,ip1jm,llm,0)
	1523	call Gather_Field(teta,ip1jmp1,llm,0)
	1524	call Gather_Field(ps,ip1jmp1,1,0)
	1525	do iq=1,nqtot
	1526	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
	1527	enddo
	1528	if (mpi_rank==0) then
	1529	#include "write_grads_dyn.h"
	1530	endif
	1531	endif ! of if (output_grads_dyn)
[774]	1532	c$OMP END MASTER
[1403]	1533	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
[1146]	1534	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
[630]	1535
	1536	IF(itau.EQ.itaufin) THEN
	1537
[774]	1538	c$OMP BARRIER
	1539	c$OMP MASTER
[630]	1540
[1454]	1541	! if (planet_type.eq."earth") then
[1146]	1542	! Write an Earth-format restart file
	1543	CALL dynredem1_p("restart.nc",0.0,
	1544	& vcov,ucov,teta,q,masse,ps)
[1454]	1545	! endif ! of if (planet_type.eq."earth")
[630]	1546
[1279]	1547	! CLOSE(99)
[774]	1548	c$OMP END MASTER
[1146]	1549	ENDIF ! of IF (itau.EQ.itaufin)
[630]	1550
	1551	c-----------------------------------------------------------------------
	1552	c gestion de l'integration temporelle:
	1553	c ------------------------------------
	1554
	1555	IF( MOD(itau,iperiod).EQ.0 ) THEN
	1556	GO TO 1
	1557	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	1558
	1559	IF( forward ) THEN
	1560	c fin du pas forward et debut du pas backward
	1561
	1562	forward = .FALSE.
	1563	leapf = .FALSE.
	1564	GO TO 2
	1565
	1566	ELSE
	1567	c fin du pas backward et debut du premier pas leapfrog
	1568
	1569	leapf = .TRUE.
	1570	dt = 2.*dtvr
	1571	GO TO 2
	1572	END IF
	1573	ELSE
	1574
	1575	c ...... pas leapfrog .....
	1576
	1577	leapf = .TRUE.
	1578	dt = 2.*dtvr
	1579	GO TO 2
[1146]	1580	END IF ! of IF (MOD(itau,iperiod).EQ.0)
	1581	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
[630]	1582
	1583
[1146]	1584	ELSE ! of IF (.not.purmats)
	1585
[630]	1586	c ........................................................
	1587	c .............. schema matsuno ...............
	1588	c ........................................................
	1589	IF( forward ) THEN
	1590
	1591	itau = itau + 1
[1279]	1592	! iday = day_ini+itau/day_step
[1403]	1593	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
[1279]	1594	!
	1595	! IF(time.GT.1.) THEN
	1596	! time = time-1.
	1597	! iday = iday+1
	1598	! ENDIF
[630]	1599
	1600	forward = .FALSE.
	1601	IF( itau. EQ. itaufinp1 ) then
[774]	1602	c$OMP MASTER
[1279]	1603	call fin_getparam
[949]	1604	call finalize_parallel
[774]	1605	c$OMP END MASTER
[630]	1606	abort_message = 'Simulation finished'
	1607	call abort_gcm(modname,abort_message,0)
[985]	1608	RETURN
[630]	1609	ENDIF
	1610	GO TO 2
	1611
[1403]	1612	ELSE ! of IF(forward) i.e. backward step
[630]	1613
[1146]	1614	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
[630]	1615	IF(itau.EQ.itaufin) THEN
	1616	iav=1
	1617	ELSE
	1618	iav=0
	1619	ENDIF
	1620	#ifdef CPP_IOIPSL
[1146]	1621	IF (ok_dynzon) THEN
[774]	1622	c$OMP BARRIER
[1146]	1623	call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
	1624	call SendRequest(TestRequest)
[985]	1625	c$OMP BARRIER
[1146]	1626	call WaitRequest(TestRequest)
[1000]	1627	c$OMP BARRIER
[985]	1628	c$OMP MASTER
[1279]	1629	! CALL writedynav_p(histaveid, itau,vcov ,
	1630	! , ucov,teta,pk,phi,q,masse,ps,phis)
[1146]	1631	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
[985]	1632	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[774]	1633	c$OMP END MASTER
[1146]	1634	END IF !ok_dynzon
[630]	1635	#endif
[1403]	1636	IF (ok_dyn_ave) THEN
	1637	!$OMP MASTER
	1638	#ifdef CPP_IOIPSL
	1639	! Ehouarn: Gather fields and make master send to output
	1640	call Gather_Field(vcov,ip1jm,llm,0)
	1641	call Gather_Field(ucov,ip1jmp1,llm,0)
	1642	call Gather_Field(teta,ip1jmp1,llm,0)
	1643	call Gather_Field(pk,ip1jmp1,llm,0)
	1644	call Gather_Field(phi,ip1jmp1,llm,0)
	1645	do iq=1,nqtot
	1646	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
	1647	enddo
	1648	call Gather_Field(masse,ip1jmp1,llm,0)
	1649	call Gather_Field(ps,ip1jmp1,1,0)
	1650	call Gather_Field(phis,ip1jmp1,1,0)
	1651	if (mpi_rank==0) then
	1652	CALL writedynav(itau,vcov,
	1653	& ucov,teta,pk,phi,q,masse,ps,phis)
	1654	endif
	1655	#endif
	1656	!$OMP END MASTER
	1657	ENDIF ! of IF (ok_dyn_ave)
	1658
[1146]	1659	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
[630]	1660
[1146]	1661
[1403]	1662	IF(MOD(itau,iecri ).EQ.0) THEN
	1663	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
[774]	1664	c$OMP BARRIER
	1665	c$OMP MASTER
[1146]	1666	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
[630]	1667
	1668	cym unat=0.
[1146]	1669	ijb=ij_begin
	1670	ije=ij_end
[949]	1671
[1146]	1672	if (pole_nord) then
	1673	ijb=ij_begin+iip1
	1674	unat(1:iip1,:)=0.
	1675	endif
[949]	1676
[1146]	1677	if (pole_sud) then
	1678	ije=ij_end-iip1
	1679	unat(ij_end-iip1+1:ij_end,:)=0.
	1680	endif
[949]	1681
[1146]	1682	do l=1,llm
	1683	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
	1684	enddo
[630]	1685
[1146]	1686	ijb=ij_begin
	1687	ije=ij_end
	1688	if (pole_sud) ije=ij_end-iip1
[949]	1689
[1146]	1690	do l=1,llm
	1691	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
	1692	enddo
[630]	1693
	1694	#ifdef CPP_IOIPSL
[1403]	1695	if (ok_dyn_ins) then
	1696	! Ehouarn: Gather fields and make master send to output
	1697	call Gather_Field(vcov,ip1jm,llm,0)
	1698	call Gather_Field(ucov,ip1jmp1,llm,0)
	1699	call Gather_Field(teta,ip1jmp1,llm,0)
	1700	call Gather_Field(phi,ip1jmp1,llm,0)
	1701	do iq=1,nqtot
	1702	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
	1703	enddo
	1704	call Gather_Field(masse,ip1jmp1,llm,0)
	1705	call Gather_Field(ps,ip1jmp1,1,0)
	1706	call Gather_Field(phis,ip1jmp1,1,0)
	1707	if (mpi_rank==0) then
	1708	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	1709	endif
[1279]	1710	! CALL writehist_p(histid, histvid, itau,vcov ,
	1711	! & ucov,teta,phi,q,masse,ps,phis)
[1403]	1712	endif ! of if (ok_dyn_ins)
[1146]	1713	#endif
	1714	! For some Grads output (but does it work?)
	1715	if (output_grads_dyn) then
	1716	call Gather_Field(unat,ip1jmp1,llm,0)
	1717	call Gather_Field(vnat,ip1jm,llm,0)
	1718	call Gather_Field(teta,ip1jmp1,llm,0)
	1719	call Gather_Field(ps,ip1jmp1,1,0)
	1720	do iq=1,nqtot
	1721	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
	1722	enddo
[630]	1723	c
[1146]	1724	if (mpi_rank==0) then
	1725	#include "write_grads_dyn.h"
	1726	endif
	1727	endif ! of if (output_grads_dyn)
[630]	1728
[774]	1729	c$OMP END MASTER
[1403]	1730	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
[630]	1731
[1146]	1732	IF(itau.EQ.itaufin) THEN
[1454]	1733	! if (planet_type.eq."earth") then
[774]	1734	c$OMP MASTER
	1735	CALL dynredem1_p("restart.nc",0.0,
[1146]	1736	. vcov,ucov,teta,q,masse,ps)
[774]	1737	c$OMP END MASTER
[1454]	1738	! endif ! of if (planet_type.eq."earth")
[1146]	1739	ENDIF ! of IF(itau.EQ.itaufin)
[630]	1740
[1146]	1741	forward = .TRUE.
	1742	GO TO 1
[630]	1743
[1146]	1744	ENDIF ! of IF (forward)
	1745
	1746	END IF ! of IF(.not.purmats)
[774]	1747	c$OMP MASTER
[1279]	1748	call fin_getparam
[1146]	1749	call finalize_parallel
[774]	1750	c$OMP END MASTER
[1146]	1751	RETURN
[630]	1752	END

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

Download in other formats:

Original Format