Context Navigation

leapfrog_p.F @ 1322

Last change on this file since 1322 was 1322, checked in by Laurent Fairhead, 14 years ago

Improvements concerning wake parametrisation (from JYG, NR, IT, with more to come).
Alp_offset is read in form physiq.def file

Améliorations à la paramétrisation des poches froides (de JYG, NR, IT, d'autres
sont à venir)
Alp_offset est rajouté à la liste des paramètres lus dans physiq.def

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

File size: 49.4 KB

Rev	Line
[764]	1	!
[1279]	2	! $Id: leapfrog_p.F 1322 2010-03-12 10:54:11Z 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
[1299]	22	USE control_mod
[630]	23
	24	IMPLICIT NONE
	25
	26	c ...... Version du 10/01/98 ..........
	27
	28	c avec coordonnees verticales hybrides
	29	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
	30
	31	c=======================================================================
	32	c
	33	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
	34	c -------
	35	c
	36	c Objet:
	37	c ------
	38	c
	39	c GCM LMD nouvelle grille
	40	c
	41	c=======================================================================
	42	c
	43	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
	44	c et possibilite d'appeler une fonction f(y) a derivee tangente
	45	c hyperbolique a la place de la fonction a derivee sinusoidale.
	46
	47	c ... Possibilite de choisir le shema pour l'advection de
	48	c q , en modifiant iadv dans traceur.def (10/02) .
	49	c
	50	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
	51	c Pour Van-Leer iadv=10
	52	c
	53	c-----------------------------------------------------------------------
	54	c Declarations:
	55	c -------------
	56
	57	#include "dimensions.h"
	58	#include "paramet.h"
	59	#include "comconst.h"
	60	#include "comdissnew.h"
	61	#include "comvert.h"
	62	#include "comgeom.h"
	63	#include "logic.h"
	64	#include "temps.h"
	65	#include "ener.h"
	66	#include "description.h"
	67	#include "serre.h"
	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
[1299]	214	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
[1279]	215	! IF(time.GT.1.) THEN
	216	! time = time-1.
	217	! iday = iday+1
	218	! ENDIF
[630]	219
[1146]	220	c Allocate variables depending on dynamic variable nqtot
	221	c$OMP MASTER
	222	IF (firstcall) THEN
	223	firstcall=.FALSE.
	224	ALLOCATE(dq(ip1jmp1,llm,nqtot))
	225	ALLOCATE(dqfi(ip1jmp1,llm,nqtot))
	226	ALLOCATE(dqfi_tmp(iip1,llm,nqtot))
	227	END IF
	228	c$OMP END MASTER
	229	c$OMP BARRIER
	230
[630]	231	c-----------------------------------------------------------------------
	232	c On initialise la pression et la fonction d'Exner :
	233	c --------------------------------------------------
	234
[774]	235	c$OMP MASTER
[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
[1322]	679	if (prt_level >= 10) then
[1146]	680	write(lunout,*)
	681	& 'leapfrog_p: Entree dans la physique : Iteration No ',true_itau
[1322]	682	endif
[764]	683	c$OMP END MASTER
	684
[630]	685	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
[764]	686
[985]	687	c$OMP BARRIER
[630]	688	CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
[985]	689	c$OMP BARRIER
[1279]	690	jD_cur = jD_ref + day_ini - day_ref
	691	$ + int (itau * dtvr / daysec)
	692	jH_cur = jH_ref + &
	693	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
	694	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
[630]	695
	696	c rajout debug
	697	c lafin = .true.
	698
	699
	700	c Inbterface avec les routines de phylmd (phymars ... )
	701	c -----------------------------------------------------
	702
	703	c+jld
	704
[764]	705	c Diagnostique de conservation de l'energie : initialisation
[630]	706	IF (ip_ebil_dyn.ge.1 ) THEN
	707	ztit='bil dyn'
[1146]	708	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
	709	IF (planet_type.eq."earth") THEN
	710	CALL diagedyn(ztit,2,1,1,dtphys
	711	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	712	ENDIF
[630]	713	ENDIF
	714	c-jld
[764]	715	c$OMP BARRIER
	716	c$OMP MASTER
[630]	717	call VTb(VThallo)
[985]	718	c$OMP END MASTER
	719
[630]	720	call SetTag(Request_physic,800)
[949]	721
[630]	722	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
	723	* jj_Nb_physic,2,2,Request_physic)
	724
	725	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
	726	* jj_Nb_physic,2,2,Request_physic)
	727
	728	call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
	729	* jj_Nb_physic,2,2,Request_physic)
	730
[1279]	731	call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
	732	* jj_Nb_physic,1,2,Request_physic)
	733
[630]	734	call Register_SwapFieldHallo(p,p,ip1jmp1,llmp1,
	735	* jj_Nb_physic,2,2,Request_physic)
	736
	737	call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
	738	* jj_Nb_physic,2,2,Request_physic)
	739
	740	call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
	741	* jj_Nb_physic,2,2,Request_physic)
	742
	743	call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
	744	* jj_Nb_physic,2,2,Request_physic)
	745
	746	call Register_SwapFieldHallo(w,w,ip1jmp1,llm,
	747	* jj_Nb_physic,2,2,Request_physic)
	748
[949]	749	c call SetDistrib(jj_nb_vanleer)
[1146]	750	do j=1,nqtot
[630]	751
	752	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
	753	* jj_Nb_physic,2,2,Request_physic)
	754	enddo
[960]	755
[630]	756	call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
	757	* jj_Nb_physic,2,2,Request_physic)
[949]	758
[630]	759	call SendRequest(Request_Physic)
[985]	760	c$OMP BARRIER
[630]	761	call WaitRequest(Request_Physic)
[985]	762
	763	c$OMP BARRIER
	764	c$OMP MASTER
	765	call SetDistrib(jj_nb_Physic)
[949]	766	call VTe(VThallo)
	767
	768	call VTb(VTphysiq)
[764]	769	c$OMP END MASTER
	770	c$OMP BARRIER
	771
[949]	772	cc$OMP MASTER
[764]	773	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
	774	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
	775	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
	776	c call WriteField_p('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
	777	c call WriteField_p('pkfi',reshape(pk,(/iip1,jmp1,llm/)))
	778	cc$OMP END MASTER
	779	cc$OMP BARRIER
[985]	780	! CALL FTRACE_REGION_BEGIN("calfis")
[1279]	781	CALL calfis_p(lafin ,jD_cur, jH_cur,
[630]	782	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
[960]	783	$ du,dv,dteta,dq,
[630]	784	$ flxw,
	785	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
[985]	786	! CALL FTRACE_REGION_END("calfis")
[630]	787	ijb=ij_begin
	788	ije=ij_end
	789	if ( .not. pole_nord) then
[764]	790	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	791	DO l=1,llm
	792	dufi_tmp(1:iip1,l) = dufi(ijb:ijb+iim,l)
	793	dvfi_tmp(1:iip1,l) = dvfi(ijb:ijb+iim,l)
	794	dtetafi_tmp(1:iip1,l)= dtetafi(ijb:ijb+iim,l)
	795	dqfi_tmp(1:iip1,l,:) = dqfi(ijb:ijb+iim,l,:)
	796	ENDDO
	797	c$OMP END DO NOWAIT
	798
	799	c$OMP MASTER
[630]	800	dpfi_tmp(1:iip1) = dpfi(ijb:ijb+iim)
[764]	801	c$OMP END MASTER
[1146]	802	endif ! of if ( .not. pole_nord)
[764]	803
	804	c$OMP BARRIER
	805	c$OMP MASTER
[949]	806	call SetDistrib(jj_nb_Physic_bis)
[630]	807
	808	call VTb(VThallo)
[985]	809	c$OMP END MASTER
	810	c$OMP BARRIER
[630]	811
	812	call Register_Hallo(dufi,ip1jmp1,llm,
	813	* 1,0,0,1,Request_physic)
	814
	815	call Register_Hallo(dvfi,ip1jm,llm,
	816	* 1,0,0,1,Request_physic)
	817
	818	call Register_Hallo(dtetafi,ip1jmp1,llm,
	819	* 1,0,0,1,Request_physic)
	820
	821	call Register_Hallo(dpfi,ip1jmp1,1,
	822	* 1,0,0,1,Request_physic)
	823
[1146]	824	do j=1,nqtot
[630]	825	call Register_Hallo(dqfi(1,1,j),ip1jmp1,llm,
	826	* 1,0,0,1,Request_physic)
	827	enddo
	828
	829	call SendRequest(Request_Physic)
[985]	830	c$OMP BARRIER
[630]	831	call WaitRequest(Request_Physic)
	832
[985]	833	c$OMP BARRIER
	834	c$OMP MASTER
[630]	835	call VTe(VThallo)
	836
	837	call SetDistrib(jj_nb_Physic)
[764]	838	c$OMP END MASTER
[949]	839	c$OMP BARRIER
	840	ijb=ij_begin
	841	if (.not. pole_nord) then
	842
[764]	843	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	844	DO l=1,llm
	845	dufi(ijb:ijb+iim,l) = dufi(ijb:ijb+iim,l)+dufi_tmp(1:iip1,l)
	846	dvfi(ijb:ijb+iim,l) = dvfi(ijb:ijb+iim,l)+dvfi_tmp(1:iip1,l)
	847	dtetafi(ijb:ijb+iim,l) = dtetafi(ijb:ijb+iim,l)
	848	& +dtetafi_tmp(1:iip1,l)
	849	dqfi(ijb:ijb+iim,l,:) = dqfi(ijb:ijb+iim,l,:)
	850	& + dqfi_tmp(1:iip1,l,:)
	851	ENDDO
	852	c$OMP END DO NOWAIT
	853
	854	c$OMP MASTER
[630]	855	dpfi(ijb:ijb+iim) = dpfi(ijb:ijb+iim)+ dpfi_tmp(1:iip1)
[764]	856	c$OMP END MASTER
[949]	857
[1146]	858	endif ! of if (.not. pole_nord)
[764]	859	c$OMP BARRIER
[949]	860	cc$OMP MASTER
[630]	861	c call WriteField_p('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
	862	c call WriteField_p('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
	863	c call WriteField_p('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
	864	c call WriteField_p('dpfi',reshape(dpfi,(/iip1,jmp1/)))
[764]	865	cc$OMP END MASTER
[630]	866	c
[1146]	867	c do j=1,nqtot
[630]	868	c call WriteField_p('dqfi'//trim(int2str(j)),
	869	c . reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
	870	c enddo
	871
	872	c ajout des tendances physiques:
	873	c ------------------------------
[1000]	874	IF (ok_strato) THEN
[1279]	875	CALL top_bound_p( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
[1000]	876	ENDIF
	877
[1146]	878	CALL addfi_p( dtphys, leapf, forward ,
[630]	879	$ ucov, vcov, teta , q ,ps ,
	880	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	881
[764]	882	c$OMP BARRIER
	883	c$OMP MASTER
[630]	884	call VTe(VTphysiq)
	885
	886	call VTb(VThallo)
[985]	887	c$OMP END MASTER
[630]	888
	889	call SetTag(Request_physic,800)
	890	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
	891	* jj_Nb_caldyn,Request_physic)
	892
	893	call Register_SwapField(vcov,vcov,ip1jm,llm,
	894	* jj_Nb_caldyn,Request_physic)
	895
	896	call Register_SwapField(teta,teta,ip1jmp1,llm,
	897	* jj_Nb_caldyn,Request_physic)
	898
[1279]	899	call Register_SwapField(masse,masse,ip1jmp1,llm,
	900	* jj_Nb_caldyn,Request_physic)
	901
[630]	902	call Register_SwapField(p,p,ip1jmp1,llmp1,
	903	* jj_Nb_caldyn,Request_physic)
	904
	905	call Register_SwapField(pk,pk,ip1jmp1,llm,
	906	* jj_Nb_caldyn,Request_physic)
	907
	908	call Register_SwapField(phis,phis,ip1jmp1,1,
	909	* jj_Nb_caldyn,Request_physic)
	910
	911	call Register_SwapField(phi,phi,ip1jmp1,llm,
	912	* jj_Nb_caldyn,Request_physic)
	913
	914	call Register_SwapField(w,w,ip1jmp1,llm,
	915	* jj_Nb_caldyn,Request_physic)
	916
[1146]	917	do j=1,nqtot
[630]	918
	919	call Register_SwapField(q(1,1,j),q(1,1,j),ip1jmp1,llm,
	920	* jj_Nb_caldyn,Request_physic)
	921
	922	enddo
	923
	924	call SendRequest(Request_Physic)
[985]	925	c$OMP BARRIER
[630]	926	call WaitRequest(Request_Physic)
	927
[985]	928	c$OMP BARRIER
	929	c$OMP MASTER
	930	call VTe(VThallo)
	931	call SetDistrib(jj_Nb_caldyn)
[764]	932	c$OMP END MASTER
	933	c$OMP BARRIER
[630]	934	c
[764]	935	c Diagnostique de conservation de l'energie : difference
[630]	936	IF (ip_ebil_dyn.ge.1 ) THEN
	937	ztit='bil phys'
	938	CALL diagedyn(ztit,2,1,1,dtphys
	939	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	940	ENDIF
[764]	941
	942	cc$OMP MASTER
	943	c if (debug) then
	944	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
	945	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
	946	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
	947	c endif
	948	cc$OMP END MASTER
	949
[630]	950
[1146]	951	c-jld
	952	c$OMP MASTER
	953	call resume_timer(timer_caldyn)
	954	if (FirstPhysic) then
	955	ok_start_timer=.TRUE.
	956	FirstPhysic=.false.
	957	endif
	958	c$OMP END MASTER
	959	ENDIF ! of IF( apphys )
	960
	961	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
[630]	962	c Calcul academique de la physique = Rappel Newtonien + fritcion
	963	c --------------------------------------------------------------
	964	cym teta(:,:)=teta(:,:)
	965	cym s -iphysiqdtvr(teta(:,:)-tetarappel(:,:))/taurappel
	966	ijb=ij_begin
	967	ije=ij_end
	968	teta(ijb:ije,:)=teta(ijb:ije,:)
	969	s -iphysiqdtvr(teta(ijb:ije,:)-tetarappel(ijb:ije,:))/taurappel
	970
	971	call Register_Hallo(ucov,ip1jmp1,llm,0,1,1,0,Request_Physic)
	972	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Physic)
	973	call SendRequest(Request_Physic)
[1279]	974	c$OMP BARRIER
[630]	975	call WaitRequest(Request_Physic)
	976
	977	call friction_p(ucov,vcov,iphysiq*dtvr)
[1146]	978	ENDIF ! of IF(iflag_phys.EQ.2)
[630]	979
	980
[985]	981	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
[774]	982	c$OMP BARRIER
[630]	983	CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
[774]	984	c$OMP BARRIER
[630]	985
[774]	986	cc$OMP END PARALLEL
[630]	987
	988	c-----------------------------------------------------------------------
	989	c dissipation horizontale et verticale des petites echelles:
	990	c ----------------------------------------------------------
	991
	992	IF(apdiss) THEN
[774]	993	cc$OMP PARALLEL DEFAULT(SHARED)
	994	cc$OMP+ PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
[764]	995	c$OMP MASTER
[630]	996	call suspend_timer(timer_caldyn)
[949]	997
	998	c print*,'Entree dans la dissipation : Iteration No ',true_itau
[630]	999	c calcul de l'energie cinetique avant dissipation
[949]	1000	c print *,'Passage dans la dissipation'
[630]	1001
	1002	call VTb(VThallo)
[764]	1003	c$OMP END MASTER
[630]	1004
[764]	1005	c$OMP BARRIER
[985]	1006
[630]	1007	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
	1008	* jj_Nb_dissip,1,1,Request_dissip)
	1009
	1010	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
	1011	* jj_Nb_dissip,1,1,Request_dissip)
	1012
	1013	call Register_SwapField(teta,teta,ip1jmp1,llm,
	1014	* jj_Nb_dissip,Request_dissip)
	1015
	1016	call Register_SwapField(p,p,ip1jmp1,llmp1,
	1017	* jj_Nb_dissip,Request_dissip)
	1018
	1019	call Register_SwapField(pk,pk,ip1jmp1,llm,
	1020	* jj_Nb_dissip,Request_dissip)
	1021
	1022	call SendRequest(Request_dissip)
[985]	1023	c$OMP BARRIER
[630]	1024	call WaitRequest(Request_dissip)
[985]	1025
	1026	c$OMP BARRIER
	1027	c$OMP MASTER
[630]	1028	call SetDistrib(jj_Nb_dissip)
[949]	1029	call VTe(VThallo)
[630]	1030	call VTb(VTdissipation)
[949]	1031	call start_timer(timer_dissip)
[764]	1032	c$OMP END MASTER
	1033	c$OMP BARRIER
	1034
[630]	1035	call covcont_p(llm,ucov,vcov,ucont,vcont)
	1036	call enercin_p(vcov,ucov,vcont,ucont,ecin0)
	1037
	1038	c dissipation
	1039
[985]	1040	! CALL FTRACE_REGION_BEGIN("dissip")
[630]	1041	CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
[985]	1042	! CALL FTRACE_REGION_END("dissip")
[764]	1043
[949]	1044	ijb=ij_begin
	1045	ije=ij_end
	1046	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1047	DO l=1,llm
	1048	ucov(ijb:ije,l)=ucov(ijb:ije,l)+dudis(ijb:ije,l)
[764]	1049	ENDDO
[949]	1050	c$OMP END DO NOWAIT
	1051	if (pole_sud) ije=ije-iip1
	1052	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1053	DO l=1,llm
	1054	vcov(ijb:ije,l)=vcov(ijb:ije,l)+dvdis(ijb:ije,l)
[764]	1055	ENDDO
[949]	1056	c$OMP END DO NOWAIT
[764]	1057
[630]	1058	c teta=teta+dtetadis
	1059
	1060
	1061	c------------------------------------------------------------------------
	1062	if (dissip_conservative) then
	1063	C On rajoute la tendance due a la transform. Ec -> E therm. cree
	1064	C lors de la dissipation
[764]	1065	c$OMP BARRIER
	1066	c$OMP MASTER
[630]	1067	call suspend_timer(timer_dissip)
[949]	1068	call VTb(VThallo)
[985]	1069	c$OMP END MASTER
[630]	1070	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Request_Dissip)
	1071	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Dissip)
	1072	call SendRequest(Request_Dissip)
[985]	1073	c$OMP BARRIER
[630]	1074	call WaitRequest(Request_Dissip)
[985]	1075	c$OMP MASTER
[630]	1076	call VTe(VThallo)
	1077	call resume_timer(timer_dissip)
[764]	1078	c$OMP END MASTER
[949]	1079	c$OMP BARRIER
[630]	1080	call covcont_p(llm,ucov,vcov,ucont,vcont)
	1081	call enercin_p(vcov,ucov,vcont,ucont,ecin)
	1082
[949]	1083	ijb=ij_begin
	1084	ije=ij_end
	1085	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1086	do l=1,llm
	1087	do ij=ijb,ije
	1088	dtetaecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
	1089	dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
[630]	1090	enddo
[949]	1091	enddo
	1092	c$OMP END DO NOWAIT
[630]	1093	endif
	1094
	1095	ijb=ij_begin
	1096	ije=ij_end
[949]	1097	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1098	do l=1,llm
	1099	do ij=ijb,ije
[630]	1100	teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
[949]	1101	enddo
	1102	enddo
	1103	c$OMP END DO NOWAIT
[630]	1104	c------------------------------------------------------------------------
	1105
	1106
	1107	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
	1108	c ... Calcul de la valeur moyenne, unique de h aux poles .....
	1109	c
	1110
	1111	ijb=ij_begin
	1112	ije=ij_end
	1113
	1114	if (pole_nord) then
[764]	1115	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	1116	DO l = 1, llm
	1117	DO ij = 1,iim
	1118	tppn(ij) = aire( ij ) * teta( ij ,l)
	1119	ENDDO
	1120	tpn = SSUM(iim,tppn,1)/apoln
	1121
	1122	DO ij = 1, iip1
	1123	teta( ij ,l) = tpn
	1124	ENDDO
	1125	ENDDO
[764]	1126	c$OMP END DO NOWAIT
	1127
	1128	c$OMP MASTER
[630]	1129	DO ij = 1,iim
	1130	tppn(ij) = aire( ij ) * ps ( ij )
	1131	ENDDO
	1132	tpn = SSUM(iim,tppn,1)/apoln
	1133
	1134	DO ij = 1, iip1
	1135	ps( ij ) = tpn
	1136	ENDDO
[764]	1137	c$OMP END MASTER
[630]	1138	endif
	1139
	1140	if (pole_sud) then
[764]	1141	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	1142	DO l = 1, llm
	1143	DO ij = 1,iim
	1144	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	1145	ENDDO
	1146	tps = SSUM(iim,tpps,1)/apols
	1147
	1148	DO ij = 1, iip1
	1149	teta(ij+ip1jm,l) = tps
	1150	ENDDO
	1151	ENDDO
[764]	1152	c$OMP END DO NOWAIT
	1153
	1154	c$OMP MASTER
[630]	1155	DO ij = 1,iim
	1156	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
	1157	ENDDO
	1158	tps = SSUM(iim,tpps,1)/apols
	1159
	1160	DO ij = 1, iip1
	1161	ps(ij+ip1jm) = tps
	1162	ENDDO
[764]	1163	c$OMP END MASTER
[630]	1164	endif
	1165
[764]	1166
	1167	c$OMP BARRIER
	1168	c$OMP MASTER
[630]	1169	call VTe(VTdissipation)
	1170
	1171	call stop_timer(timer_dissip)
[949]	1172
[630]	1173	call VTb(VThallo)
[985]	1174	c$OMP END MASTER
[630]	1175	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
	1176	* jj_Nb_caldyn,Request_dissip)
	1177
	1178	call Register_SwapField(vcov,vcov,ip1jm,llm,
	1179	* jj_Nb_caldyn,Request_dissip)
	1180
	1181	call Register_SwapField(teta,teta,ip1jmp1,llm,
	1182	* jj_Nb_caldyn,Request_dissip)
	1183
	1184	call Register_SwapField(p,p,ip1jmp1,llmp1,
	1185	* jj_Nb_caldyn,Request_dissip)
	1186
	1187	call Register_SwapField(pk,pk,ip1jmp1,llm,
	1188	* jj_Nb_caldyn,Request_dissip)
	1189
	1190	call SendRequest(Request_dissip)
[985]	1191	c$OMP BARRIER
[630]	1192	call WaitRequest(Request_dissip)
[985]	1193
	1194	c$OMP BARRIER
	1195	c$OMP MASTER
[630]	1196	call SetDistrib(jj_Nb_caldyn)
	1197	call VTe(VThallo)
	1198	call resume_timer(timer_caldyn)
[961]	1199	c print *,'fin dissipation'
[764]	1200	c$OMP END MASTER
[985]	1201	c$OMP BARRIER
[630]	1202	END IF
	1203
[774]	1204	cc$OMP END PARALLEL
	1205
[630]	1206	c ajout debug
	1207	c IF( lafin ) then
	1208	c abort_message = 'Simulation finished'
	1209	c call abort_gcm(modname,abort_message,0)
	1210	c ENDIF
	1211
	1212	c ********************************************************************
	1213	c ********************************************************************
	1214	c .... fin de l'integration dynamique et physique pour le pas itau ..
	1215	c ********************************************************************
	1216	c ********************************************************************
	1217
	1218	c preparation du pas d'integration suivant ......
	1219	cym call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	1220	cym call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
[774]	1221	c$OMP MASTER
[630]	1222	call stop_timer(timer_caldyn)
[774]	1223	c$OMP END MASTER
[630]	1224	IF (itau==itaumax) then
[774]	1225	c$OMP MASTER
[630]	1226	call allgather_timer_average
	1227
	1228	if (mpi_rank==0) then
	1229
	1230	print ,'********************************'
[949]	1231	print ,'*** TIMER CALDYN ****'
	1232	do i=0,mpi_size-1
	1233	print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
[630]	1234	& ' : temps moyen :',
	1235	& timer_average(jj_nb_caldyn(i),timer_caldyn,i)
	1236	enddo
	1237
	1238	print ,'********************************'
[949]	1239	print ,'*** TIMER VANLEER ****'
	1240	do i=0,mpi_size-1
	1241	print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
[630]	1242	& ' : temps moyen :',
	1243	& timer_average(jj_nb_vanleer(i),timer_vanleer,i)
	1244	enddo
	1245
	1246	print ,'********************************'
[949]	1247	print ,'*** TIMER DISSIP ****'
	1248	do i=0,mpi_size-1
	1249	print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
[630]	1250	& ' : temps moyen :',
	1251	& timer_average(jj_nb_dissip(i),timer_dissip,i)
	1252	enddo
[949]	1253
	1254	print ,'********************************'
	1255	print ,'*** TIMER PHYSIC ****'
	1256	do i=0,mpi_size-1
	1257	print *,'proc',i,' : Nb Bandes :',jj_nb_physic(i),
[630]	1258	& ' : temps moyen :',
	1259	& timer_average(jj_nb_physic(i),timer_physic,i)
	1260	enddo
[949]	1261
[630]	1262	endif
	1263
	1264	print ,'Taille du Buffer MPI (REAL8)',MaxBufferSize
	1265	print ,'Taille du Buffer MPI utilise (REAL8)',MaxBufferSize_Used
	1266	print *, 'Temps total ecoule sur la parallelisation :',DiffTime()
	1267	print *, 'Temps CPU ecoule sur la parallelisation :',DiffCpuTime()
[985]	1268	CALL print_filtre_timer
[1279]	1269	call fin_getparam
[630]	1270	call finalize_parallel
[774]	1271	c$OMP END MASTER
[985]	1272	c$OMP BARRIER
[949]	1273	RETURN
[630]	1274	ENDIF
	1275
	1276	IF ( .NOT.purmats ) THEN
	1277	c ........................................................
	1278	c .............. schema matsuno + leapfrog ..............
	1279	c ........................................................
	1280
	1281	IF(forward. OR. leapf) THEN
	1282	itau= itau + 1
[1279]	1283	! iday= day_ini+itau/day_step
[1299]	1284	! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
[1279]	1285	! IF(time.GT.1.) THEN
	1286	! time = time-1.
	1287	! iday = iday+1
	1288	! ENDIF
[630]	1289	ENDIF
	1290
	1291
[1286]	1292	IF( itau. EQ. itaufinp1 ) then
[630]	1293
[1286]	1294	if (flag_verif) then
	1295	write(79,*) 'ucov',ucov
	1296	write(80,*) 'vcov',vcov
	1297	write(81,*) 'teta',teta
	1298	write(82,*) 'ps',ps
	1299	write(83,*) 'q',q
	1300	WRITE(85,*) 'q1 = ',q(:,:,1)
	1301	WRITE(86,*) 'q3 = ',q(:,:,3)
	1302	endif
	1303
	1304
[774]	1305	c$OMP MASTER
[1279]	1306	call fin_getparam
[764]	1307	call finalize_parallel
[774]	1308	c$OMP END MASTER
[630]	1309	abort_message = 'Simulation finished'
	1310	call abort_gcm(modname,abort_message,0)
[985]	1311	RETURN
[630]	1312	ENDIF
	1313	c-----------------------------------------------------------------------
	1314	c ecriture du fichier histoire moyenne:
	1315	c -------------------------------------
	1316
	1317	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
[774]	1318	c$OMP BARRIER
[630]	1319	IF(itau.EQ.itaufin) THEN
	1320	iav=1
	1321	ELSE
	1322	iav=0
	1323	ENDIF
	1324	#ifdef CPP_IOIPSL
[1146]	1325	IF (ok_dynzon) THEN
[774]	1326	call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
[985]	1327	call SendRequest(TestRequest)
	1328	c$OMP BARRIER
[630]	1329	call WaitRequest(TestRequest)
[1000]	1330	c$OMP BARRIER
[985]	1331	c$OMP MASTER
[1279]	1332	! CALL writedynav_p(histaveid, itau,vcov ,
	1333	! , ucov,teta,pk,phi,q,masse,ps,phis)
[1146]	1334
	1335	c ATTENTION!!! bilan_dyn_p ne marche probablement pas avec OpenMP
	1336	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
	1337	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[985]	1338	c$OMP END MASTER
[1146]	1339	ENDIF !ok_dynzon
[630]	1340	#endif
	1341	ENDIF
	1342
	1343	c-----------------------------------------------------------------------
	1344	c ecriture de la bande histoire:
	1345	c ------------------------------
	1346
	1347	c IF( MOD(itau,iecri ).EQ.0) THEN
[774]	1348
[1146]	1349	IF( MOD(itau,iecri*day_step).EQ.0) THEN
[774]	1350	c$OMP BARRIER
	1351	c$OMP MASTER
[1146]	1352	nbetat = nbetatdem
	1353	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
[774]	1354
[949]	1355	cym unat=0.
	1356
[1146]	1357	ijb=ij_begin
	1358	ije=ij_end
[949]	1359
[1146]	1360	if (pole_nord) then
	1361	ijb=ij_begin+iip1
	1362	unat(1:iip1,:)=0.
	1363	endif
[949]	1364
[1146]	1365	if (pole_sud) then
	1366	ije=ij_end-iip1
	1367	unat(ij_end-iip1+1:ij_end,:)=0.
	1368	endif
[949]	1369
[1146]	1370	do l=1,llm
	1371	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
	1372	enddo
[630]	1373
[1146]	1374	ijb=ij_begin
	1375	ije=ij_end
	1376	if (pole_sud) ije=ij_end-iip1
[949]	1377
[1146]	1378	do l=1,llm
	1379	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
	1380	enddo
[949]	1381
[630]	1382	#ifdef CPP_IOIPSL
	1383
[1279]	1384	! CALL writehist_p(histid,histvid, itau,vcov,
	1385	! & ucov,teta,phi,q,masse,ps,phis)
[1000]	1386
[630]	1387	#endif
[1146]	1388	! For some Grads outputs of fields
	1389	if (output_grads_dyn) then
	1390	! Ehouarn: hope this works the way I think it does:
	1391	call Gather_Field(unat,ip1jmp1,llm,0)
	1392	call Gather_Field(vnat,ip1jm,llm,0)
	1393	call Gather_Field(teta,ip1jmp1,llm,0)
	1394	call Gather_Field(ps,ip1jmp1,1,0)
	1395	do iq=1,nqtot
	1396	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
	1397	enddo
	1398	if (mpi_rank==0) then
	1399	#include "write_grads_dyn.h"
	1400	endif
	1401	endif ! of if (output_grads_dyn)
[774]	1402	c$OMP END MASTER
[1146]	1403	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
[630]	1404
	1405	IF(itau.EQ.itaufin) THEN
	1406
[774]	1407	c$OMP BARRIER
	1408	c$OMP MASTER
[630]	1409
[1146]	1410	if (planet_type.eq."earth") then
	1411	! Write an Earth-format restart file
	1412	CALL dynredem1_p("restart.nc",0.0,
	1413	& vcov,ucov,teta,q,masse,ps)
	1414	endif ! of if (planet_type.eq."earth")
[630]	1415
[1279]	1416	! CLOSE(99)
[774]	1417	c$OMP END MASTER
[1146]	1418	ENDIF ! of IF (itau.EQ.itaufin)
[630]	1419
	1420	c-----------------------------------------------------------------------
	1421	c gestion de l'integration temporelle:
	1422	c ------------------------------------
	1423
	1424	IF( MOD(itau,iperiod).EQ.0 ) THEN
	1425	GO TO 1
	1426	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	1427
	1428	IF( forward ) THEN
	1429	c fin du pas forward et debut du pas backward
	1430
	1431	forward = .FALSE.
	1432	leapf = .FALSE.
	1433	GO TO 2
	1434
	1435	ELSE
	1436	c fin du pas backward et debut du premier pas leapfrog
	1437
	1438	leapf = .TRUE.
	1439	dt = 2.*dtvr
	1440	GO TO 2
	1441	END IF
	1442	ELSE
	1443
	1444	c ...... pas leapfrog .....
	1445
	1446	leapf = .TRUE.
	1447	dt = 2.*dtvr
	1448	GO TO 2
[1146]	1449	END IF ! of IF (MOD(itau,iperiod).EQ.0)
	1450	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
[630]	1451
	1452
[1146]	1453	ELSE ! of IF (.not.purmats)
	1454
[630]	1455	c ........................................................
	1456	c .............. schema matsuno ...............
	1457	c ........................................................
	1458	IF( forward ) THEN
	1459
	1460	itau = itau + 1
[1279]	1461	! iday = day_ini+itau/day_step
[1299]	1462	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
[1279]	1463	!
	1464	! IF(time.GT.1.) THEN
	1465	! time = time-1.
	1466	! iday = iday+1
	1467	! ENDIF
[630]	1468
	1469	forward = .FALSE.
	1470	IF( itau. EQ. itaufinp1 ) then
[774]	1471	c$OMP MASTER
[1279]	1472	call fin_getparam
[949]	1473	call finalize_parallel
[774]	1474	c$OMP END MASTER
[630]	1475	abort_message = 'Simulation finished'
	1476	call abort_gcm(modname,abort_message,0)
[985]	1477	RETURN
[630]	1478	ENDIF
	1479	GO TO 2
	1480
[1146]	1481	ELSE ! of IF(forward)
[630]	1482
[1146]	1483	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
[630]	1484	IF(itau.EQ.itaufin) THEN
	1485	iav=1
	1486	ELSE
	1487	iav=0
	1488	ENDIF
	1489	#ifdef CPP_IOIPSL
[1146]	1490	IF (ok_dynzon) THEN
[774]	1491	c$OMP BARRIER
	1492
[1146]	1493	call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
	1494	call SendRequest(TestRequest)
[985]	1495	c$OMP BARRIER
[1146]	1496	call WaitRequest(TestRequest)
[630]	1497
[1000]	1498	c$OMP BARRIER
[985]	1499	c$OMP MASTER
[1279]	1500	! CALL writedynav_p(histaveid, itau,vcov ,
	1501	! , ucov,teta,pk,phi,q,masse,ps,phis)
[1146]	1502	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
[985]	1503	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[774]	1504	c$OMP END MASTER
[1146]	1505	END IF !ok_dynzon
[630]	1506	#endif
[1146]	1507	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
[630]	1508
[1146]	1509
[630]	1510	c IF(MOD(itau,iecri ).EQ.0) THEN
	1511	IF(MOD(itau,iecri*day_step).EQ.0) THEN
[774]	1512	c$OMP BARRIER
	1513	c$OMP MASTER
[1146]	1514	nbetat = nbetatdem
	1515	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
[630]	1516
	1517	cym unat=0.
[1146]	1518	ijb=ij_begin
	1519	ije=ij_end
[949]	1520
[1146]	1521	if (pole_nord) then
	1522	ijb=ij_begin+iip1
	1523	unat(1:iip1,:)=0.
	1524	endif
[949]	1525
[1146]	1526	if (pole_sud) then
	1527	ije=ij_end-iip1
	1528	unat(ij_end-iip1+1:ij_end,:)=0.
	1529	endif
[949]	1530
[1146]	1531	do l=1,llm
	1532	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
	1533	enddo
[630]	1534
[1146]	1535	ijb=ij_begin
	1536	ije=ij_end
	1537	if (pole_sud) ije=ij_end-iip1
[949]	1538
[1146]	1539	do l=1,llm
	1540	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
	1541	enddo
[630]	1542
	1543	#ifdef CPP_IOIPSL
	1544
[1279]	1545	! CALL writehist_p(histid, histvid, itau,vcov ,
	1546	! & ucov,teta,phi,q,masse,ps,phis)
[1146]	1547	#endif
	1548	! For some Grads output (but does it work?)
	1549	if (output_grads_dyn) then
	1550	call Gather_Field(unat,ip1jmp1,llm,0)
	1551	call Gather_Field(vnat,ip1jm,llm,0)
	1552	call Gather_Field(teta,ip1jmp1,llm,0)
	1553	call Gather_Field(ps,ip1jmp1,1,0)
	1554	do iq=1,nqtot
	1555	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
	1556	enddo
[630]	1557	c
[1146]	1558	if (mpi_rank==0) then
	1559	#include "write_grads_dyn.h"
	1560	endif
	1561	endif ! of if (output_grads_dyn)
[630]	1562
[774]	1563	c$OMP END MASTER
[1146]	1564	ENDIF ! of IF(MOD(itau,iecri*day_step).EQ.0)
[630]	1565
[1146]	1566	IF(itau.EQ.itaufin) THEN
	1567	if (planet_type.eq."earth") then
[774]	1568	c$OMP MASTER
	1569	CALL dynredem1_p("restart.nc",0.0,
[1146]	1570	. vcov,ucov,teta,q,masse,ps)
[774]	1571	c$OMP END MASTER
[1146]	1572	endif ! of if (planet_type.eq."earth")
	1573	ENDIF ! of IF(itau.EQ.itaufin)
[630]	1574
[1146]	1575	forward = .TRUE.
	1576	GO TO 1
[630]	1577
[1146]	1578	ENDIF ! of IF (forward)
	1579
	1580	END IF ! of IF(.not.purmats)
[774]	1581	c$OMP MASTER
[1279]	1582	call fin_getparam
[1146]	1583	call finalize_parallel
[774]	1584	c$OMP END MASTER
[1146]	1585	RETURN
[630]	1586	END

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

Download in other formats:

Original Format