Context Navigation

leapfrog.F @ 5404

Last change on this file since 5404 was 935, checked in by Laurent Fairhead, 17 years ago
Rajout ok_dynzon pour faire ou non les sorties bilan 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
[692]	65	#include "clesphys.h"
	66
[524]	67	integer nq
	68
	69	INTEGER longcles
	70	PARAMETER ( longcles = 20 )
	71	REAL clesphy0( longcles )
	72
	73	real zqmin,zqmax
	74	INTEGER nbetatmoy, nbetatdem,nbetat
	75
	76	c variables dynamiques
	77	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
	78	REAL teta(ip1jmp1,llm) ! temperature potentielle
	79	REAL q(ip1jmp1,llm,nqmx) ! champs advectes
	80	REAL ps(ip1jmp1) ! pression au sol
	81	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
	82	REAL pks(ip1jmp1) ! exner au sol
	83	REAL pk(ip1jmp1,llm) ! exner au milieu des couches
	84	REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
	85	REAL masse(ip1jmp1,llm) ! masse d'air
	86	REAL phis(ip1jmp1) ! geopotentiel au sol
	87	REAL phi(ip1jmp1,llm) ! geopotentiel
	88	REAL w(ip1jmp1,llm) ! vitesse verticale
	89
	90	c variables dynamiques intermediaire pour le transport
	91	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
	92
	93	c variables dynamiques au pas -1
	94	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
	95	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
	96	REAL massem1(ip1jmp1,llm)
	97
	98	c tendances dynamiques
	99	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
	100	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
	101
	102	c tendances de la dissipation
	103	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
	104	REAL dtetadis(ip1jmp1,llm)
	105
	106	c tendances physiques
	107	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
	108	REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
	109
	110	c variables pour le fichier histoire
	111	REAL dtav ! intervalle de temps elementaire
	112
	113	REAL tppn(iim),tpps(iim),tpn,tps
	114	c
	115	INTEGER itau,itaufinp1,iav
	116	INTEGER*4 iday ! jour julien
	117	REAL time ! Heure de la journee en fraction d'1 jour
	118
	119	REAL SSUM
	120	REAL time_0 , finvmaold(ip1jmp1,llm)
	121
[566]	122	cym LOGICAL lafin
	123	LOGICAL :: lafin=.false.
[524]	124	INTEGER ij,iq,l
	125	INTEGER ik
	126
	127	real time_step, t_wrt, t_ops
	128
	129	REAL rdayvrai,rdaym_ini
	130	LOGICAL first,callinigrads
[692]	131	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
	132	save first
	133	data first/.true./
	134	real dt_cum, zjulian
	135	character*10 infile
	136	integer zan, tau0, thoriid
	137	integer nid_ctesGCM
	138	save nid_ctesGCM
	139	real degres
	140	real rlong(iip1), rlatg(jjp1)
	141	real zx_tmp_2d(iip1,jjp1)
	142	integer ndex2d(iip1*jjp1)
	143	logical ok_sync
	144	parameter (ok_sync = .true.)
[524]	145
	146	data callinigrads/.true./
	147	character*10 string10
	148
	149	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
[616]	150	#ifdef INCA
[524]	151	REAL :: flxw(ip1jmp1,llm)
	152	#endif
	153
	154	c+jld variables test conservation energie
	155	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
	156	C Tendance de la temp. potentiel d (theta)/ d t due a la
	157	C tansformation d'energie cinetique en energie thermique
	158	C cree par la dissipation
	159	REAL dtetaecdt(ip1jmp1,llm)
	160	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
	161	REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
	162	REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
	163	CHARACTER*15 ztit
[692]	164	!IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
	165	!IM SAVE ip_ebil_dyn
	166	!IM DATA ip_ebil_dyn/0/
[524]	167	c-jld
	168
	169	character*80 dynhist_file, dynhistave_file
	170	character*20 modname
	171	character*80 abort_message
	172
	173	C Calendrier
	174	LOGICAL true_calendar
	175	PARAMETER (true_calendar = .false.)
	176
	177	logical dissip_conservative
	178	save dissip_conservative
	179	data dissip_conservative/.true./
	180
	181	LOGICAL prem
	182	save prem
	183	DATA prem/.true./
	184	INTEGER testita
	185	PARAMETER (testita = 9)
	186
[845]	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
	271	s .and. iflag_phys.NE.0 ) apphys = .TRUE.
	272	ELSE
	273	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
	274	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
	275	IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.NE.0) apphys=.TRUE.
	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
	294	c
[616]	295	#ifdef INCA
[524]	296	CALL caladvtrac(q,pbaru,pbarv,
	297	* p, masse, dq, teta,
	298	. flxw,
[762]	299	. pk)
[524]	300	#else
	301	CALL caladvtrac(q,pbaru,pbarv,
	302	* p, masse, dq, teta,
	303	. pk)
	304	#endif
	305
	306	IF (offline) THEN
	307	Cmaf stokage du flux de masse pour traceurs OFF-LINE
	308
	309	#ifdef CPP_IOIPSL
[541]	310	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
	311	. dtvr, itau)
[524]	312	#endif
	313
	314
	315	ENDIF
	316	c
	317	ENDIF
	318
	319
	320	c-----------------------------------------------------------------------
	321	c integrations dynamique et traceurs:
	322	c ----------------------------------
	323
	324
	325	CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
	326	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
	327	$ finvmaold )
	328
	329
	330	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
	331	c
	332	c-----------------------------------------------------------------------
	333	c calcul des tendances physiques:
	334	c -------------------------------
	335	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
	336	c
	337	IF( purmats ) THEN
	338	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
	339	ELSE
	340	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
	341	ENDIF
	342	c
	343	c
	344	IF( apphys ) THEN
	345	c
	346	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
	347	c
	348
	349	CALL pression ( ip1jmp1, ap, bp, ps, p )
	350	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
	351
	352	rdaym_ini = itau * dtvr / daysec
	353	rdayvrai = rdaym_ini + day_ini
	354
	355
	356	c rajout debug
	357	c lafin = .true.
	358
	359
	360	c Inbterface avec les routines de phylmd (phymars ... )
	361	c -----------------------------------------------------
	362
	363	#ifdef CPP_PHYS
	364	c+jld
	365
	366	c Diagnostique de conservation de l'énergie : initialisation
	367	IF (ip_ebil_dyn.ge.1 ) THEN
	368	ztit='bil dyn'
	369	CALL diagedyn(ztit,2,1,1,dtphys
	370	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	371	ENDIF
	372	c-jld
[692]	373	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
	374	IF (first) THEN
	375	first=.false.
	376	#include "ini_paramLMDZ_dyn.h"
	377	ENDIF
	378	c
	379	#include "write_paramLMDZ_dyn.h"
	380	c
[524]	381
	382	CALL calfis( nq, lafin ,rdayvrai,time ,
	383	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
	384	$ du,dv,dteta,dq,w,
[616]	385	#ifdef INCA
[524]	386	$ flxw,
	387	#endif
	388	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
	389
	390	c ajout des tendances physiques:
	391	c ------------------------------
	392	CALL addfi( nqmx, dtphys, leapf, forward ,
	393	$ ucov, vcov, teta , q ,ps ,
	394	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	395	c
	396	c Diagnostique de conservation de l'énergie : difference
	397	IF (ip_ebil_dyn.ge.1 ) THEN
	398	ztit='bil phys'
	399	CALL diagedyn(ztit,2,1,1,dtphys
	400	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	401	ENDIF
	402	#else
	403
	404	c Calcul academique de la physique = Rappel Newtonien + fritcion
	405	c --------------------------------------------------------------
	406	teta(:,:)=teta(:,:)
	407	s -iphysiqdtvr(teta(:,:)-tetarappel(:,:))/taurappel
	408	call friction(ucov,vcov,iphysiq*dtvr)
	409
	410	#endif
	411
	412	c-jld
	413	ENDIF
	414
	415	CALL pression ( ip1jmp1, ap, bp, ps, p )
	416	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	417
	418
	419	c-----------------------------------------------------------------------
	420	c dissipation horizontale et verticale des petites echelles:
	421	c ----------------------------------------------------------
	422
	423	IF(apdiss) THEN
	424
	425
	426	c calcul de l'energie cinetique avant dissipation
	427	call covcont(llm,ucov,vcov,ucont,vcont)
	428	call enercin(vcov,ucov,vcont,ucont,ecin0)
	429
	430	c dissipation
	431	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
	432	ucov=ucov+dudis
	433	vcov=vcov+dvdis
	434	c teta=teta+dtetadis
	435
	436
	437	c------------------------------------------------------------------------
	438	if (dissip_conservative) then
	439	C On rajoute la tendance due a la transform. Ec -> E therm. cree
	440	C lors de la dissipation
	441	call covcont(llm,ucov,vcov,ucont,vcont)
	442	call enercin(vcov,ucov,vcont,ucont,ecin)
	443	dtetaecdt= (ecin0-ecin)/ pk
	444	c teta=teta+dtetaecdt
	445	dtetadis=dtetadis+dtetaecdt
	446	endif
	447	teta=teta+dtetadis
	448	c------------------------------------------------------------------------
	449
	450
	451	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
	452	c ... Calcul de la valeur moyenne, unique de h aux poles .....
	453	c
	454
	455	DO l = 1, llm
	456	DO ij = 1,iim
	457	tppn(ij) = aire( ij ) * teta( ij ,l)
	458	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	459	ENDDO
	460	tpn = SSUM(iim,tppn,1)/apoln
	461	tps = SSUM(iim,tpps,1)/apols
	462
	463	DO ij = 1, iip1
	464	teta( ij ,l) = tpn
	465	teta(ij+ip1jm,l) = tps
	466	ENDDO
	467	ENDDO
	468
	469	DO ij = 1,iim
	470	tppn(ij) = aire( ij ) * ps ( ij )
	471	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
	472	ENDDO
	473	tpn = SSUM(iim,tppn,1)/apoln
	474	tps = SSUM(iim,tpps,1)/apols
	475
	476	DO ij = 1, iip1
	477	ps( ij ) = tpn
	478	ps(ij+ip1jm) = tps
	479	ENDDO
	480
	481
	482	END IF
	483
	484	c ajout debug
	485	c IF( lafin ) then
	486	c abort_message = 'Simulation finished'
	487	c call abort_gcm(modname,abort_message,0)
	488	c ENDIF
	489
	490	c ********************************************************************
	491	c ********************************************************************
	492	c .... fin de l'integration dynamique et physique pour le pas itau ..
	493	c ********************************************************************
	494	c ********************************************************************
	495
	496	c preparation du pas d'integration suivant ......
	497
	498	IF ( .NOT.purmats ) THEN
	499	c ........................................................
	500	c .............. schema matsuno + leapfrog ..............
	501	c ........................................................
	502
	503	IF(forward. OR. leapf) THEN
	504	itau= itau + 1
	505	iday= day_ini+itau/day_step
	506	time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	507	IF(time.GT.1.) THEN
	508	time = time-1.
	509	iday = iday+1
	510	ENDIF
	511	ENDIF
	512
	513
	514	IF( itau. EQ. itaufinp1 ) then
[773]	515	c write(79,*) 'ucov',ucov
	516	c write(80,*) 'vcov',vcov
	517	c write(81,*) 'teta',teta
	518	c write(82,*) 'ps',ps
	519	c write(83,*) 'q',q
	520	c WRITE(85,*) 'q1 = ',q(:,:,1)
	521	c WRITE(86,*) 'q3 = ',q(:,:,3)
[524]	522
	523	abort_message = 'Simulation finished'
	524
	525	call abort_gcm(modname,abort_message,0)
	526	ENDIF
	527	c-----------------------------------------------------------------------
	528	c ecriture du fichier histoire moyenne:
	529	c -------------------------------------
	530
	531	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	532	IF(itau.EQ.itaufin) THEN
	533	iav=1
	534	ELSE
	535	iav=0
	536	ENDIF
	537	#ifdef CPP_IOIPSL
[935]	538	IF (ok_dynzon) THEN
[524]	539	CALL writedynav(histaveid, nqmx, itau,vcov ,
	540	, ucov,teta,pk,phi,q,masse,ps,phis)
	541	call bilan_dyn (2,dtvriperiod,dtvrday_step*periodav,
	542	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[935]	543	ENDIF
[524]	544	#endif
	545
	546	ENDIF
	547
	548	c-----------------------------------------------------------------------
	549	c ecriture de la bande histoire:
	550	c ------------------------------
	551
	552	IF( MOD(itau,iecri ).EQ.0) THEN
	553	c IF( MOD(itau,iecri*day_step).EQ.0) THEN
	554
	555	nbetat = nbetatdem
	556	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	557	unat=0.
	558	do l=1,llm
	559	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	560	vnat(:,l)=vcov(:,l)/cv(:)
	561	enddo
	562	#ifdef CPP_IOIPSL
	563	c CALL writehist(histid,histvid, nqmx,itau,vcov,
	564	c s ucov,teta,phi,q,masse,ps,phis)
	565	#else
	566	#include "write_grads_dyn.h"
	567	#endif
	568
	569
	570	ENDIF
	571
	572	IF(itau.EQ.itaufin) THEN
	573
	574
	575	#ifdef CPP_IOIPSL
	576	CALL dynredem1("restart.nc",0.0,
	577	, vcov,ucov,teta,q,nqmx,masse,ps)
	578	#endif
	579
	580	CLOSE(99)
	581	ENDIF
	582
	583	c-----------------------------------------------------------------------
	584	c gestion de l'integration temporelle:
	585	c ------------------------------------
	586
	587	IF( MOD(itau,iperiod).EQ.0 ) THEN
	588	GO TO 1
	589	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	590
	591	IF( forward ) THEN
	592	c fin du pas forward et debut du pas backward
	593
	594	forward = .FALSE.
	595	leapf = .FALSE.
	596	GO TO 2
	597
	598	ELSE
	599	c fin du pas backward et debut du premier pas leapfrog
	600
	601	leapf = .TRUE.
	602	dt = 2.*dtvr
	603	GO TO 2
	604	END IF
	605	ELSE
	606
	607	c ...... pas leapfrog .....
	608
	609	leapf = .TRUE.
	610	dt = 2.*dtvr
	611	GO TO 2
	612	END IF
	613
	614	ELSE
	615
	616	c ........................................................
	617	c .............. schema matsuno ...............
	618	c ........................................................
	619	IF( forward ) THEN
	620
	621	itau = itau + 1
	622	iday = day_ini+itau/day_step
	623	time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	624
	625	IF(time.GT.1.) THEN
	626	time = time-1.
	627	iday = iday+1
	628	ENDIF
	629
	630	forward = .FALSE.
	631	IF( itau. EQ. itaufinp1 ) then
	632	abort_message = 'Simulation finished'
	633	call abort_gcm(modname,abort_message,0)
	634	ENDIF
	635	GO TO 2
	636
	637	ELSE
	638
	639	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	640	IF(itau.EQ.itaufin) THEN
	641	iav=1
	642	ELSE
	643	iav=0
	644	ENDIF
	645	#ifdef CPP_IOIPSL
[935]	646	IF (ok_dynzon) THEN
	647	CALL writedynav(histaveid, nqmx, itau,vcov ,
[524]	648	, ucov,teta,pk,phi,q,masse,ps,phis)
[935]	649	call bilan_dyn (2,dtvriperiod,dtvrday_step*periodav,
[524]	650	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[935]	651	ENDIF
[524]	652	#endif
	653
	654	ENDIF
	655
	656	IF(MOD(itau,iecri ).EQ.0) THEN
	657	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
	658	nbetat = nbetatdem
	659	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	660	unat=0.
	661	do l=1,llm
	662	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	663	vnat(:,l)=vcov(:,l)/cv(:)
	664	enddo
	665	#ifdef CPP_IOIPSL
	666	c CALL writehist( histid, histvid, nqmx, itau,vcov ,
	667	c , ucov,teta,phi,q,masse,ps,phis)
	668	#else
	669	#include "write_grads_dyn.h"
	670	#endif
	671
	672
	673	ENDIF
	674
	675	#ifdef CPP_IOIPSL
	676	IF(itau.EQ.itaufin)
	677	. CALL dynredem1("restart.nc",0.0,
	678	. vcov,ucov,teta,q,nqmx,masse,ps)
	679	#endif
	680
	681	forward = .TRUE.
	682	GO TO 1
	683
	684	ENDIF
	685
	686	END IF
	687
	688	STOP
	689	END

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

Download in other formats:

Original Format