Context Navigation

leapfrog.F @ 1060

Last change on this file since 1060 was 1060, checked in by Laurent Fairhead, 16 years ago
Correction bug cas physique Newtonnienne (iflag_phys=2) pour ne plus passer par la physique EM/FH/IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 19.9 KB

Rev	Line
[524]	1	!
	2	c
	3	c
[559]	4	SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0,
	5	& time_0)
[524]	6
	7
[692]	8	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
	9	USE IOIPSL
	10
[524]	11	IMPLICIT NONE
	12
	13	c ...... Version du 10/01/98 ..........
	14
	15	c avec coordonnees verticales hybrides
	16	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
	17
	18	c=======================================================================
	19	c
	20	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
	21	c -------
	22	c
	23	c Objet:
	24	c ------
	25	c
	26	c GCM LMD nouvelle grille
	27	c
	28	c=======================================================================
	29	c
	30	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
	31	c et possibilite d'appeler une fonction f(y) a derivee tangente
	32	c hyperbolique a la place de la fonction a derivee sinusoidale.
	33
	34	c ... Possibilite de choisir le shema pour l'advection de
	35	c q , en modifiant iadv dans traceur.def (10/02) .
	36	c
	37	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
	38	c Pour Van-Leer iadv=10
	39	c
	40	c-----------------------------------------------------------------------
	41	c Declarations:
	42	c -------------
	43
	44	#include "dimensions.h"
	45	#include "paramet.h"
	46	#include "comconst.h"
	47	#include "comdissnew.h"
	48	#include "comvert.h"
	49	#include "comgeom.h"
	50	#include "logic.h"
	51	#include "temps.h"
	52	#include "control.h"
	53	#include "ener.h"
	54	#include "description.h"
	55	#include "serre.h"
	56	#include "com_io_dyn.h"
	57	#include "iniprint.h"
[762]	58	#include "advtrac.h"
[524]	59	c#include "tracstoke.h"
	60
	61	#include "academic.h"
	62
[956]	63	! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
	64	! #include "clesphys.h"
[692]	65
[524]	66	integer nq
	67
	68	INTEGER longcles
	69	PARAMETER ( longcles = 20 )
	70	REAL clesphy0( longcles )
	71
	72	real zqmin,zqmax
	73	INTEGER nbetatmoy, nbetatdem,nbetat
	74
	75	c variables dynamiques
	76	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
	77	REAL teta(ip1jmp1,llm) ! temperature potentielle
	78	REAL q(ip1jmp1,llm,nqmx) ! champs advectes
	79	REAL ps(ip1jmp1) ! pression au sol
	80	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
	81	REAL pks(ip1jmp1) ! exner au sol
	82	REAL pk(ip1jmp1,llm) ! exner au milieu des couches
	83	REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
	84	REAL masse(ip1jmp1,llm) ! masse d'air
	85	REAL phis(ip1jmp1) ! geopotentiel au sol
	86	REAL phi(ip1jmp1,llm) ! geopotentiel
	87	REAL w(ip1jmp1,llm) ! vitesse verticale
	88
	89	c variables dynamiques intermediaire pour le transport
	90	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
	91
	92	c variables dynamiques au pas -1
	93	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
	94	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
	95	REAL massem1(ip1jmp1,llm)
	96
	97	c tendances dynamiques
	98	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
	99	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
	100
	101	c tendances de la dissipation
	102	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
	103	REAL dtetadis(ip1jmp1,llm)
	104
	105	c tendances physiques
	106	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
	107	REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
	108
	109	c variables pour le fichier histoire
	110	REAL dtav ! intervalle de temps elementaire
	111
	112	REAL tppn(iim),tpps(iim),tpn,tps
	113	c
	114	INTEGER itau,itaufinp1,iav
	115	INTEGER*4 iday ! jour julien
	116	REAL time ! Heure de la journee en fraction d'1 jour
	117
	118	REAL SSUM
	119	REAL time_0 , finvmaold(ip1jmp1,llm)
	120
[566]	121	cym LOGICAL lafin
	122	LOGICAL :: lafin=.false.
[524]	123	INTEGER ij,iq,l
	124	INTEGER ik
	125
	126	real time_step, t_wrt, t_ops
	127
	128	REAL rdayvrai,rdaym_ini
	129	LOGICAL first,callinigrads
[692]	130	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
	131	save first
	132	data first/.true./
	133	real dt_cum, zjulian
	134	character*10 infile
	135	integer zan, tau0, thoriid
	136	integer nid_ctesGCM
	137	save nid_ctesGCM
	138	real degres
	139	real rlong(iip1), rlatg(jjp1)
	140	real zx_tmp_2d(iip1,jjp1)
	141	integer ndex2d(iip1*jjp1)
	142	logical ok_sync
	143	parameter (ok_sync = .true.)
[524]	144
	145	data callinigrads/.true./
	146	character*10 string10
	147
	148	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
[960]	149	REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
[524]	150
	151	c+jld variables test conservation energie
	152	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
	153	C Tendance de la temp. potentiel d (theta)/ d t due a la
	154	C tansformation d'energie cinetique en energie thermique
	155	C cree par la dissipation
	156	REAL dtetaecdt(ip1jmp1,llm)
	157	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
	158	REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
	159	REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
	160	CHARACTER*15 ztit
[692]	161	!IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
	162	!IM SAVE ip_ebil_dyn
	163	!IM DATA ip_ebil_dyn/0/
[524]	164	c-jld
	165
	166	character*80 dynhist_file, dynhistave_file
	167	character*20 modname
	168	character*80 abort_message
	169
	170	C Calendrier
	171	LOGICAL true_calendar
	172	PARAMETER (true_calendar = .false.)
	173
	174	logical dissip_conservative
	175	save dissip_conservative
	176	data dissip_conservative/.true./
	177
	178	LOGICAL prem
	179	save prem
	180	DATA prem/.true./
	181	INTEGER testita
	182	PARAMETER (testita = 9)
	183
[1007]	184	logical , parameter :: flag_verif = .false.
[999]	185
	186
[825]	187	integer itau_w ! pas de temps ecriture = itap + itau_phy
	188
	189
[524]	190	itaufin = nday*day_step
	191	itaufinp1 = itaufin +1
	192
	193
	194	itau = 0
	195	iday = day_ini+itau/day_step
	196	time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	197	IF(time.GT.1.) THEN
	198	time = time-1.
	199	iday = iday+1
	200	ENDIF
	201
	202
	203	c-----------------------------------------------------------------------
	204	c On initialise la pression et la fonction d'Exner :
	205	c --------------------------------------------------
	206
	207	dq=0.
	208	CALL pression ( ip1jmp1, ap, bp, ps, p )
	209	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	210
	211	c-----------------------------------------------------------------------
	212	c Debut de l'integration temporelle:
	213	c ----------------------------------
	214
	215	1 CONTINUE
	216
	217
	218	#ifdef CPP_IOIPSL
[541]	219	if (ok_guide.and.(itaufin-itau-1)*dtvr.gt.21600) then
[524]	220	call guide(itau,ucov,vcov,teta,q,masse,ps)
[541]	221	else
[550]	222	IF(prt_level>9)WRITE(,)'attention on ne guide pas les ',
	223	. '6 dernieres heures'
[524]	224	endif
	225	#endif
	226	c
	227	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
	228	c CALL test_period ( ucov,vcov,teta,q,p,phis )
	229	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
	230	c ENDIF
	231	c
	232	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
	233	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
	234	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
	235	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
	236	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
	237
	238	forward = .TRUE.
	239	leapf = .FALSE.
	240	dt = dtvr
	241
	242	c ... P.Le Van .26/04/94 ....
	243
	244	CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
	245	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
	246
	247	call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
	248
	249	2 CONTINUE
	250
	251	c-----------------------------------------------------------------------
	252
	253	c date:
	254	c -----
	255
	256
	257	c gestion des appels de la physique et des dissipations:
	258	c ------------------------------------------------------
	259	c
	260	c ... P.Le Van ( 6/02/95 ) ....
	261
	262	apphys = .FALSE.
	263	statcl = .FALSE.
	264	conser = .FALSE.
	265	apdiss = .FALSE.
	266
	267	IF( purmats ) THEN
	268	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
	269	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
	270	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
[1060]	271	s .and. iflag_phys.EQ.1 ) apphys = .TRUE.
[524]	272	ELSE
	273	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
	274	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
[1060]	275	IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
[524]	276	END IF
	277
	278	c-----------------------------------------------------------------------
	279	c calcul des tendances dynamiques:
	280	c --------------------------------
	281
	282	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	283
	284	CALL caldyn
	285	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
	286	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
	287
	288	c-----------------------------------------------------------------------
	289	c calcul des tendances advection des traceurs (dont l'humidite)
	290	c -------------------------------------------------------------
	291
	292	IF( forward. OR . leapf ) THEN
	293
[960]	294	CALL caladvtrac(q,pbaru,pbarv,
	295	* p, masse, dq, teta,
	296	. flxw, pk)
	297
[524]	298	IF (offline) THEN
	299	Cmaf stokage du flux de masse pour traceurs OFF-LINE
	300
	301	#ifdef CPP_IOIPSL
[541]	302	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
	303	. dtvr, itau)
[524]	304	#endif
	305
	306
	307	ENDIF
	308	c
	309	ENDIF
	310
	311
	312	c-----------------------------------------------------------------------
	313	c integrations dynamique et traceurs:
	314	c ----------------------------------
	315
	316
	317	CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
	318	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
	319	$ finvmaold )
	320
	321
	322	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
	323	c
	324	c-----------------------------------------------------------------------
	325	c calcul des tendances physiques:
	326	c -------------------------------
	327	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
	328	c
	329	IF( purmats ) THEN
	330	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
	331	ELSE
	332	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
	333	ENDIF
	334	c
	335	c
	336	IF( apphys ) THEN
	337	c
	338	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
	339	c
	340
	341	CALL pression ( ip1jmp1, ap, bp, ps, p )
	342	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
	343
	344	rdaym_ini = itau * dtvr / daysec
	345	rdayvrai = rdaym_ini + day_ini
	346
	347
	348	c rajout debug
	349	c lafin = .true.
	350
	351
	352	c Inbterface avec les routines de phylmd (phymars ... )
	353	c -----------------------------------------------------
	354
	355	#ifdef CPP_PHYS
	356	c+jld
	357
	358	c Diagnostique de conservation de l'énergie : initialisation
	359	IF (ip_ebil_dyn.ge.1 ) THEN
	360	ztit='bil dyn'
	361	CALL diagedyn(ztit,2,1,1,dtphys
	362	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	363	ENDIF
	364	c-jld
[692]	365	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
	366	IF (first) THEN
	367	first=.false.
	368	#include "ini_paramLMDZ_dyn.h"
	369	ENDIF
	370	c
	371	#include "write_paramLMDZ_dyn.h"
	372	c
[524]	373
	374	CALL calfis( nq, lafin ,rdayvrai,time ,
	375	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
[960]	376	$ du,dv,dteta,dq,
[524]	377	$ flxw,
	378	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
	379
[999]	380	IF (ok_strato) THEN
	381	CALL top_bound( vcov,ucov,teta, dufi,dvfi,dtetafi)
	382	ENDIF
	383
[524]	384	c ajout des tendances physiques:
	385	c ------------------------------
	386	CALL addfi( nqmx, dtphys, leapf, forward ,
	387	$ ucov, vcov, teta , q ,ps ,
	388	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	389	c
	390	c Diagnostique de conservation de l'énergie : difference
	391	IF (ip_ebil_dyn.ge.1 ) THEN
	392	ztit='bil phys'
	393	CALL diagedyn(ztit,2,1,1,dtphys
	394	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	395	ENDIF
[1060]	396	#endif
	397	ENDIF ! of IF( apphys )
[524]	398
[1060]	399	IF(iflag_phys.EQ.2) THEN ! "Newtonian physics" case
	400	c Calcul academique de la physique = Rappel Newtonien + friction
[524]	401	c --------------------------------------------------------------
	402	teta(:,:)=teta(:,:)
	403	s -iphysiqdtvr(teta(:,:)-tetarappel(:,:))/taurappel
	404	call friction(ucov,vcov,iphysiq*dtvr)
[1060]	405	ENDIF
[524]	406
	407
	408	c-jld
	409
	410	CALL pression ( ip1jmp1, ap, bp, ps, p )
	411	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	412
	413
	414	c-----------------------------------------------------------------------
	415	c dissipation horizontale et verticale des petites echelles:
	416	c ----------------------------------------------------------
	417
	418	IF(apdiss) THEN
	419
	420
	421	c calcul de l'energie cinetique avant dissipation
	422	call covcont(llm,ucov,vcov,ucont,vcont)
	423	call enercin(vcov,ucov,vcont,ucont,ecin0)
	424
	425	c dissipation
	426	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
	427	ucov=ucov+dudis
	428	vcov=vcov+dvdis
	429	c teta=teta+dtetadis
	430
	431
	432	c------------------------------------------------------------------------
	433	if (dissip_conservative) then
	434	C On rajoute la tendance due a la transform. Ec -> E therm. cree
	435	C lors de la dissipation
	436	call covcont(llm,ucov,vcov,ucont,vcont)
	437	call enercin(vcov,ucov,vcont,ucont,ecin)
	438	dtetaecdt= (ecin0-ecin)/ pk
	439	c teta=teta+dtetaecdt
	440	dtetadis=dtetadis+dtetaecdt
	441	endif
	442	teta=teta+dtetadis
	443	c------------------------------------------------------------------------
	444
	445
	446	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
	447	c ... Calcul de la valeur moyenne, unique de h aux poles .....
	448	c
	449
	450	DO l = 1, llm
	451	DO ij = 1,iim
	452	tppn(ij) = aire( ij ) * teta( ij ,l)
	453	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	454	ENDDO
	455	tpn = SSUM(iim,tppn,1)/apoln
	456	tps = SSUM(iim,tpps,1)/apols
	457
	458	DO ij = 1, iip1
	459	teta( ij ,l) = tpn
	460	teta(ij+ip1jm,l) = tps
	461	ENDDO
	462	ENDDO
	463
	464	DO ij = 1,iim
	465	tppn(ij) = aire( ij ) * ps ( ij )
	466	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
	467	ENDDO
	468	tpn = SSUM(iim,tppn,1)/apoln
	469	tps = SSUM(iim,tpps,1)/apols
	470
	471	DO ij = 1, iip1
	472	ps( ij ) = tpn
	473	ps(ij+ip1jm) = tps
	474	ENDDO
	475
	476
	477	END IF
	478
	479	c ajout debug
	480	c IF( lafin ) then
	481	c abort_message = 'Simulation finished'
	482	c call abort_gcm(modname,abort_message,0)
	483	c ENDIF
	484
	485	c ********************************************************************
	486	c ********************************************************************
	487	c .... fin de l'integration dynamique et physique pour le pas itau ..
	488	c ********************************************************************
	489	c ********************************************************************
	490
	491	c preparation du pas d'integration suivant ......
	492
	493	IF ( .NOT.purmats ) THEN
	494	c ........................................................
	495	c .............. schema matsuno + leapfrog ..............
	496	c ........................................................
	497
	498	IF(forward. OR. leapf) THEN
	499	itau= itau + 1
	500	iday= day_ini+itau/day_step
	501	time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	502	IF(time.GT.1.) THEN
	503	time = time-1.
	504	iday = iday+1
	505	ENDIF
	506	ENDIF
	507
	508
	509	IF( itau. EQ. itaufinp1 ) then
[999]	510	if (flag_verif) then
	511	write(79,*) 'ucov',ucov
	512	write(80,*) 'vcov',vcov
	513	write(81,*) 'teta',teta
	514	write(82,*) 'ps',ps
	515	write(83,*) 'q',q
	516	WRITE(85,*) 'q1 = ',q(:,:,1)
	517	WRITE(86,*) 'q3 = ',q(:,:,3)
	518	endif
[524]	519
	520	abort_message = 'Simulation finished'
	521
	522	call abort_gcm(modname,abort_message,0)
	523	ENDIF
	524	c-----------------------------------------------------------------------
	525	c ecriture du fichier histoire moyenne:
	526	c -------------------------------------
	527
	528	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	529	IF(itau.EQ.itaufin) THEN
	530	iav=1
	531	ELSE
	532	iav=0
	533	ENDIF
	534	#ifdef CPP_IOIPSL
	535	CALL writedynav(histaveid, nqmx, itau,vcov ,
	536	, ucov,teta,pk,phi,q,masse,ps,phis)
	537	call bilan_dyn (2,dtvriperiod,dtvrday_step*periodav,
	538	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
	539	#endif
	540
	541	ENDIF
	542
	543	c-----------------------------------------------------------------------
	544	c ecriture de la bande histoire:
	545	c ------------------------------
	546
	547	IF( MOD(itau,iecri ).EQ.0) THEN
	548	c IF( MOD(itau,iecri*day_step).EQ.0) THEN
	549
	550	nbetat = nbetatdem
	551	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	552	unat=0.
	553	do l=1,llm
	554	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	555	vnat(:,l)=vcov(:,l)/cv(:)
	556	enddo
	557	#ifdef CPP_IOIPSL
	558	c CALL writehist(histid,histvid, nqmx,itau,vcov,
	559	c s ucov,teta,phi,q,masse,ps,phis)
	560	#else
	561	#include "write_grads_dyn.h"
	562	#endif
	563
	564
	565	ENDIF
	566
	567	IF(itau.EQ.itaufin) THEN
	568
	569
	570	#ifdef CPP_IOIPSL
	571	CALL dynredem1("restart.nc",0.0,
	572	, vcov,ucov,teta,q,nqmx,masse,ps)
	573	#endif
	574
	575	CLOSE(99)
	576	ENDIF
	577
	578	c-----------------------------------------------------------------------
	579	c gestion de l'integration temporelle:
	580	c ------------------------------------
	581
	582	IF( MOD(itau,iperiod).EQ.0 ) THEN
	583	GO TO 1
	584	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	585
	586	IF( forward ) THEN
	587	c fin du pas forward et debut du pas backward
	588
	589	forward = .FALSE.
	590	leapf = .FALSE.
	591	GO TO 2
	592
	593	ELSE
	594	c fin du pas backward et debut du premier pas leapfrog
	595
	596	leapf = .TRUE.
	597	dt = 2.*dtvr
	598	GO TO 2
	599	END IF
	600	ELSE
	601
	602	c ...... pas leapfrog .....
	603
	604	leapf = .TRUE.
	605	dt = 2.*dtvr
	606	GO TO 2
	607	END IF
	608
	609	ELSE
	610
	611	c ........................................................
	612	c .............. schema matsuno ...............
	613	c ........................................................
	614	IF( forward ) THEN
	615
	616	itau = itau + 1
	617	iday = day_ini+itau/day_step
	618	time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	619
	620	IF(time.GT.1.) THEN
	621	time = time-1.
	622	iday = iday+1
	623	ENDIF
	624
	625	forward = .FALSE.
	626	IF( itau. EQ. itaufinp1 ) then
	627	abort_message = 'Simulation finished'
	628	call abort_gcm(modname,abort_message,0)
	629	ENDIF
	630	GO TO 2
	631
	632	ELSE
	633
	634	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	635	IF(itau.EQ.itaufin) THEN
	636	iav=1
	637	ELSE
	638	iav=0
	639	ENDIF
	640	#ifdef CPP_IOIPSL
	641	CALL writedynav(histaveid, nqmx, itau,vcov ,
	642	, ucov,teta,pk,phi,q,masse,ps,phis)
	643	call bilan_dyn (2,dtvriperiod,dtvrday_step*periodav,
	644	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
	645	#endif
	646
	647	ENDIF
	648
	649	IF(MOD(itau,iecri ).EQ.0) THEN
	650	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
	651	nbetat = nbetatdem
	652	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	653	unat=0.
	654	do l=1,llm
	655	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	656	vnat(:,l)=vcov(:,l)/cv(:)
	657	enddo
	658	#ifdef CPP_IOIPSL
	659	c CALL writehist( histid, histvid, nqmx, itau,vcov ,
	660	c , ucov,teta,phi,q,masse,ps,phis)
	661	#else
	662	#include "write_grads_dyn.h"
	663	#endif
	664
	665
	666	ENDIF
	667
	668	#ifdef CPP_IOIPSL
	669	IF(itau.EQ.itaufin)
	670	. CALL dynredem1("restart.nc",0.0,
	671	. vcov,ucov,teta,q,nqmx,masse,ps)
	672	#endif
	673
	674	forward = .TRUE.
	675	GO TO 1
	676
	677	ENDIF
	678
	679	END IF
	680
	681	STOP
	682	END

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

Download in other formats:

Original Format