Context Navigation

gcm.F @ 1036

Last change on this file since 1036 was 1036, checked in by emillour, 11 years ago

Mars GCM: (a first step towards using parallel dynamics)

IMPORTANT CHANGE: Implemented dynamic tracers. It is no longer necessary to compile the model with the '-t #' option, number of tracers is simply read from tracer.def file (as before). Adapted makegcm_* scripts (and co.) accordingly. Technical aspects of the switch to dynamic tracers are:
- advtrac.h (in dyn3d) removed and replaced by module infotrac.F
- tracer.h (in phymars) removed and replaced by module tracer_mod.F90 (which contains nqmx, the number of tracers, etc. and can be used anywhere in the physics).
Included some side cleanups: removed unused files (in dyn3d) anldoppler2.F, anl_mcdstats.F and anl_stats-diag.F, and all the unecessary dimensions.* files in grid/dimension.
Checked that changes are clean and that GCM yields identical results (in debug mode) to previous svn version.

File size: 19.5 KB

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

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

Download in other formats:

Original Format