Context Navigation

leapfrog.F @ 2520

Last change on this file since 2520 was 726, checked in by Laurent Fairhead, 18 years ago
Modifications pour rendre INCA plus independant de LMDZ ACo LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 19.7 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"
[726]	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
	187	itaufin = nday*day_step
	188	itaufinp1 = itaufin +1
	189
	190
	191	itau = 0
	192	iday = day_ini+itau/day_step
	193	time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	194	IF(time.GT.1.) THEN
	195	time = time-1.
	196	iday = iday+1
	197	ENDIF
	198
	199
	200	c-----------------------------------------------------------------------
	201	c On initialise la pression et la fonction d'Exner :
	202	c --------------------------------------------------
	203
	204	dq=0.
	205	CALL pression ( ip1jmp1, ap, bp, ps, p )
	206	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	207
	208	c-----------------------------------------------------------------------
	209	c Debut de l'integration temporelle:
	210	c ----------------------------------
	211
	212	1 CONTINUE
	213
	214
	215	#ifdef CPP_IOIPSL
[541]	216	if (ok_guide.and.(itaufin-itau-1)*dtvr.gt.21600) then
[524]	217	call guide(itau,ucov,vcov,teta,q,masse,ps)
[541]	218	else
[550]	219	IF(prt_level>9)WRITE(,)'attention on ne guide pas les ',
	220	. '6 dernieres heures'
[524]	221	endif
	222	#endif
	223	c
	224	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
	225	c CALL test_period ( ucov,vcov,teta,q,p,phis )
	226	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
	227	c ENDIF
	228	c
	229	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
	230	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
	231	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
	232	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
	233	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
	234
	235	forward = .TRUE.
	236	leapf = .FALSE.
	237	dt = dtvr
	238
	239	c ... P.Le Van .26/04/94 ....
	240
	241	CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
	242	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
	243
	244	call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
	245
	246	2 CONTINUE
	247
	248	c-----------------------------------------------------------------------
	249
	250	c date:
	251	c -----
	252
	253
	254	c gestion des appels de la physique et des dissipations:
	255	c ------------------------------------------------------
	256	c
	257	c ... P.Le Van ( 6/02/95 ) ....
	258
	259	apphys = .FALSE.
	260	statcl = .FALSE.
	261	conser = .FALSE.
	262	apdiss = .FALSE.
	263
	264	IF( purmats ) THEN
	265	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
	266	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
	267	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
	268	s .and. iflag_phys.NE.0 ) apphys = .TRUE.
	269	ELSE
	270	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
	271	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
	272	IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.NE.0) apphys=.TRUE.
	273	END IF
	274
	275	c-----------------------------------------------------------------------
	276	c calcul des tendances dynamiques:
	277	c --------------------------------
	278
	279	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	280
	281	CALL caldyn
	282	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
	283	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
	284
	285	c-----------------------------------------------------------------------
	286	c calcul des tendances advection des traceurs (dont l'humidite)
	287	c -------------------------------------------------------------
	288
	289	IF( forward. OR . leapf ) THEN
	290
	291	c
[616]	292	#ifdef INCA
[524]	293	CALL caladvtrac(q,pbaru,pbarv,
	294	* p, masse, dq, teta,
	295	. flxw,
[726]	296	. pk)
[524]	297	#else
	298	CALL caladvtrac(q,pbaru,pbarv,
	299	* p, masse, dq, teta,
	300	. pk)
	301	#endif
	302
	303	IF (offline) THEN
	304	Cmaf stokage du flux de masse pour traceurs OFF-LINE
	305
	306	#ifdef CPP_IOIPSL
[541]	307	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
	308	. dtvr, itau)
[524]	309	#endif
	310
	311
	312	ENDIF
	313	c
	314	ENDIF
	315
	316
	317	c-----------------------------------------------------------------------
	318	c integrations dynamique et traceurs:
	319	c ----------------------------------
	320
	321
	322	CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
	323	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
	324	$ finvmaold )
	325
	326
	327	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
	328	c
	329	c-----------------------------------------------------------------------
	330	c calcul des tendances physiques:
	331	c -------------------------------
	332	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
	333	c
	334	IF( purmats ) THEN
	335	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
	336	ELSE
	337	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
	338	ENDIF
	339	c
	340	c
	341	IF( apphys ) THEN
	342	c
	343	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
	344	c
	345
	346	CALL pression ( ip1jmp1, ap, bp, ps, p )
	347	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
	348
	349	rdaym_ini = itau * dtvr / daysec
	350	rdayvrai = rdaym_ini + day_ini
	351
	352
	353	c rajout debug
	354	c lafin = .true.
	355
	356
	357	c Inbterface avec les routines de phylmd (phymars ... )
	358	c -----------------------------------------------------
	359
	360	#ifdef CPP_PHYS
	361	c+jld
	362
	363	c Diagnostique de conservation de l'énergie : initialisation
	364	IF (ip_ebil_dyn.ge.1 ) THEN
	365	ztit='bil dyn'
	366	CALL diagedyn(ztit,2,1,1,dtphys
	367	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	368	ENDIF
	369	c-jld
[692]	370	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
	371	IF (first) THEN
	372	first=.false.
	373	#include "ini_paramLMDZ_dyn.h"
	374	ENDIF
	375	c
	376	#include "write_paramLMDZ_dyn.h"
	377	c
[524]	378
	379	CALL calfis( nq, lafin ,rdayvrai,time ,
	380	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
	381	$ du,dv,dteta,dq,w,
[616]	382	#ifdef INCA
[524]	383	$ flxw,
	384	#endif
	385	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
	386
	387	c ajout des tendances physiques:
	388	c ------------------------------
	389	CALL addfi( nqmx, dtphys, leapf, forward ,
	390	$ ucov, vcov, teta , q ,ps ,
	391	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	392	c
	393	c Diagnostique de conservation de l'énergie : difference
	394	IF (ip_ebil_dyn.ge.1 ) THEN
	395	ztit='bil phys'
	396	CALL diagedyn(ztit,2,1,1,dtphys
	397	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	398	ENDIF
	399	#else
	400
	401	c Calcul academique de la physique = Rappel Newtonien + fritcion
	402	c --------------------------------------------------------------
	403	teta(:,:)=teta(:,:)
	404	s -iphysiqdtvr(teta(:,:)-tetarappel(:,:))/taurappel
	405	call friction(ucov,vcov,iphysiq*dtvr)
	406
	407	#endif
	408
	409	c-jld
	410	ENDIF
	411
	412	CALL pression ( ip1jmp1, ap, bp, ps, p )
	413	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	414
	415
	416	c-----------------------------------------------------------------------
	417	c dissipation horizontale et verticale des petites echelles:
	418	c ----------------------------------------------------------
	419
	420	IF(apdiss) THEN
	421
	422
	423	c calcul de l'energie cinetique avant dissipation
	424	call covcont(llm,ucov,vcov,ucont,vcont)
	425	call enercin(vcov,ucov,vcont,ucont,ecin0)
	426
	427	c dissipation
	428	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
	429	ucov=ucov+dudis
	430	vcov=vcov+dvdis
	431	c teta=teta+dtetadis
	432
	433
	434	c------------------------------------------------------------------------
	435	if (dissip_conservative) then
	436	C On rajoute la tendance due a la transform. Ec -> E therm. cree
	437	C lors de la dissipation
	438	call covcont(llm,ucov,vcov,ucont,vcont)
	439	call enercin(vcov,ucov,vcont,ucont,ecin)
	440	dtetaecdt= (ecin0-ecin)/ pk
	441	c teta=teta+dtetaecdt
	442	dtetadis=dtetadis+dtetaecdt
	443	endif
	444	teta=teta+dtetadis
	445	c------------------------------------------------------------------------
	446
	447
	448	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
	449	c ... Calcul de la valeur moyenne, unique de h aux poles .....
	450	c
	451
	452	DO l = 1, llm
	453	DO ij = 1,iim
	454	tppn(ij) = aire( ij ) * teta( ij ,l)
	455	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	456	ENDDO
	457	tpn = SSUM(iim,tppn,1)/apoln
	458	tps = SSUM(iim,tpps,1)/apols
	459
	460	DO ij = 1, iip1
	461	teta( ij ,l) = tpn
	462	teta(ij+ip1jm,l) = tps
	463	ENDDO
	464	ENDDO
	465
	466	DO ij = 1,iim
	467	tppn(ij) = aire( ij ) * ps ( ij )
	468	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
	469	ENDDO
	470	tpn = SSUM(iim,tppn,1)/apoln
	471	tps = SSUM(iim,tpps,1)/apols
	472
	473	DO ij = 1, iip1
	474	ps( ij ) = tpn
	475	ps(ij+ip1jm) = tps
	476	ENDDO
	477
	478
	479	END IF
	480
	481	c ajout debug
	482	c IF( lafin ) then
	483	c abort_message = 'Simulation finished'
	484	c call abort_gcm(modname,abort_message,0)
	485	c ENDIF
	486
	487	c ********************************************************************
	488	c ********************************************************************
	489	c .... fin de l'integration dynamique et physique pour le pas itau ..
	490	c ********************************************************************
	491	c ********************************************************************
	492
	493	c preparation du pas d'integration suivant ......
	494
	495	IF ( .NOT.purmats ) THEN
	496	c ........................................................
	497	c .............. schema matsuno + leapfrog ..............
	498	c ........................................................
	499
	500	IF(forward. OR. leapf) THEN
	501	itau= itau + 1
	502	iday= day_ini+itau/day_step
	503	time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	504	IF(time.GT.1.) THEN
	505	time = time-1.
	506	iday = iday+1
	507	ENDIF
	508	ENDIF
	509
	510
	511	IF( itau. EQ. itaufinp1 ) then
[550]	512	c$$$ write(79,*) 'ucov',ucov
	513	c$$$ write(80,*) 'vcov',vcov
	514	c$$$ write(81,*) 'teta',teta
	515	c$$$ write(82,*) 'ps',ps
	516	c$$$ write(83,*) 'q',q
	517	c$$$ WRITE(85,*) 'q1 = ',q(:,:,1)
	518	c$$$ WRITE(86,*) 'q3 = ',q(:,:,3)
[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