Context Navigation

leapfrog.F @ 3396

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

Some cleanup and fixing the possibility to output fields in the dynamics, on the dynamical grids.

CLEANUPS:

arch-PW6_VARGAS.fcm : add potentially benefic compiling options
removed obsolete "control.h" in dyn3d/dyn3dpar (module control_mod.F90 is used instead)

OUTPUTS in the dynamics (3 sets of files, one for each grid: scalar, u, v):

removed "com_io_dyn.h" common; use module "com_io_dyn_mod.F90" instead
updated "initdynav.F","inithist.F","writehist.F" and "writedynav.F" in bibio: which field will be written is hard coded there.
flags "ok_dyn_ins" and "ok_dyn_ave" (loaded via conf_gcm.F) trigger output of fields in the dynamics: if ok_dyn_ins is true, then files "dyn_hist.nc", "dyn_histu.nc" and "dyn_histv.nc" are written (the frequency of the outputs is given by 'iecri' in run.def; values are written every 'iecri' dynamical step). if ok_dyn_ave is true then files "dyn_hist_ave.nc", "dyn_histu_ave.nc" and "dyn_histv_ave.nc" are written (the rate at which averages and made/written, in days, is given by 'periodav' in run.def).

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

File size: 23.1 KB

Rev	Line
[524]	1	!
[1279]	2	! $Id: leapfrog.F 1357 2010-04-14 14:03:19Z idelkadi $
	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
[1299]	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"
[1357]	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
[1357]	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
	211	dq=0.
	212	CALL pression ( ip1jmp1, ap, bp, ps, p )
	213	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	214
	215	c-----------------------------------------------------------------------
	216	c Debut de l'integration temporelle:
	217	c ----------------------------------
	218
	219	1 CONTINUE
	220
[1279]	221	jD_cur = jD_ref + day_ini - day_ref + int (itau * dtvr / daysec)
	222	jH_cur = jH_ref + &
	223	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
[524]	224
[1279]	225
[524]	226	#ifdef CPP_IOIPSL
[1279]	227	if (ok_guide) then
	228	call guide_main(itau,ucov,vcov,teta,q,masse,ps)
[524]	229	endif
	230	#endif
[1279]	231
	232
[524]	233	c
	234	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
	235	c CALL test_period ( ucov,vcov,teta,q,p,phis )
	236	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
	237	c ENDIF
	238	c
[1146]	239
	240	! Save fields obtained at previous time step as '...m1'
[524]	241	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
	242	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
	243	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
	244	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
	245	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
	246
	247	forward = .TRUE.
	248	leapf = .FALSE.
	249	dt = dtvr
	250
	251	c ... P.Le Van .26/04/94 ....
	252
	253	CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
	254	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
	255
[1146]	256	! Ehouarn: what is this for? zqmin & zqmax are not used anyway ...
	257	! call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
[524]	258
	259	2 CONTINUE
	260
	261	c-----------------------------------------------------------------------
	262
	263	c date:
	264	c -----
	265
	266
	267	c gestion des appels de la physique et des dissipations:
	268	c ------------------------------------------------------
	269	c
	270	c ... P.Le Van ( 6/02/95 ) ....
	271
	272	apphys = .FALSE.
	273	statcl = .FALSE.
	274	conser = .FALSE.
	275	apdiss = .FALSE.
	276
	277	IF( purmats ) THEN
[1357]	278	! Purely Matsuno time stepping
[524]	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
[1060]	282	s .and. iflag_phys.EQ.1 ) apphys = .TRUE.
[524]	283	ELSE
[1357]	284	! Leapfrog/Matsuno time stepping
[524]	285	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
	286	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
[1060]	287	IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
[524]	288	END IF
	289
[1357]	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
[524]	296	c-----------------------------------------------------------------------
	297	c calcul des tendances dynamiques:
	298	c --------------------------------
	299
	300	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	301
[1279]	302	time = jD_cur + jH_cur
[524]	303	CALL caldyn
	304	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
[1279]	305	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
[524]	306
[1279]	307
[524]	308	c-----------------------------------------------------------------------
	309	c calcul des tendances advection des traceurs (dont l'humidite)
	310	c -------------------------------------------------------------
	311
	312	IF( forward. OR . leapf ) THEN
	313
[960]	314	CALL caladvtrac(q,pbaru,pbarv,
	315	* p, masse, dq, teta,
	316	. flxw, pk)
	317
[524]	318	IF (offline) THEN
	319	Cmaf stokage du flux de masse pour traceurs OFF-LINE
	320
	321	#ifdef CPP_IOIPSL
[541]	322	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
	323	. dtvr, itau)
[524]	324	#endif
	325
	326
[1146]	327	ENDIF ! of IF (offline)
[524]	328	c
[1146]	329	ENDIF ! of IF( forward. OR . leapf )
[524]	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
[1279]	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
[524]	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
[1146]	384	IF (ip_ebil_dyn.ge.1 ) THEN
[524]	385	ztit='bil dyn'
[1146]	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 )
[524]	392	c-jld
[1146]	393	#ifdef CPP_IOIPSL
[692]	394	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
[1146]	395	IF (first) THEN
	396	first=.false.
[692]	397	#include "ini_paramLMDZ_dyn.h"
[1146]	398	ENDIF
[692]	399	c
	400	#include "write_paramLMDZ_dyn.h"
	401	c
[1146]	402	#endif
	403	! #endif of #ifdef CPP_IOIPSL
[1279]	404	CALL calfis( lafin , jD_cur, jH_cur,
[524]	405	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
[960]	406	$ du,dv,dteta,dq,
[524]	407	$ flxw,
	408	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
	409
[1146]	410	IF (ok_strato) THEN
[1279]	411	CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
[1146]	412	ENDIF
[999]	413
[524]	414	c ajout des tendances physiques:
	415	c ------------------------------
[1146]	416	CALL addfi( dtphys, leapf, forward ,
[524]	417	$ ucov, vcov, teta , q ,ps ,
	418	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	419	c
	420	c Diagnostique de conservation de l'énergie : difference
[1146]	421	IF (ip_ebil_dyn.ge.1 ) THEN
[524]	422	ztit='bil phys'
[1146]	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
[1060]	429	ENDIF ! of IF( apphys )
[524]	430
[1146]	431	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
[1060]	432	c Calcul academique de la physique = Rappel Newtonien + friction
[524]	433	c --------------------------------------------------------------
	434	teta(:,:)=teta(:,:)
	435	s -iphysiqdtvr(teta(:,:)-tetarappel(:,:))/taurappel
	436	call friction(ucov,vcov,iphysiq*dtvr)
[1060]	437	ENDIF
[524]	438
	439
	440	c-jld
	441
	442	CALL pression ( ip1jmp1, ap, bp, ps, p )
	443	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	444
	445
	446	c-----------------------------------------------------------------------
	447	c dissipation horizontale et verticale des petites echelles:
	448	c ----------------------------------------------------------
	449
	450	IF(apdiss) THEN
	451
	452
	453	c calcul de l'energie cinetique avant dissipation
	454	call covcont(llm,ucov,vcov,ucont,vcont)
	455	call enercin(vcov,ucov,vcont,ucont,ecin0)
	456
	457	c dissipation
	458	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
	459	ucov=ucov+dudis
	460	vcov=vcov+dvdis
	461	c teta=teta+dtetadis
	462
	463
	464	c------------------------------------------------------------------------
	465	if (dissip_conservative) then
	466	C On rajoute la tendance due a la transform. Ec -> E therm. cree
	467	C lors de la dissipation
	468	call covcont(llm,ucov,vcov,ucont,vcont)
	469	call enercin(vcov,ucov,vcont,ucont,ecin)
	470	dtetaecdt= (ecin0-ecin)/ pk
	471	c teta=teta+dtetaecdt
	472	dtetadis=dtetadis+dtetaecdt
	473	endif
	474	teta=teta+dtetadis
	475	c------------------------------------------------------------------------
	476
	477
	478	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
	479	c ... Calcul de la valeur moyenne, unique de h aux poles .....
	480	c
	481
	482	DO l = 1, llm
	483	DO ij = 1,iim
	484	tppn(ij) = aire( ij ) * teta( ij ,l)
	485	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	486	ENDDO
	487	tpn = SSUM(iim,tppn,1)/apoln
	488	tps = SSUM(iim,tpps,1)/apols
	489
	490	DO ij = 1, iip1
	491	teta( ij ,l) = tpn
	492	teta(ij+ip1jm,l) = tps
	493	ENDDO
	494	ENDDO
	495
	496	DO ij = 1,iim
	497	tppn(ij) = aire( ij ) * ps ( ij )
	498	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
	499	ENDDO
	500	tpn = SSUM(iim,tppn,1)/apoln
	501	tps = SSUM(iim,tpps,1)/apols
	502
	503	DO ij = 1, iip1
	504	ps( ij ) = tpn
	505	ps(ij+ip1jm) = tps
	506	ENDDO
	507
	508
[1146]	509	END IF ! of IF(apdiss)
[524]	510
	511	c ajout debug
	512	c IF( lafin ) then
	513	c abort_message = 'Simulation finished'
	514	c call abort_gcm(modname,abort_message,0)
	515	c ENDIF
	516
	517	c ********************************************************************
	518	c ********************************************************************
	519	c .... fin de l'integration dynamique et physique pour le pas itau ..
	520	c ********************************************************************
	521	c ********************************************************************
	522
	523	c preparation du pas d'integration suivant ......
	524
	525	IF ( .NOT.purmats ) THEN
	526	c ........................................................
	527	c .............. schema matsuno + leapfrog ..............
	528	c ........................................................
	529
	530	IF(forward. OR. leapf) THEN
	531	itau= itau + 1
[1357]	532	c iday= day_ini+itau/day_step
	533	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	534	c IF(time.GT.1.) THEN
	535	c time = time-1.
	536	c iday = iday+1
	537	c ENDIF
[524]	538	ENDIF
	539
	540
	541	IF( itau. EQ. itaufinp1 ) then
[999]	542	if (flag_verif) then
[1279]	543	write(79,*) 'ucov',ucov
	544	write(80,*) 'vcov',vcov
	545	write(81,*) 'teta',teta
	546	write(82,*) 'ps',ps
	547	write(83,*) 'q',q
[999]	548	WRITE(85,*) 'q1 = ',q(:,:,1)
	549	WRITE(86,*) 'q3 = ',q(:,:,3)
	550	endif
[524]	551
	552	abort_message = 'Simulation finished'
	553
	554	call abort_gcm(modname,abort_message,0)
	555	ENDIF
	556	c-----------------------------------------------------------------------
	557	c ecriture du fichier histoire moyenne:
	558	c -------------------------------------
	559
	560	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	561	IF(itau.EQ.itaufin) THEN
	562	iav=1
	563	ELSE
	564	iav=0
	565	ENDIF
[1146]	566
	567	IF (ok_dynzon) THEN
[524]	568	#ifdef CPP_IOIPSL
[1357]	569	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
	570	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[524]	571	#endif
[1146]	572	END IF
[1357]	573	IF (ok_dyn_ave) THEN
	574	#ifdef CPP_IOIPSL
	575	CALL writedynav(itau,vcov,
	576	& ucov,teta,pk,phi,q,masse,ps,phis)
	577	#endif
	578	ENDIF
[524]	579
[1357]	580	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
[524]	581
	582	c-----------------------------------------------------------------------
	583	c ecriture de la bande histoire:
	584	c ------------------------------
	585
[1357]	586	IF( MOD(itau,iecri).EQ.0) THEN
	587	! Ehouarn: output only during LF or Backward Matsuno
	588	if (leapf.or.(.not.leapf.and.(.not.forward))) then
[1146]	589	nbetat = nbetatdem
	590	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	591	unat=0.
	592	do l=1,llm
	593	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	594	vnat(:,l)=vcov(:,l)/cv(:)
	595	enddo
[524]	596	#ifdef CPP_IOIPSL
[1357]	597	if (ok_dyn_ins) then
	598	! write(lunout,*) "leapfrog: call writehist, itau=",itau
	599	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	600	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	601	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
	602	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
	603	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
	604	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
	605	endif ! of if (ok_dyn_ins)
[1146]	606	#endif
	607	! For some Grads outputs of fields
[1357]	608	if (output_grads_dyn) then
[524]	609	#include "write_grads_dyn.h"
[1357]	610	endif
	611	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
[1146]	612	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
[524]	613
	614	IF(itau.EQ.itaufin) THEN
	615
	616
[1146]	617	if (planet_type.eq."earth") then
	618	! Write an Earth-format restart file
	619	CALL dynredem1("restart.nc",0.0,
	620	& vcov,ucov,teta,q,masse,ps)
	621	endif ! of if (planet_type.eq."earth")
[524]	622
	623	CLOSE(99)
[1146]	624	ENDIF ! of IF (itau.EQ.itaufin)
[524]	625
	626	c-----------------------------------------------------------------------
	627	c gestion de l'integration temporelle:
	628	c ------------------------------------
	629
	630	IF( MOD(itau,iperiod).EQ.0 ) THEN
	631	GO TO 1
	632	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	633
	634	IF( forward ) THEN
	635	c fin du pas forward et debut du pas backward
	636
	637	forward = .FALSE.
	638	leapf = .FALSE.
	639	GO TO 2
	640
	641	ELSE
	642	c fin du pas backward et debut du premier pas leapfrog
	643
	644	leapf = .TRUE.
	645	dt = 2.*dtvr
[1146]	646	GO TO 2
	647	END IF ! of IF (forward)
[524]	648	ELSE
	649
	650	c ...... pas leapfrog .....
	651
	652	leapf = .TRUE.
	653	dt = 2.*dtvr
	654	GO TO 2
[1146]	655	END IF ! of IF (MOD(itau,iperiod).EQ.0)
	656	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
[524]	657
[1146]	658	ELSE ! of IF (.not.purmats)
[524]	659
	660	c ........................................................
	661	c .............. schema matsuno ...............
	662	c ........................................................
	663	IF( forward ) THEN
	664
	665	itau = itau + 1
[1357]	666	c iday = day_ini+itau/day_step
	667	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	668	c
	669	c IF(time.GT.1.) THEN
	670	c time = time-1.
	671	c iday = iday+1
	672	c ENDIF
[524]	673
	674	forward = .FALSE.
	675	IF( itau. EQ. itaufinp1 ) then
	676	abort_message = 'Simulation finished'
	677	call abort_gcm(modname,abort_message,0)
	678	ENDIF
	679	GO TO 2
	680
[1357]	681	ELSE ! of IF(forward) i.e. backward step
[524]	682
[1146]	683	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
[524]	684	IF(itau.EQ.itaufin) THEN
	685	iav=1
	686	ELSE
	687	iav=0
	688	ENDIF
[1146]	689
	690	IF (ok_dynzon) THEN
[524]	691	#ifdef CPP_IOIPSL
[1357]	692	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
	693	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[524]	694	#endif
[1357]	695	ENDIF
	696	IF (ok_dyn_ave) THEN
	697	#ifdef CPP_IOIPSL
	698	CALL writedynav(itau,vcov,
	699	& ucov,teta,pk,phi,q,masse,ps,phis)
	700	#endif
	701	ENDIF
[524]	702
[1146]	703	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
[524]	704
[1146]	705	IF(MOD(itau,iecri ).EQ.0) THEN
[524]	706	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
[1146]	707	nbetat = nbetatdem
	708	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	709	unat=0.
	710	do l=1,llm
	711	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	712	vnat(:,l)=vcov(:,l)/cv(:)
	713	enddo
[524]	714	#ifdef CPP_IOIPSL
[1357]	715	if (ok_dyn_ins) then
	716	! write(lunout,*) "leapfrog: call writehist (b)",
	717	! & itau,iecri
	718	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	719	endif ! of if (ok_dyn_ins)
[1146]	720	#endif
	721	! For some Grads outputs
	722	if (output_grads_dyn) then
[524]	723	#include "write_grads_dyn.h"
[1146]	724	endif
[524]	725
[1146]	726	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
[524]	727
[1146]	728	IF(itau.EQ.itaufin) THEN
	729	if (planet_type.eq."earth") then
	730	CALL dynredem1("restart.nc",0.0,
	731	& vcov,ucov,teta,q,masse,ps)
	732	endif ! of if (planet_type.eq."earth")
	733	ENDIF ! of IF(itau.EQ.itaufin)
[524]	734
[1146]	735	forward = .TRUE.
	736	GO TO 1
[524]	737
[1146]	738	ENDIF ! of IF (forward)
[524]	739
[1146]	740	END IF ! of IF(.not.purmats)
[524]	741
	742	STOP
	743	END

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

Download in other formats:

Original Format