Context Navigation

gcm.F @ 331

Last change on this file since 331 was 195, checked in by lmdz, 24 years ago
Chgt pour faire tourner le offline FH LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 20.8 KB

Rev	Line
[2]	1
	2	PROGRAM gcm
	3
[57]	4	USE IOIPSL
	5
[2]	6	IMPLICIT NONE
	7
	8	c ...... Version du 10/01/98 ..........
	9
	10	c avec coordonnees verticales hybrides
	11	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
	12
	13	c=======================================================================
	14	c
	15	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
	16	c -------
	17	c
	18	c Objet:
	19	c ------
	20	c
	21	c GCM LMD nouvelle grille
	22	c
	23	c=======================================================================
[78]	24	c
	25	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
	26	c et possibilite d'appeler une fonction f(y) a derivee tangente
	27	c hyperbolique a la place de la fonction a derivee sinusoidale.
[2]	28
[78]	29	c ... Possibilite de choisir le shema de Van-leer pour l'advection de
	30	c q , en faisant iadv = 3 dans traceur (29/04/97) .
[2]	31	c
[78]	32	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
	33	c
[2]	34	c-----------------------------------------------------------------------
	35	c Declarations:
	36	c -------------
	37
	38	#include "dimensions.h"
	39	#include "paramet.h"
	40	#include "comconst.h"
	41	#include "comdissnew.h"
	42	#include "comvert.h"
	43	#include "comgeom.h"
	44	#include "logic.h"
	45	#include "temps.h"
	46	#include "control.h"
	47	#include "ener.h"
	48	#include "netcdf.inc"
	49	#include "description.h"
	50	#include "serre.h"
	51	#include "tracstoke.h"
	52
	53
	54	INTEGER longcles
	55	PARAMETER ( longcles = 20 )
	56	REAL clesphy0( longcles )
	57	SAVE clesphy0
	58
	59	INTEGER*4 iday ! jour julien
	60	REAL time ! Heure de la journee en fraction d'1 jour
	61	REAL zdtvr
	62	INTEGER nbetatmoy, nbetatdem,nbetat
	63
	64	c variables dynamiques
	65	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
	66	REAL teta(ip1jmp1,llm) ! temperature potentielle
	67	REAL q(ip1jmp1,llm,nqmx) ! champs advectes
	68	REAL ps(ip1jmp1) ! pression au sol
	69	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
	70	REAL pks(ip1jmp1) ! exner au sol
	71	REAL pk(ip1jmp1,llm) ! exner au milieu des couches
	72	REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
	73	REAL masse(ip1jmp1,llm) ! masse d'air
	74	REAL phis(ip1jmp1) ! geopotentiel au sol
	75	REAL phi(ip1jmp1,llm) ! geopotentiel
	76	REAL w(ip1jmp1,llm) ! vitesse verticale
	77
	78	c variables dynamiques intermediaire pour le transport
	79	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
	80
	81	c variables dynamiques au pas -1
	82	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
	83	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
	84	REAL massem1(ip1jmp1,llm)
	85
	86	c tendances dynamiques
	87	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
	88	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
	89
	90	c tendances de la dissipation
	91	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
	92	REAL dhdis(ip1jmp1,llm)
	93
	94	c tendances physiques
	95	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
	96	REAL dhfi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
	97
	98	c variables pour le fichier histoire
	99	REAL dtav ! intervalle de temps elementaire
	100
	101	REAL tppn(iim),tpps(iim),tpn,tps
	102	c
	103	INTEGER iadv(nqmx) ! indice schema de transport pour le traceur iq
	104
	105	INTEGER itau,itaufinp1,iav
	106
	107	EXTERNAL caldyn, traceur
	108	EXTERNAL dissip,geopot,iniconst,inifilr
	109	EXTERNAL integrd,SCOPY
	110	EXTERNAL iniav,writeav,writehist
	111	EXTERNAL inigeom
	112	EXTERNAL exner_hyb,addit
	113	EXTERNAL defrun_new, test_period
	114	REAL SSUM
	115	REAL time_0 , finvmaold(ip1jmp1,llm)
	116
	117	LOGICAL lafin
	118	INTEGER ij,iq,l,numvanle,iapp_tracvl
	119
	120
[57]	121	INTEGER fluxid, fluxvid,fluxdid
[2]	122	integer histid, histvid, histaveid
	123	real time_step, t_wrt, t_ops
	124
	125	REAL rdayvrai,rdaym_ini,rday_ecri
	126	LOGICAL first
	127	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
	128
	129	LOGICAL offline ! Controle du stockage ds "fluxmass"
[57]	130	PARAMETER (offline=.true.)
[2]	131
	132	character*80 dynhist_file, dynhistave_file
	133	character*20 modname
	134	character*80 abort_message
	135
[57]	136	C Calendrier
	137	LOGICAL true_calendar
	138	PARAMETER (true_calendar = .false.)
[190]	139	C Run nudge
	140	LOGICAL ok_nudge
	141	PARAMETER (ok_nudge = .false.)
[57]	142
[2]	143	c-----------------------------------------------------------------------
	144	c Initialisations:
	145	c ----------------
	146
	147	modname = 'gcm'
	148	descript = 'Run GCM LMDZ'
	149	lafin = .FALSE.
	150	dynhist_file = 'dyn_hist.nc'
	151	dynhistave_file = 'dyn_hist_ave.nc'
[57]	152
	153	if (true_calendar) then
	154	call ioconf_calendar('gregorian')
	155	else
	156	call ioconf_calendar('360d')
	157	endif
[2]	158	c-----------------------------------------------------------------------
	159	c
	160	c .... Choix des shemas d'advection pour l'eau et les traceurs ...
	161	c ...................................................................
	162	c
	163	c iadv = 1 shema transport type "humidite specifique LMD"
	164	c iadv = 2 shema amont
	165	c iadv = 3 shema Van-leer
[78]	166	c iadv = 4 schema Van-leer + humidite specifique
	167	c Modif F.Codron
[2]	168	c
	169	c dans le tableau q(ij,l,iq) , iq = 1 pour l'eau vapeur
	170	c , iq = 2 pour l'eau liquide
	171	c et eventuellement , iq = 3, nqmx pour les autres traceurs
	172	c
	173	c iadv(1): choix pour l'eau vap. et iadv(2) : choix pour l'eau liq.
	174	c
	175	DO iq = 1, nqmx
	176	iadv( iq ) = 3
	177	ENDDO
	178	c
[78]	179	IF( iadv(1).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
	180	* ,' pour la vapeur d''eau'
	181	IF( iadv(1).EQ.2 ) PRINT *,' Choix du shema amont',' pour la'
	182	* ,' vapeur d''eau '
	183	IF( iadv(1).EQ.3 ) PRINT *,' Choix du shema Van-Leer ',' pour'
	184	* ,'la vapeur d''eau'
	185	IF( iadv(1).EQ.4 ) PRINT *,' Choix du shema Van-Leer + humidite'
	186	* ,' specifique pour la vapeur d''eau'
	187	c
	188	IF( iadv(1).LE.0.OR.iadv(1).GT.4 ) THEN
	189	PRINT *,' Erreur dans le choix de iadv (1).Corriger et repasser
	190	* , car iadv(1) = ', iadv(1)
	191	CALL ABORT
	192	ENDIF
	193
	194	DO iq = 2, nqmx
[2]	195	IF( iadv(iq).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
	196	* ,' pour le traceur no ', iq
	197	IF( iadv(iq).EQ.2 ) PRINT *,' Choix du shema amont',' pour le'
	198	* ,' traceur no ', iq
	199	IF( iadv(iq).EQ.3 ) PRINT *,' Choix du shema Van-Leer ',' pour'
	200	* ,'le traceur no ', iq
[78]	201
	202	IF( iadv(iq).EQ.4 ) THEN
	203	PRINT *,' Le shema Van-Leer + humidite specifique ',
	204	* ' est uniquement pour la vapeur d eau .'
	205	PRINT *,' Corriger iadv( ',iq, ') et repasser ! '
	206	CALL ABORT
	207	ENDIF
	208
	209	IF( iadv(iq).LE.0.OR.iadv(iq).GT.4 ) THEN
[2]	210	PRINT *,' Erreur dans le choix de iadv . Corriger et repasser
	211	* . '
[78]	212	CALL ABORT
[2]	213	ENDIF
	214	ENDDO
	215	c
	216	first = .TRUE.
	217	numvanle = nqmx + 1
	218	DO iq = 1, nqmx
[78]	219	IF((iadv(iq).EQ.3.OR.iadv(iq).EQ.4).AND.first ) THEN
[2]	220	numvanle = iq
	221	first = .FALSE.
	222	ENDIF
	223	ENDDO
	224	c
	225	DO iq = 1, nqmx
[78]	226	IF( (iadv(iq).NE.3.AND.iadv(iq).NE.4).AND.iq.GT.numvanle ) THEN
[2]	227	PRINT *,' Il y a discontinuite dans le choix du shema de ',
	228	* 'Van-leer pour les traceurs . Corriger et repasser . '
[78]	229	CALL ABORT
	230	ENDIF
	231	IF( iadv(iq).LT.1.OR.iadv(iq).GT.4 ) THEN
[2]	232	PRINT *,' Le choix de iadv est errone pour le traceur ',
	233	* iq
	234	STOP
	235	ENDIF
	236	ENDDO
	237	c
	238
	239	CALL dynetat0("start.nc",nqmx,vcov,ucov,
	240	. teta,q,masse,ps,phis, time_0)
	241
	242
	243	c OPEN( 99,file ='run.def',status='old',form='formatted')
	244	CALL defrun_new( 99, .TRUE. , clesphy0 )
	245	c CLOSE(99)
	246
	247	c on recalcule eventuellement le pas de temps
	248
	249	IF(MOD(day_step,iperiod).NE.0)
	250	* STOP'Il faut choisir un nb de pas par jour multiple de iperiod'
	251
	252	IF(MOD(day_step,iphysiq).NE.0)
	253	* STOP'Il faut choisir un nb de pas par jour multiple de iphysiq'
	254
	255	zdtvr = daysec/FLOAT(day_step)
	256	IF(dtvr.NE.zdtvr) THEN
	257	PRINT*,'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
	258	ENDIF
	259
	260	c nombre d'etats dans les fichiers demarrage et histoire
	261	nbetatdem = nday / iecri
	262	nbetatmoy = nday / periodav + 1
	263
	264	dtvr = zdtvr
	265	CALL iniconst
	266	CALL inigeom
	267
	268	CALL inifilr
	269	c
	270	c ...... P.Le Van ( modif le 29/04/97 ) .........
	271	c
	272	CALL inidissip( lstardis, nitergdiv, nitergrot, niterh ,
	273	* tetagdiv, tetagrot , tetatemp )
	274	c
	275
	276	CALL pression ( ip1jmp1, ap, bp, ps, p )
	277	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	278	c
	279
	280	c numero de stockage pour les fichiers de redemarrage:
	281
	282	c-----------------------------------------------------------------------
	283	c temps de depart et de fin:
	284	c --------------------------
	285
	286	itau = 0
	287	iday = day_ini+itau/day_step
	288	time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	289	IF(time.GT.1.) THEN
	290	time = time-1.
	291	iday = iday+1
	292	ENDIF
	293	itaufin = nday*day_step
	294	itaufinp1 = itaufin +1
	295
	296	day_end = day_ini + nday
	297	PRINT 300, itau,itaufin,day_ini,day_end
	298
	299	CALL dynredem0("restart.nc",day_end,anne_ini,phis,nqmx)
	300
	301	ecripar = .TRUE.
	302
	303	time_step = zdtvr
	304	t_ops = iecri * daysec
	305	t_wrt = iecri * daysec
	306	CALL inithist(dynhist_file,day_ini,anne_ini,time_step,
	307	. t_ops, t_wrt, nqmx, histid, histvid)
	308
	309	t_ops = iperiod * time_step
	310	t_wrt = periodav * daysec
	311	CALL initdynav(dynhistave_file,day_ini,anne_ini,time_step,
	312	. t_ops, t_wrt, nqmx, histaveid)
	313
	314	dtav = iperiod*dtvr/daysec
	315
	316
	317	c Quelques initialisations pour les traceurs
	318	call initial0(ijp1llm*nqmx,dq)
	319	istdyn=day_step/4 ! stockage toutes les 6h=1jour/4
	320	istphy=istdyn/iphysiq
	321
	322	c-----------------------------------------------------------------------
	323	c Debut de l'integration temporelle:
	324	c ----------------------------------
	325
	326	1 CONTINUE
[190]	327
	328	if (ok_nudge) then
	329	call nudge(itau,ucov,vcov,teta,masse,ps)
	330	endif
[2]	331	c
	332	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
	333	CALL test_period ( ucov,vcov,teta,q,p,phis )
	334	PRINT *,' ---- Test_period apres continue OK ! -----', itau
	335	c ENDIF
	336	c
	337	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
	338	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
	339	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
	340	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
	341	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
	342
	343	forward = .TRUE.
	344	leapf = .FALSE.
	345	dt = dtvr
	346
	347	c ... P.Le Van .26/04/94 ....
	348
	349	CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
	350	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
	351
	352
	353	2 CONTINUE
	354
	355	c-----------------------------------------------------------------------
	356
	357	c date:
	358	c -----
	359
	360
	361	c gestion des appels de la physique et des dissipations:
	362	c ------------------------------------------------------
	363	c
	364	c ... P.Le Van ( 6/02/95 ) ....
	365
	366	apphys = .FALSE.
	367	statcl = .FALSE.
	368	conser = .FALSE.
	369	apdiss = .FALSE.
	370
	371	IF( purmats ) THEN
	372	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
	373	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
	374	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
	375	$ .AND. physic ) apphys = .TRUE.
	376	ELSE
	377	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
	378	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
	379	IF( MOD(itau+1,iphysiq).EQ.0. AND. physic ) apphys = .TRUE.
	380	END IF
	381
	382	c-----------------------------------------------------------------------
	383	c calcul des tendances dynamiques:
	384	c --------------------------------
	385
	386	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	387	c
	388	c
	389	CALL caldyn
	390	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
	391	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
	392
	393	c-----------------------------------------------------------------------
	394	c calcul des tendances advection des traceurs (dont l'humidite)
	395	c -------------------------------------------------------------
	396
	397	IF( forward. OR . leapf ) THEN
	398	c
	399	DO 15 iq = 1, nqmx
	400	c
	401	IF( iadv(iq).EQ.1.OR.iadv(iq).EQ.2 ) THEN
	402	CALL traceur( iq,iadv,q,teta,pk,w, pbaru, pbarv, dq )
	403
	404	ELSE IF( iq.EQ. nqmx ) THEN
	405	c
	406	iapp_tracvl = 5
	407	c
	408	cccc iapp_tracvl est la frequence en pas du groupement des flux
	409	cccc de masse pour Van-Leer dans la routine tracvl .
	410	c
	411	CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv,
[78]	412	* p, masse, dq, iadv(1), teta, pk )
[2]	413	c
[78]	414	c ... Modif F.Codron ....
	415	c
[2]	416	ENDIF
	417	c
	418	15 CONTINUE
	419	C
	420	IF (offline) THEN
	421	Cmaf stokage du flux de masse pour traceurs OFF-LINE
[57]	422
	423	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
[195]	424	. time_step, itau)
[57]	425
[195]	426
[2]	427	ENDIF
	428	c
	429	ENDIF
	430
	431
	432	c-----------------------------------------------------------------------
	433	c integrations dynamique et traceurs:
	434	c ----------------------------------
	435
	436
	437	CALL integrd ( iadv, 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
	438	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
	439	$ finvmaold )
	440
	441	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
	442	c
	443	c-----------------------------------------------------------------------
	444	c calcul des tendances physiques:
	445	c -------------------------------
	446	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
	447	c
	448	IF( purmats ) THEN
	449	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
	450	ELSE
	451	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
	452	ENDIF
	453	c
	454	c
	455	IF( apphys ) THEN
	456	c
	457	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
	458	c
	459
	460	CALL pression ( ip1jmp1, ap, bp, ps, p )
	461	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
	462
	463	rdaym_ini = itau * dtvr / daysec
	464	rdayvrai = rdaym_ini + day_ini
	465
	466	IF ( ecritphy.LT.1. ) THEN
	467	rday_ecri = rdaym_ini
	468	ELSE
	469	rday_ecri = NINT( rdayvrai )
	470	ENDIF
	471	c
	472	CALL calfis( nqmx, lafin ,rdayvrai,rday_ecri,time ,
	473	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
	474	$ du,dv,dteta,dq,w,clesphy0, dufi,dvfi,dhfi,dqfi,dpfi )
	475
	476	c ajout des tendances physiques:
	477	c ------------------------------
	478	CALL addfi( nqmx, dtphys, leapf, forward ,
	479	$ ucov, vcov, teta , q ,ps ,
	480	$ dufi, dvfi, dhfi , dqfi ,dpfi )
	481	c
	482	ENDIF
	483
	484	CALL pression ( ip1jmp1, ap, bp, ps, p )
	485	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	486
	487	c-----------------------------------------------------------------------
	488	c
	489	c
	490	c dissipation horizontale et verticale des petites echelles:
	491	c ----------------------------------------------------------
	492
	493	IF(apdiss) THEN
	494
	495	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dhdis)
	496
	497	CALL addit( ijp1llm,ucov ,dudis,ucov )
	498	CALL addit( ijmllm ,vcov ,dvdis,vcov )
	499	CALL addit( ijp1llm,teta ,dhdis,teta )
	500
	501
	502	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
	503	c ... Calcul de la valeur moyenne, unique de h aux poles .....
	504	c
	505
	506	DO l = 1, llm
	507	DO ij = 1,iim
	508	tppn(ij) = aire( ij ) * teta( ij ,l)
	509	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	510	ENDDO
	511	tpn = SSUM(iim,tppn,1)/apoln
	512	tps = SSUM(iim,tpps,1)/apols
	513
	514	DO ij = 1, iip1
	515	teta( ij ,l) = tpn
	516	teta(ij+ip1jm,l) = tps
	517	ENDDO
	518	ENDDO
	519
[49]	520	DO ij = 1,iim
	521	tppn(ij) = aire( ij ) * ps ( ij )
	522	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
	523	ENDDO
	524	tpn = SSUM(iim,tppn,1)/apoln
	525	tps = SSUM(iim,tpps,1)/apols
[2]	526
[49]	527	DO ij = 1, iip1
	528	ps( ij ) = tpn
	529	ps(ij+ip1jm) = tps
	530	ENDDO
	531
	532
[2]	533	END IF
	534
	535	c ********************************************************************
	536	c ********************************************************************
	537	c .... fin de l'integration dynamique et physique pour le pas itau ..
	538	c ********************************************************************
	539	c ********************************************************************
	540
	541	c preparation du pas d'integration suivant ......
	542
	543	IF ( .NOT.purmats ) THEN
	544	c ........................................................
	545	c .............. schema matsuno + leapfrog ..............
	546	c ........................................................
	547
	548	IF(forward. OR. leapf) THEN
	549	itau= itau + 1
	550	iday= day_ini+itau/day_step
	551	time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	552	IF(time.GT.1.) THEN
	553	time = time-1.
	554	iday = iday+1
	555	ENDIF
	556	ENDIF
	557
	558
	559	IF( itau. EQ. itaufinp1 ) then
	560	abort_message = 'Simulation finished'
	561	call abort_gcm(modname,abort_message,0)
	562	ENDIF
	563	c-----------------------------------------------------------------------
	564	c ecriture du fichier histoire moyenne:
	565	c -------------------------------------
	566
	567	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	568	IF(itau.EQ.itaufin) THEN
	569	iav=1
	570	ELSE
	571	iav=0
	572	ENDIF
	573	CALL writedynav(histaveid, nqmx, itau,vcov ,
	574	, ucov,teta,pk,phi,q,masse,ps,phis)
	575	ENDIF
	576
	577	c-----------------------------------------------------------------------
	578	c ecriture de la bande histoire:
	579	c ------------------------------
	580
	581	IF( MOD(itau,iecri*day_step).EQ.0) THEN
	582
	583	nbetat = nbetatdem
	584	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	585	CALL writehist( histid, histvid, nqmx, itau,vcov ,
	586	, ucov,teta,phi,q,masse,ps,phis)
	587
	588
	589	ENDIF
	590
	591	IF(itau.EQ.itaufin) THEN
	592
	593
	594	PRINT *,' Appel test_period avant redem ', itau
	595	CALL test_period ( ucov,vcov,teta,q,p,phis )
	596	CALL dynredem1("restart.nc",0.0,
	597	. vcov,ucov,teta,q,nqmx,masse,ps)
	598
	599	CLOSE(99)
	600	ENDIF
	601
	602	c-----------------------------------------------------------------------
	603	c gestion de l'integration temporelle:
	604	c ------------------------------------
	605
	606	IF( MOD(itau,iperiod).EQ.0 ) THEN
	607	GO TO 1
	608	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	609
	610	IF( forward ) THEN
	611	c fin du pas forward et debut du pas backward
	612
	613	forward = .FALSE.
	614	leapf = .FALSE.
	615	GO TO 2
	616
	617	ELSE
	618	c fin du pas backward et debut du premier pas leapfrog
	619
	620	leapf = .TRUE.
	621	dt = 2.*dtvr
	622	GO TO 2
	623	END IF
	624	ELSE
	625
	626	c ...... pas leapfrog .....
	627
	628	leapf = .TRUE.
	629	dt = 2.*dtvr
	630	GO TO 2
	631	END IF
	632
	633	ELSE
	634
	635	c ........................................................
	636	c .............. schema matsuno ...............
	637	c ........................................................
	638	IF( forward ) THEN
	639
	640	itau = itau + 1
	641	iday = day_ini+itau/day_step
	642	time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	643
	644	IF(time.GT.1.) THEN
	645	time = time-1.
	646	iday = iday+1
	647	ENDIF
	648
	649	forward = .FALSE.
	650	IF( itau. EQ. itaufinp1 ) then
	651	abort_message = 'Simulation finished'
	652	call abort_gcm(modname,abort_message,0)
	653	ENDIF
	654	GO TO 2
	655
	656	ELSE
	657
	658	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	659	IF(itau.EQ.itaufin) THEN
	660	iav=1
	661	ELSE
	662	iav=0
	663	ENDIF
	664	CALL writedynav(histaveid, nqmx, itau,vcov ,
	665	, ucov,teta,pk,phi,q,masse,ps,phis)
	666
	667	ENDIF
	668
	669	IF(MOD(itau,iecri*day_step).EQ.0) THEN
	670	nbetat = nbetatdem
	671	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	672	CALL writehist( histid, histvid, nqmx, itau,vcov ,
	673	, ucov,teta,phi,q,masse,ps,phis)
	674	ENDIF
	675
	676	IF(itau.EQ.itaufin)
	677	. CALL dynredem1("restart.nc",0.0,
	678	. vcov,ucov,teta,q,nqmx,masse,ps)
	679
	680	forward = .TRUE.
	681	GO TO 1
	682
	683	ENDIF
	684
	685	END IF
	686
	687	300 FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x,
	688	. 'c''est a dire du jour',i7,3x,'au jour',i7//)
	689
	690	STOP
	691	END

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

Download in other formats:

Original Format