Context Navigation

leapfrog.F @ 3998

Last change on this file since 3998 was 1446, checked in by Ehouarn Millour, 14 years ago

Implemented modifications to enable running with only one tracer for planet types different from "earth". Rem: If flag 'planet_type' is set to "earth" (default behaviour) then there must be at least 2 tracers for the dynamics to function properly.

These updates do not induce any changes in model outputs with respect to previous revisions.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 23.5 KB

Rev	Line
[524]	1	!
[1279]	2	! $Id: leapfrog.F 1446 2010-10-22 09:27:25Z 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
[1403]	17	USE control_mod
[524]	18	IMPLICIT NONE
	19
	20	c ...... Version du 10/01/98 ..........
	21
	22	c avec coordonnees verticales hybrides
	23	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
	24
	25	c=======================================================================
	26	c
	27	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
	28	c -------
	29	c
	30	c Objet:
	31	c ------
	32	c
	33	c GCM LMD nouvelle grille
	34	c
	35	c=======================================================================
	36	c
	37	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
	38	c et possibilite d'appeler une fonction f(y) a derivee tangente
	39	c hyperbolique a la place de la fonction a derivee sinusoidale.
	40
	41	c ... Possibilite de choisir le shema pour l'advection de
	42	c q , en modifiant iadv dans traceur.def (10/02) .
	43	c
	44	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
	45	c Pour Van-Leer iadv=10
	46	c
	47	c-----------------------------------------------------------------------
	48	c Declarations:
	49	c -------------
	50
	51	#include "dimensions.h"
	52	#include "paramet.h"
	53	#include "comconst.h"
	54	#include "comdissnew.h"
	55	#include "comvert.h"
	56	#include "comgeom.h"
	57	#include "logic.h"
	58	#include "temps.h"
	59	#include "ener.h"
	60	#include "description.h"
	61	#include "serre.h"
[1403]	62	!#include "com_io_dyn.h"
[524]	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
[1403]	199	c iday = day_ini+itau/day_step
	200	c time = REAL(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
[1446]	211	dq(:,:,:)=0.
[524]	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
	256	2 CONTINUE
	257
	258	c-----------------------------------------------------------------------
	259
	260	c date:
	261	c -----
	262
	263
	264	c gestion des appels de la physique et des dissipations:
	265	c ------------------------------------------------------
	266	c
	267	c ... P.Le Van ( 6/02/95 ) ....
	268
	269	apphys = .FALSE.
	270	statcl = .FALSE.
	271	conser = .FALSE.
	272	apdiss = .FALSE.
	273
	274	IF( purmats ) THEN
[1403]	275	! Purely Matsuno time stepping
[524]	276	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
	277	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
	278	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
[1060]	279	s .and. iflag_phys.EQ.1 ) apphys = .TRUE.
[524]	280	ELSE
[1403]	281	! Leapfrog/Matsuno time stepping
[524]	282	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
	283	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
[1060]	284	IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
[524]	285	END IF
	286
[1403]	287	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
	288	! supress dissipation step
	289	if (llm.eq.1) then
	290	apdiss=.false.
	291	endif
	292
[524]	293	c-----------------------------------------------------------------------
	294	c calcul des tendances dynamiques:
	295	c --------------------------------
	296
	297	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	298
[1279]	299	time = jD_cur + jH_cur
[524]	300	CALL caldyn
	301	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
[1279]	302	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
[524]	303
[1279]	304
[524]	305	c-----------------------------------------------------------------------
	306	c calcul des tendances advection des traceurs (dont l'humidite)
	307	c -------------------------------------------------------------
	308
	309	IF( forward. OR . leapf ) THEN
	310
[960]	311	CALL caladvtrac(q,pbaru,pbarv,
	312	* p, masse, dq, teta,
	313	. flxw, pk)
	314
[524]	315	IF (offline) THEN
	316	Cmaf stokage du flux de masse pour traceurs OFF-LINE
	317
	318	#ifdef CPP_IOIPSL
[541]	319	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
	320	. dtvr, itau)
[524]	321	#endif
	322
	323
[1146]	324	ENDIF ! of IF (offline)
[524]	325	c
[1146]	326	ENDIF ! of IF( forward. OR . leapf )
[524]	327
	328
	329	c-----------------------------------------------------------------------
	330	c integrations dynamique et traceurs:
	331	c ----------------------------------
	332
	333
	334	CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
	335	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
	336	$ finvmaold )
	337
	338
	339	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
	340	c
	341	c-----------------------------------------------------------------------
	342	c calcul des tendances physiques:
	343	c -------------------------------
	344	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
	345	c
	346	IF( purmats ) THEN
	347	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
	348	ELSE
	349	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
	350	ENDIF
	351	c
	352	c
	353	IF( apphys ) THEN
	354	c
	355	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
	356	c
	357
	358	CALL pression ( ip1jmp1, ap, bp, ps, p )
	359	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
	360
[1279]	361	! rdaym_ini = itau * dtvr / daysec
	362	! rdayvrai = rdaym_ini + day_ini
	363	jD_cur = jD_ref + day_ini - day_ref
	364	$ + int (itau * dtvr / daysec)
	365	jH_cur = jH_ref + &
	366	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
	367	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
	368	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
	369	! write(lunout,*)'current date = ',an, mois, jour, secondes
[524]	370
	371	c rajout debug
	372	c lafin = .true.
	373
	374
	375	c Inbterface avec les routines de phylmd (phymars ... )
	376	c -----------------------------------------------------
	377
	378	c+jld
	379
	380	c Diagnostique de conservation de l'énergie : initialisation
[1146]	381	IF (ip_ebil_dyn.ge.1 ) THEN
[524]	382	ztit='bil dyn'
[1146]	383	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
	384	IF (planet_type.eq."earth") THEN
	385	CALL diagedyn(ztit,2,1,1,dtphys
	386	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	387	ENDIF
	388	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
[524]	389	c-jld
[1146]	390	#ifdef CPP_IOIPSL
[692]	391	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
[1146]	392	IF (first) THEN
	393	first=.false.
[692]	394	#include "ini_paramLMDZ_dyn.h"
[1146]	395	ENDIF
[692]	396	c
	397	#include "write_paramLMDZ_dyn.h"
	398	c
[1146]	399	#endif
	400	! #endif of #ifdef CPP_IOIPSL
[1279]	401	CALL calfis( lafin , jD_cur, jH_cur,
[524]	402	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
[960]	403	$ du,dv,dteta,dq,
[524]	404	$ flxw,
	405	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
	406
[1146]	407	IF (ok_strato) THEN
[1279]	408	CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
[1146]	409	ENDIF
[999]	410
[524]	411	c ajout des tendances physiques:
	412	c ------------------------------
[1146]	413	CALL addfi( dtphys, leapf, forward ,
[524]	414	$ ucov, vcov, teta , q ,ps ,
	415	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	416	c
	417	c Diagnostique de conservation de l'énergie : difference
[1146]	418	IF (ip_ebil_dyn.ge.1 ) THEN
[524]	419	ztit='bil phys'
[1146]	420	IF (planet_type.eq."earth") THEN
	421	CALL diagedyn(ztit,2,1,1,dtphys
	422	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	423	ENDIF
	424	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
	425
[1060]	426	ENDIF ! of IF( apphys )
[524]	427
[1146]	428	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
[1437]	429	! Academic case : Simple friction and Newtonan relaxation
	430	! -------------------------------------------------------
	431	DO l=1,llm
	432	DO ij=1,ip1jmp1
	433	teta(ij,l)=teta(ij,l)-dtvr*
	434	& (teta(ij,l)-tetarappel(ij,l))(knewt_g+knewt_t(l)clat4(ij))
	435	ENDDO
	436	ENDDO ! of DO l=1,llm
	437	call friction(ucov,vcov,dtvr)
	438
	439	! Sponge layer (if any)
	440	IF (ok_strato) THEN
	441	dufi(:,:)=0.
	442	dvfi(:,:)=0.
	443	dtetafi(:,:)=0.
	444	dqfi(:,:,:)=0.
	445	dpfi(:)=0.
	446	CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
	447	CALL addfi( dtvr, leapf, forward ,
	448	$ ucov, vcov, teta , q ,ps ,
	449	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	450	ENDIF ! of IF (ok_strato)
	451	ENDIF ! of IF (iflag_phys.EQ.2)
[524]	452
	453
	454	c-jld
	455
	456	CALL pression ( ip1jmp1, ap, bp, ps, p )
	457	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	458
	459
	460	c-----------------------------------------------------------------------
	461	c dissipation horizontale et verticale des petites echelles:
	462	c ----------------------------------------------------------
	463
	464	IF(apdiss) THEN
	465
	466
	467	c calcul de l'energie cinetique avant dissipation
	468	call covcont(llm,ucov,vcov,ucont,vcont)
	469	call enercin(vcov,ucov,vcont,ucont,ecin0)
	470
	471	c dissipation
	472	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
	473	ucov=ucov+dudis
	474	vcov=vcov+dvdis
	475	c teta=teta+dtetadis
	476
	477
	478	c------------------------------------------------------------------------
	479	if (dissip_conservative) then
	480	C On rajoute la tendance due a la transform. Ec -> E therm. cree
	481	C lors de la dissipation
	482	call covcont(llm,ucov,vcov,ucont,vcont)
	483	call enercin(vcov,ucov,vcont,ucont,ecin)
	484	dtetaecdt= (ecin0-ecin)/ pk
	485	c teta=teta+dtetaecdt
	486	dtetadis=dtetadis+dtetaecdt
	487	endif
	488	teta=teta+dtetadis
	489	c------------------------------------------------------------------------
	490
	491
	492	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
	493	c ... Calcul de la valeur moyenne, unique de h aux poles .....
	494	c
	495
	496	DO l = 1, llm
	497	DO ij = 1,iim
	498	tppn(ij) = aire( ij ) * teta( ij ,l)
	499	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	500	ENDDO
	501	tpn = SSUM(iim,tppn,1)/apoln
	502	tps = SSUM(iim,tpps,1)/apols
	503
	504	DO ij = 1, iip1
	505	teta( ij ,l) = tpn
	506	teta(ij+ip1jm,l) = tps
	507	ENDDO
	508	ENDDO
	509
	510	DO ij = 1,iim
	511	tppn(ij) = aire( ij ) * ps ( ij )
	512	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
	513	ENDDO
	514	tpn = SSUM(iim,tppn,1)/apoln
	515	tps = SSUM(iim,tpps,1)/apols
	516
	517	DO ij = 1, iip1
	518	ps( ij ) = tpn
	519	ps(ij+ip1jm) = tps
	520	ENDDO
	521
	522
[1146]	523	END IF ! of IF(apdiss)
[524]	524
	525	c ajout debug
	526	c IF( lafin ) then
	527	c abort_message = 'Simulation finished'
	528	c call abort_gcm(modname,abort_message,0)
	529	c ENDIF
	530
	531	c ********************************************************************
	532	c ********************************************************************
	533	c .... fin de l'integration dynamique et physique pour le pas itau ..
	534	c ********************************************************************
	535	c ********************************************************************
	536
	537	c preparation du pas d'integration suivant ......
	538
	539	IF ( .NOT.purmats ) THEN
	540	c ........................................................
	541	c .............. schema matsuno + leapfrog ..............
	542	c ........................................................
	543
	544	IF(forward. OR. leapf) THEN
	545	itau= itau + 1
[1403]	546	c iday= day_ini+itau/day_step
	547	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	548	c IF(time.GT.1.) THEN
	549	c time = time-1.
	550	c iday = iday+1
	551	c ENDIF
[524]	552	ENDIF
	553
	554
	555	IF( itau. EQ. itaufinp1 ) then
[999]	556	if (flag_verif) then
[1279]	557	write(79,*) 'ucov',ucov
	558	write(80,*) 'vcov',vcov
	559	write(81,*) 'teta',teta
	560	write(82,*) 'ps',ps
	561	write(83,*) 'q',q
[999]	562	WRITE(85,*) 'q1 = ',q(:,:,1)
	563	WRITE(86,*) 'q3 = ',q(:,:,3)
	564	endif
[524]	565
	566	abort_message = 'Simulation finished'
	567
	568	call abort_gcm(modname,abort_message,0)
	569	ENDIF
	570	c-----------------------------------------------------------------------
	571	c ecriture du fichier histoire moyenne:
	572	c -------------------------------------
	573
	574	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	575	IF(itau.EQ.itaufin) THEN
	576	iav=1
	577	ELSE
	578	iav=0
	579	ENDIF
[1146]	580
	581	IF (ok_dynzon) THEN
[524]	582	#ifdef CPP_IOIPSL
[1403]	583	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
	584	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[524]	585	#endif
[1146]	586	END IF
[1403]	587	IF (ok_dyn_ave) THEN
	588	#ifdef CPP_IOIPSL
	589	CALL writedynav(itau,vcov,
	590	& ucov,teta,pk,phi,q,masse,ps,phis)
	591	#endif
	592	ENDIF
[524]	593
[1403]	594	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
[524]	595
	596	c-----------------------------------------------------------------------
	597	c ecriture de la bande histoire:
	598	c ------------------------------
	599
[1403]	600	IF( MOD(itau,iecri).EQ.0) THEN
	601	! Ehouarn: output only during LF or Backward Matsuno
	602	if (leapf.or.(.not.leapf.and.(.not.forward))) then
[1146]	603	nbetat = nbetatdem
	604	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	605	unat=0.
	606	do l=1,llm
	607	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	608	vnat(:,l)=vcov(:,l)/cv(:)
	609	enddo
[524]	610	#ifdef CPP_IOIPSL
[1403]	611	if (ok_dyn_ins) then
	612	! write(lunout,*) "leapfrog: call writehist, itau=",itau
	613	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	614	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	615	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
	616	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
	617	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
	618	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
	619	endif ! of if (ok_dyn_ins)
[1146]	620	#endif
	621	! For some Grads outputs of fields
[1403]	622	if (output_grads_dyn) then
[524]	623	#include "write_grads_dyn.h"
[1403]	624	endif
	625	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
[1146]	626	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
[524]	627
	628	IF(itau.EQ.itaufin) THEN
	629
	630
[1446]	631	! if (planet_type.eq."earth") then
[1146]	632	! Write an Earth-format restart file
	633	CALL dynredem1("restart.nc",0.0,
	634	& vcov,ucov,teta,q,masse,ps)
[1446]	635	! endif ! of if (planet_type.eq."earth")
[524]	636
	637	CLOSE(99)
[1146]	638	ENDIF ! of IF (itau.EQ.itaufin)
[524]	639
	640	c-----------------------------------------------------------------------
	641	c gestion de l'integration temporelle:
	642	c ------------------------------------
	643
	644	IF( MOD(itau,iperiod).EQ.0 ) THEN
	645	GO TO 1
	646	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	647
	648	IF( forward ) THEN
	649	c fin du pas forward et debut du pas backward
	650
	651	forward = .FALSE.
	652	leapf = .FALSE.
	653	GO TO 2
	654
	655	ELSE
	656	c fin du pas backward et debut du premier pas leapfrog
	657
	658	leapf = .TRUE.
	659	dt = 2.*dtvr
[1146]	660	GO TO 2
	661	END IF ! of IF (forward)
[524]	662	ELSE
	663
	664	c ...... pas leapfrog .....
	665
	666	leapf = .TRUE.
	667	dt = 2.*dtvr
	668	GO TO 2
[1146]	669	END IF ! of IF (MOD(itau,iperiod).EQ.0)
	670	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
[524]	671
[1146]	672	ELSE ! of IF (.not.purmats)
[524]	673
	674	c ........................................................
	675	c .............. schema matsuno ...............
	676	c ........................................................
	677	IF( forward ) THEN
	678
	679	itau = itau + 1
[1403]	680	c iday = day_ini+itau/day_step
	681	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	682	c
	683	c IF(time.GT.1.) THEN
	684	c time = time-1.
	685	c iday = iday+1
	686	c ENDIF
[524]	687
	688	forward = .FALSE.
	689	IF( itau. EQ. itaufinp1 ) then
	690	abort_message = 'Simulation finished'
	691	call abort_gcm(modname,abort_message,0)
	692	ENDIF
	693	GO TO 2
	694
[1403]	695	ELSE ! of IF(forward) i.e. backward step
[524]	696
[1146]	697	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
[524]	698	IF(itau.EQ.itaufin) THEN
	699	iav=1
	700	ELSE
	701	iav=0
	702	ENDIF
[1146]	703
	704	IF (ok_dynzon) THEN
[524]	705	#ifdef CPP_IOIPSL
[1403]	706	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
	707	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[524]	708	#endif
[1403]	709	ENDIF
	710	IF (ok_dyn_ave) THEN
	711	#ifdef CPP_IOIPSL
	712	CALL writedynav(itau,vcov,
	713	& ucov,teta,pk,phi,q,masse,ps,phis)
	714	#endif
	715	ENDIF
[524]	716
[1146]	717	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
[524]	718
[1146]	719	IF(MOD(itau,iecri ).EQ.0) THEN
[524]	720	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
[1146]	721	nbetat = nbetatdem
	722	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	723	unat=0.
	724	do l=1,llm
	725	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	726	vnat(:,l)=vcov(:,l)/cv(:)
	727	enddo
[524]	728	#ifdef CPP_IOIPSL
[1403]	729	if (ok_dyn_ins) then
	730	! write(lunout,*) "leapfrog: call writehist (b)",
	731	! & itau,iecri
	732	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	733	endif ! of if (ok_dyn_ins)
[1146]	734	#endif
	735	! For some Grads outputs
	736	if (output_grads_dyn) then
[524]	737	#include "write_grads_dyn.h"
[1146]	738	endif
[524]	739
[1146]	740	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
[524]	741
[1146]	742	IF(itau.EQ.itaufin) THEN
[1446]	743	! if (planet_type.eq."earth") then
[1146]	744	CALL dynredem1("restart.nc",0.0,
	745	& vcov,ucov,teta,q,masse,ps)
[1446]	746	! endif ! of if (planet_type.eq."earth")
[1146]	747	ENDIF ! of IF(itau.EQ.itaufin)
[524]	748
[1146]	749	forward = .TRUE.
	750	GO TO 1
[524]	751
[1146]	752	ENDIF ! of IF (forward)
[524]	753
[1146]	754	END IF ! of IF(.not.purmats)
[524]	755
	756	STOP
	757	END

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

Download in other formats:

Original Format