Context Navigation

leapfrog.F @ 1578

Last change on this file since 1578 was 1577, checked in by Laurent Fairhead, 13 years ago

Modifications au code qui permettent de commencer une simulation à n'importe
quelle heure de la journée. On fait toujours un nombre entier de jours de
simulation.
On spécifie cette heure de départ dans la variable starttime du run.def (la
valeur est en jour et elle est à zéro par défaut).
La valeur est sauvegardée dans le fichier restart.nc. Les valeurs lues dans
le fichier start et le run.def sont comparées en début de simulation. La
simulation s'arrête si elles ne sont pas égales sauf si une remise à zéro de
la date a été demandée.
Par ailleurs, la fréquence de lecture des conditions aux limites a été modifiée
pour qu'à chaque changement de jour, celles-ci soient mises à jour (jusqu'à
maintenant elles étaient mises à jour à une fréquence donnée qui, en cas de
départ de simulation à une heure différente de minuit, ne correspondait pas
forcèment à un changement dans la date).
Validation effectuée en traçant le flux solaire descendant au sommet de
l'atmosphère à différentes heures de la journée, après un redémarrage, en
s'assurant que le maximum est bien là où il est sensé être.

Modifications to the code to enable it to be started at any time of the day.
The code still runs for an integer number of days.
The start time is specified using variable starttime in the run.def file (the
value is in days and is zero by default).
The start time is saved in the restart.nc file at the end of the simulation.
The values read in from the start.nc file and the run.def file are compared
at the start of the simulation. If they differ, the simulation is aborted
unless the raz_date variable has been set.
Furthermore, the frequency at which boundary conditions are read in has been
modified so that they are updated everyday at midnight (until now, they were
updated at a certain frequency that, in case of a simulation starting at a time
other than midnight, did not ensure that those conditions would be updated each
day at midnight)
The modifications were validated by plotting the downward solaf flux at TOA at
different times of the day (and after having restarted the simulation) and
ensuring that the maximum of flux was at the right place according to local
time.

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

File size: 24.8 KB

Rev	Line
[524]	1	!
[1279]	2	! $Id: leapfrog.F 1577 2011-10-20 15:06:47Z jghattas $
	3	!
[524]	4	c
	5	c
[1146]	6	SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
[559]	7	& time_0)
[524]	8
	9
[692]	10	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
[1146]	11	#ifdef CPP_IOIPSL
	12	use IOIPSL
	13	#endif
	14	USE infotrac
[1279]	15	USE guide_mod, ONLY : guide_main
	16	USE write_field
[1403]	17	USE control_mod
[524]	18	IMPLICIT NONE
	19
	20	c ...... Version du 10/01/98 ..........
	21
	22	c avec coordonnees verticales hybrides
	23	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
	24
	25	c=======================================================================
	26	c
	27	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
	28	c -------
	29	c
	30	c Objet:
	31	c ------
	32	c
	33	c GCM LMD nouvelle grille
	34	c
	35	c=======================================================================
	36	c
	37	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
	38	c et possibilite d'appeler une fonction f(y) a derivee tangente
	39	c hyperbolique a la place de la fonction a derivee sinusoidale.
	40
	41	c ... Possibilite de choisir le shema pour l'advection de
	42	c q , en modifiant iadv dans traceur.def (10/02) .
	43	c
	44	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
	45	c Pour Van-Leer iadv=10
	46	c
	47	c-----------------------------------------------------------------------
	48	c Declarations:
	49	c -------------
	50
	51	#include "dimensions.h"
	52	#include "paramet.h"
	53	#include "comconst.h"
	54	#include "comdissnew.h"
	55	#include "comvert.h"
	56	#include "comgeom.h"
	57	#include "logic.h"
	58	#include "temps.h"
	59	#include "ener.h"
	60	#include "description.h"
	61	#include "serre.h"
[1403]	62	!#include "com_io_dyn.h"
[524]	63	#include "iniprint.h"
	64	#include "academic.h"
	65
[956]	66	! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
	67	! #include "clesphys.h"
[692]	68
[524]	69	INTEGER longcles
	70	PARAMETER ( longcles = 20 )
	71	REAL clesphy0( longcles )
	72
	73	real zqmin,zqmax
	74	INTEGER nbetatmoy, nbetatdem,nbetat
	75
	76	c variables dynamiques
	77	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
	78	REAL teta(ip1jmp1,llm) ! temperature potentielle
[1146]	79	REAL q(ip1jmp1,llm,nqtot) ! champs advectes
[524]	80	REAL ps(ip1jmp1) ! pression au sol
	81	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
	82	REAL pks(ip1jmp1) ! exner au sol
	83	REAL pk(ip1jmp1,llm) ! exner au milieu des couches
	84	REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
	85	REAL masse(ip1jmp1,llm) ! masse d'air
	86	REAL phis(ip1jmp1) ! geopotentiel au sol
	87	REAL phi(ip1jmp1,llm) ! geopotentiel
	88	REAL w(ip1jmp1,llm) ! vitesse verticale
	89
	90	c variables dynamiques intermediaire pour le transport
	91	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
	92
	93	c variables dynamiques au pas -1
	94	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
	95	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
	96	REAL massem1(ip1jmp1,llm)
	97
	98	c tendances dynamiques
	99	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
[1146]	100	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
[524]	101
	102	c tendances de la dissipation
	103	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
	104	REAL dtetadis(ip1jmp1,llm)
	105
	106	c tendances physiques
	107	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
[1146]	108	REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
[524]	109
	110	c variables pour le fichier histoire
	111	REAL dtav ! intervalle de temps elementaire
	112
	113	REAL tppn(iim),tpps(iim),tpn,tps
	114	c
	115	INTEGER itau,itaufinp1,iav
[1279]	116	! INTEGER iday ! jour julien
	117	REAL time
[524]	118
	119	REAL SSUM
	120	REAL time_0 , finvmaold(ip1jmp1,llm)
	121
[566]	122	cym LOGICAL lafin
	123	LOGICAL :: lafin=.false.
[524]	124	INTEGER ij,iq,l
	125	INTEGER ik
	126
	127	real time_step, t_wrt, t_ops
	128
[1279]	129	! REAL rdayvrai,rdaym_ini
	130	! jD_cur: jour julien courant
	131	! jH_cur: heure julienne courante
	132	REAL :: jD_cur, jH_cur
	133	INTEGER :: an, mois, jour
	134	REAL :: secondes
	135
[524]	136	LOGICAL first,callinigrads
[692]	137	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
	138	save first
	139	data first/.true./
[1279]	140	real dt_cum
[692]	141	character*10 infile
	142	integer zan, tau0, thoriid
	143	integer nid_ctesGCM
	144	save nid_ctesGCM
	145	real degres
	146	real rlong(iip1), rlatg(jjp1)
	147	real zx_tmp_2d(iip1,jjp1)
	148	integer ndex2d(iip1*jjp1)
	149	logical ok_sync
	150	parameter (ok_sync = .true.)
[1529]	151	logical physic
[524]	152
	153	data callinigrads/.true./
	154	character*10 string10
	155
	156	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
[960]	157	REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
[524]	158
	159	c+jld variables test conservation energie
	160	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
	161	C Tendance de la temp. potentiel d (theta)/ d t due a la
	162	C tansformation d'energie cinetique en energie thermique
	163	C cree par la dissipation
	164	REAL dtetaecdt(ip1jmp1,llm)
	165	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
	166	REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
	167	REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
	168	CHARACTER*15 ztit
[692]	169	!IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
	170	!IM SAVE ip_ebil_dyn
	171	!IM DATA ip_ebil_dyn/0/
[524]	172	c-jld
	173
	174	character*80 dynhist_file, dynhistave_file
[1520]	175	character(len=*),parameter :: modname="leapfrog"
[524]	176	character*80 abort_message
	177
	178	logical dissip_conservative
	179	save dissip_conservative
	180	data dissip_conservative/.true./
	181
	182	LOGICAL prem
	183	save prem
	184	DATA prem/.true./
	185	INTEGER testita
	186	PARAMETER (testita = 9)
	187
[1286]	188	logical , parameter :: flag_verif = .false.
[999]	189
	190
[825]	191	integer itau_w ! pas de temps ecriture = itap + itau_phy
	192
	193
[524]	194	itaufin = nday*day_step
	195	itaufinp1 = itaufin +1
[1529]	196	itau = 0
	197	physic=.true.
	198	if (iflag_phys==0.or.iflag_phys==2) physic=.false.
[524]	199
[1403]	200	c iday = day_ini+itau/day_step
	201	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	202	c IF(time.GT.1.) THEN
	203	c time = time-1.
	204	c iday = iday+1
	205	c ENDIF
[524]	206
	207
	208	c-----------------------------------------------------------------------
	209	c On initialise la pression et la fonction d'Exner :
	210	c --------------------------------------------------
	211
[1454]	212	dq(:,:,:)=0.
[524]	213	CALL pression ( ip1jmp1, ap, bp, ps, p )
[1520]	214	if (disvert_type==1) then
	215	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	216	else ! we assume that we are in the disvert_type==2 case
	217	CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
	218	endif
[524]	219
	220	c-----------------------------------------------------------------------
	221	c Debut de l'integration temporelle:
	222	c ----------------------------------
	223
	224	1 CONTINUE
	225
[1577]	226	jD_cur = jD_ref + day_ini - day_ref + &
	227	& int (itau * dtvr / daysec)
	228	jH_cur = jH_ref + start_time + &
[1279]	229	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
[1577]	230	jD_cur = jD_cur + int(jH_cur)
	231	jH_cur = jH_cur - int(jH_cur)
[524]	232
[1279]	233
[524]	234	#ifdef CPP_IOIPSL
[1279]	235	if (ok_guide) then
	236	call guide_main(itau,ucov,vcov,teta,q,masse,ps)
[524]	237	endif
	238	#endif
[1279]	239
	240
[524]	241	c
	242	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
	243	c CALL test_period ( ucov,vcov,teta,q,p,phis )
	244	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
	245	c ENDIF
	246	c
[1146]	247
	248	! Save fields obtained at previous time step as '...m1'
[524]	249	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
	250	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
	251	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
	252	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
	253	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
	254
	255	forward = .TRUE.
	256	leapf = .FALSE.
	257	dt = dtvr
	258
	259	c ... P.Le Van .26/04/94 ....
	260
	261	CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
	262	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
	263
	264	2 CONTINUE
	265
	266	c-----------------------------------------------------------------------
	267
	268	c date:
	269	c -----
	270
	271
	272	c gestion des appels de la physique et des dissipations:
	273	c ------------------------------------------------------
	274	c
	275	c ... P.Le Van ( 6/02/95 ) ....
	276
	277	apphys = .FALSE.
	278	statcl = .FALSE.
	279	conser = .FALSE.
	280	apdiss = .FALSE.
	281
	282	IF( purmats ) THEN
[1403]	283	! Purely Matsuno time stepping
[524]	284	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
[1502]	285	IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward )
	286	s apdiss = .TRUE.
[524]	287	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
[1529]	288	s .and. physic ) apphys = .TRUE.
[524]	289	ELSE
[1403]	290	! Leapfrog/Matsuno time stepping
[524]	291	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
[1502]	292	IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
	293	s apdiss = .TRUE.
[1529]	294	IF( MOD(itau+1,iphysiq).EQ.0.AND.physic ) apphys=.TRUE.
[524]	295	END IF
	296
[1403]	297	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
	298	! supress dissipation step
	299	if (llm.eq.1) then
	300	apdiss=.false.
	301	endif
	302
[524]	303	c-----------------------------------------------------------------------
	304	c calcul des tendances dynamiques:
	305	c --------------------------------
	306
	307	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	308
[1279]	309	time = jD_cur + jH_cur
[524]	310	CALL caldyn
	311	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
[1279]	312	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
[524]	313
[1279]	314
[524]	315	c-----------------------------------------------------------------------
	316	c calcul des tendances advection des traceurs (dont l'humidite)
	317	c -------------------------------------------------------------
	318
	319	IF( forward. OR . leapf ) THEN
	320
[960]	321	CALL caladvtrac(q,pbaru,pbarv,
	322	* p, masse, dq, teta,
	323	. flxw, pk)
	324
[524]	325	IF (offline) THEN
	326	Cmaf stokage du flux de masse pour traceurs OFF-LINE
	327
	328	#ifdef CPP_IOIPSL
[541]	329	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
	330	. dtvr, itau)
[524]	331	#endif
	332
	333
[1146]	334	ENDIF ! of IF (offline)
[524]	335	c
[1146]	336	ENDIF ! of IF( forward. OR . leapf )
[524]	337
	338
	339	c-----------------------------------------------------------------------
	340	c integrations dynamique et traceurs:
	341	c ----------------------------------
	342
	343
	344	CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
	345	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
	346	$ finvmaold )
	347
	348
	349	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
	350	c
	351	c-----------------------------------------------------------------------
	352	c calcul des tendances physiques:
	353	c -------------------------------
	354	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
	355	c
	356	IF( purmats ) THEN
	357	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
	358	ELSE
	359	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
	360	ENDIF
	361	c
	362	c
	363	IF( apphys ) THEN
	364	c
	365	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
	366	c
	367
	368	CALL pression ( ip1jmp1, ap, bp, ps, p )
[1520]	369	if (disvert_type==1) then
	370	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
	371	else ! we assume that we are in the disvert_type==2 case
	372	CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
	373	endif
[524]	374
[1279]	375	! rdaym_ini = itau * dtvr / daysec
	376	! rdayvrai = rdaym_ini + day_ini
[1577]	377	! jD_cur = jD_ref + day_ini - day_ref
	378	! $ + int (itau * dtvr / daysec)
	379	! jH_cur = jH_ref + &
	380	! & (itau * dtvr / daysec - int(itau * dtvr / daysec))
	381	jD_cur = jD_ref + day_ini - day_ref + &
	382	& int (itau * dtvr / daysec)
	383	jH_cur = jH_ref + start_time + &
	384	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
	385	jD_cur = jD_cur + int(jH_cur)
	386	jH_cur = jH_cur - int(jH_cur)
[1279]	387	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
	388	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
	389	! write(lunout,*)'current date = ',an, mois, jour, secondes
[524]	390
	391	c rajout debug
	392	c lafin = .true.
	393
	394
	395	c Inbterface avec les routines de phylmd (phymars ... )
	396	c -----------------------------------------------------
	397
	398	c+jld
	399
	400	c Diagnostique de conservation de l'énergie : initialisation
[1146]	401	IF (ip_ebil_dyn.ge.1 ) THEN
[524]	402	ztit='bil dyn'
[1146]	403	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
	404	IF (planet_type.eq."earth") THEN
	405	CALL diagedyn(ztit,2,1,1,dtphys
	406	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	407	ENDIF
	408	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
[524]	409	c-jld
[1146]	410	#ifdef CPP_IOIPSL
[692]	411	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
[1146]	412	IF (first) THEN
	413	first=.false.
[692]	414	#include "ini_paramLMDZ_dyn.h"
[1146]	415	ENDIF
[692]	416	c
	417	#include "write_paramLMDZ_dyn.h"
	418	c
[1146]	419	#endif
	420	! #endif of #ifdef CPP_IOIPSL
[1279]	421	CALL calfis( lafin , jD_cur, jH_cur,
[524]	422	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
[960]	423	$ du,dv,dteta,dq,
[524]	424	$ flxw,
	425	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
	426
[1146]	427	IF (ok_strato) THEN
[1279]	428	CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
[1146]	429	ENDIF
[999]	430
[524]	431	c ajout des tendances physiques:
	432	c ------------------------------
[1146]	433	CALL addfi( dtphys, leapf, forward ,
[524]	434	$ ucov, vcov, teta , q ,ps ,
	435	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	436	c
	437	c Diagnostique de conservation de l'énergie : difference
[1146]	438	IF (ip_ebil_dyn.ge.1 ) THEN
[524]	439	ztit='bil phys'
[1146]	440	IF (planet_type.eq."earth") THEN
	441	CALL diagedyn(ztit,2,1,1,dtphys
	442	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	443	ENDIF
	444	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
	445
[1060]	446	ENDIF ! of IF( apphys )
[524]	447
[1146]	448	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
[1454]	449	! Academic case : Simple friction and Newtonan relaxation
	450	! -------------------------------------------------------
	451	DO l=1,llm
	452	DO ij=1,ip1jmp1
	453	teta(ij,l)=teta(ij,l)-dtvr*
	454	& (teta(ij,l)-tetarappel(ij,l))(knewt_g+knewt_t(l)clat4(ij))
	455	ENDDO
	456	ENDDO ! of DO l=1,llm
	457
[1505]	458	if (planet_type.eq."giant") then
	459	! add an intrinsic heat flux at the base of the atmosphere
	460	teta(:,1)=teta(:,1)+dtvraire(:)ihf/cpp/masse(:,1)
	461	endif
	462
	463	call friction(ucov,vcov,dtvr)
	464
[1454]	465	! Sponge layer (if any)
	466	IF (ok_strato) THEN
	467	dufi(:,:)=0.
	468	dvfi(:,:)=0.
	469	dtetafi(:,:)=0.
	470	dqfi(:,:,:)=0.
	471	dpfi(:)=0.
	472	CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
	473	CALL addfi( dtvr, leapf, forward ,
	474	$ ucov, vcov, teta , q ,ps ,
	475	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	476	ENDIF ! of IF (ok_strato)
	477	ENDIF ! of IF (iflag_phys.EQ.2)
[524]	478
	479
	480	c-jld
	481
	482	CALL pression ( ip1jmp1, ap, bp, ps, p )
[1520]	483	if (disvert_type==1) then
	484	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	485	else ! we assume that we are in the disvert_type==2 case
	486	CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
	487	endif
[524]	488
	489
	490	c-----------------------------------------------------------------------
	491	c dissipation horizontale et verticale des petites echelles:
	492	c ----------------------------------------------------------
	493
	494	IF(apdiss) THEN
	495
	496
	497	c calcul de l'energie cinetique avant dissipation
	498	call covcont(llm,ucov,vcov,ucont,vcont)
	499	call enercin(vcov,ucov,vcont,ucont,ecin0)
	500
	501	c dissipation
	502	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
	503	ucov=ucov+dudis
	504	vcov=vcov+dvdis
	505	c teta=teta+dtetadis
	506
	507
	508	c------------------------------------------------------------------------
	509	if (dissip_conservative) then
	510	C On rajoute la tendance due a la transform. Ec -> E therm. cree
	511	C lors de la dissipation
	512	call covcont(llm,ucov,vcov,ucont,vcont)
	513	call enercin(vcov,ucov,vcont,ucont,ecin)
	514	dtetaecdt= (ecin0-ecin)/ pk
	515	c teta=teta+dtetaecdt
	516	dtetadis=dtetadis+dtetaecdt
	517	endif
	518	teta=teta+dtetadis
	519	c------------------------------------------------------------------------
	520
	521
	522	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
	523	c ... Calcul de la valeur moyenne, unique de h aux poles .....
	524	c
	525
	526	DO l = 1, llm
	527	DO ij = 1,iim
	528	tppn(ij) = aire( ij ) * teta( ij ,l)
	529	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	530	ENDDO
	531	tpn = SSUM(iim,tppn,1)/apoln
	532	tps = SSUM(iim,tpps,1)/apols
	533
	534	DO ij = 1, iip1
	535	teta( ij ,l) = tpn
	536	teta(ij+ip1jm,l) = tps
	537	ENDDO
	538	ENDDO
	539
	540	DO ij = 1,iim
	541	tppn(ij) = aire( ij ) * ps ( ij )
	542	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
	543	ENDDO
	544	tpn = SSUM(iim,tppn,1)/apoln
	545	tps = SSUM(iim,tpps,1)/apols
	546
	547	DO ij = 1, iip1
	548	ps( ij ) = tpn
	549	ps(ij+ip1jm) = tps
	550	ENDDO
	551
	552
[1146]	553	END IF ! of IF(apdiss)
[524]	554
	555	c ajout debug
	556	c IF( lafin ) then
	557	c abort_message = 'Simulation finished'
	558	c call abort_gcm(modname,abort_message,0)
	559	c ENDIF
	560
	561	c ********************************************************************
	562	c ********************************************************************
	563	c .... fin de l'integration dynamique et physique pour le pas itau ..
	564	c ********************************************************************
	565	c ********************************************************************
	566
	567	c preparation du pas d'integration suivant ......
	568
	569	IF ( .NOT.purmats ) THEN
	570	c ........................................................
	571	c .............. schema matsuno + leapfrog ..............
	572	c ........................................................
	573
	574	IF(forward. OR. leapf) THEN
	575	itau= itau + 1
[1403]	576	c iday= day_ini+itau/day_step
	577	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	578	c IF(time.GT.1.) THEN
	579	c time = time-1.
	580	c iday = iday+1
	581	c ENDIF
[524]	582	ENDIF
	583
	584
	585	IF( itau. EQ. itaufinp1 ) then
[999]	586	if (flag_verif) then
[1279]	587	write(79,*) 'ucov',ucov
	588	write(80,*) 'vcov',vcov
	589	write(81,*) 'teta',teta
	590	write(82,*) 'ps',ps
	591	write(83,*) 'q',q
[999]	592	WRITE(85,*) 'q1 = ',q(:,:,1)
	593	WRITE(86,*) 'q3 = ',q(:,:,3)
	594	endif
[524]	595
	596	abort_message = 'Simulation finished'
	597
	598	call abort_gcm(modname,abort_message,0)
	599	ENDIF
	600	c-----------------------------------------------------------------------
	601	c ecriture du fichier histoire moyenne:
	602	c -------------------------------------
	603
	604	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	605	IF(itau.EQ.itaufin) THEN
	606	iav=1
	607	ELSE
	608	iav=0
	609	ENDIF
[1146]	610
	611	IF (ok_dynzon) THEN
[524]	612	#ifdef CPP_IOIPSL
[1403]	613	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
	614	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[524]	615	#endif
[1146]	616	END IF
[1403]	617	IF (ok_dyn_ave) THEN
	618	#ifdef CPP_IOIPSL
	619	CALL writedynav(itau,vcov,
	620	& ucov,teta,pk,phi,q,masse,ps,phis)
	621	#endif
	622	ENDIF
[524]	623
[1403]	624	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
[524]	625
	626	c-----------------------------------------------------------------------
	627	c ecriture de la bande histoire:
	628	c ------------------------------
	629
[1403]	630	IF( MOD(itau,iecri).EQ.0) THEN
	631	! Ehouarn: output only during LF or Backward Matsuno
	632	if (leapf.or.(.not.leapf.and.(.not.forward))) then
[1146]	633	nbetat = nbetatdem
	634	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	635	unat=0.
	636	do l=1,llm
	637	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	638	vnat(:,l)=vcov(:,l)/cv(:)
	639	enddo
[524]	640	#ifdef CPP_IOIPSL
[1403]	641	if (ok_dyn_ins) then
	642	! write(lunout,*) "leapfrog: call writehist, itau=",itau
	643	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	644	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	645	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
	646	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
	647	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
	648	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
	649	endif ! of if (ok_dyn_ins)
[1146]	650	#endif
	651	! For some Grads outputs of fields
[1403]	652	if (output_grads_dyn) then
[524]	653	#include "write_grads_dyn.h"
[1403]	654	endif
	655	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
[1146]	656	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
[524]	657
	658	IF(itau.EQ.itaufin) THEN
	659
	660
[1454]	661	! if (planet_type.eq."earth") then
[1146]	662	! Write an Earth-format restart file
[1577]	663	CALL dynredem1("restart.nc",start_time,
[1146]	664	& vcov,ucov,teta,q,masse,ps)
[1454]	665	! endif ! of if (planet_type.eq."earth")
[524]	666
	667	CLOSE(99)
[1146]	668	ENDIF ! of IF (itau.EQ.itaufin)
[524]	669
	670	c-----------------------------------------------------------------------
	671	c gestion de l'integration temporelle:
	672	c ------------------------------------
	673
	674	IF( MOD(itau,iperiod).EQ.0 ) THEN
	675	GO TO 1
	676	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	677
	678	IF( forward ) THEN
	679	c fin du pas forward et debut du pas backward
	680
	681	forward = .FALSE.
	682	leapf = .FALSE.
	683	GO TO 2
	684
	685	ELSE
	686	c fin du pas backward et debut du premier pas leapfrog
	687
	688	leapf = .TRUE.
	689	dt = 2.*dtvr
[1146]	690	GO TO 2
	691	END IF ! of IF (forward)
[524]	692	ELSE
	693
	694	c ...... pas leapfrog .....
	695
	696	leapf = .TRUE.
	697	dt = 2.*dtvr
	698	GO TO 2
[1146]	699	END IF ! of IF (MOD(itau,iperiod).EQ.0)
	700	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
[524]	701
[1146]	702	ELSE ! of IF (.not.purmats)
[524]	703
	704	c ........................................................
	705	c .............. schema matsuno ...............
	706	c ........................................................
	707	IF( forward ) THEN
	708
	709	itau = itau + 1
[1403]	710	c iday = day_ini+itau/day_step
	711	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	712	c
	713	c IF(time.GT.1.) THEN
	714	c time = time-1.
	715	c iday = iday+1
	716	c ENDIF
[524]	717
	718	forward = .FALSE.
	719	IF( itau. EQ. itaufinp1 ) then
	720	abort_message = 'Simulation finished'
	721	call abort_gcm(modname,abort_message,0)
	722	ENDIF
	723	GO TO 2
	724
[1403]	725	ELSE ! of IF(forward) i.e. backward step
[524]	726
[1146]	727	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
[524]	728	IF(itau.EQ.itaufin) THEN
	729	iav=1
	730	ELSE
	731	iav=0
	732	ENDIF
[1146]	733
	734	IF (ok_dynzon) THEN
[524]	735	#ifdef CPP_IOIPSL
[1403]	736	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
	737	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[524]	738	#endif
[1403]	739	ENDIF
	740	IF (ok_dyn_ave) THEN
	741	#ifdef CPP_IOIPSL
	742	CALL writedynav(itau,vcov,
	743	& ucov,teta,pk,phi,q,masse,ps,phis)
	744	#endif
	745	ENDIF
[524]	746
[1146]	747	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
[524]	748
[1146]	749	IF(MOD(itau,iecri ).EQ.0) THEN
[524]	750	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
[1146]	751	nbetat = nbetatdem
	752	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	753	unat=0.
	754	do l=1,llm
	755	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	756	vnat(:,l)=vcov(:,l)/cv(:)
	757	enddo
[524]	758	#ifdef CPP_IOIPSL
[1403]	759	if (ok_dyn_ins) then
	760	! write(lunout,*) "leapfrog: call writehist (b)",
	761	! & itau,iecri
	762	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	763	endif ! of if (ok_dyn_ins)
[1146]	764	#endif
	765	! For some Grads outputs
	766	if (output_grads_dyn) then
[524]	767	#include "write_grads_dyn.h"
[1146]	768	endif
[524]	769
[1146]	770	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
[524]	771
[1146]	772	IF(itau.EQ.itaufin) THEN
[1454]	773	! if (planet_type.eq."earth") then
[1577]	774	CALL dynredem1("restart.nc",start_time,
[1146]	775	& vcov,ucov,teta,q,masse,ps)
[1454]	776	! endif ! of if (planet_type.eq."earth")
[1146]	777	ENDIF ! of IF(itau.EQ.itaufin)
[524]	778
[1146]	779	forward = .TRUE.
	780	GO TO 1
[524]	781
[1146]	782	ENDIF ! of IF (forward)
[524]	783
[1146]	784	END IF ! of IF(.not.purmats)
[524]	785
	786	STOP
	787	END

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

Download in other formats:

Original Format