Context Navigation

gcm.F @ 876

Last change on this file since 876 was 344, checked in by lmdz, 24 years ago
Inclusion des modifs de D. Hauglustaine pour la version 1 de INCA LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 21.8 KB

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

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

Download in other formats:

Original Format