Context Navigation

leapfrog_p.F @ 1286

Last change on this file since 1286 was 1286, checked in by Laurent Fairhead, 15 years ago
Les parametres definissant la transition eau glace dans les nuages sont maintenant lus dans physiq.def. Les valeurs par defaut donnent les memes resultats que precedemment. JLD
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 49.3 KB

Rev	Line
[764]	1	!
[1279]	2	! $Id: leapfrog_p.F 1286 2009-12-17 13:47:10Z fairhead $
[630]	3	!
	4	c
	5	c
	6
[1146]	7	SUBROUTINE leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
[630]	8	& time_0)
	9
	10	USE misc_mod
	11	USE parallel
	12	USE times
	13	USE mod_hallo
	14	USE Bands
	15	USE Write_Field
	16	USE Write_Field_p
	17	USE vampir
[985]	18	USE timer_filtre, ONLY : print_filtre_timer
[1146]	19	USE infotrac
[1279]	20	USE guide_p_mod, ONLY : guide_main
	21	USE getparam
[630]	22
	23	IMPLICIT NONE
	24
	25	c ...... Version du 10/01/98 ..........
	26
	27	c avec coordonnees verticales hybrides
	28	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
	29
	30	c=======================================================================
	31	c
	32	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
	33	c -------
	34	c
	35	c Objet:
	36	c ------
	37	c
	38	c GCM LMD nouvelle grille
	39	c
	40	c=======================================================================
	41	c
	42	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
	43	c et possibilite d'appeler une fonction f(y) a derivee tangente
	44	c hyperbolique a la place de la fonction a derivee sinusoidale.
	45
	46	c ... Possibilite de choisir le shema pour l'advection de
	47	c q , en modifiant iadv dans traceur.def (10/02) .
	48	c
	49	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
	50	c Pour Van-Leer iadv=10
	51	c
	52	c-----------------------------------------------------------------------
	53	c Declarations:
	54	c -------------
	55
	56	#include "dimensions.h"
	57	#include "paramet.h"
	58	#include "comconst.h"
	59	#include "comdissnew.h"
	60	#include "comvert.h"
	61	#include "comgeom.h"
	62	#include "logic.h"
	63	#include "temps.h"
	64	#include "control.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	INTEGER longcles
	73	PARAMETER ( longcles = 20 )
	74	REAL clesphy0( longcles )
	75
	76	real zqmin,zqmax
	77	INTEGER nbetatmoy, nbetatdem,nbetat
	78
	79	c variables dynamiques
[774]	80	REAL :: vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
	81	REAL :: teta(ip1jmp1,llm) ! temperature potentielle
[1146]	82	REAL :: q(ip1jmp1,llm,nqtot) ! champs advectes
[774]	83	REAL :: ps(ip1jmp1) ! pression au sol
	84	REAL,SAVE :: p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
	85	REAL,SAVE :: pks(ip1jmp1) ! exner au sol
	86	REAL,SAVE :: pk(ip1jmp1,llm) ! exner au milieu des couches
	87	REAL,SAVE :: pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
	88	REAL :: masse(ip1jmp1,llm) ! masse d'air
	89	REAL :: phis(ip1jmp1) ! geopotentiel au sol
	90	REAL,SAVE :: phi(ip1jmp1,llm) ! geopotentiel
	91	REAL,SAVE :: w(ip1jmp1,llm) ! vitesse verticale
[630]	92
	93	c variables dynamiques intermediaire pour le transport
[774]	94	REAL,SAVE :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
[630]	95
	96	c variables dynamiques au pas -1
[774]	97	REAL,SAVE :: vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
	98	REAL,SAVE :: tetam1(ip1jmp1,llm),psm1(ip1jmp1)
	99	REAL,SAVE :: massem1(ip1jmp1,llm)
[630]	100
	101	c tendances dynamiques
[774]	102	REAL,SAVE :: dv(ip1jm,llm),du(ip1jmp1,llm)
[1146]	103	REAL,SAVE :: dteta(ip1jmp1,llm),dp(ip1jmp1)
	104	REAL,DIMENSION(:,:,:), ALLOCATABLE, SAVE :: dq
[630]	105
	106	c tendances de la dissipation
[774]	107	REAL,SAVE :: dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
	108	REAL,SAVE :: dtetadis(ip1jmp1,llm)
[630]	109
	110	c tendances physiques
[774]	111	REAL,SAVE :: dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
	112	REAL,SAVE :: dtetafi(ip1jmp1,llm)
[1146]	113	REAL,SAVE :: dpfi(ip1jmp1)
	114	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi
[630]	115
	116	c variables pour le fichier histoire
	117	REAL dtav ! intervalle de temps elementaire
	118
	119	REAL tppn(iim),tpps(iim),tpn,tps
	120	c
	121	INTEGER itau,itaufinp1,iav
[1279]	122	! INTEGER iday ! jour julien
	123	REAL time
[630]	124
	125	REAL SSUM
[774]	126	REAL time_0
	127	REAL,SAVE :: finvmaold(ip1jmp1,llm)
[630]	128
	129	cym LOGICAL lafin
[774]	130	LOGICAL :: lafin
[630]	131	INTEGER ij,iq,l
	132	INTEGER ik
	133
	134	real time_step, t_wrt, t_ops
	135
[1279]	136	! jD_cur: jour julien courant
	137	! jH_cur: heure julienne courante
	138	REAL :: jD_cur, jH_cur
	139	INTEGER :: an, mois, jour
	140	REAL :: secondes
	141
[630]	142	LOGICAL first,callinigrads
	143
	144	data callinigrads/.true./
	145	character*10 string10
	146
[774]	147	REAL,SAVE :: alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
[960]	148	REAL,SAVE :: flxw(ip1jmp1,llm) ! flux de masse verticale
[630]	149
	150	c+jld variables test conservation energie
[774]	151	REAL,SAVE :: ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
[630]	152	C Tendance de la temp. potentiel d (theta)/ d t due a la
	153	C tansformation d'energie cinetique en energie thermique
	154	C cree par la dissipation
[774]	155	REAL,SAVE :: dtetaecdt(ip1jmp1,llm)
	156	REAL,SAVE :: vcont(ip1jm,llm),ucont(ip1jmp1,llm)
	157	REAL,SAVE :: vnat(ip1jm,llm),unat(ip1jmp1,llm)
[630]	158	REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
	159	CHARACTER*15 ztit
[764]	160	! INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
	161	! SAVE ip_ebil_dyn
	162	! DATA ip_ebil_dyn/0/
[630]	163	c-jld
	164
	165	character*80 dynhist_file, dynhistave_file
	166	character*20 modname
	167	character*80 abort_message
	168
	169
[774]	170	logical,PARAMETER :: dissip_conservative=.TRUE.
	171
[630]	172	INTEGER testita
	173	PARAMETER (testita = 9)
[1286]	174
	175	logical , parameter :: flag_verif = .false.
[630]	176
[764]	177	c declaration liees au parallelisme
[630]	178	INTEGER :: ierr
[774]	179	LOGICAL :: FirstCaldyn
	180	LOGICAL :: FirstPhysic
[630]	181	INTEGER :: ijb,ije,j,i
	182	type(Request) :: TestRequest
	183	type(Request) :: Request_Dissip
	184	type(Request) :: Request_physic
[774]	185	REAL,SAVE :: dvfi_tmp(iip1,llm),dufi_tmp(iip1,llm)
[1146]	186	REAL,SAVE :: dtetafi_tmp(iip1,llm)
	187	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi_tmp
[774]	188	REAL,SAVE :: dpfi_tmp(iip1)
	189
	190	INTEGER :: true_itau
[630]	191	LOGICAL :: verbose=.true.
[774]	192	INTEGER :: iapptrac
	193	INTEGER :: AdjustCount
[1279]	194	! INTEGER :: var_time
[774]	195	LOGICAL :: ok_start_timer=.FALSE.
[1146]	196	LOGICAL, SAVE :: firstcall=.TRUE.
[774]	197
	198	c$OMP MASTER
[630]	199	ItCount=0
[774]	200	c$OMP END MASTER
	201	true_itau=0
	202	FirstCaldyn=.TRUE.
	203	FirstPhysic=.TRUE.
	204	iapptrac=0
	205	AdjustCount = 0
	206	lafin=.false.
[630]	207
	208	itaufin = nday*day_step
	209	itaufinp1 = itaufin +1
[1146]	210	modname="leapfrog_p"
[630]	211
	212	itau = 0
[1279]	213	! iday = day_ini+itau/day_step
	214	! time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	215	! IF(time.GT.1.) THEN
	216	! time = time-1.
	217	! iday = iday+1
	218	! ENDIF
[630]	219
[1146]	220	c Allocate variables depending on dynamic variable nqtot
	221	c$OMP MASTER
	222	IF (firstcall) THEN
	223	firstcall=.FALSE.
	224	ALLOCATE(dq(ip1jmp1,llm,nqtot))
	225	ALLOCATE(dqfi(ip1jmp1,llm,nqtot))
	226	ALLOCATE(dqfi_tmp(iip1,llm,nqtot))
	227	END IF
	228	c$OMP END MASTER
	229	c$OMP BARRIER
	230
[630]	231	c-----------------------------------------------------------------------
	232	c On initialise la pression et la fonction d'Exner :
	233	c --------------------------------------------------
	234
[774]	235	c$OMP MASTER
[630]	236	dq=0.
	237	CALL pression ( ip1jmp1, ap, bp, ps, p )
	238	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
[774]	239	c$OMP END MASTER
[630]	240	c-----------------------------------------------------------------------
	241	c Debut de l'integration temporelle:
	242	c ----------------------------------
[764]	243	c et du parallelisme !!
[630]	244
	245	1 CONTINUE
	246
[1279]	247	jD_cur = jD_ref + day_ini - day_ref + int (itau * dtvr / daysec)
	248	jH_cur = jH_ref + &
	249	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
[985]	250
[630]	251
	252	#ifdef CPP_IOIPSL
[1279]	253	if (ok_guide) then
	254	!$OMP MASTER
	255	call guide_main(itau,ucov,vcov,teta,q,masse,ps)
	256	!$OMP END MASTER
	257	!$OMP BARRIER
[630]	258	endif
	259	#endif
[1279]	260
[630]	261	c
	262	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
	263	c CALL test_period ( ucov,vcov,teta,q,p,phis )
	264	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
	265	c ENDIF
	266	c
	267	cym CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
	268	cym CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
	269	cym CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
	270	cym CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
	271	cym CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
	272
	273	if (FirstCaldyn) then
[774]	274	c$OMP MASTER
[630]	275	ucovm1=ucov
	276	vcovm1=vcov
	277	tetam1= teta
	278	massem1= masse
	279	psm1= ps
	280
	281	finvmaold = masse
	282	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
[774]	283	c$OMP END MASTER
	284	c$OMP BARRIER
[630]	285	else
[1146]	286	! Save fields obtained at previous time step as '...m1'
[630]	287	ijb=ij_begin
	288	ije=ij_end
[774]	289
	290	c$OMP MASTER
[630]	291	psm1 (ijb:ije) = ps (ijb:ije)
[774]	292	c$OMP END MASTER
	293
	294	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[949]	295	DO l=1,llm
[774]	296	ije=ij_end
[949]	297	ucovm1 (ijb:ije,l) = ucov (ijb:ije,l)
[774]	298	tetam1 (ijb:ije,l) = teta (ijb:ije,l)
	299	massem1 (ijb:ije,l) = masse (ijb:ije,l)
	300	finvmaold(ijb:ije,l)=masse(ijb:ije,l)
[949]	301
[774]	302	if (pole_sud) ije=ij_end-iip1
	303	vcovm1(ijb:ije,l) = vcov (ijb:ije,l)
[630]	304
[774]	305
	306	ENDDO
[985]	307	c$OMP ENDDO
[774]	308
	309
	310	CALL filtreg_p ( finvmaold ,jj_begin,jj_end,jjp1,
[630]	311	. llm, -2,2, .TRUE., 1 )
	312
[1146]	313	endif ! of if (FirstCaldyn)
[630]	314
	315	forward = .TRUE.
	316	leapf = .FALSE.
	317	dt = dtvr
	318
	319	c ... P.Le Van .26/04/94 ....
	320
	321	cym CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
	322	cym CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
	323
[764]	324	cym ne sert a rien
[630]	325	cym call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
	326
	327	2 CONTINUE
	328
[774]	329	c$OMP MASTER
[630]	330	ItCount=ItCount+1
[995]	331	if (MOD(ItCount,1)==1) then
[630]	332	debug=.true.
	333	else
	334	debug=.false.
	335	endif
[774]	336	c$OMP END MASTER
[630]	337	c-----------------------------------------------------------------------
	338
	339	c date:
	340	c -----
	341
	342
	343	c gestion des appels de la physique et des dissipations:
	344	c ------------------------------------------------------
	345	c
	346	c ... P.Le Van ( 6/02/95 ) ....
	347
	348	apphys = .FALSE.
	349	statcl = .FALSE.
	350	conser = .FALSE.
	351	apdiss = .FALSE.
[764]	352	c idissip=1
[630]	353	IF( purmats ) THEN
	354	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
	355	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
	356	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
[1146]	357	s .and. iflag_phys.EQ.1 ) apphys = .TRUE.
[630]	358	ELSE
	359	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
	360	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
[1146]	361	IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
[630]	362	END IF
	363
	364	cym ---> Pour le moment
	365	cym apphys = .FALSE.
	366	statcl = .FALSE.
	367	conser = .FALSE.
	368
	369	if (firstCaldyn) then
[774]	370	c$OMP MASTER
[630]	371	call SetDistrib(jj_Nb_Caldyn)
[774]	372	c$OMP END MASTER
	373	c$OMP BARRIER
[630]	374	firstCaldyn=.FALSE.
	375	cym call InitTime
[774]	376	c$OMP MASTER
[630]	377	call Init_timer
[774]	378	c$OMP END MASTER
[630]	379	endif
[774]	380
	381	c$OMP MASTER
	382	IF (ok_start_timer) THEN
	383	CALL InitTime
[949]	384	ok_start_timer=.FALSE.
[774]	385	ENDIF
	386	c$OMP END MASTER
	387
[630]	388	if (Adjust) then
[774]	389	c$OMP MASTER
	390	AdjustCount=AdjustCount+1
	391	if (iapptrac==iapp_tracvl .and. (forward. OR . leapf)
	392	& .and. itau/iphysiq>2 .and. Adjustcount>30) then
[630]	393	AdjustCount=0
	394	call allgather_timer_average
	395
	396	if (Verbose) then
	397
[949]	398	print ,'********************************'
	399	print ,'*** TIMER CALDYN ****'
	400	do i=0,mpi_size-1
	401	print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
[630]	402	& ' : temps moyen :',
	403	& timer_average(jj_nb_caldyn(i),timer_caldyn,i),
	404	& '+-',timer_delta(jj_nb_caldyn(i),timer_caldyn,i)
	405	enddo
	406
	407	print ,'********************************'
[949]	408	print ,'*** TIMER VANLEER ****'
	409	do i=0,mpi_size-1
	410	print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
[630]	411	& ' : temps moyen :',
	412	& timer_average(jj_nb_vanleer(i),timer_vanleer,i),
	413	& '+-',timer_delta(jj_nb_vanleer(i),timer_vanleer,i)
	414	enddo
	415
	416	print ,'********************************'
[949]	417	print ,'*** TIMER DISSIP ****'
	418	do i=0,mpi_size-1
	419	print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
[630]	420	& ' : temps moyen :',
	421	& timer_average(jj_nb_dissip(i),timer_dissip,i),
	422	& '+-',timer_delta(jj_nb_dissip(i),timer_dissip,i)
	423	enddo
	424
[949]	425	if (mpi_rank==0) call WriteBands
	426
[630]	427	endif
	428
[949]	429	call AdjustBands_caldyn
[630]	430	if (mpi_rank==0) call WriteBands
[949]	431
	432	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
[630]	433	& jj_Nb_caldyn,0,0,TestRequest)
	434	call Register_SwapFieldHallo(ucovm1,ucovm1,ip1jmp1,llm,
	435	& jj_Nb_caldyn,0,0,TestRequest)
	436	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
	437	& jj_Nb_caldyn,0,0,TestRequest)
	438	call Register_SwapFieldHallo(vcovm1,vcovm1,ip1jm,llm,
	439	& jj_Nb_caldyn,0,0,TestRequest)
	440	call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
	441	& jj_Nb_caldyn,0,0,TestRequest)
	442	call Register_SwapFieldHallo(tetam1,tetam1,ip1jmp1,llm,
	443	& jj_Nb_caldyn,0,0,TestRequest)
	444	call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
	445	& jj_Nb_caldyn,0,0,TestRequest)
	446	call Register_SwapFieldHallo(massem1,massem1,ip1jmp1,llm,
	447	& jj_Nb_caldyn,0,0,TestRequest)
	448	call Register_SwapFieldHallo(ps,ps,ip1jmp1,1,
	449	& jj_Nb_caldyn,0,0,TestRequest)
	450	call Register_SwapFieldHallo(psm1,psm1,ip1jmp1,1,
	451	& jj_Nb_caldyn,0,0,TestRequest)
	452	call Register_SwapFieldHallo(pkf,pkf,ip1jmp1,llm,
	453	& jj_Nb_caldyn,0,0,TestRequest)
	454	call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
	455	& jj_Nb_caldyn,0,0,TestRequest)
	456	call Register_SwapFieldHallo(pks,pks,ip1jmp1,1,
	457	& jj_Nb_caldyn,0,0,TestRequest)
	458	call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
	459	& jj_Nb_caldyn,0,0,TestRequest)
[774]	460	call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
[630]	461	& jj_Nb_caldyn,0,0,TestRequest)
	462	call Register_SwapFieldHallo(finvmaold,finvmaold,ip1jmp1,llm,
	463	& jj_Nb_caldyn,0,0,TestRequest)
	464
[1146]	465	do j=1,nqtot
[764]	466	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
	467	& jj_nb_caldyn,0,0,TestRequest)
	468	enddo
	469
[630]	470	call SetDistrib(jj_nb_caldyn)
	471	call SendRequest(TestRequest)
	472	call WaitRequest(TestRequest)
	473
[949]	474	call AdjustBands_dissip
	475	call AdjustBands_physic
[774]	476
[630]	477	endif
[774]	478	c$OMP END MASTER
[630]	479	endif
[774]	480
[630]	481
	482
	483	c-----------------------------------------------------------------------
	484	c calcul des tendances dynamiques:
	485	c --------------------------------
[774]	486	c$OMP BARRIER
	487	c$OMP MASTER
[630]	488	call VTb(VThallo)
[985]	489	c$OMP END MASTER
	490
[630]	491	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,TestRequest)
	492	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,TestRequest)
	493	call Register_Hallo(teta,ip1jmp1,llm,1,1,1,1,TestRequest)
	494	call Register_Hallo(ps,ip1jmp1,1,1,2,2,1,TestRequest)
	495	call Register_Hallo(pkf,ip1jmp1,llm,1,1,1,1,TestRequest)
	496	call Register_Hallo(pk,ip1jmp1,llm,1,1,1,1,TestRequest)
	497	call Register_Hallo(pks,ip1jmp1,1,1,1,1,1,TestRequest)
	498	call Register_Hallo(p,ip1jmp1,llmp1,1,1,1,1,TestRequest)
	499
[1146]	500	c do j=1,nqtot
[630]	501	c call Register_Hallo(q(1,1,j),ip1jmp1,llm,1,1,1,1,
	502	c * TestRequest)
	503	c enddo
	504
	505	call SendRequest(TestRequest)
[985]	506	c$OMP BARRIER
[630]	507	call WaitRequest(TestRequest)
[985]	508
	509	c$OMP MASTER
[630]	510	call VTe(VThallo)
[985]	511	c$OMP END MASTER
	512	c$OMP BARRIER
[764]	513
[949]	514	if (debug) then
[985]	515	!$OMP BARRIER
	516	!$OMP MASTER
[630]	517	call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	518	call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
	519	call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
	520	call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
	521	call WriteField_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
[949]	522	call WriteField_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
[630]	523	call WriteField_p('pks',reshape(pks,(/iip1,jmp1/)))
	524	call WriteField_p('pkf',reshape(pkf,(/iip1,jmp1,llm/)))
	525	call WriteField_p('phis',reshape(phis,(/iip1,jmp1/)))
[1146]	526	do j=1,nqtot
[764]	527	call WriteField_p('q'//trim(int2str(j)),
	528	. reshape(q(:,:,j),(/iip1,jmp1,llm/)))
[985]	529	enddo
	530	!$OMP END MASTER
	531	c$OMP BARRIER
[630]	532	endif
[985]	533
[630]	534
	535	True_itau=True_itau+1
[774]	536
	537	c$OMP MASTER
[1146]	538	IF (prt_level>9) THEN
	539	WRITE(lunout,*)"leapfrog_p: Iteration No",True_itau
	540	ENDIF
[630]	541
	542
	543	call start_timer(timer_caldyn)
	544
	545	CALL geopot_p ( ip1jmp1, teta , pk , pks, phis , phi )
	546
	547
	548	call VTb(VTcaldyn)
[774]	549	c$OMP END MASTER
[1279]	550	! var_time=time+iday-day_ini
[764]	551
[985]	552	c$OMP BARRIER
	553	! CALL FTRACE_REGION_BEGIN("caldyn")
[1279]	554	time = jD_cur + jH_cur
[630]	555	CALL caldyn_p
	556	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
[1279]	557	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
[764]	558
[985]	559	! CALL FTRACE_REGION_END("caldyn")
[1279]	560
[774]	561	c$OMP MASTER
[630]	562	call VTe(VTcaldyn)
[774]	563	c$OMP END MASTER
[985]	564
	565	cc$OMP BARRIER
	566	cc$OMP MASTER
[1279]	567	! call WriteField_p('du',reshape(du,(/iip1,jmp1,llm/)))
	568	! call WriteField_p('dv',reshape(dv,(/iip1,jjm,llm/)))
	569	! call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
	570	! call WriteField_p('dp',reshape(dp,(/iip1,jmp1/)))
	571	! call WriteField_p('w',reshape(w,(/iip1,jmp1,llm/)))
	572	! call WriteField_p('pbaru',reshape(pbaru,(/iip1,jmp1,llm/)))
	573	! call WriteField_p('pbarv',reshape(pbarv,(/iip1,jjm,llm/)))
	574	! call WriteField_p('p',reshape(p,(/iip1,jmp1,llmp1/)))
	575	! call WriteField_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
	576	! call WriteField_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
[985]	577	cc$OMP END MASTER
[630]	578
	579	c-----------------------------------------------------------------------
	580	c calcul des tendances advection des traceurs (dont l'humidite)
	581	c -------------------------------------------------------------
	582
	583	IF( forward. OR . leapf ) THEN
[774]	584	cc$OMP PARALLEL DEFAULT(SHARED)
[630]	585	c
[960]	586	CALL caladvtrac_p(q,pbaru,pbarv,
	587	* p, masse, dq, teta,
	588	. flxw,pk, iapptrac)
[985]	589
[764]	590	IF (offline) THEN
[630]	591	Cmaf stokage du flux de masse pour traceurs OFF-LINE
[985]	592
[630]	593	#ifdef CPP_IOIPSL
[985]	594	CALL fluxstokenc_p(pbaru,pbarv,masse,teta,phi,phis,
[630]	595	. dtvr, itau)
	596	#endif
	597
	598
[1146]	599	ENDIF ! of IF (offline)
[630]	600	c
[1146]	601	ENDIF ! of IF( forward. OR . leapf )
[774]	602	cc$OMP END PARALLEL
[630]	603
	604	c-----------------------------------------------------------------------
	605	c integrations dynamique et traceurs:
	606	c ----------------------------------
[774]	607
	608	c$OMP MASTER
[630]	609	call VTb(VTintegre)
[774]	610	c$OMP END MASTER
[764]	611	c call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
	612	c call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
	613	c call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
	614	c call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
	615	c call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	616	c call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
	617	c call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
	618	c call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
[774]	619	cc$OMP PARALLEL DEFAULT(SHARED)
[985]	620	c$OMP BARRIER
	621	! CALL FTRACE_REGION_BEGIN("integrd")
[1146]	622
[630]	623	CALL integrd_p ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
	624	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
	625	$ finvmaold )
	626
[985]	627	! CALL FTRACE_REGION_END("integrd")
	628	c$OMP BARRIER
	629	cc$OMP MASTER
[764]	630	c call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
	631	c call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
	632	c call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
	633	c call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
	634	c call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	635	c call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
	636	c call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
	637	c call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
[985]	638	c
	639	c call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
[1146]	640	c do j=1,nqtot
[985]	641	c call WriteField_p('q'//trim(int2str(j)),
	642	c . reshape(q(:,:,j),(/iip1,jmp1,llm/)))
	643	c call WriteField_p('dq'//trim(int2str(j)),
	644	c . reshape(dq(:,:,j),(/iip1,jmp1,llm/)))
	645	c enddo
	646	cc$OMP END MASTER
[764]	647
[985]	648
[774]	649	c$OMP MASTER
[630]	650	call VTe(VTintegre)
[774]	651	c$OMP END MASTER
[630]	652	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
	653	c
	654	c-----------------------------------------------------------------------
	655	c calcul des tendances physiques:
	656	c -------------------------------
	657	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
	658	c
	659	IF( purmats ) THEN
	660	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
	661	ELSE
	662	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
	663	ENDIF
[774]	664
	665	cc$OMP END PARALLEL
	666
[630]	667	c
	668	c
	669	IF( apphys ) THEN
	670	c
	671	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
	672	c
[774]	673	cc$OMP PARALLEL DEFAULT(SHARED)
	674	cc$OMP+ PRIVATE(rdaym_ini,rdayvrai,ijb,ije)
[764]	675
	676	c$OMP MASTER
[630]	677	call suspend_timer(timer_caldyn)
[1146]	678
	679	write(lunout,*)
	680	& 'leapfrog_p: Entree dans la physique : Iteration No ',true_itau
[764]	681	c$OMP END MASTER
	682
[630]	683	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
[764]	684
[985]	685	c$OMP BARRIER
[630]	686	CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
[985]	687	c$OMP BARRIER
[1279]	688	jD_cur = jD_ref + day_ini - day_ref
	689	$ + int (itau * dtvr / daysec)
	690	jH_cur = jH_ref + &
	691	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
	692	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
[630]	693
	694	c rajout debug
	695	c lafin = .true.
	696
	697
	698	c Inbterface avec les routines de phylmd (phymars ... )
	699	c -----------------------------------------------------
	700
	701	c+jld
	702
[764]	703	c Diagnostique de conservation de l'energie : initialisation
[630]	704	IF (ip_ebil_dyn.ge.1 ) THEN
	705	ztit='bil dyn'
[1146]	706	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
	707	IF (planet_type.eq."earth") THEN
	708	CALL diagedyn(ztit,2,1,1,dtphys
	709	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	710	ENDIF
[630]	711	ENDIF
	712	c-jld
[764]	713	c$OMP BARRIER
	714	c$OMP MASTER
[630]	715	call VTb(VThallo)
[985]	716	c$OMP END MASTER
	717
[630]	718	call SetTag(Request_physic,800)
[949]	719
[630]	720	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
	721	* jj_Nb_physic,2,2,Request_physic)
	722
	723	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
	724	* jj_Nb_physic,2,2,Request_physic)
	725
	726	call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
	727	* jj_Nb_physic,2,2,Request_physic)
	728
[1279]	729	call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
	730	* jj_Nb_physic,1,2,Request_physic)
	731
[630]	732	call Register_SwapFieldHallo(p,p,ip1jmp1,llmp1,
	733	* jj_Nb_physic,2,2,Request_physic)
	734
	735	call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
	736	* jj_Nb_physic,2,2,Request_physic)
	737
	738	call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
	739	* jj_Nb_physic,2,2,Request_physic)
	740
	741	call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
	742	* jj_Nb_physic,2,2,Request_physic)
	743
	744	call Register_SwapFieldHallo(w,w,ip1jmp1,llm,
	745	* jj_Nb_physic,2,2,Request_physic)
	746
[949]	747	c call SetDistrib(jj_nb_vanleer)
[1146]	748	do j=1,nqtot
[630]	749
	750	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
	751	* jj_Nb_physic,2,2,Request_physic)
	752	enddo
[960]	753
[630]	754	call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
	755	* jj_Nb_physic,2,2,Request_physic)
[949]	756
[630]	757	call SendRequest(Request_Physic)
[985]	758	c$OMP BARRIER
[630]	759	call WaitRequest(Request_Physic)
[985]	760
	761	c$OMP BARRIER
	762	c$OMP MASTER
	763	call SetDistrib(jj_nb_Physic)
[949]	764	call VTe(VThallo)
	765
	766	call VTb(VTphysiq)
[764]	767	c$OMP END MASTER
	768	c$OMP BARRIER
	769
[949]	770	cc$OMP MASTER
[764]	771	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
	772	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
	773	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
	774	c call WriteField_p('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
	775	c call WriteField_p('pkfi',reshape(pk,(/iip1,jmp1,llm/)))
	776	cc$OMP END MASTER
	777	cc$OMP BARRIER
[985]	778	! CALL FTRACE_REGION_BEGIN("calfis")
[1279]	779	CALL calfis_p(lafin ,jD_cur, jH_cur,
[630]	780	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
[960]	781	$ du,dv,dteta,dq,
[630]	782	$ flxw,
	783	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
[985]	784	! CALL FTRACE_REGION_END("calfis")
[630]	785	ijb=ij_begin
	786	ije=ij_end
	787	if ( .not. pole_nord) then
[764]	788	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	789	DO l=1,llm
	790	dufi_tmp(1:iip1,l) = dufi(ijb:ijb+iim,l)
	791	dvfi_tmp(1:iip1,l) = dvfi(ijb:ijb+iim,l)
	792	dtetafi_tmp(1:iip1,l)= dtetafi(ijb:ijb+iim,l)
	793	dqfi_tmp(1:iip1,l,:) = dqfi(ijb:ijb+iim,l,:)
	794	ENDDO
	795	c$OMP END DO NOWAIT
	796
	797	c$OMP MASTER
[630]	798	dpfi_tmp(1:iip1) = dpfi(ijb:ijb+iim)
[764]	799	c$OMP END MASTER
[1146]	800	endif ! of if ( .not. pole_nord)
[764]	801
	802	c$OMP BARRIER
	803	c$OMP MASTER
[949]	804	call SetDistrib(jj_nb_Physic_bis)
[630]	805
	806	call VTb(VThallo)
[985]	807	c$OMP END MASTER
	808	c$OMP BARRIER
[630]	809
	810	call Register_Hallo(dufi,ip1jmp1,llm,
	811	* 1,0,0,1,Request_physic)
	812
	813	call Register_Hallo(dvfi,ip1jm,llm,
	814	* 1,0,0,1,Request_physic)
	815
	816	call Register_Hallo(dtetafi,ip1jmp1,llm,
	817	* 1,0,0,1,Request_physic)
	818
	819	call Register_Hallo(dpfi,ip1jmp1,1,
	820	* 1,0,0,1,Request_physic)
	821
[1146]	822	do j=1,nqtot
[630]	823	call Register_Hallo(dqfi(1,1,j),ip1jmp1,llm,
	824	* 1,0,0,1,Request_physic)
	825	enddo
	826
	827	call SendRequest(Request_Physic)
[985]	828	c$OMP BARRIER
[630]	829	call WaitRequest(Request_Physic)
	830
[985]	831	c$OMP BARRIER
	832	c$OMP MASTER
[630]	833	call VTe(VThallo)
	834
	835	call SetDistrib(jj_nb_Physic)
[764]	836	c$OMP END MASTER
[949]	837	c$OMP BARRIER
	838	ijb=ij_begin
	839	if (.not. pole_nord) then
	840
[764]	841	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	842	DO l=1,llm
	843	dufi(ijb:ijb+iim,l) = dufi(ijb:ijb+iim,l)+dufi_tmp(1:iip1,l)
	844	dvfi(ijb:ijb+iim,l) = dvfi(ijb:ijb+iim,l)+dvfi_tmp(1:iip1,l)
	845	dtetafi(ijb:ijb+iim,l) = dtetafi(ijb:ijb+iim,l)
	846	& +dtetafi_tmp(1:iip1,l)
	847	dqfi(ijb:ijb+iim,l,:) = dqfi(ijb:ijb+iim,l,:)
	848	& + dqfi_tmp(1:iip1,l,:)
	849	ENDDO
	850	c$OMP END DO NOWAIT
	851
	852	c$OMP MASTER
[630]	853	dpfi(ijb:ijb+iim) = dpfi(ijb:ijb+iim)+ dpfi_tmp(1:iip1)
[764]	854	c$OMP END MASTER
[949]	855
[1146]	856	endif ! of if (.not. pole_nord)
[764]	857	c$OMP BARRIER
[949]	858	cc$OMP MASTER
[630]	859	c call WriteField_p('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
	860	c call WriteField_p('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
	861	c call WriteField_p('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
	862	c call WriteField_p('dpfi',reshape(dpfi,(/iip1,jmp1/)))
[764]	863	cc$OMP END MASTER
[630]	864	c
[1146]	865	c do j=1,nqtot
[630]	866	c call WriteField_p('dqfi'//trim(int2str(j)),
	867	c . reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
	868	c enddo
	869
	870	c ajout des tendances physiques:
	871	c ------------------------------
[1000]	872	IF (ok_strato) THEN
[1279]	873	CALL top_bound_p( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
[1000]	874	ENDIF
	875
[1146]	876	CALL addfi_p( dtphys, leapf, forward ,
[630]	877	$ ucov, vcov, teta , q ,ps ,
	878	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	879
[764]	880	c$OMP BARRIER
	881	c$OMP MASTER
[630]	882	call VTe(VTphysiq)
	883
	884	call VTb(VThallo)
[985]	885	c$OMP END MASTER
[630]	886
	887	call SetTag(Request_physic,800)
	888	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
	889	* jj_Nb_caldyn,Request_physic)
	890
	891	call Register_SwapField(vcov,vcov,ip1jm,llm,
	892	* jj_Nb_caldyn,Request_physic)
	893
	894	call Register_SwapField(teta,teta,ip1jmp1,llm,
	895	* jj_Nb_caldyn,Request_physic)
	896
[1279]	897	call Register_SwapField(masse,masse,ip1jmp1,llm,
	898	* jj_Nb_caldyn,Request_physic)
	899
[630]	900	call Register_SwapField(p,p,ip1jmp1,llmp1,
	901	* jj_Nb_caldyn,Request_physic)
	902
	903	call Register_SwapField(pk,pk,ip1jmp1,llm,
	904	* jj_Nb_caldyn,Request_physic)
	905
	906	call Register_SwapField(phis,phis,ip1jmp1,1,
	907	* jj_Nb_caldyn,Request_physic)
	908
	909	call Register_SwapField(phi,phi,ip1jmp1,llm,
	910	* jj_Nb_caldyn,Request_physic)
	911
	912	call Register_SwapField(w,w,ip1jmp1,llm,
	913	* jj_Nb_caldyn,Request_physic)
	914
[1146]	915	do j=1,nqtot
[630]	916
	917	call Register_SwapField(q(1,1,j),q(1,1,j),ip1jmp1,llm,
	918	* jj_Nb_caldyn,Request_physic)
	919
	920	enddo
	921
	922	call SendRequest(Request_Physic)
[985]	923	c$OMP BARRIER
[630]	924	call WaitRequest(Request_Physic)
	925
[985]	926	c$OMP BARRIER
	927	c$OMP MASTER
	928	call VTe(VThallo)
	929	call SetDistrib(jj_Nb_caldyn)
[764]	930	c$OMP END MASTER
	931	c$OMP BARRIER
[630]	932	c
[764]	933	c Diagnostique de conservation de l'energie : difference
[630]	934	IF (ip_ebil_dyn.ge.1 ) THEN
	935	ztit='bil phys'
	936	CALL diagedyn(ztit,2,1,1,dtphys
	937	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	938	ENDIF
[764]	939
	940	cc$OMP MASTER
	941	c if (debug) then
	942	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
	943	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
	944	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
	945	c endif
	946	cc$OMP END MASTER
	947
[630]	948
[1146]	949	c-jld
	950	c$OMP MASTER
	951	call resume_timer(timer_caldyn)
	952	if (FirstPhysic) then
	953	ok_start_timer=.TRUE.
	954	FirstPhysic=.false.
	955	endif
	956	c$OMP END MASTER
	957	ENDIF ! of IF( apphys )
	958
	959	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
[630]	960	c Calcul academique de la physique = Rappel Newtonien + fritcion
	961	c --------------------------------------------------------------
	962	cym teta(:,:)=teta(:,:)
	963	cym s -iphysiqdtvr(teta(:,:)-tetarappel(:,:))/taurappel
	964	ijb=ij_begin
	965	ije=ij_end
	966	teta(ijb:ije,:)=teta(ijb:ije,:)
	967	s -iphysiqdtvr(teta(ijb:ije,:)-tetarappel(ijb:ije,:))/taurappel
	968
	969	call Register_Hallo(ucov,ip1jmp1,llm,0,1,1,0,Request_Physic)
	970	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Physic)
	971	call SendRequest(Request_Physic)
[1279]	972	c$OMP BARRIER
[630]	973	call WaitRequest(Request_Physic)
	974
	975	call friction_p(ucov,vcov,iphysiq*dtvr)
[1146]	976	ENDIF ! of IF(iflag_phys.EQ.2)
[630]	977
	978
[985]	979	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
[774]	980	c$OMP BARRIER
[630]	981	CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
[774]	982	c$OMP BARRIER
[630]	983
[774]	984	cc$OMP END PARALLEL
[630]	985
	986	c-----------------------------------------------------------------------
	987	c dissipation horizontale et verticale des petites echelles:
	988	c ----------------------------------------------------------
	989
	990	IF(apdiss) THEN
[774]	991	cc$OMP PARALLEL DEFAULT(SHARED)
	992	cc$OMP+ PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
[764]	993	c$OMP MASTER
[630]	994	call suspend_timer(timer_caldyn)
[949]	995
	996	c print*,'Entree dans la dissipation : Iteration No ',true_itau
[630]	997	c calcul de l'energie cinetique avant dissipation
[949]	998	c print *,'Passage dans la dissipation'
[630]	999
	1000	call VTb(VThallo)
[764]	1001	c$OMP END MASTER
[630]	1002
[764]	1003	c$OMP BARRIER
[985]	1004
[630]	1005	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
	1006	* jj_Nb_dissip,1,1,Request_dissip)
	1007
	1008	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
	1009	* jj_Nb_dissip,1,1,Request_dissip)
	1010
	1011	call Register_SwapField(teta,teta,ip1jmp1,llm,
	1012	* jj_Nb_dissip,Request_dissip)
	1013
	1014	call Register_SwapField(p,p,ip1jmp1,llmp1,
	1015	* jj_Nb_dissip,Request_dissip)
	1016
	1017	call Register_SwapField(pk,pk,ip1jmp1,llm,
	1018	* jj_Nb_dissip,Request_dissip)
	1019
	1020	call SendRequest(Request_dissip)
[985]	1021	c$OMP BARRIER
[630]	1022	call WaitRequest(Request_dissip)
[985]	1023
	1024	c$OMP BARRIER
	1025	c$OMP MASTER
[630]	1026	call SetDistrib(jj_Nb_dissip)
[949]	1027	call VTe(VThallo)
[630]	1028	call VTb(VTdissipation)
[949]	1029	call start_timer(timer_dissip)
[764]	1030	c$OMP END MASTER
	1031	c$OMP BARRIER
	1032
[630]	1033	call covcont_p(llm,ucov,vcov,ucont,vcont)
	1034	call enercin_p(vcov,ucov,vcont,ucont,ecin0)
	1035
	1036	c dissipation
	1037
[985]	1038	! CALL FTRACE_REGION_BEGIN("dissip")
[630]	1039	CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
[985]	1040	! CALL FTRACE_REGION_END("dissip")
[764]	1041
[949]	1042	ijb=ij_begin
	1043	ije=ij_end
	1044	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1045	DO l=1,llm
	1046	ucov(ijb:ije,l)=ucov(ijb:ije,l)+dudis(ijb:ije,l)
[764]	1047	ENDDO
[949]	1048	c$OMP END DO NOWAIT
	1049	if (pole_sud) ije=ije-iip1
	1050	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1051	DO l=1,llm
	1052	vcov(ijb:ije,l)=vcov(ijb:ije,l)+dvdis(ijb:ije,l)
[764]	1053	ENDDO
[949]	1054	c$OMP END DO NOWAIT
[764]	1055
[630]	1056	c teta=teta+dtetadis
	1057
	1058
	1059	c------------------------------------------------------------------------
	1060	if (dissip_conservative) then
	1061	C On rajoute la tendance due a la transform. Ec -> E therm. cree
	1062	C lors de la dissipation
[764]	1063	c$OMP BARRIER
	1064	c$OMP MASTER
[630]	1065	call suspend_timer(timer_dissip)
[949]	1066	call VTb(VThallo)
[985]	1067	c$OMP END MASTER
[630]	1068	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Request_Dissip)
	1069	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Dissip)
	1070	call SendRequest(Request_Dissip)
[985]	1071	c$OMP BARRIER
[630]	1072	call WaitRequest(Request_Dissip)
[985]	1073	c$OMP MASTER
[630]	1074	call VTe(VThallo)
	1075	call resume_timer(timer_dissip)
[764]	1076	c$OMP END MASTER
[949]	1077	c$OMP BARRIER
[630]	1078	call covcont_p(llm,ucov,vcov,ucont,vcont)
	1079	call enercin_p(vcov,ucov,vcont,ucont,ecin)
	1080
[949]	1081	ijb=ij_begin
	1082	ije=ij_end
	1083	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1084	do l=1,llm
	1085	do ij=ijb,ije
	1086	dtetaecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
	1087	dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
[630]	1088	enddo
[949]	1089	enddo
	1090	c$OMP END DO NOWAIT
[630]	1091	endif
	1092
	1093	ijb=ij_begin
	1094	ije=ij_end
[949]	1095	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1096	do l=1,llm
	1097	do ij=ijb,ije
[630]	1098	teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
[949]	1099	enddo
	1100	enddo
	1101	c$OMP END DO NOWAIT
[630]	1102	c------------------------------------------------------------------------
	1103
	1104
	1105	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
	1106	c ... Calcul de la valeur moyenne, unique de h aux poles .....
	1107	c
	1108
	1109	ijb=ij_begin
	1110	ije=ij_end
	1111
	1112	if (pole_nord) then
[764]	1113	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	1114	DO l = 1, llm
	1115	DO ij = 1,iim
	1116	tppn(ij) = aire( ij ) * teta( ij ,l)
	1117	ENDDO
	1118	tpn = SSUM(iim,tppn,1)/apoln
	1119
	1120	DO ij = 1, iip1
	1121	teta( ij ,l) = tpn
	1122	ENDDO
	1123	ENDDO
[764]	1124	c$OMP END DO NOWAIT
	1125
	1126	c$OMP MASTER
[630]	1127	DO ij = 1,iim
	1128	tppn(ij) = aire( ij ) * ps ( ij )
	1129	ENDDO
	1130	tpn = SSUM(iim,tppn,1)/apoln
	1131
	1132	DO ij = 1, iip1
	1133	ps( ij ) = tpn
	1134	ENDDO
[764]	1135	c$OMP END MASTER
[630]	1136	endif
	1137
	1138	if (pole_sud) then
[764]	1139	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	1140	DO l = 1, llm
	1141	DO ij = 1,iim
	1142	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	1143	ENDDO
	1144	tps = SSUM(iim,tpps,1)/apols
	1145
	1146	DO ij = 1, iip1
	1147	teta(ij+ip1jm,l) = tps
	1148	ENDDO
	1149	ENDDO
[764]	1150	c$OMP END DO NOWAIT
	1151
	1152	c$OMP MASTER
[630]	1153	DO ij = 1,iim
	1154	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
	1155	ENDDO
	1156	tps = SSUM(iim,tpps,1)/apols
	1157
	1158	DO ij = 1, iip1
	1159	ps(ij+ip1jm) = tps
	1160	ENDDO
[764]	1161	c$OMP END MASTER
[630]	1162	endif
	1163
[764]	1164
	1165	c$OMP BARRIER
	1166	c$OMP MASTER
[630]	1167	call VTe(VTdissipation)
	1168
	1169	call stop_timer(timer_dissip)
[949]	1170
[630]	1171	call VTb(VThallo)
[985]	1172	c$OMP END MASTER
[630]	1173	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
	1174	* jj_Nb_caldyn,Request_dissip)
	1175
	1176	call Register_SwapField(vcov,vcov,ip1jm,llm,
	1177	* jj_Nb_caldyn,Request_dissip)
	1178
	1179	call Register_SwapField(teta,teta,ip1jmp1,llm,
	1180	* jj_Nb_caldyn,Request_dissip)
	1181
	1182	call Register_SwapField(p,p,ip1jmp1,llmp1,
	1183	* jj_Nb_caldyn,Request_dissip)
	1184
	1185	call Register_SwapField(pk,pk,ip1jmp1,llm,
	1186	* jj_Nb_caldyn,Request_dissip)
	1187
	1188	call SendRequest(Request_dissip)
[985]	1189	c$OMP BARRIER
[630]	1190	call WaitRequest(Request_dissip)
[985]	1191
	1192	c$OMP BARRIER
	1193	c$OMP MASTER
[630]	1194	call SetDistrib(jj_Nb_caldyn)
	1195	call VTe(VThallo)
	1196	call resume_timer(timer_caldyn)
[961]	1197	c print *,'fin dissipation'
[764]	1198	c$OMP END MASTER
[985]	1199	c$OMP BARRIER
[630]	1200	END IF
	1201
[774]	1202	cc$OMP END PARALLEL
	1203
[630]	1204	c ajout debug
	1205	c IF( lafin ) then
	1206	c abort_message = 'Simulation finished'
	1207	c call abort_gcm(modname,abort_message,0)
	1208	c ENDIF
	1209
	1210	c ********************************************************************
	1211	c ********************************************************************
	1212	c .... fin de l'integration dynamique et physique pour le pas itau ..
	1213	c ********************************************************************
	1214	c ********************************************************************
	1215
	1216	c preparation du pas d'integration suivant ......
	1217	cym call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	1218	cym call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
[774]	1219	c$OMP MASTER
[630]	1220	call stop_timer(timer_caldyn)
[774]	1221	c$OMP END MASTER
[630]	1222	IF (itau==itaumax) then
[774]	1223	c$OMP MASTER
[630]	1224	call allgather_timer_average
	1225
	1226	if (mpi_rank==0) then
	1227
	1228	print ,'********************************'
[949]	1229	print ,'*** TIMER CALDYN ****'
	1230	do i=0,mpi_size-1
	1231	print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
[630]	1232	& ' : temps moyen :',
	1233	& timer_average(jj_nb_caldyn(i),timer_caldyn,i)
	1234	enddo
	1235
	1236	print ,'********************************'
[949]	1237	print ,'*** TIMER VANLEER ****'
	1238	do i=0,mpi_size-1
	1239	print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
[630]	1240	& ' : temps moyen :',
	1241	& timer_average(jj_nb_vanleer(i),timer_vanleer,i)
	1242	enddo
	1243
	1244	print ,'********************************'
[949]	1245	print ,'*** TIMER DISSIP ****'
	1246	do i=0,mpi_size-1
	1247	print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
[630]	1248	& ' : temps moyen :',
	1249	& timer_average(jj_nb_dissip(i),timer_dissip,i)
	1250	enddo
[949]	1251
	1252	print ,'********************************'
	1253	print ,'*** TIMER PHYSIC ****'
	1254	do i=0,mpi_size-1
	1255	print *,'proc',i,' : Nb Bandes :',jj_nb_physic(i),
[630]	1256	& ' : temps moyen :',
	1257	& timer_average(jj_nb_physic(i),timer_physic,i)
	1258	enddo
[949]	1259
[630]	1260	endif
	1261
	1262	print ,'Taille du Buffer MPI (REAL8)',MaxBufferSize
	1263	print ,'Taille du Buffer MPI utilise (REAL8)',MaxBufferSize_Used
	1264	print *, 'Temps total ecoule sur la parallelisation :',DiffTime()
	1265	print *, 'Temps CPU ecoule sur la parallelisation :',DiffCpuTime()
[985]	1266	CALL print_filtre_timer
[1279]	1267	call fin_getparam
[630]	1268	call finalize_parallel
[774]	1269	c$OMP END MASTER
[985]	1270	c$OMP BARRIER
[949]	1271	RETURN
[630]	1272	ENDIF
	1273
	1274	IF ( .NOT.purmats ) THEN
	1275	c ........................................................
	1276	c .............. schema matsuno + leapfrog ..............
	1277	c ........................................................
	1278
	1279	IF(forward. OR. leapf) THEN
	1280	itau= itau + 1
[1279]	1281	! iday= day_ini+itau/day_step
	1282	! time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	1283	! IF(time.GT.1.) THEN
	1284	! time = time-1.
	1285	! iday = iday+1
	1286	! ENDIF
[630]	1287	ENDIF
	1288
	1289
[1286]	1290	IF( itau. EQ. itaufinp1 ) then
[630]	1291
[1286]	1292	if (flag_verif) then
	1293	write(79,*) 'ucov',ucov
	1294	write(80,*) 'vcov',vcov
	1295	write(81,*) 'teta',teta
	1296	write(82,*) 'ps',ps
	1297	write(83,*) 'q',q
	1298	WRITE(85,*) 'q1 = ',q(:,:,1)
	1299	WRITE(86,*) 'q3 = ',q(:,:,3)
	1300	endif
	1301
	1302
[774]	1303	c$OMP MASTER
[1279]	1304	call fin_getparam
[764]	1305	call finalize_parallel
[774]	1306	c$OMP END MASTER
[630]	1307	abort_message = 'Simulation finished'
	1308	call abort_gcm(modname,abort_message,0)
[985]	1309	RETURN
[630]	1310	ENDIF
	1311	c-----------------------------------------------------------------------
	1312	c ecriture du fichier histoire moyenne:
	1313	c -------------------------------------
	1314
	1315	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
[774]	1316	c$OMP BARRIER
[630]	1317	IF(itau.EQ.itaufin) THEN
	1318	iav=1
	1319	ELSE
	1320	iav=0
	1321	ENDIF
	1322	#ifdef CPP_IOIPSL
[1146]	1323	IF (ok_dynzon) THEN
[774]	1324	call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
[985]	1325	call SendRequest(TestRequest)
	1326	c$OMP BARRIER
[630]	1327	call WaitRequest(TestRequest)
[1000]	1328	c$OMP BARRIER
[985]	1329	c$OMP MASTER
[1279]	1330	! CALL writedynav_p(histaveid, itau,vcov ,
	1331	! , ucov,teta,pk,phi,q,masse,ps,phis)
[1146]	1332
	1333	c ATTENTION!!! bilan_dyn_p ne marche probablement pas avec OpenMP
	1334	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
	1335	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[985]	1336	c$OMP END MASTER
[1146]	1337	ENDIF !ok_dynzon
[630]	1338	#endif
	1339	ENDIF
	1340
	1341	c-----------------------------------------------------------------------
	1342	c ecriture de la bande histoire:
	1343	c ------------------------------
	1344
	1345	c IF( MOD(itau,iecri ).EQ.0) THEN
[774]	1346
[1146]	1347	IF( MOD(itau,iecri*day_step).EQ.0) THEN
[774]	1348	c$OMP BARRIER
	1349	c$OMP MASTER
[1146]	1350	nbetat = nbetatdem
	1351	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
[774]	1352
[949]	1353	cym unat=0.
	1354
[1146]	1355	ijb=ij_begin
	1356	ije=ij_end
[949]	1357
[1146]	1358	if (pole_nord) then
	1359	ijb=ij_begin+iip1
	1360	unat(1:iip1,:)=0.
	1361	endif
[949]	1362
[1146]	1363	if (pole_sud) then
	1364	ije=ij_end-iip1
	1365	unat(ij_end-iip1+1:ij_end,:)=0.
	1366	endif
[949]	1367
[1146]	1368	do l=1,llm
	1369	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
	1370	enddo
[630]	1371
[1146]	1372	ijb=ij_begin
	1373	ije=ij_end
	1374	if (pole_sud) ije=ij_end-iip1
[949]	1375
[1146]	1376	do l=1,llm
	1377	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
	1378	enddo
[949]	1379
[630]	1380	#ifdef CPP_IOIPSL
	1381
[1279]	1382	! CALL writehist_p(histid,histvid, itau,vcov,
	1383	! & ucov,teta,phi,q,masse,ps,phis)
[1000]	1384
[630]	1385	#endif
[1146]	1386	! For some Grads outputs of fields
	1387	if (output_grads_dyn) then
	1388	! Ehouarn: hope this works the way I think it does:
	1389	call Gather_Field(unat,ip1jmp1,llm,0)
	1390	call Gather_Field(vnat,ip1jm,llm,0)
	1391	call Gather_Field(teta,ip1jmp1,llm,0)
	1392	call Gather_Field(ps,ip1jmp1,1,0)
	1393	do iq=1,nqtot
	1394	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
	1395	enddo
	1396	if (mpi_rank==0) then
	1397	#include "write_grads_dyn.h"
	1398	endif
	1399	endif ! of if (output_grads_dyn)
[774]	1400	c$OMP END MASTER
[1146]	1401	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
[630]	1402
	1403	IF(itau.EQ.itaufin) THEN
	1404
[774]	1405	c$OMP BARRIER
	1406	c$OMP MASTER
[630]	1407
[1146]	1408	if (planet_type.eq."earth") then
	1409	! Write an Earth-format restart file
	1410	CALL dynredem1_p("restart.nc",0.0,
	1411	& vcov,ucov,teta,q,masse,ps)
	1412	endif ! of if (planet_type.eq."earth")
[630]	1413
[1279]	1414	! CLOSE(99)
[774]	1415	c$OMP END MASTER
[1146]	1416	ENDIF ! of IF (itau.EQ.itaufin)
[630]	1417
	1418	c-----------------------------------------------------------------------
	1419	c gestion de l'integration temporelle:
	1420	c ------------------------------------
	1421
	1422	IF( MOD(itau,iperiod).EQ.0 ) THEN
	1423	GO TO 1
	1424	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	1425
	1426	IF( forward ) THEN
	1427	c fin du pas forward et debut du pas backward
	1428
	1429	forward = .FALSE.
	1430	leapf = .FALSE.
	1431	GO TO 2
	1432
	1433	ELSE
	1434	c fin du pas backward et debut du premier pas leapfrog
	1435
	1436	leapf = .TRUE.
	1437	dt = 2.*dtvr
	1438	GO TO 2
	1439	END IF
	1440	ELSE
	1441
	1442	c ...... pas leapfrog .....
	1443
	1444	leapf = .TRUE.
	1445	dt = 2.*dtvr
	1446	GO TO 2
[1146]	1447	END IF ! of IF (MOD(itau,iperiod).EQ.0)
	1448	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
[630]	1449
	1450
[1146]	1451	ELSE ! of IF (.not.purmats)
	1452
[630]	1453	c ........................................................
	1454	c .............. schema matsuno ...............
	1455	c ........................................................
	1456	IF( forward ) THEN
	1457
	1458	itau = itau + 1
[1279]	1459	! iday = day_ini+itau/day_step
	1460	! time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	1461	!
	1462	! IF(time.GT.1.) THEN
	1463	! time = time-1.
	1464	! iday = iday+1
	1465	! ENDIF
[630]	1466
	1467	forward = .FALSE.
	1468	IF( itau. EQ. itaufinp1 ) then
[774]	1469	c$OMP MASTER
[1279]	1470	call fin_getparam
[949]	1471	call finalize_parallel
[774]	1472	c$OMP END MASTER
[630]	1473	abort_message = 'Simulation finished'
	1474	call abort_gcm(modname,abort_message,0)
[985]	1475	RETURN
[630]	1476	ENDIF
	1477	GO TO 2
	1478
[1146]	1479	ELSE ! of IF(forward)
[630]	1480
[1146]	1481	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
[630]	1482	IF(itau.EQ.itaufin) THEN
	1483	iav=1
	1484	ELSE
	1485	iav=0
	1486	ENDIF
	1487	#ifdef CPP_IOIPSL
[1146]	1488	IF (ok_dynzon) THEN
[774]	1489	c$OMP BARRIER
	1490
[1146]	1491	call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
	1492	call SendRequest(TestRequest)
[985]	1493	c$OMP BARRIER
[1146]	1494	call WaitRequest(TestRequest)
[630]	1495
[1000]	1496	c$OMP BARRIER
[985]	1497	c$OMP MASTER
[1279]	1498	! CALL writedynav_p(histaveid, itau,vcov ,
	1499	! , ucov,teta,pk,phi,q,masse,ps,phis)
[1146]	1500	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
[985]	1501	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[774]	1502	c$OMP END MASTER
[1146]	1503	END IF !ok_dynzon
[630]	1504	#endif
[1146]	1505	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
[630]	1506
[1146]	1507
[630]	1508	c IF(MOD(itau,iecri ).EQ.0) THEN
	1509	IF(MOD(itau,iecri*day_step).EQ.0) THEN
[774]	1510	c$OMP BARRIER
	1511	c$OMP MASTER
[1146]	1512	nbetat = nbetatdem
	1513	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
[630]	1514
	1515	cym unat=0.
[1146]	1516	ijb=ij_begin
	1517	ije=ij_end
[949]	1518
[1146]	1519	if (pole_nord) then
	1520	ijb=ij_begin+iip1
	1521	unat(1:iip1,:)=0.
	1522	endif
[949]	1523
[1146]	1524	if (pole_sud) then
	1525	ije=ij_end-iip1
	1526	unat(ij_end-iip1+1:ij_end,:)=0.
	1527	endif
[949]	1528
[1146]	1529	do l=1,llm
	1530	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
	1531	enddo
[630]	1532
[1146]	1533	ijb=ij_begin
	1534	ije=ij_end
	1535	if (pole_sud) ije=ij_end-iip1
[949]	1536
[1146]	1537	do l=1,llm
	1538	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
	1539	enddo
[630]	1540
	1541	#ifdef CPP_IOIPSL
	1542
[1279]	1543	! CALL writehist_p(histid, histvid, itau,vcov ,
	1544	! & ucov,teta,phi,q,masse,ps,phis)
[1146]	1545	#endif
	1546	! For some Grads output (but does it work?)
	1547	if (output_grads_dyn) then
	1548	call Gather_Field(unat,ip1jmp1,llm,0)
	1549	call Gather_Field(vnat,ip1jm,llm,0)
	1550	call Gather_Field(teta,ip1jmp1,llm,0)
	1551	call Gather_Field(ps,ip1jmp1,1,0)
	1552	do iq=1,nqtot
	1553	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
	1554	enddo
[630]	1555	c
[1146]	1556	if (mpi_rank==0) then
	1557	#include "write_grads_dyn.h"
	1558	endif
	1559	endif ! of if (output_grads_dyn)
[630]	1560
[774]	1561	c$OMP END MASTER
[1146]	1562	ENDIF ! of IF(MOD(itau,iecri*day_step).EQ.0)
[630]	1563
[1146]	1564	IF(itau.EQ.itaufin) THEN
	1565	if (planet_type.eq."earth") then
[774]	1566	c$OMP MASTER
	1567	CALL dynredem1_p("restart.nc",0.0,
[1146]	1568	. vcov,ucov,teta,q,masse,ps)
[774]	1569	c$OMP END MASTER
[1146]	1570	endif ! of if (planet_type.eq."earth")
	1571	ENDIF ! of IF(itau.EQ.itaufin)
[630]	1572
[1146]	1573	forward = .TRUE.
	1574	GO TO 1
[630]	1575
[1146]	1576	ENDIF ! of IF (forward)
	1577
	1578	END IF ! of IF(.not.purmats)
[774]	1579	c$OMP MASTER
[1279]	1580	call fin_getparam
[1146]	1581	call finalize_parallel
[774]	1582	c$OMP END MASTER
[1146]	1583	RETURN
[630]	1584	END

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

Download in other formats:

Original Format