Context Navigation

gcm.F @ 242

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

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

Download in other formats:

Original Format