Context Navigation

leapfrog_p.F @ 1000

Last change on this file since 1000 was 953, checked in by slebonnois, 12 years ago
SL: optimisation pour le parallèle suite à tests Venus / petite correction appels routines secondaires dans Venus et Titan
File size: 62.1 KB

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

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

Download in other formats:

Original Format