Context Navigation

leapfrog.F @ 1366

Last change on this file since 1366 was 1286, checked in by Laurent Fairhead, 15 years ago
Les parametres definissant la transition eau glace dans les nuages sont maintenant lus dans physiq.def. Les valeurs par defaut donnent les memes resultats que precedemment. JLD
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 22.0 KB

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

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

Download in other formats:

Original Format