Context Navigation

leapfrog.F @ 1750

Last change on this file since 1750 was 1676, checked in by aslmd, 12 years ago

Repare le probleme "un jour sans fin" des versions precedentes pour le modele generique. Pas d'impact sur modele terrestre.

NB: changements egalement faits sur dyn3dmem en prevision de son utilisation

Fix the bug "Groundhog Day" which was a problem for previous versions using generic physics. No impact whatsoever on terrestrial model.

NB: changes also done on dyn3dmem which is planned to be used soon

Author : AS equipe planeto

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

File size: 25.4 KB

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

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

Download in other formats:

Original Format