Context Navigation

gcm.F @ 182

Last change on this file since 182 was 78, checked in by lmdz, 25 years ago
Modification de F.Codron pour le schema d'advection de l'eau
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 20.7 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.)
	139
[2]	140	c-----------------------------------------------------------------------
	141	c Initialisations:
	142	c ----------------
	143
	144	modname = 'gcm'
	145	descript = 'Run GCM LMDZ'
	146	lafin = .FALSE.
	147	dynhist_file = 'dyn_hist.nc'
	148	dynhistave_file = 'dyn_hist_ave.nc'
[57]	149
	150	if (true_calendar) then
	151	call ioconf_calendar('gregorian')
	152	else
	153	call ioconf_calendar('360d')
	154	endif
[2]	155	c-----------------------------------------------------------------------
	156	c
	157	c .... Choix des shemas d'advection pour l'eau et les traceurs ...
	158	c ...................................................................
	159	c
	160	c iadv = 1 shema transport type "humidite specifique LMD"
	161	c iadv = 2 shema amont
	162	c iadv = 3 shema Van-leer
[78]	163	c iadv = 4 schema Van-leer + humidite specifique
	164	c Modif F.Codron
[2]	165	c
	166	c dans le tableau q(ij,l,iq) , iq = 1 pour l'eau vapeur
	167	c , iq = 2 pour l'eau liquide
	168	c et eventuellement , iq = 3, nqmx pour les autres traceurs
	169	c
	170	c iadv(1): choix pour l'eau vap. et iadv(2) : choix pour l'eau liq.
	171	c
	172	DO iq = 1, nqmx
	173	iadv( iq ) = 3
	174	ENDDO
	175	c
[78]	176	IF( iadv(1).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
	177	* ,' pour la vapeur d''eau'
	178	IF( iadv(1).EQ.2 ) PRINT *,' Choix du shema amont',' pour la'
	179	* ,' vapeur d''eau '
	180	IF( iadv(1).EQ.3 ) PRINT *,' Choix du shema Van-Leer ',' pour'
	181	* ,'la vapeur d''eau'
	182	IF( iadv(1).EQ.4 ) PRINT *,' Choix du shema Van-Leer + humidite'
	183	* ,' specifique pour la vapeur d''eau'
	184	c
	185	IF( iadv(1).LE.0.OR.iadv(1).GT.4 ) THEN
	186	PRINT *,' Erreur dans le choix de iadv (1).Corriger et repasser
	187	* , car iadv(1) = ', iadv(1)
	188	CALL ABORT
	189	ENDIF
	190
	191	DO iq = 2, nqmx
[2]	192	IF( iadv(iq).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
	193	* ,' pour le traceur no ', iq
	194	IF( iadv(iq).EQ.2 ) PRINT *,' Choix du shema amont',' pour le'
	195	* ,' traceur no ', iq
	196	IF( iadv(iq).EQ.3 ) PRINT *,' Choix du shema Van-Leer ',' pour'
	197	* ,'le traceur no ', iq
[78]	198
	199	IF( iadv(iq).EQ.4 ) THEN
	200	PRINT *,' Le shema Van-Leer + humidite specifique ',
	201	* ' est uniquement pour la vapeur d eau .'
	202	PRINT *,' Corriger iadv( ',iq, ') et repasser ! '
	203	CALL ABORT
	204	ENDIF
	205
	206	IF( iadv(iq).LE.0.OR.iadv(iq).GT.4 ) THEN
[2]	207	PRINT *,' Erreur dans le choix de iadv . Corriger et repasser
	208	* . '
[78]	209	CALL ABORT
[2]	210	ENDIF
	211	ENDDO
	212	c
	213	first = .TRUE.
	214	numvanle = nqmx + 1
	215	DO iq = 1, nqmx
[78]	216	IF((iadv(iq).EQ.3.OR.iadv(iq).EQ.4).AND.first ) THEN
[2]	217	numvanle = iq
	218	first = .FALSE.
	219	ENDIF
	220	ENDDO
	221	c
	222	DO iq = 1, nqmx
[78]	223	IF( (iadv(iq).NE.3.AND.iadv(iq).NE.4).AND.iq.GT.numvanle ) THEN
[2]	224	PRINT *,' Il y a discontinuite dans le choix du shema de ',
	225	* 'Van-leer pour les traceurs . Corriger et repasser . '
[78]	226	CALL ABORT
	227	ENDIF
	228	IF( iadv(iq).LT.1.OR.iadv(iq).GT.4 ) THEN
[2]	229	PRINT *,' Le choix de iadv est errone pour le traceur ',
	230	* iq
	231	STOP
	232	ENDIF
	233	ENDDO
	234	c
	235
	236	CALL dynetat0("start.nc",nqmx,vcov,ucov,
	237	. teta,q,masse,ps,phis, time_0)
	238
	239
	240	c OPEN( 99,file ='run.def',status='old',form='formatted')
	241	CALL defrun_new( 99, .TRUE. , clesphy0 )
	242	c CLOSE(99)
	243
	244	c on recalcule eventuellement le pas de temps
	245
	246	IF(MOD(day_step,iperiod).NE.0)
	247	* STOP'Il faut choisir un nb de pas par jour multiple de iperiod'
	248
	249	IF(MOD(day_step,iphysiq).NE.0)
	250	* STOP'Il faut choisir un nb de pas par jour multiple de iphysiq'
	251
	252	zdtvr = daysec/FLOAT(day_step)
	253	IF(dtvr.NE.zdtvr) THEN
	254	PRINT*,'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
	255	ENDIF
	256
	257	c nombre d'etats dans les fichiers demarrage et histoire
	258	nbetatdem = nday / iecri
	259	nbetatmoy = nday / periodav + 1
	260
	261	dtvr = zdtvr
	262	CALL iniconst
	263	CALL inigeom
	264
	265	CALL inifilr
	266	c
	267	c ...... P.Le Van ( modif le 29/04/97 ) .........
	268	c
	269	CALL inidissip( lstardis, nitergdiv, nitergrot, niterh ,
	270	* tetagdiv, tetagrot , tetatemp )
	271	c
	272
	273	CALL pression ( ip1jmp1, ap, bp, ps, p )
	274	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	275	c
	276
	277	c numero de stockage pour les fichiers de redemarrage:
	278
	279	c-----------------------------------------------------------------------
	280	c temps de depart et de fin:
	281	c --------------------------
	282
	283	itau = 0
	284	iday = day_ini+itau/day_step
	285	time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	286	IF(time.GT.1.) THEN
	287	time = time-1.
	288	iday = iday+1
	289	ENDIF
	290	itaufin = nday*day_step
	291	itaufinp1 = itaufin +1
	292
	293	day_end = day_ini + nday
	294	PRINT 300, itau,itaufin,day_ini,day_end
	295
	296	CALL dynredem0("restart.nc",day_end,anne_ini,phis,nqmx)
	297
	298	ecripar = .TRUE.
	299
	300	time_step = zdtvr
	301	t_ops = iecri * daysec
	302	t_wrt = iecri * daysec
	303	CALL inithist(dynhist_file,day_ini,anne_ini,time_step,
	304	. t_ops, t_wrt, nqmx, histid, histvid)
	305
	306	t_ops = iperiod * time_step
	307	t_wrt = periodav * daysec
	308	CALL initdynav(dynhistave_file,day_ini,anne_ini,time_step,
	309	. t_ops, t_wrt, nqmx, histaveid)
	310
	311	dtav = iperiod*dtvr/daysec
	312
	313
	314	c Quelques initialisations pour les traceurs
	315	call initial0(ijp1llm*nqmx,dq)
	316	istdyn=day_step/4 ! stockage toutes les 6h=1jour/4
	317	istphy=istdyn/iphysiq
	318
	319	c-----------------------------------------------------------------------
	320	c Debut de l'integration temporelle:
	321	c ----------------------------------
	322
	323	1 CONTINUE
	324	c call nudge(itau,ucov,vcov,teta,masse,ps)
	325	c
	326	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
	327	CALL test_period ( ucov,vcov,teta,q,p,phis )
	328	PRINT *,' ---- Test_period apres continue OK ! -----', itau
	329	c ENDIF
	330	c
	331	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
	332	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
	333	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
	334	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
	335	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
	336
	337	forward = .TRUE.
	338	leapf = .FALSE.
	339	dt = dtvr
	340
	341	c ... P.Le Van .26/04/94 ....
	342
	343	CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
	344	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
	345
	346
	347	2 CONTINUE
	348
	349	c-----------------------------------------------------------------------
	350
	351	c date:
	352	c -----
	353
	354
	355	c gestion des appels de la physique et des dissipations:
	356	c ------------------------------------------------------
	357	c
	358	c ... P.Le Van ( 6/02/95 ) ....
	359
	360	apphys = .FALSE.
	361	statcl = .FALSE.
	362	conser = .FALSE.
	363	apdiss = .FALSE.
	364
	365	IF( purmats ) THEN
	366	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
	367	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
	368	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
	369	$ .AND. physic ) apphys = .TRUE.
	370	ELSE
	371	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
	372	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
	373	IF( MOD(itau+1,iphysiq).EQ.0. AND. physic ) apphys = .TRUE.
	374	END IF
	375
	376	c-----------------------------------------------------------------------
	377	c calcul des tendances dynamiques:
	378	c --------------------------------
	379
	380	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	381	c
	382	c
	383	CALL caldyn
	384	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
	385	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
	386
	387	c-----------------------------------------------------------------------
	388	c calcul des tendances advection des traceurs (dont l'humidite)
	389	c -------------------------------------------------------------
	390
	391	IF( forward. OR . leapf ) THEN
	392	c
	393	DO 15 iq = 1, nqmx
	394	c
	395	IF( iadv(iq).EQ.1.OR.iadv(iq).EQ.2 ) THEN
	396	CALL traceur( iq,iadv,q,teta,pk,w, pbaru, pbarv, dq )
	397
	398	ELSE IF( iq.EQ. nqmx ) THEN
	399	c
	400	iapp_tracvl = 5
	401	c
	402	cccc iapp_tracvl est la frequence en pas du groupement des flux
	403	cccc de masse pour Van-Leer dans la routine tracvl .
	404	c
	405	CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv,
[78]	406	* p, masse, dq, iadv(1), teta, pk )
[2]	407	c
[78]	408	c ... Modif F.Codron ....
	409	c
[2]	410	ENDIF
	411	c
	412	15 CONTINUE
	413	C
	414	IF (offline) THEN
	415	Cmaf stokage du flux de masse pour traceurs OFF-LINE
[57]	416
	417	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
	418	. time_step, itau,fluxid, fluxvid,fluxdid )
	419
[2]	420	ENDIF
	421	c
	422	ENDIF
	423
	424
	425	c-----------------------------------------------------------------------
	426	c integrations dynamique et traceurs:
	427	c ----------------------------------
	428
	429
	430	CALL integrd ( iadv, 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
	431	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
	432	$ finvmaold )
	433
	434	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
	435	c
	436	c-----------------------------------------------------------------------
	437	c calcul des tendances physiques:
	438	c -------------------------------
	439	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
	440	c
	441	IF( purmats ) THEN
	442	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
	443	ELSE
	444	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
	445	ENDIF
	446	c
	447	c
	448	IF( apphys ) THEN
	449	c
	450	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
	451	c
	452
	453	CALL pression ( ip1jmp1, ap, bp, ps, p )
	454	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
	455
	456	rdaym_ini = itau * dtvr / daysec
	457	rdayvrai = rdaym_ini + day_ini
	458
	459	IF ( ecritphy.LT.1. ) THEN
	460	rday_ecri = rdaym_ini
	461	ELSE
	462	rday_ecri = NINT( rdayvrai )
	463	ENDIF
	464	c
	465	CALL calfis( nqmx, lafin ,rdayvrai,rday_ecri,time ,
	466	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
	467	$ du,dv,dteta,dq,w,clesphy0, dufi,dvfi,dhfi,dqfi,dpfi )
	468
	469	c ajout des tendances physiques:
	470	c ------------------------------
	471	CALL addfi( nqmx, dtphys, leapf, forward ,
	472	$ ucov, vcov, teta , q ,ps ,
	473	$ dufi, dvfi, dhfi , dqfi ,dpfi )
	474	c
	475	ENDIF
	476
	477	CALL pression ( ip1jmp1, ap, bp, ps, p )
	478	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	479
	480	c-----------------------------------------------------------------------
	481	c
	482	c
	483	c dissipation horizontale et verticale des petites echelles:
	484	c ----------------------------------------------------------
	485
	486	IF(apdiss) THEN
	487
	488	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dhdis)
	489
	490	CALL addit( ijp1llm,ucov ,dudis,ucov )
	491	CALL addit( ijmllm ,vcov ,dvdis,vcov )
	492	CALL addit( ijp1llm,teta ,dhdis,teta )
	493
	494
	495	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
	496	c ... Calcul de la valeur moyenne, unique de h aux poles .....
	497	c
	498
	499	DO l = 1, llm
	500	DO ij = 1,iim
	501	tppn(ij) = aire( ij ) * teta( ij ,l)
	502	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	503	ENDDO
	504	tpn = SSUM(iim,tppn,1)/apoln
	505	tps = SSUM(iim,tpps,1)/apols
	506
	507	DO ij = 1, iip1
	508	teta( ij ,l) = tpn
	509	teta(ij+ip1jm,l) = tps
	510	ENDDO
	511	ENDDO
	512
[49]	513	DO ij = 1,iim
	514	tppn(ij) = aire( ij ) * ps ( ij )
	515	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
	516	ENDDO
	517	tpn = SSUM(iim,tppn,1)/apoln
	518	tps = SSUM(iim,tpps,1)/apols
[2]	519
[49]	520	DO ij = 1, iip1
	521	ps( ij ) = tpn
	522	ps(ij+ip1jm) = tps
	523	ENDDO
	524
	525
[2]	526	END IF
	527
	528	c ********************************************************************
	529	c ********************************************************************
	530	c .... fin de l'integration dynamique et physique pour le pas itau ..
	531	c ********************************************************************
	532	c ********************************************************************
	533
	534	c preparation du pas d'integration suivant ......
	535
	536	IF ( .NOT.purmats ) THEN
	537	c ........................................................
	538	c .............. schema matsuno + leapfrog ..............
	539	c ........................................................
	540
	541	IF(forward. OR. leapf) THEN
	542	itau= itau + 1
	543	iday= day_ini+itau/day_step
	544	time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	545	IF(time.GT.1.) THEN
	546	time = time-1.
	547	iday = iday+1
	548	ENDIF
	549	ENDIF
	550
	551
	552	IF( itau. EQ. itaufinp1 ) then
	553	abort_message = 'Simulation finished'
	554	call abort_gcm(modname,abort_message,0)
	555	ENDIF
	556	c-----------------------------------------------------------------------
	557	c ecriture du fichier histoire moyenne:
	558	c -------------------------------------
	559
	560	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	561	IF(itau.EQ.itaufin) THEN
	562	iav=1
	563	ELSE
	564	iav=0
	565	ENDIF
	566	CALL writedynav(histaveid, nqmx, itau,vcov ,
	567	, ucov,teta,pk,phi,q,masse,ps,phis)
	568	ENDIF
	569
	570	c-----------------------------------------------------------------------
	571	c ecriture de la bande histoire:
	572	c ------------------------------
	573
	574	IF( MOD(itau,iecri*day_step).EQ.0) THEN
	575
	576	nbetat = nbetatdem
	577	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	578	CALL writehist( histid, histvid, nqmx, itau,vcov ,
	579	, ucov,teta,phi,q,masse,ps,phis)
	580
	581
	582	ENDIF
	583
	584	IF(itau.EQ.itaufin) THEN
	585
	586
	587	PRINT *,' Appel test_period avant redem ', itau
	588	CALL test_period ( ucov,vcov,teta,q,p,phis )
	589	CALL dynredem1("restart.nc",0.0,
	590	. vcov,ucov,teta,q,nqmx,masse,ps)
	591
	592	CLOSE(99)
	593	ENDIF
	594
	595	c-----------------------------------------------------------------------
	596	c gestion de l'integration temporelle:
	597	c ------------------------------------
	598
	599	IF( MOD(itau,iperiod).EQ.0 ) THEN
	600	GO TO 1
	601	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	602
	603	IF( forward ) THEN
	604	c fin du pas forward et debut du pas backward
	605
	606	forward = .FALSE.
	607	leapf = .FALSE.
	608	GO TO 2
	609
	610	ELSE
	611	c fin du pas backward et debut du premier pas leapfrog
	612
	613	leapf = .TRUE.
	614	dt = 2.*dtvr
	615	GO TO 2
	616	END IF
	617	ELSE
	618
	619	c ...... pas leapfrog .....
	620
	621	leapf = .TRUE.
	622	dt = 2.*dtvr
	623	GO TO 2
	624	END IF
	625
	626	ELSE
	627
	628	c ........................................................
	629	c .............. schema matsuno ...............
	630	c ........................................................
	631	IF( forward ) THEN
	632
	633	itau = itau + 1
	634	iday = day_ini+itau/day_step
	635	time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	636
	637	IF(time.GT.1.) THEN
	638	time = time-1.
	639	iday = iday+1
	640	ENDIF
	641
	642	forward = .FALSE.
	643	IF( itau. EQ. itaufinp1 ) then
	644	abort_message = 'Simulation finished'
	645	call abort_gcm(modname,abort_message,0)
	646	ENDIF
	647	GO TO 2
	648
	649	ELSE
	650
	651	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	652	IF(itau.EQ.itaufin) THEN
	653	iav=1
	654	ELSE
	655	iav=0
	656	ENDIF
	657	CALL writedynav(histaveid, nqmx, itau,vcov ,
	658	, ucov,teta,pk,phi,q,masse,ps,phis)
	659
	660	ENDIF
	661
	662	IF(MOD(itau,iecri*day_step).EQ.0) THEN
	663	nbetat = nbetatdem
	664	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	665	CALL writehist( histid, histvid, nqmx, itau,vcov ,
	666	, ucov,teta,phi,q,masse,ps,phis)
	667	ENDIF
	668
	669	IF(itau.EQ.itaufin)
	670	. CALL dynredem1("restart.nc",0.0,
	671	. vcov,ucov,teta,q,nqmx,masse,ps)
	672
	673	forward = .TRUE.
	674	GO TO 1
	675
	676	ENDIF
	677
	678	END IF
	679
	680	300 FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x,
	681	. 'c''est a dire du jour',i7,3x,'au jour',i7//)
	682
	683	STOP
	684	END

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

Download in other formats:

Original Format