Context Navigation

leapfrog.F @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago
Import initial LMDZ5
File size: 23.7 KB

Rev	Line
[1]	1	!
	2	! $Id: leapfrog.F 1437 2010-09-30 08:29:10Z emillour $
	3	!
	4	c
	5	c
	6	SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
	7	& time_0)
	8
	9
	10	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
	11	#ifdef CPP_IOIPSL
	12	use IOIPSL
	13	#endif
	14	USE infotrac
	15	USE guide_mod, ONLY : guide_main
	16	USE write_field
	17	USE control_mod
	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"
	62	!#include "com_io_dyn.h"
	63	#include "iniprint.h"
	64	#include "academic.h"
	65
	66	! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
	67	! #include "clesphys.h"
	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,nqtot) ! 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,nqtot),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,nqtot),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 iday ! jour julien
	117	REAL time
	118
	119	REAL SSUM
	120	REAL time_0 , finvmaold(ip1jmp1,llm)
	121
	122	cym LOGICAL lafin
	123	LOGICAL :: lafin=.false.
	124	INTEGER ij,iq,l
	125	INTEGER ik
	126
	127	real time_step, t_wrt, t_ops
	128
	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
	136	LOGICAL first,callinigrads
	137	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
	138	save first
	139	data first/.true./
	140	real dt_cum
	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.)
	151
	152	data callinigrads/.true./
	153	character*10 string10
	154
	155	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
	156	REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
	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
	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/
	171	c-jld
	172
	173	character*80 dynhist_file, dynhistave_file
	174	character(len=20) :: modname
	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
	187	logical , parameter :: flag_verif = .false.
	188
	189
	190	integer itau_w ! pas de temps ecriture = itap + itau_phy
	191
	192
	193	itaufin = nday*day_step
	194	itaufinp1 = itaufin +1
	195	modname="leapfrog"
	196
	197
	198	itau = 0
	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
	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
	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))
	224
	225
	226	#ifdef CPP_IOIPSL
	227	if (ok_guide) then
	228	call guide_main(itau,ucov,vcov,teta,q,masse,ps)
	229	endif
	230	#endif
	231
	232
	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
	239
	240	! Save fields obtained at previous time step as '...m1'
	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	! Ehouarn: what is this for? zqmin & zqmax are not used anyway ...
	257	! call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
	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	! Purely Matsuno time stepping
	279	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
	280	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
	281	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
	282	s .and. iflag_phys.EQ.1 ) apphys = .TRUE.
	283	ELSE
	284	! Leapfrog/Matsuno time stepping
	285	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
	286	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
	287	IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
	288	END IF
	289
	290	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
	291	! supress dissipation step
	292	if (llm.eq.1) then
	293	apdiss=.false.
	294	endif
	295
	296	c-----------------------------------------------------------------------
	297	c calcul des tendances dynamiques:
	298	c --------------------------------
	299
	300	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	301
	302	time = jD_cur + jH_cur
	303	CALL caldyn
	304	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
	305	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
	306
	307
	308	c-----------------------------------------------------------------------
	309	c calcul des tendances advection des traceurs (dont l'humidite)
	310	c -------------------------------------------------------------
	311
	312	IF( forward. OR . leapf ) THEN
	313
	314	CALL caladvtrac(q,pbaru,pbarv,
	315	* p, masse, dq, teta,
	316	. flxw, pk)
	317
	318	IF (offline) THEN
	319	Cmaf stokage du flux de masse pour traceurs OFF-LINE
	320
	321	#ifdef CPP_IOIPSL
	322	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
	323	. dtvr, itau)
	324	#endif
	325
	326
	327	ENDIF ! of IF (offline)
	328	c
	329	ENDIF ! of IF( forward. OR . leapf )
	330
	331
	332	c-----------------------------------------------------------------------
	333	c integrations dynamique et traceurs:
	334	c ----------------------------------
	335
	336
	337	CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
	338	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
	339	$ finvmaold )
	340
	341
	342	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
	343	c
	344	c-----------------------------------------------------------------------
	345	c calcul des tendances physiques:
	346	c -------------------------------
	347	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
	348	c
	349	IF( purmats ) THEN
	350	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
	351	ELSE
	352	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
	353	ENDIF
	354	c
	355	c
	356	IF( apphys ) THEN
	357	c
	358	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
	359	c
	360
	361	CALL pression ( ip1jmp1, ap, bp, ps, p )
	362	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
	363
	364	! rdaym_ini = itau * dtvr / daysec
	365	! rdayvrai = rdaym_ini + day_ini
	366	jD_cur = jD_ref + day_ini - day_ref
	367	$ + int (itau * dtvr / daysec)
	368	jH_cur = jH_ref + &
	369	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
	370	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
	371	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
	372	! write(lunout,*)'current date = ',an, mois, jour, secondes
	373
	374	c rajout debug
	375	c lafin = .true.
	376
	377
	378	c Inbterface avec les routines de phylmd (phymars ... )
	379	c -----------------------------------------------------
	380
	381	c+jld
	382
	383	c Diagnostique de conservation de l'énergie : initialisation
	384	IF (ip_ebil_dyn.ge.1 ) THEN
	385	ztit='bil dyn'
	386	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
	387	IF (planet_type.eq."earth") THEN
	388	CALL diagedyn(ztit,2,1,1,dtphys
	389	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	390	ENDIF
	391	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
	392	c-jld
	393	#ifdef CPP_IOIPSL
	394	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
	395	IF (first) THEN
	396	first=.false.
	397	#include "ini_paramLMDZ_dyn.h"
	398	ENDIF
	399	c
	400	#include "write_paramLMDZ_dyn.h"
	401	c
	402	#endif
	403	! #endif of #ifdef CPP_IOIPSL
	404	CALL calfis( lafin , jD_cur, jH_cur,
	405	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
	406	$ du,dv,dteta,dq,
	407	$ flxw,
	408	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
	409
	410	IF (ok_strato) THEN
	411	CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
	412	ENDIF
	413
	414	c ajout des tendances physiques:
	415	c ------------------------------
	416	CALL addfi( dtphys, leapf, forward ,
	417	$ ucov, vcov, teta , q ,ps ,
	418	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	419	c
	420	c Diagnostique de conservation de l'énergie : difference
	421	IF (ip_ebil_dyn.ge.1 ) THEN
	422	ztit='bil phys'
	423	IF (planet_type.eq."earth") THEN
	424	CALL diagedyn(ztit,2,1,1,dtphys
	425	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	426	ENDIF
	427	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
	428
	429	ENDIF ! of IF( apphys )
	430
	431	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
	432	! Academic case : Simple friction and Newtonan relaxation
	433	! -------------------------------------------------------
	434	DO l=1,llm
	435	DO ij=1,ip1jmp1
	436	teta(ij,l)=teta(ij,l)-dtvr*
	437	& (teta(ij,l)-tetarappel(ij,l))(knewt_g+knewt_t(l)clat4(ij))
	438	ENDDO
	439	ENDDO ! of DO l=1,llm
	440	call friction(ucov,vcov,dtvr)
	441
	442	! Sponge layer (if any)
	443	IF (ok_strato) THEN
	444	dufi(:,:)=0.
	445	dvfi(:,:)=0.
	446	dtetafi(:,:)=0.
	447	dqfi(:,:,:)=0.
	448	dpfi(:)=0.
	449	CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
	450	CALL addfi( dtvr, leapf, forward ,
	451	$ ucov, vcov, teta , q ,ps ,
	452	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	453	ENDIF ! of IF (ok_strato)
	454	ENDIF ! of IF (iflag_phys.EQ.2)
	455
	456
	457	c-jld
	458
	459	CALL pression ( ip1jmp1, ap, bp, ps, p )
	460	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	461
	462
	463	c-----------------------------------------------------------------------
	464	c dissipation horizontale et verticale des petites echelles:
	465	c ----------------------------------------------------------
	466
	467	IF(apdiss) THEN
	468
	469
	470	c calcul de l'energie cinetique avant dissipation
	471	call covcont(llm,ucov,vcov,ucont,vcont)
	472	call enercin(vcov,ucov,vcont,ucont,ecin0)
	473
	474	c dissipation
	475	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
	476	ucov=ucov+dudis
	477	vcov=vcov+dvdis
	478	c teta=teta+dtetadis
	479
	480
	481	c------------------------------------------------------------------------
	482	if (dissip_conservative) then
	483	C On rajoute la tendance due a la transform. Ec -> E therm. cree
	484	C lors de la dissipation
	485	call covcont(llm,ucov,vcov,ucont,vcont)
	486	call enercin(vcov,ucov,vcont,ucont,ecin)
	487	dtetaecdt= (ecin0-ecin)/ pk
	488	c teta=teta+dtetaecdt
	489	dtetadis=dtetadis+dtetaecdt
	490	endif
	491	teta=teta+dtetadis
	492	c------------------------------------------------------------------------
	493
	494
	495	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
	496	c ... Calcul de la valeur moyenne, unique de h aux poles .....
	497	c
	498
	499	DO l = 1, llm
	500	DO ij = 1,iim
	501	tppn(ij) = aire( ij ) * teta( ij ,l)
	502	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	503	ENDDO
	504	tpn = SSUM(iim,tppn,1)/apoln
	505	tps = SSUM(iim,tpps,1)/apols
	506
	507	DO ij = 1, iip1
	508	teta( ij ,l) = tpn
	509	teta(ij+ip1jm,l) = tps
	510	ENDDO
	511	ENDDO
	512
	513	DO ij = 1,iim
	514	tppn(ij) = aire( ij ) * ps ( ij )
	515	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
	516	ENDDO
	517	tpn = SSUM(iim,tppn,1)/apoln
	518	tps = SSUM(iim,tpps,1)/apols
	519
	520	DO ij = 1, iip1
	521	ps( ij ) = tpn
	522	ps(ij+ip1jm) = tps
	523	ENDDO
	524
	525
	526	END IF ! of IF(apdiss)
	527
	528	c ajout debug
	529	c IF( lafin ) then
	530	c abort_message = 'Simulation finished'
	531	c call abort_gcm(modname,abort_message,0)
	532	c ENDIF
	533
	534	c ********************************************************************
	535	c ********************************************************************
	536	c .... fin de l'integration dynamique et physique pour le pas itau ..
	537	c ********************************************************************
	538	c ********************************************************************
	539
	540	c preparation du pas d'integration suivant ......
	541
	542	IF ( .NOT.purmats ) THEN
	543	c ........................................................
	544	c .............. schema matsuno + leapfrog ..............
	545	c ........................................................
	546
	547	IF(forward. OR. leapf) THEN
	548	itau= itau + 1
	549	c iday= day_ini+itau/day_step
	550	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	551	c IF(time.GT.1.) THEN
	552	c time = time-1.
	553	c iday = iday+1
	554	c ENDIF
	555	ENDIF
	556
	557
	558	IF( itau. EQ. itaufinp1 ) then
	559	if (flag_verif) then
	560	write(79,*) 'ucov',ucov
	561	write(80,*) 'vcov',vcov
	562	write(81,*) 'teta',teta
	563	write(82,*) 'ps',ps
	564	write(83,*) 'q',q
	565	WRITE(85,*) 'q1 = ',q(:,:,1)
	566	WRITE(86,*) 'q3 = ',q(:,:,3)
	567	endif
	568
	569	abort_message = 'Simulation finished'
	570
	571	call abort_gcm(modname,abort_message,0)
	572	ENDIF
	573	c-----------------------------------------------------------------------
	574	c ecriture du fichier histoire moyenne:
	575	c -------------------------------------
	576
	577	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	578	IF(itau.EQ.itaufin) THEN
	579	iav=1
	580	ELSE
	581	iav=0
	582	ENDIF
	583
	584	IF (ok_dynzon) THEN
	585	#ifdef CPP_IOIPSL
	586	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
	587	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
	588	#endif
	589	END IF
	590	IF (ok_dyn_ave) THEN
	591	#ifdef CPP_IOIPSL
	592	CALL writedynav(itau,vcov,
	593	& ucov,teta,pk,phi,q,masse,ps,phis)
	594	#endif
	595	ENDIF
	596
	597	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
	598
	599	c-----------------------------------------------------------------------
	600	c ecriture de la bande histoire:
	601	c ------------------------------
	602
	603	IF( MOD(itau,iecri).EQ.0) THEN
	604	! Ehouarn: output only during LF or Backward Matsuno
	605	if (leapf.or.(.not.leapf.and.(.not.forward))) then
	606	nbetat = nbetatdem
	607	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	608	unat=0.
	609	do l=1,llm
	610	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	611	vnat(:,l)=vcov(:,l)/cv(:)
	612	enddo
	613	#ifdef CPP_IOIPSL
	614	if (ok_dyn_ins) then
	615	! write(lunout,*) "leapfrog: call writehist, itau=",itau
	616	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	617	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	618	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
	619	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
	620	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
	621	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
	622	endif ! of if (ok_dyn_ins)
	623	#endif
	624	! For some Grads outputs of fields
	625	if (output_grads_dyn) then
	626	#include "write_grads_dyn.h"
	627	endif
	628	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
	629	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
	630
	631	IF(itau.EQ.itaufin) THEN
	632
	633
	634	if (planet_type.eq."earth") then
	635	! Write an Earth-format restart file
	636	CALL dynredem1("restart.nc",0.0,
	637	& vcov,ucov,teta,q,masse,ps)
	638	endif ! of if (planet_type.eq."earth")
	639
	640	CLOSE(99)
	641	ENDIF ! of IF (itau.EQ.itaufin)
	642
	643	c-----------------------------------------------------------------------
	644	c gestion de l'integration temporelle:
	645	c ------------------------------------
	646
	647	IF( MOD(itau,iperiod).EQ.0 ) THEN
	648	GO TO 1
	649	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	650
	651	IF( forward ) THEN
	652	c fin du pas forward et debut du pas backward
	653
	654	forward = .FALSE.
	655	leapf = .FALSE.
	656	GO TO 2
	657
	658	ELSE
	659	c fin du pas backward et debut du premier pas leapfrog
	660
	661	leapf = .TRUE.
	662	dt = 2.*dtvr
	663	GO TO 2
	664	END IF ! of IF (forward)
	665	ELSE
	666
	667	c ...... pas leapfrog .....
	668
	669	leapf = .TRUE.
	670	dt = 2.*dtvr
	671	GO TO 2
	672	END IF ! of IF (MOD(itau,iperiod).EQ.0)
	673	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
	674
	675	ELSE ! of IF (.not.purmats)
	676
	677	c ........................................................
	678	c .............. schema matsuno ...............
	679	c ........................................................
	680	IF( forward ) THEN
	681
	682	itau = itau + 1
	683	c iday = day_ini+itau/day_step
	684	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	685	c
	686	c IF(time.GT.1.) THEN
	687	c time = time-1.
	688	c iday = iday+1
	689	c ENDIF
	690
	691	forward = .FALSE.
	692	IF( itau. EQ. itaufinp1 ) then
	693	abort_message = 'Simulation finished'
	694	call abort_gcm(modname,abort_message,0)
	695	ENDIF
	696	GO TO 2
	697
	698	ELSE ! of IF(forward) i.e. backward step
	699
	700	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	701	IF(itau.EQ.itaufin) THEN
	702	iav=1
	703	ELSE
	704	iav=0
	705	ENDIF
	706
	707	IF (ok_dynzon) THEN
	708	#ifdef CPP_IOIPSL
	709	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
	710	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
	711	#endif
	712	ENDIF
	713	IF (ok_dyn_ave) THEN
	714	#ifdef CPP_IOIPSL
	715	CALL writedynav(itau,vcov,
	716	& ucov,teta,pk,phi,q,masse,ps,phis)
	717	#endif
	718	ENDIF
	719
	720	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
	721
	722	IF(MOD(itau,iecri ).EQ.0) THEN
	723	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
	724	nbetat = nbetatdem
	725	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	726	unat=0.
	727	do l=1,llm
	728	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	729	vnat(:,l)=vcov(:,l)/cv(:)
	730	enddo
	731	#ifdef CPP_IOIPSL
	732	if (ok_dyn_ins) then
	733	! write(lunout,*) "leapfrog: call writehist (b)",
	734	! & itau,iecri
	735	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	736	endif ! of if (ok_dyn_ins)
	737	#endif
	738	! For some Grads outputs
	739	if (output_grads_dyn) then
	740	#include "write_grads_dyn.h"
	741	endif
	742
	743	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
	744
	745	IF(itau.EQ.itaufin) THEN
	746	if (planet_type.eq."earth") then
	747	CALL dynredem1("restart.nc",0.0,
	748	& vcov,ucov,teta,q,masse,ps)
	749	endif ! of if (planet_type.eq."earth")
	750	ENDIF ! of IF(itau.EQ.itaufin)
	751
	752	forward = .TRUE.
	753	GO TO 1
	754
	755	ENDIF ! of IF (forward)
	756
	757	END IF ! of IF(.not.purmats)
	758
	759	STOP
	760	END

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

Download in other formats:

Original Format