Context Navigation

leapfrog_p.F @ 282

Last change on this file since 282 was 270, checked in by emillour, 14 years ago

Ehouarn: Mise a jour des dynamiques (seq et ) pour suivre
les changements et developpements de LMDZ5 terrestre
(mise a niveau avec LMDZ5 trunk, rev 1560). Ce qui ne devrais pas changer grand chose au fonctionnement de base du GCM).
Modification importante: correction du bug sur le cumul des flux de masse pour le transport des traceurs (mais dans les faits semble avoir peu d'impact).
(voir "commit_importants.log" pour les details des ajouts et modifications).

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

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

Download in other formats:

Original Format