Context Navigation

leapfrog_p.F @ 524

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

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

Download in other formats:

Original Format