Context Navigation

leapfrog.F @ 1032

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

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

Download in other formats:

Original Format