Context Navigation

leapfrog.F @ 1249

Last change on this file since 1249 was 1146, checked in by Laurent Fairhead, 16 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 21.5 KB

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

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

Download in other formats:

Original Format