Context Navigation

leapfrog.F @ 1995

Last change on this file since 1995 was 1987, checked in by Ehouarn Millour, 11 years ago
Add updating pressure, mass and Exner function (ie: all variables which depend on surface pressure) after adding physics tendencies (which include a surface pressure tendency). Note that this change induces slight changes in GCM results with respect to previous svn version of the code, even if surface pressure tendency is zero (because of recomputation of polar values as an average over polar points on the dynamics grid). EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 26.2 KB

Rev	Line
[524]	1	!
[1279]	2	! $Id: leapfrog.F 1987 2014-02-24 15:05:47Z fairhead $
	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
[1987]	14	USE infotrac, ONLY: nqtot
[1279]	15	USE guide_mod, ONLY : guide_main
[1987]	16	USE write_field, ONLY: writefield
	17	USE control_mod, ONLY: nday, day_step, planet_type, offline,
	18	& iconser, iphysiq, iperiod, dissip_period,
	19	& iecri, ip_ebil_dyn, ok_dynzon, ok_dyn_ins,
	20	& periodav, ok_dyn_ave, output_grads_dyn
[524]	21	IMPLICIT NONE
	22
	23	c ...... Version du 10/01/98 ..........
	24
	25	c avec coordonnees verticales hybrides
	26	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
	27
	28	c=======================================================================
	29	c
	30	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
	31	c -------
	32	c
	33	c Objet:
	34	c ------
	35	c
	36	c GCM LMD nouvelle grille
	37	c
	38	c=======================================================================
	39	c
	40	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
	41	c et possibilite d'appeler une fonction f(y) a derivee tangente
	42	c hyperbolique a la place de la fonction a derivee sinusoidale.
	43
	44	c ... Possibilite de choisir le shema pour l'advection de
	45	c q , en modifiant iadv dans traceur.def (10/02) .
	46	c
	47	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
	48	c Pour Van-Leer iadv=10
	49	c
	50	c-----------------------------------------------------------------------
	51	c Declarations:
	52	c -------------
	53
	54	#include "dimensions.h"
	55	#include "paramet.h"
	56	#include "comconst.h"
	57	#include "comdissnew.h"
	58	#include "comvert.h"
	59	#include "comgeom.h"
	60	#include "logic.h"
	61	#include "temps.h"
	62	#include "ener.h"
	63	#include "description.h"
	64	#include "serre.h"
[1403]	65	!#include "com_io_dyn.h"
[524]	66	#include "iniprint.h"
	67	#include "academic.h"
	68
[956]	69	! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
	70	! #include "clesphys.h"
[692]	71
[1987]	72	INTEGER,PARAMETER :: longcles = 20
	73	REAL,INTENT(IN) :: clesphy0( longcles ) ! not used
	74	REAL,INTENT(IN) :: time_0 ! not used
[524]	75
[1987]	76	c dynamical variables:
	77	REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
	78	REAL,INTENT(INOUT) :: vcov(ip1jm,llm) ! meridional covariant wind
	79	REAL,INTENT(INOUT) :: teta(ip1jmp1,llm) ! potential temperature
	80	REAL,INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
	81	REAL,INTENT(INOUT) :: masse(ip1jmp1,llm) ! air mass
	82	REAL,INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
	83	REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
	84
	85	REAL p (ip1jmp1,llmp1 ) ! interlayer pressure
	86	REAL pks(ip1jmp1) ! exner at the surface
	87	REAL pk(ip1jmp1,llm) ! exner at mid-layer
	88	REAL pkf(ip1jmp1,llm) ! filtered exner at mid-layer
	89	REAL phi(ip1jmp1,llm) ! geopotential
	90	REAL w(ip1jmp1,llm) ! vertical velocity
	91
[524]	92	real zqmin,zqmax
	93
	94	c variables dynamiques intermediaire pour le transport
	95	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
	96
	97	c variables dynamiques au pas -1
	98	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
	99	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
	100	REAL massem1(ip1jmp1,llm)
	101
	102	c tendances dynamiques
	103	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
[1146]	104	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
[524]	105
	106	c tendances de la dissipation
	107	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
	108	REAL dtetadis(ip1jmp1,llm)
	109
	110	c tendances physiques
	111	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
[1146]	112	REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
[524]	113
	114	c variables pour le fichier histoire
	115	REAL dtav ! intervalle de temps elementaire
	116
	117	REAL tppn(iim),tpps(iim),tpn,tps
	118	c
	119	INTEGER itau,itaufinp1,iav
[1279]	120	! INTEGER iday ! jour julien
	121	REAL time
[524]	122
	123	REAL SSUM
[1616]	124	! REAL finvmaold(ip1jmp1,llm)
[524]	125
[566]	126	cym LOGICAL lafin
	127	LOGICAL :: lafin=.false.
[524]	128	INTEGER ij,iq,l
	129	INTEGER ik
	130
	131	real time_step, t_wrt, t_ops
	132
[1279]	133	! REAL rdayvrai,rdaym_ini
	134	! jD_cur: jour julien courant
	135	! jH_cur: heure julienne courante
	136	REAL :: jD_cur, jH_cur
	137	INTEGER :: an, mois, jour
	138	REAL :: secondes
	139
[524]	140	LOGICAL first,callinigrads
[692]	141	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
	142	save first
	143	data first/.true./
[1279]	144	real dt_cum
[692]	145	character*10 infile
	146	integer zan, tau0, thoriid
	147	integer nid_ctesGCM
	148	save nid_ctesGCM
	149	real degres
	150	real rlong(iip1), rlatg(jjp1)
	151	real zx_tmp_2d(iip1,jjp1)
	152	integer ndex2d(iip1*jjp1)
	153	logical ok_sync
	154	parameter (ok_sync = .true.)
[1529]	155	logical physic
[524]	156
	157	data callinigrads/.true./
	158	character*10 string10
	159
	160	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
[960]	161	REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
[524]	162
	163	c+jld variables test conservation energie
	164	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
	165	C Tendance de la temp. potentiel d (theta)/ d t due a la
	166	C tansformation d'energie cinetique en energie thermique
	167	C cree par la dissipation
	168	REAL dtetaecdt(ip1jmp1,llm)
	169	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
	170	REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
	171	REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
	172	CHARACTER*15 ztit
[692]	173	!IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
	174	!IM SAVE ip_ebil_dyn
	175	!IM DATA ip_ebil_dyn/0/
[524]	176	c-jld
	177
	178	character*80 dynhist_file, dynhistave_file
[1520]	179	character(len=*),parameter :: modname="leapfrog"
[524]	180	character*80 abort_message
	181
	182	logical dissip_conservative
	183	save dissip_conservative
	184	data dissip_conservative/.true./
	185
	186	LOGICAL prem
	187	save prem
	188	DATA prem/.true./
	189	INTEGER testita
	190	PARAMETER (testita = 9)
	191
[1286]	192	logical , parameter :: flag_verif = .false.
[999]	193
	194
[825]	195	integer itau_w ! pas de temps ecriture = itap + itau_phy
	196
	197
[524]	198	itaufin = nday*day_step
	199	itaufinp1 = itaufin +1
[1529]	200	itau = 0
	201	physic=.true.
	202	if (iflag_phys==0.or.iflag_phys==2) physic=.false.
[524]	203
[1403]	204	c iday = day_ini+itau/day_step
	205	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	206	c IF(time.GT.1.) THEN
	207	c time = time-1.
	208	c iday = iday+1
	209	c ENDIF
[524]	210
	211
	212	c-----------------------------------------------------------------------
	213	c On initialise la pression et la fonction d'Exner :
	214	c --------------------------------------------------
	215
[1454]	216	dq(:,:,:)=0.
[524]	217	CALL pression ( ip1jmp1, ap, bp, ps, p )
[1625]	218	if (pressure_exner) then
[1520]	219	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
[1625]	220	else
[1520]	221	CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
	222	endif
[524]	223
	224	c-----------------------------------------------------------------------
	225	c Debut de l'integration temporelle:
	226	c ----------------------------------
	227
[1614]	228	1 CONTINUE ! Matsuno Forward step begins here
[524]	229
[1577]	230	jD_cur = jD_ref + day_ini - day_ref + &
[1626]	231	& itau/day_step
[1577]	232	jH_cur = jH_ref + start_time + &
[1626]	233	& mod(itau,day_step)/float(day_step)
[1577]	234	jD_cur = jD_cur + int(jH_cur)
	235	jH_cur = jH_cur - int(jH_cur)
[524]	236
[1279]	237
[524]	238	#ifdef CPP_IOIPSL
[1279]	239	if (ok_guide) then
	240	call guide_main(itau,ucov,vcov,teta,q,masse,ps)
[524]	241	endif
	242	#endif
[1279]	243
	244
[524]	245	c
	246	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
	247	c CALL test_period ( ucov,vcov,teta,q,p,phis )
	248	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
	249	c ENDIF
	250	c
[1146]	251
	252	! Save fields obtained at previous time step as '...m1'
[524]	253	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
	254	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
	255	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
	256	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
	257	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
	258
	259	forward = .TRUE.
	260	leapf = .FALSE.
	261	dt = dtvr
	262
	263	c ... P.Le Van .26/04/94 ....
[1616]	264	! Ehouarn: finvmaold is actually not used
	265	! CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
	266	! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
[524]	267
[1614]	268	2 CONTINUE ! Matsuno backward or leapfrog step begins here
[524]	269
	270	c-----------------------------------------------------------------------
	271
	272	c date:
	273	c -----
	274
	275
	276	c gestion des appels de la physique et des dissipations:
	277	c ------------------------------------------------------
	278	c
	279	c ... P.Le Van ( 6/02/95 ) ....
	280
	281	apphys = .FALSE.
	282	statcl = .FALSE.
	283	conser = .FALSE.
	284	apdiss = .FALSE.
	285
	286	IF( purmats ) THEN
[1403]	287	! Purely Matsuno time stepping
[524]	288	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
[1502]	289	IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward )
	290	s apdiss = .TRUE.
[524]	291	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
[1529]	292	s .and. physic ) apphys = .TRUE.
[524]	293	ELSE
[1403]	294	! Leapfrog/Matsuno time stepping
[524]	295	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
[1502]	296	IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
	297	s apdiss = .TRUE.
[1529]	298	IF( MOD(itau+1,iphysiq).EQ.0.AND.physic ) apphys=.TRUE.
[524]	299	END IF
	300
[1403]	301	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
	302	! supress dissipation step
	303	if (llm.eq.1) then
	304	apdiss=.false.
	305	endif
	306
[524]	307	c-----------------------------------------------------------------------
	308	c calcul des tendances dynamiques:
	309	c --------------------------------
	310
[1614]	311	! compute geopotential phi()
[524]	312	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	313
[1279]	314	time = jD_cur + jH_cur
[524]	315	CALL caldyn
	316	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
[1279]	317	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
[524]	318
[1279]	319
[524]	320	c-----------------------------------------------------------------------
	321	c calcul des tendances advection des traceurs (dont l'humidite)
	322	c -------------------------------------------------------------
	323
	324	IF( forward. OR . leapf ) THEN
[1987]	325	! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
[960]	326	CALL caladvtrac(q,pbaru,pbarv,
	327	* p, masse, dq, teta,
	328	. flxw, pk)
	329
[524]	330	IF (offline) THEN
	331	Cmaf stokage du flux de masse pour traceurs OFF-LINE
	332
	333	#ifdef CPP_IOIPSL
[541]	334	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
	335	. dtvr, itau)
[524]	336	#endif
	337
	338
[1146]	339	ENDIF ! of IF (offline)
[524]	340	c
[1146]	341	ENDIF ! of IF( forward. OR . leapf )
[524]	342
	343
	344	c-----------------------------------------------------------------------
	345	c integrations dynamique et traceurs:
	346	c ----------------------------------
	347
	348
	349	CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
[1616]	350	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
	351	! $ finvmaold )
[524]	352
	353
	354	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
	355	c
	356	c-----------------------------------------------------------------------
	357	c calcul des tendances physiques:
	358	c -------------------------------
	359	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
	360	c
	361	IF( purmats ) THEN
	362	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
	363	ELSE
	364	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
	365	ENDIF
	366	c
	367	c
	368	IF( apphys ) THEN
	369	c
	370	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
	371	c
	372
	373	CALL pression ( ip1jmp1, ap, bp, ps, p )
[1625]	374	if (pressure_exner) then
[1520]	375	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
[1625]	376	else
[1520]	377	CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
	378	endif
[524]	379
[1279]	380	! rdaym_ini = itau * dtvr / daysec
	381	! rdayvrai = rdaym_ini + day_ini
[1577]	382	! jD_cur = jD_ref + day_ini - day_ref
	383	! $ + int (itau * dtvr / daysec)
	384	! jH_cur = jH_ref + &
	385	! & (itau * dtvr / daysec - int(itau * dtvr / daysec))
	386	jD_cur = jD_ref + day_ini - day_ref + &
[1626]	387	& itau/day_step
[1676]	388
	389	IF (planet_type .eq."generic") THEN
	390	! AS: we make jD_cur to be pday
	391	jD_cur = int(day_ini + itau/day_step)
	392	ENDIF
	393
[1577]	394	jH_cur = jH_ref + start_time + &
[1626]	395	& mod(itau,day_step)/float(day_step)
[1577]	396	jD_cur = jD_cur + int(jH_cur)
	397	jH_cur = jH_cur - int(jH_cur)
[1279]	398	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
	399	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
	400	! write(lunout,*)'current date = ',an, mois, jour, secondes
[524]	401
	402	c rajout debug
	403	c lafin = .true.
	404
	405
	406	c Inbterface avec les routines de phylmd (phymars ... )
	407	c -----------------------------------------------------
	408
	409	c+jld
	410
	411	c Diagnostique de conservation de l'énergie : initialisation
[1146]	412	IF (ip_ebil_dyn.ge.1 ) THEN
[524]	413	ztit='bil dyn'
[1146]	414	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
	415	IF (planet_type.eq."earth") THEN
[1615]	416	#ifdef CPP_EARTH
[1146]	417	CALL diagedyn(ztit,2,1,1,dtphys
	418	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
[1615]	419	#endif
[1146]	420	ENDIF
	421	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
[524]	422	c-jld
[1146]	423	#ifdef CPP_IOIPSL
[1656]	424	cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
	425	cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
	426	c IF (first) THEN
	427	c first=.false.
	428	c#include "ini_paramLMDZ_dyn.h"
	429	c ENDIF
[692]	430	c
[1656]	431	c#include "write_paramLMDZ_dyn.h"
[692]	432	c
[1146]	433	#endif
	434	! #endif of #ifdef CPP_IOIPSL
[1279]	435	CALL calfis( lafin , jD_cur, jH_cur,
[524]	436	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
[960]	437	$ du,dv,dteta,dq,
[524]	438	$ flxw,
	439	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
	440
	441	c ajout des tendances physiques:
	442	c ------------------------------
[1146]	443	CALL addfi( dtphys, leapf, forward ,
[524]	444	$ ucov, vcov, teta , q ,ps ,
	445	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
[1987]	446	! since addfi updates ps(), also update p(), masse() and pk()
	447	CALL pression (ip1jmp1,ap,bp,ps,p)
	448	CALL massdair(p,masse)
	449	if (pressure_exner) then
	450	CALL exner_hyb(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
	451	else
	452	CALL exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
	453	endif
[1793]	454
	455	IF (ok_strato) THEN
	456	CALL top_bound( vcov,ucov,teta,masse,dtphys)
	457	ENDIF
	458
[524]	459	c
	460	c Diagnostique de conservation de l'énergie : difference
[1146]	461	IF (ip_ebil_dyn.ge.1 ) THEN
[524]	462	ztit='bil phys'
[1146]	463	IF (planet_type.eq."earth") THEN
	464	CALL diagedyn(ztit,2,1,1,dtphys
	465	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	466	ENDIF
	467	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
	468
[1060]	469	ENDIF ! of IF( apphys )
[524]	470
[1146]	471	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
[1454]	472	! Academic case : Simple friction and Newtonan relaxation
	473	! -------------------------------------------------------
	474	DO l=1,llm
	475	DO ij=1,ip1jmp1
	476	teta(ij,l)=teta(ij,l)-dtvr*
	477	& (teta(ij,l)-tetarappel(ij,l))(knewt_g+knewt_t(l)clat4(ij))
	478	ENDDO
	479	ENDDO ! of DO l=1,llm
	480
[1505]	481	if (planet_type.eq."giant") then
	482	! add an intrinsic heat flux at the base of the atmosphere
	483	teta(:,1)=teta(:,1)+dtvraire(:)ihf/cpp/masse(:,1)
	484	endif
	485
	486	call friction(ucov,vcov,dtvr)
	487
[1454]	488	! Sponge layer (if any)
	489	IF (ok_strato) THEN
[1793]	490	! dufi(:,:)=0.
	491	! dvfi(:,:)=0.
	492	! dtetafi(:,:)=0.
	493	! dqfi(:,:,:)=0.
	494	! dpfi(:)=0.
	495	! CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
	496	CALL top_bound( vcov,ucov,teta,masse,dtvr)
	497	! CALL addfi( dtvr, leapf, forward ,
	498	! $ ucov, vcov, teta , q ,ps ,
	499	! $ dufi, dvfi, dtetafi , dqfi ,dpfi )
[1454]	500	ENDIF ! of IF (ok_strato)
	501	ENDIF ! of IF (iflag_phys.EQ.2)
[524]	502
	503
	504	c-jld
	505
	506	CALL pression ( ip1jmp1, ap, bp, ps, p )
[1625]	507	if (pressure_exner) then
[1520]	508	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
[1625]	509	else
[1520]	510	CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
	511	endif
[1987]	512	CALL massdair(p,masse)
[524]	513
	514
	515	c-----------------------------------------------------------------------
	516	c dissipation horizontale et verticale des petites echelles:
	517	c ----------------------------------------------------------
	518
	519	IF(apdiss) THEN
	520
	521
	522	c calcul de l'energie cinetique avant dissipation
	523	call covcont(llm,ucov,vcov,ucont,vcont)
	524	call enercin(vcov,ucov,vcont,ucont,ecin0)
	525
	526	c dissipation
	527	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
	528	ucov=ucov+dudis
	529	vcov=vcov+dvdis
	530	c teta=teta+dtetadis
	531
	532
	533	c------------------------------------------------------------------------
	534	if (dissip_conservative) then
	535	C On rajoute la tendance due a la transform. Ec -> E therm. cree
	536	C lors de la dissipation
	537	call covcont(llm,ucov,vcov,ucont,vcont)
	538	call enercin(vcov,ucov,vcont,ucont,ecin)
	539	dtetaecdt= (ecin0-ecin)/ pk
	540	c teta=teta+dtetaecdt
	541	dtetadis=dtetadis+dtetaecdt
	542	endif
	543	teta=teta+dtetadis
	544	c------------------------------------------------------------------------
	545
	546
	547	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
	548	c ... Calcul de la valeur moyenne, unique de h aux poles .....
	549	c
	550
	551	DO l = 1, llm
	552	DO ij = 1,iim
	553	tppn(ij) = aire( ij ) * teta( ij ,l)
	554	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	555	ENDDO
	556	tpn = SSUM(iim,tppn,1)/apoln
	557	tps = SSUM(iim,tpps,1)/apols
	558
	559	DO ij = 1, iip1
	560	teta( ij ,l) = tpn
	561	teta(ij+ip1jm,l) = tps
	562	ENDDO
	563	ENDDO
	564
[1614]	565	if (1 == 0) then
	566	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
	567	!!! 2) should probably not be here anyway
	568	!!! but are kept for those who would want to revert to previous behaviour
	569	DO ij = 1,iim
	570	tppn(ij) = aire( ij ) * ps ( ij )
	571	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
	572	ENDDO
	573	tpn = SSUM(iim,tppn,1)/apoln
	574	tps = SSUM(iim,tpps,1)/apols
[524]	575
[1614]	576	DO ij = 1, iip1
	577	ps( ij ) = tpn
	578	ps(ij+ip1jm) = tps
	579	ENDDO
	580	endif ! of if (1 == 0)
[524]	581
[1146]	582	END IF ! of IF(apdiss)
[524]	583
	584	c ajout debug
	585	c IF( lafin ) then
	586	c abort_message = 'Simulation finished'
	587	c call abort_gcm(modname,abort_message,0)
	588	c ENDIF
	589
	590	c ********************************************************************
	591	c ********************************************************************
	592	c .... fin de l'integration dynamique et physique pour le pas itau ..
	593	c ********************************************************************
	594	c ********************************************************************
	595
	596	c preparation du pas d'integration suivant ......
	597
	598	IF ( .NOT.purmats ) THEN
	599	c ........................................................
	600	c .............. schema matsuno + leapfrog ..............
	601	c ........................................................
	602
	603	IF(forward. OR. leapf) THEN
	604	itau= itau + 1
[1403]	605	c iday= day_ini+itau/day_step
	606	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	607	c IF(time.GT.1.) THEN
	608	c time = time-1.
	609	c iday = iday+1
	610	c ENDIF
[524]	611	ENDIF
	612
	613
	614	IF( itau. EQ. itaufinp1 ) then
[999]	615	if (flag_verif) then
[1279]	616	write(79,*) 'ucov',ucov
	617	write(80,*) 'vcov',vcov
	618	write(81,*) 'teta',teta
	619	write(82,*) 'ps',ps
	620	write(83,*) 'q',q
[999]	621	WRITE(85,*) 'q1 = ',q(:,:,1)
	622	WRITE(86,*) 'q3 = ',q(:,:,3)
	623	endif
[524]	624
	625	abort_message = 'Simulation finished'
	626
	627	call abort_gcm(modname,abort_message,0)
	628	ENDIF
	629	c-----------------------------------------------------------------------
	630	c ecriture du fichier histoire moyenne:
	631	c -------------------------------------
	632
	633	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	634	IF(itau.EQ.itaufin) THEN
	635	iav=1
	636	ELSE
	637	iav=0
	638	ENDIF
[1146]	639
	640	IF (ok_dynzon) THEN
[524]	641	#ifdef CPP_IOIPSL
[1403]	642	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
	643	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[524]	644	#endif
[1146]	645	END IF
[1403]	646	IF (ok_dyn_ave) THEN
	647	#ifdef CPP_IOIPSL
	648	CALL writedynav(itau,vcov,
	649	& ucov,teta,pk,phi,q,masse,ps,phis)
	650	#endif
	651	ENDIF
[524]	652
[1403]	653	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
[524]	654
	655	c-----------------------------------------------------------------------
	656	c ecriture de la bande histoire:
	657	c ------------------------------
	658
[1403]	659	IF( MOD(itau,iecri).EQ.0) THEN
	660	! Ehouarn: output only during LF or Backward Matsuno
	661	if (leapf.or.(.not.leapf.and.(.not.forward))) then
[1146]	662	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	663	unat=0.
	664	do l=1,llm
	665	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	666	vnat(:,l)=vcov(:,l)/cv(:)
	667	enddo
[524]	668	#ifdef CPP_IOIPSL
[1403]	669	if (ok_dyn_ins) then
	670	! write(lunout,*) "leapfrog: call writehist, itau=",itau
	671	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	672	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	673	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
	674	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
	675	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
	676	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
	677	endif ! of if (ok_dyn_ins)
[1146]	678	#endif
	679	! For some Grads outputs of fields
[1403]	680	if (output_grads_dyn) then
[524]	681	#include "write_grads_dyn.h"
[1403]	682	endif
	683	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
[1146]	684	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
[524]	685
	686	IF(itau.EQ.itaufin) THEN
	687
	688
[1454]	689	! if (planet_type.eq."earth") then
[1146]	690	! Write an Earth-format restart file
[1577]	691	CALL dynredem1("restart.nc",start_time,
[1146]	692	& vcov,ucov,teta,q,masse,ps)
[1454]	693	! endif ! of if (planet_type.eq."earth")
[524]	694
	695	CLOSE(99)
[1614]	696	!!! Ehouarn: Why not stop here and now?
[1146]	697	ENDIF ! of IF (itau.EQ.itaufin)
[524]	698
	699	c-----------------------------------------------------------------------
	700	c gestion de l'integration temporelle:
	701	c ------------------------------------
	702
	703	IF( MOD(itau,iperiod).EQ.0 ) THEN
	704	GO TO 1
	705	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	706
	707	IF( forward ) THEN
	708	c fin du pas forward et debut du pas backward
	709
	710	forward = .FALSE.
	711	leapf = .FALSE.
	712	GO TO 2
	713
	714	ELSE
	715	c fin du pas backward et debut du premier pas leapfrog
	716
	717	leapf = .TRUE.
	718	dt = 2.*dtvr
[1146]	719	GO TO 2
	720	END IF ! of IF (forward)
[524]	721	ELSE
	722
	723	c ...... pas leapfrog .....
	724
	725	leapf = .TRUE.
	726	dt = 2.*dtvr
	727	GO TO 2
[1146]	728	END IF ! of IF (MOD(itau,iperiod).EQ.0)
	729	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
[524]	730
[1146]	731	ELSE ! of IF (.not.purmats)
[524]	732
	733	c ........................................................
	734	c .............. schema matsuno ...............
	735	c ........................................................
	736	IF( forward ) THEN
	737
	738	itau = itau + 1
[1403]	739	c iday = day_ini+itau/day_step
	740	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	741	c
	742	c IF(time.GT.1.) THEN
	743	c time = time-1.
	744	c iday = iday+1
	745	c ENDIF
[524]	746
	747	forward = .FALSE.
	748	IF( itau. EQ. itaufinp1 ) then
	749	abort_message = 'Simulation finished'
	750	call abort_gcm(modname,abort_message,0)
	751	ENDIF
	752	GO TO 2
	753
[1403]	754	ELSE ! of IF(forward) i.e. backward step
[524]	755
[1146]	756	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
[524]	757	IF(itau.EQ.itaufin) THEN
	758	iav=1
	759	ELSE
	760	iav=0
	761	ENDIF
[1146]	762
	763	IF (ok_dynzon) THEN
[524]	764	#ifdef CPP_IOIPSL
[1403]	765	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
	766	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[524]	767	#endif
[1403]	768	ENDIF
	769	IF (ok_dyn_ave) THEN
	770	#ifdef CPP_IOIPSL
	771	CALL writedynav(itau,vcov,
	772	& ucov,teta,pk,phi,q,masse,ps,phis)
	773	#endif
	774	ENDIF
[524]	775
[1146]	776	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
[524]	777
[1146]	778	IF(MOD(itau,iecri ).EQ.0) THEN
[524]	779	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
[1146]	780	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	781	unat=0.
	782	do l=1,llm
	783	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	784	vnat(:,l)=vcov(:,l)/cv(:)
	785	enddo
[524]	786	#ifdef CPP_IOIPSL
[1403]	787	if (ok_dyn_ins) then
	788	! write(lunout,*) "leapfrog: call writehist (b)",
	789	! & itau,iecri
	790	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	791	endif ! of if (ok_dyn_ins)
[1146]	792	#endif
	793	! For some Grads outputs
	794	if (output_grads_dyn) then
[524]	795	#include "write_grads_dyn.h"
[1146]	796	endif
[524]	797
[1146]	798	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
[524]	799
[1146]	800	IF(itau.EQ.itaufin) THEN
[1454]	801	! if (planet_type.eq."earth") then
[1577]	802	CALL dynredem1("restart.nc",start_time,
[1146]	803	& vcov,ucov,teta,q,masse,ps)
[1454]	804	! endif ! of if (planet_type.eq."earth")
[1146]	805	ENDIF ! of IF(itau.EQ.itaufin)
[524]	806
[1146]	807	forward = .TRUE.
	808	GO TO 1
[524]	809
[1146]	810	ENDIF ! of IF (forward)
[524]	811
[1146]	812	END IF ! of IF(.not.purmats)
[524]	813
	814	STOP
	815	END

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

Download in other formats:

Original Format