Context Navigation

leapfrog.F @ 2125

Last change on this file since 2125 was 1824, checked in by emillour, 8 years ago

Common dynamics:

enable possiblity to store multiple time steps in the restart.nc file (flag "ecritstart" gives the frequency, in dynamical steps).
fixed dynredem_mod.F90 to correctly write multiple time steps.
fixed computation of JH_cur in the mars case where "hour_ini" contains the initial time of day read from the start.nc file
minor fix in dynetat0.F90

File size: 34.6 KB

Rev	Line
[1]	1	!
[7]	2	! $Id: leapfrog.F 1446 2010-10-22 09:27:25Z emillour $
[1]	3	!
	4	c
	5	c
[97]	6	SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,
[1]	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
[1508]	14	USE infotrac, ONLY: nqtot,ok_iso_verif
[1]	15	USE guide_mod, ONLY : guide_main
[1189]	16	USE write_field, ONLY: writefield
	17	USE control_mod, ONLY: planet_type,nday,day_step,iperiod,iphysiq,
[1022]	18	& less1day,fractday,ndynstep,iconser,
	19	& dissip_period,offline,ip_ebil_dyn,
	20	& ok_dynzon,periodav,ok_dyn_ave,iecri,
[1824]	21	& ok_dyn_ins,output_grads_dyn,ecritstart
[1302]	22	use exner_hyb_m, only: exner_hyb
	23	use exner_milieu_m, only: exner_milieu
[1017]	24	use cpdet_mod, only: cpdet,tpot2t,t2tpot
	25	use sponge_mod, only: callsponge,mode_sponge,sponge
[1024]	26	use comuforc_h
[1422]	27	USE comvert_mod, ONLY: ap,bp,pressure_exner,presnivs
	28	USE comconst_mod, ONLY: daysec,dtvr,dtphys,dtdiss,
	29	. cpp,ihf,iflag_top_bound,pi
	30	USE logic_mod, ONLY: iflag_phys,ok_guide,forward,leapf,apphys,
	31	. statcl,conser,apdiss,purmats,tidal,ok_strato
	32	USE temps_mod, ONLY: jD_ref,jH_ref,itaufin,day_ini,day_ref,
[1824]	33	. start_time,dt,hour_ini
[1024]	34
[1]	35	IMPLICIT NONE
	36
	37	c ...... Version du 10/01/98 ..........
	38
	39	c avec coordonnees verticales hybrides
	40	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
	41
	42	c=======================================================================
	43	c
	44	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
	45	c -------
	46	c
	47	c Objet:
	48	c ------
	49	c
	50	c GCM LMD nouvelle grille
	51	c
	52	c=======================================================================
	53	c
	54	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
	55	c et possibilite d'appeler une fonction f(y) a derivee tangente
	56	c hyperbolique a la place de la fonction a derivee sinusoidale.
	57
	58	c ... Possibilite de choisir le shema pour l'advection de
	59	c q , en modifiant iadv dans traceur.def (10/02) .
	60	c
	61	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
	62	c Pour Van-Leer iadv=10
	63	c
	64	c-----------------------------------------------------------------------
	65	c Declarations:
	66	c -------------
	67
	68	#include "dimensions.h"
	69	#include "paramet.h"
	70	#include "comdissnew.h"
	71	#include "comgeom.h"
	72	!#include "com_io_dyn.h"
	73	#include "iniprint.h"
	74	#include "academic.h"
	75
	76	! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
	77	! #include "clesphys.h"
	78
[1189]	79	REAL,INTENT(IN) :: time_0 ! not used
[1]	80
[1189]	81	c dynamical variables:
	82	REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
	83	REAL,INTENT(INOUT) :: vcov(ip1jm,llm) ! meridional covariant wind
	84	REAL,INTENT(INOUT) :: teta(ip1jmp1,llm) ! potential temperature
	85	REAL,INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
	86	REAL,INTENT(INOUT) :: masse(ip1jmp1,llm) ! air mass
	87	REAL,INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
	88	REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
	89
	90	REAL p (ip1jmp1,llmp1 ) ! interlayer pressure
	91	REAL pks(ip1jmp1) ! exner at the surface
	92	REAL pk(ip1jmp1,llm) ! exner at mid-layer
	93	REAL pkf(ip1jmp1,llm) ! filtered exner at mid-layer
	94	REAL phi(ip1jmp1,llm) ! geopotential
	95	REAL w(ip1jmp1,llm) ! vertical velocity
[5]	96	! ADAPTATION GCM POUR CP(T)
	97	REAL temp(ip1jmp1,llm) ! temperature
	98	REAL tsurpk(ip1jmp1,llm) ! cpp*T/pk
[1]	99
[1572]	100	! real zqmin,zqmax
[1189]	101
[1]	102	c variables dynamiques intermediaire pour le transport
	103	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
	104
	105	c variables dynamiques au pas -1
	106	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
	107	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
	108	REAL massem1(ip1jmp1,llm)
	109
[6]	110	c tendances dynamiques en */s
[1]	111	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
	112	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
	113
[6]	114	c tendances de la dissipation en */s
[1]	115	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
	116	REAL dtetadis(ip1jmp1,llm)
	117
[6]	118	c tendances de la couche superieure */s
[1010]	119	c REAL dvtop(ip1jm,llm)
	120	REAL dutop(ip1jmp1,llm)
	121	c REAL dtetatop(ip1jmp1,llm)
	122	c REAL dqtop(ip1jmp1,llm,nqtot),dptop(ip1jmp1)
[6]	123
[495]	124	c TITAN : tendances due au forces de marees */s
	125	REAL dvtidal(ip1jm,llm),dutidal(ip1jmp1,llm)
	126
[6]	127	c tendances physiques */s
[1]	128	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
	129	REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
	130
	131	c variables pour le fichier histoire
[1572]	132	! REAL dtav ! intervalle de temps elementaire
[1824]	133	LOGICAL lrestart
[1]	134
	135	REAL tppn(iim),tpps(iim),tpn,tps
	136	c
	137	INTEGER itau,itaufinp1,iav
	138	! INTEGER iday ! jour julien
	139	REAL time
	140
	141	REAL SSUM
[776]	142	! REAL finvmaold(ip1jmp1,llm)
[1]	143
	144	cym LOGICAL lafin
	145	LOGICAL :: lafin=.false.
	146	INTEGER ij,iq,l
[1572]	147	! INTEGER ik
[1]	148
[1572]	149	! real time_step, t_wrt, t_ops
[1]	150
[497]	151	REAL rdaym_ini
[1]	152	! jD_cur: jour julien courant
	153	! jH_cur: heure julienne courante
	154	REAL :: jD_cur, jH_cur
[1572]	155	! INTEGER :: an, mois, jour
	156	! REAL :: secondes
[1]	157
	158	LOGICAL first,callinigrads
	159	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
	160	save first
	161	data first/.true./
[1572]	162	! real dt_cum
	163	! character*10 infile
	164	! integer zan, tau0, thoriid
	165	! integer nid_ctesGCM
	166	! save nid_ctesGCM
	167	! real degres
	168	! real rlong(iip1), rlatg(jjp1)
	169	! real zx_tmp_2d(iip1,jjp1)
	170	! integer ndex2d(iip1*jjp1)
[1]	171	logical ok_sync
	172	parameter (ok_sync = .true.)
[1403]	173	logical physics
[1]	174
	175	data callinigrads/.true./
	176	character*10 string10
	177
[1572]	178	! REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
[1]	179	REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
	180
	181	c+jld variables test conservation energie
	182	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
	183	C Tendance de la temp. potentiel d (theta)/ d t due a la
	184	C tansformation d'energie cinetique en energie thermique
	185	C cree par la dissipation
	186	REAL dtetaecdt(ip1jmp1,llm)
	187	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
	188	REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
[1572]	189	! REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
[1]	190	CHARACTER*15 ztit
	191	!IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
	192	!IM SAVE ip_ebil_dyn
	193	!IM DATA ip_ebil_dyn/0/
	194	c-jld
	195
[1572]	196	! integer :: itau_w ! for write_paramLMDZ_dyn.h
[37]	197
[1572]	198	! character*80 dynhist_file, dynhistave_file
[127]	199	character(len=*),parameter :: modname="leapfrog"
[1]	200	character*80 abort_message
	201
	202	logical dissip_conservative
	203	save dissip_conservative
	204	data dissip_conservative/.true./
	205
	206	INTEGER testita
	207	PARAMETER (testita = 9)
	208
	209	logical , parameter :: flag_verif = .false.
	210
[37]	211	! for CP(T)
	212	real :: dtec
	213	real :: ztetaec(ip1jmp1,llm)
[1]	214
[1302]	215	if (nday>=0) then
	216	itaufin = nday*day_step
	217	else
	218	! to run a given (-nday) number of dynamical steps
	219	itaufin = -nday
	220	endif
[97]	221	if (less1day) then
	222	c MODIF VENUS: to run less than one day:
	223	itaufin = int(fractday*day_step)
	224	endif
[1022]	225	if (ndynstep.gt.0) then
	226	! running a given number of dynamical steps
	227	itaufin=ndynstep
	228	endif
[1]	229	itaufinp1 = itaufin +1
	230
[6]	231	c INITIALISATIONS
	232	dudis(:,:) =0.
	233	dvdis(:,:) =0.
	234	dtetadis(:,:)=0.
	235	dutop(:,:) =0.
[1010]	236	c dvtop(:,:) =0.
	237	c dtetatop(:,:)=0.
	238	c dqtop(:,:,:) =0.
	239	c dptop(:) =0.
[6]	240	dufi(:,:) =0.
	241	dvfi(:,:) =0.
	242	dtetafi(:,:)=0.
	243	dqfi(:,:,:) =0.
	244	dpfi(:) =0.
[1]	245
	246	itau = 0
[1403]	247	physics=.true.
	248	if (iflag_phys==0.or.iflag_phys==2) physics=.false.
[270]	249
[1]	250	c iday = day_ini+itau/day_step
	251	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	252	c IF(time.GT.1.) THEN
	253	c time = time-1.
	254	c iday = iday+1
	255	c ENDIF
	256
	257
	258	c-----------------------------------------------------------------------
	259	c On initialise la pression et la fonction d'Exner :
	260	c --------------------------------------------------
	261
[7]	262	dq(:,:,:)=0.
[1]	263	CALL pression ( ip1jmp1, ap, bp, ps, p )
[776]	264	if (pressure_exner) then
[1302]	265	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
[776]	266	else
[1302]	267	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
[109]	268	endif
	269
[97]	270	c------------------
	271	c TEST PK MONOTONE
	272	c------------------
	273	write(,) "Test PK"
	274	do ij=1,ip1jmp1
	275	do l=2,llm
	276	if(pk(ij,l).gt.pk(ij,l-1)) then
	277	c write(,) ij,l,pk(ij,l)
	278	abort_message = 'PK non strictement decroissante'
	279	call abort_gcm(modname,abort_message,1)
	280	c write(,) "ATTENTION, Test PK deconnecté..."
	281	endif
	282	enddo
	283	enddo
	284	write(,) "Fin Test PK"
	285	c stop
	286	c------------------
[1]	287
	288	c-----------------------------------------------------------------------
	289	c Debut de l'integration temporelle:
	290	c ----------------------------------
	291
[1824]	292	c RMBY: check that hour_ini and start_time are not both non-zero
	293	if ((hour_ini.ne.0.0).and.(start_time.ne.0.0)) then
	294	write(,) "ERROR: hour_ini = ", hour_ini,
	295	& "start_time = ", start_time
	296	abort_message = 'hour_ini and start_time both nonzero'
	297	call abort_gcm(modname,abort_message,1)
	298	endif
	299
[776]	300	1 CONTINUE ! Matsuno Forward step begins here
[1]	301
[1549]	302	c date: (NB: date remains unchanged for Backward step)
	303	c -----
	304
[492]	305	jD_cur = jD_ref + day_ini - day_ref + &
[1549]	306	& (itau+1)/day_step
[1824]	307	IF (planet_type .eq. "mars") THEN
	308	jH_cur = jH_ref + hour_ini + &
	309	& mod(itau+1,day_step)/float(day_step)
	310	ELSE
	311	jH_cur = jH_ref + start_time + &
	312	& mod(itau+1,day_step)/float(day_step)
	313	ENDIF
[492]	314	jD_cur = jD_cur + int(jH_cur)
	315	jH_cur = jH_cur - int(jH_cur)
[1]	316
[1508]	317	if (ok_iso_verif) then
	318	call check_isotopes_seq(q,ip1jmp1,'leapfrog 321')
	319	endif !if (ok_iso_verif) then
[1]	320
	321	#ifdef CPP_IOIPSL
[1024]	322	IF (planet_type.eq."earth") THEN
[1]	323	if (ok_guide) then
	324	call guide_main(itau,ucov,vcov,teta,q,masse,ps)
	325	endif
[1024]	326	ENDIF
[1]	327	#endif
	328
	329
	330	c
	331	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
	332	c CALL test_period ( ucov,vcov,teta,q,p,phis )
	333	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
	334	c ENDIF
	335	c
	336
	337	! Save fields obtained at previous time step as '...m1'
	338	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
	339	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
	340	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
	341	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
	342	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
	343
	344	forward = .TRUE.
	345	leapf = .FALSE.
	346	dt = dtvr
	347
	348	c ... P.Le Van .26/04/94 ....
[776]	349	! Ehouarn: finvmaold is actually not used
	350	! CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
	351	! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
[1]	352
[1508]	353	if (ok_iso_verif) then
	354	call check_isotopes_seq(q,ip1jmp1,'leapfrog 400')
	355	endif !if (ok_iso_verif) then
	356
[776]	357	2 CONTINUE ! Matsuno backward or leapfrog step begins here
[1]	358
	359	c-----------------------------------------------------------------------
	360
[1549]	361	c date: (NB: only leapfrog step requires recomputing date)
[1]	362	c -----
	363
[1549]	364	IF (leapf) THEN
	365	jD_cur = jD_ref + day_ini - day_ref +
	366	& (itau+1)/day_step
[1824]	367	IF (planet_type .eq. "mars") THEN
	368	jH_cur = jH_ref + hour_ini +
	369	& mod(itau+1,day_step)/float(day_step)
	370	ELSE
	371	jH_cur = jH_ref + start_time +
	372	& mod(itau+1,day_step)/float(day_step)
	373	ENDIF
[1549]	374	jD_cur = jD_cur + int(jH_cur)
	375	jH_cur = jH_cur - int(jH_cur)
	376	ENDIF
[1]	377
[1549]	378
[1]	379	c gestion des appels de la physique et des dissipations:
	380	c ------------------------------------------------------
	381	c
	382	c ... P.Le Van ( 6/02/95 ) ....
	383
	384	apphys = .FALSE.
	385	statcl = .FALSE.
	386	conser = .FALSE.
	387	apdiss = .FALSE.
	388
	389	IF( purmats ) THEN
	390	! Purely Matsuno time stepping
	391	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
[270]	392	IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward )
	393	s apdiss = .TRUE.
[1]	394	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
[1403]	395	s .and. physics ) apphys = .TRUE.
[1]	396	ELSE
	397	! Leapfrog/Matsuno time stepping
	398	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
[270]	399	IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
	400	s apdiss = .TRUE.
[1403]	401	IF( MOD(itau+1,iphysiq).EQ.0.AND.physics ) apphys=.TRUE.
[1]	402	END IF
	403
	404	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
	405	! supress dissipation step
	406	if (llm.eq.1) then
	407	apdiss=.false.
	408	endif
	409
[1024]	410	#ifdef NODYN
	411	apdiss=.false.
	412	#endif
	413
[1508]	414	if (ok_iso_verif) then
	415	call check_isotopes_seq(q,ip1jmp1,'leapfrog 589')
	416	endif !if (ok_iso_verif) then
	417
[1]	418	c-----------------------------------------------------------------------
	419	c calcul des tendances dynamiques:
	420	c --------------------------------
	421
[5]	422	! ADAPTATION GCM POUR CP(T)
	423	call tpot2t(ijp1llm,teta,temp,pk)
	424	tsurpk = cpp*temp/pk
[776]	425	! compute geopotential phi()
[5]	426	CALL geopot ( ip1jmp1, tsurpk , pk , pks, phis , phi )
[1]	427
[497]	428	rdaym_ini = itau * dtvr / daysec
	429
[1]	430	time = jD_cur + jH_cur
[1024]	431
	432	#ifdef NODYN
	433	dv(:,:) = 0.D+0
	434	du(:,:) = 0.D+0
	435	dteta(:,:) = 0.D+0
	436	dq(:,:,:) = 0.D+0
	437	dp(:) = 0.D+0
	438	#else
[1]	439	CALL caldyn
[5]	440	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,tsurpk,phis ,
[1]	441	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
	442
[1024]	443	! Simple zonal wind nudging for generic planetary model
	444	! AS 09/2013
	445	! ---------------------------------------------------
	446	if (planet_type.eq."generic") then
	447	if (ok_guide) then
	448	du(:,:) = du(:,:) + ((uforc(:,:)-ucov(:,:)) / facwind)
	449	endif
	450	endif
[1]	451
	452	c-----------------------------------------------------------------------
	453	c calcul des tendances advection des traceurs (dont l'humidite)
	454	c -------------------------------------------------------------
	455
[1508]	456	if (ok_iso_verif) then
	457	call check_isotopes_seq(q,ip1jmp1,
	458	& 'leapfrog 686: avant caladvtrac')
	459	endif !if (ok_iso_verif) then
	460
[1189]	461	IF( forward. OR . leapf ) THEN
	462	! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
[1]	463	CALL caladvtrac(q,pbaru,pbarv,
	464	* p, masse, dq, teta,
	465	. flxw, pk)
	466
	467	IF (offline) THEN
	468	Cmaf stokage du flux de masse pour traceurs OFF-LINE
	469
	470	#ifdef CPP_IOIPSL
	471	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
	472	. dtvr, itau)
	473	#endif
	474
	475
	476	ENDIF ! of IF (offline)
	477	c
	478	ENDIF ! of IF( forward. OR . leapf )
	479
	480
	481	c-----------------------------------------------------------------------
	482	c integrations dynamique et traceurs:
	483	c ----------------------------------
	484
[1508]	485	if (ok_iso_verif) then
	486	write(,) 'leapfrog 720'
	487	call check_isotopes_seq(q,ip1jmp1,'leapfrog 756')
	488	endif !if (ok_iso_verif) then
	489
	490	CALL integrd ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 ,
[776]	491	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
	492	! $ finvmaold )
[1]	493
[1508]	494	if (ok_iso_verif) then
	495	write(,) 'leapfrog 724'
	496	call check_isotopes_seq(q,ip1jmp1,'leapfrog 762')
	497	endif !if (ok_iso_verif) then
	498
[500]	499	IF ((planet_type.eq."titan").and.(tidal)) then
[495]	500	c-----------------------------------------------------------------------
	501	c Marées gravitationnelles causées par Saturne
	502	c B. Charnay (28/10/2010)
	503	c ----------------------------------------------------------
	504	CALL tidal_forces(rdaym_ini, dutidal, dvtidal)
	505	ucov=ucov+dutidal*dt
	506	vcov=vcov+dvtidal*dt
	507	ENDIF
[1]	508
[1024]	509	! NODYN precompiling flag
	510	#endif
	511
[1]	512	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
	513	c
	514	c-----------------------------------------------------------------------
	515	c calcul des tendances physiques:
	516	c -------------------------------
	517	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
	518	c
	519	IF( purmats ) THEN
	520	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
	521	ELSE
	522	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
	523	ENDIF
	524	c
	525	c
	526	IF( apphys ) THEN
	527	c
	528	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
	529	c
	530
	531	CALL pression ( ip1jmp1, ap, bp, ps, p )
[776]	532	if (pressure_exner) then
[1302]	533	CALL exner_hyb( ip1jmp1, ps, p,pks, pk, pkf )
[776]	534	else
[1302]	535	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
[109]	536	endif
[1]	537
[1302]	538	! Compute geopotential (physics might need it)
[1703]	539	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
[1302]	540
[492]	541	jD_cur = jD_ref + day_ini - day_ref + &
[1549]	542	& (itau+1)/day_step
[841]	543
[1107]	544	IF ((planet_type .eq."generic").or.
	545	& (planet_type .eq."mars")) THEN
[841]	546	! AS: we make jD_cur to be pday
	547	jD_cur = int(day_ini + itau/day_step)
	548	ENDIF
	549
[1824]	550	IF (planet_type .eq. "mars") THEN
	551	jH_cur = jH_ref + hour_ini + &
	552	& mod(itau,day_step)/float(day_step)
	553	ELSE IF (planet_type .eq. "generic") THEN
[1549]	554	jH_cur = jH_ref + start_time + &
[1824]	555	& mod(itau,day_step)/float(day_step)
	556	ELSE
	557	jH_cur = jH_ref + start_time + &
	558	& mod(itau+1,day_step)/float(day_step)
[1549]	559	ENDIF
[492]	560	jD_cur = jD_cur + int(jH_cur)
	561	jH_cur = jH_cur - int(jH_cur)
[1]	562	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
	563	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
	564	! write(lunout,*)'current date = ',an, mois, jour, secondes
	565
	566	c rajout debug
	567	c lafin = .true.
	568
	569
[6]	570	c Interface avec les routines de phylmd (phymars ... )
[1]	571	c -----------------------------------------------------
	572
	573	c+jld
	574
	575	c Diagnostique de conservation de l'énergie : initialisation
	576	IF (ip_ebil_dyn.ge.1 ) THEN
	577	ztit='bil dyn'
[1508]	578	! Ehouarn: be careful, diagedyn is Earth-specific!
[1]	579	IF (planet_type.eq."earth") THEN
	580	CALL diagedyn(ztit,2,1,1,dtphys
	581	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	582	ENDIF
	583	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
	584	c-jld
	585	#ifdef CPP_IOIPSL
[841]	586	cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
	587	cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
	588	c IF (first) THEN
	589	c first=.false.
	590	c#include "ini_paramLMDZ_dyn.h"
	591	c ENDIF
[1]	592	c
[841]	593	c#include "write_paramLMDZ_dyn.h"
[1]	594	c
	595	#endif
	596	! #endif of #ifdef CPP_IOIPSL
[6]	597
[1056]	598	c call WriteField('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
	599
[1]	600	CALL calfis( lafin , jD_cur, jH_cur,
	601	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
	602	$ du,dv,dteta,dq,
	603	$ flxw,
[97]	604	$ dufi,dvfi,dtetafi,dqfi,dpfi )
[1]	605
[1056]	606	c call WriteField('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
	607	c call WriteField('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
	608	c call WriteField('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
	609
[6]	610	c ajout des tendances physiques:
	611	c ------------------------------
	612	CALL addfi( dtphys, leapf, forward ,
	613	$ ucov, vcov, teta , q ,ps ,
	614	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
[1189]	615	! since addfi updates ps(), also update p(), masse() and pk()
	616	CALL pression (ip1jmp1,ap,bp,ps,p)
	617	CALL massdair(p,masse)
	618	if (pressure_exner) then
[1302]	619	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
[1189]	620	else
[1302]	621	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
[1189]	622	endif
	623
[6]	624	c Couche superieure :
	625	c -------------------
[1012]	626	IF (iflag_top_bound > 0) THEN
[1010]	627	CALL top_bound(vcov,ucov,teta,masse,dtphys,dutop)
	628	dutop(:,:)=dutop(:,:)/dtphys ! convert to a tendency in (m/s)/s
[108]	629	ENDIF
	630
[1]	631	c Diagnostique de conservation de l'énergie : difference
	632	IF (ip_ebil_dyn.ge.1 ) THEN
	633	ztit='bil phys'
	634	IF (planet_type.eq."earth") THEN
	635	CALL diagedyn(ztit,2,1,1,dtphys
	636	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	637	ENDIF
	638	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
	639
	640	ENDIF ! of IF( apphys )
	641
	642	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
	643	! Academic case : Simple friction and Newtonan relaxation
	644	! -------------------------------------------------------
	645	DO l=1,llm
	646	DO ij=1,ip1jmp1
	647	teta(ij,l)=teta(ij,l)-dtvr*
	648	& (teta(ij,l)-tetarappel(ij,l))(knewt_g+knewt_t(l)clat4(ij))
	649	ENDDO
	650	ENDDO ! of DO l=1,llm
[53]	651
[270]	652	if (planet_type.eq."giant") then
	653	! add an intrinsic heat flux at the base of the atmosphere
	654	teta(:,1)=teta(:,1)+dtvraire(:)ihf/cpp/masse(:,1)
	655	endif
[53]	656
[270]	657	call friction(ucov,vcov,dtvr)
	658
[53]	659
[1]	660	! Sponge layer (if any)
	661	IF (ok_strato) THEN
[1010]	662	CALL top_bound(vcov,ucov,teta,masse,dtvr,dutop)
	663	dutop(:,:)=dutop(:,:)/dtvr ! convert to a tendency in (m/s)/s
[1]	664	ENDIF ! of IF (ok_strato)
	665	ENDIF ! of IF (iflag_phys.EQ.2)
	666
	667
	668	c-jld
	669
	670	CALL pression ( ip1jmp1, ap, bp, ps, p )
[776]	671	if (pressure_exner) then
[1302]	672	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
[776]	673	else
[1302]	674	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
[109]	675	endif
[1189]	676	CALL massdair(p,masse)
[1]	677
[1508]	678	if (ok_iso_verif) then
	679	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196')
	680	endif !if (ok_iso_verif) then
	681
[1]	682	c-----------------------------------------------------------------------
	683	c dissipation horizontale et verticale des petites echelles:
	684	c ----------------------------------------------------------
	685
	686	IF(apdiss) THEN
	687
[1017]	688	! sponge layer
	689	if (callsponge) then
	690	CALL sponge(ucov,vcov,teta,ps,dtdiss,mode_sponge)
	691	endif
[1]	692
	693	c calcul de l'energie cinetique avant dissipation
	694	call covcont(llm,ucov,vcov,ucont,vcont)
	695	call enercin(vcov,ucov,vcont,ucont,ecin0)
	696
	697	c dissipation
[5]	698	! ADAPTATION GCM POUR CP(T)
	699	call tpot2t(ijp1llm,teta,temp,pk)
	700
[1]	701	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
	702	ucov=ucov+dudis
	703	vcov=vcov+dvdis
[6]	704	dudis=dudis/dtdiss ! passage en (m/s)/s
	705	dvdis=dvdis/dtdiss ! passage en (m/s)/s
[1]	706
	707	c------------------------------------------------------------------------
	708	if (dissip_conservative) then
	709	C On rajoute la tendance due a la transform. Ec -> E therm. cree
	710	C lors de la dissipation
	711	call covcont(llm,ucov,vcov,ucont,vcont)
	712	call enercin(vcov,ucov,vcont,ucont,ecin)
[5]	713	! ADAPTATION GCM POUR CP(T)
	714	do ij=1,ip1jmp1
	715	do l=1,llm
	716	dtec = (ecin0(ij,l)-ecin(ij,l))/cpdet(temp(ij,l))
	717	temp(ij,l) = temp(ij,l) + dtec
	718	enddo
	719	enddo
	720	call t2tpot(ijp1llm,temp,ztetaec,pk)
	721	dtetaecdt=ztetaec-teta
[1]	722	dtetadis=dtetadis+dtetaecdt
	723	endif
	724	teta=teta+dtetadis
[6]	725	dtetadis=dtetadis/dtdiss ! passage en K/s
[1]	726	c------------------------------------------------------------------------
	727
	728
	729	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
	730	c ... Calcul de la valeur moyenne, unique de h aux poles .....
	731	c
	732
	733	DO l = 1, llm
	734	DO ij = 1,iim
	735	tppn(ij) = aire( ij ) * teta( ij ,l)
	736	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	737	ENDDO
	738	tpn = SSUM(iim,tppn,1)/apoln
	739	tps = SSUM(iim,tpps,1)/apols
	740
	741	DO ij = 1, iip1
	742	teta( ij ,l) = tpn
	743	teta(ij+ip1jm,l) = tps
	744	ENDDO
	745	ENDDO
	746
[776]	747	if (1 == 0) then
	748	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
	749	!!! 2) should probably not be here anyway
	750	!!! but are kept for those who would want to revert to previous behaviour
	751	DO ij = 1,iim
	752	tppn(ij) = aire( ij ) * ps ( ij )
	753	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
	754	ENDDO
	755	tpn = SSUM(iim,tppn,1)/apoln
	756	tps = SSUM(iim,tpps,1)/apols
[1]	757
[776]	758	DO ij = 1, iip1
	759	ps( ij ) = tpn
	760	ps(ij+ip1jm) = tps
	761	ENDDO
	762	endif ! of if (1 == 0)
[1]	763
	764	END IF ! of IF(apdiss)
	765
	766	c ajout debug
	767	c IF( lafin ) then
	768	c abort_message = 'Simulation finished'
	769	c call abort_gcm(modname,abort_message,0)
	770	c ENDIF
	771
	772	c ********************************************************************
	773	c ********************************************************************
	774	c .... fin de l'integration dynamique et physique pour le pas itau ..
	775	c ********************************************************************
	776	c ********************************************************************
	777
	778	c preparation du pas d'integration suivant ......
	779
[1508]	780	if (ok_iso_verif) then
	781	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509')
	782	endif !if (ok_iso_verif) then
	783
[1]	784	IF ( .NOT.purmats ) THEN
	785	c ........................................................
	786	c .............. schema matsuno + leapfrog ..............
	787	c ........................................................
	788
	789	IF(forward. OR. leapf) THEN
	790	itau= itau + 1
	791	c iday= day_ini+itau/day_step
	792	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	793	c IF(time.GT.1.) THEN
	794	c time = time-1.
	795	c iday = iday+1
	796	c ENDIF
	797	ENDIF
	798
	799
	800	IF( itau. EQ. itaufinp1 ) then
	801	if (flag_verif) then
	802	write(79,*) 'ucov',ucov
	803	write(80,*) 'vcov',vcov
	804	write(81,*) 'teta',teta
	805	write(82,*) 'ps',ps
	806	write(83,*) 'q',q
	807	WRITE(85,*) 'q1 = ',q(:,:,1)
	808	WRITE(86,*) 'q3 = ',q(:,:,3)
	809	endif
	810
	811	abort_message = 'Simulation finished'
	812
	813	call abort_gcm(modname,abort_message,0)
	814	ENDIF
	815	c-----------------------------------------------------------------------
	816	c ecriture du fichier histoire moyenne:
	817	c -------------------------------------
	818
	819	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	820	IF(itau.EQ.itaufin) THEN
	821	iav=1
	822	ELSE
	823	iav=0
	824	ENDIF
	825
[1572]	826	! ! Ehouarn: re-compute geopotential for outputs
[1703]	827	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
[1572]	828
[1]	829	IF (ok_dynzon) THEN
	830	#ifdef CPP_IOIPSL
[6]	831	c les traceurs ne sont pas sortis, trop lourd.
	832	c Peut changer eventuellement si besoin.
	833	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
	834	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
[108]	835	& du,dudis,dutop,dufi)
[1]	836	#endif
	837	END IF
	838	IF (ok_dyn_ave) THEN
	839	#ifdef CPP_IOIPSL
	840	CALL writedynav(itau,vcov,
	841	& ucov,teta,pk,phi,q,masse,ps,phis)
	842	#endif
	843	ENDIF
	844
	845	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
	846
[1508]	847	if (ok_iso_verif) then
	848	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
	849	endif !if (ok_iso_verif) then
	850
[1]	851	c-----------------------------------------------------------------------
	852	c ecriture de la bande histoire:
	853	c ------------------------------
	854
	855	IF( MOD(itau,iecri).EQ.0) THEN
	856	! Ehouarn: output only during LF or Backward Matsuno
	857	if (leapf.or.(.not.leapf.and.(.not.forward))) then
[5]	858	! ADAPTATION GCM POUR CP(T)
	859	call tpot2t(ijp1llm,teta,temp,pk)
	860	tsurpk = cpp*temp/pk
	861	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
[1]	862	unat=0.
	863	do l=1,llm
	864	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	865	vnat(:,l)=vcov(:,l)/cv(:)
	866	enddo
	867	#ifdef CPP_IOIPSL
	868	if (ok_dyn_ins) then
	869	! write(lunout,*) "leapfrog: call writehist, itau=",itau
	870	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	871	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	872	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
	873	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
	874	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
	875	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
	876	endif ! of if (ok_dyn_ins)
	877	#endif
	878	! For some Grads outputs of fields
	879	if (output_grads_dyn) then
	880	#include "write_grads_dyn.h"
	881	endif
	882	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
	883	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
	884
[1824]	885	c Determine whether to write to the restart.nc file
	886	c Decision can't be made in one IF statement as if
	887	c ecritstart==0 there will be a divide-by-zero error
	888	lrestart = .false.
	889	IF (itau.EQ.itaufin) THEN
	890	lrestart = .true.
	891	ELSE IF (ecritstart.GT.0) THEN
	892	IF (MOD(itau,ecritstart).EQ.0) lrestart = .true.
	893	ENDIF
[1]	894
[1824]	895	c Write to the restart.nc file
	896	IF (lrestart) THEN
[1107]	897	if (planet_type=="mars") then
	898	CALL dynredem1("restart.nc",REAL(itau)/REAL(day_step),
	899	& vcov,ucov,teta,q,masse,ps)
[6]	900	else
[492]	901	CALL dynredem1("restart.nc",start_time,
[1]	902	& vcov,ucov,teta,q,masse,ps)
[1107]	903	endif
[1]	904	CLOSE(99)
[776]	905	!!! Ehouarn: Why not stop here and now?
[1824]	906	ENDIF ! of IF (lrestart)
[1]	907
	908	c-----------------------------------------------------------------------
	909	c gestion de l'integration temporelle:
	910	c ------------------------------------
	911
	912	IF( MOD(itau,iperiod).EQ.0 ) THEN
	913	GO TO 1
	914	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	915
	916	IF( forward ) THEN
	917	c fin du pas forward et debut du pas backward
	918
	919	forward = .FALSE.
	920	leapf = .FALSE.
	921	GO TO 2
	922
	923	ELSE
	924	c fin du pas backward et debut du premier pas leapfrog
	925
	926	leapf = .TRUE.
	927	dt = 2.*dtvr
	928	GO TO 2
	929	END IF ! of IF (forward)
	930	ELSE
	931
	932	c ...... pas leapfrog .....
	933
	934	leapf = .TRUE.
	935	dt = 2.*dtvr
	936	GO TO 2
	937	END IF ! of IF (MOD(itau,iperiod).EQ.0)
	938	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
	939
	940	ELSE ! of IF (.not.purmats)
	941
[1508]	942	if (ok_iso_verif) then
	943	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
	944	endif !if (ok_iso_verif) then
	945
[1]	946	c ........................................................
	947	c .............. schema matsuno ...............
	948	c ........................................................
	949	IF( forward ) THEN
	950
	951	itau = itau + 1
	952	c iday = day_ini+itau/day_step
	953	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	954	c
	955	c IF(time.GT.1.) THEN
	956	c time = time-1.
	957	c iday = iday+1
	958	c ENDIF
	959
	960	forward = .FALSE.
	961	IF( itau. EQ. itaufinp1 ) then
	962	abort_message = 'Simulation finished'
	963	call abort_gcm(modname,abort_message,0)
	964	ENDIF
	965	GO TO 2
	966
	967	ELSE ! of IF(forward) i.e. backward step
	968
[1508]	969	if (ok_iso_verif) then
	970	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
	971	endif !if (ok_iso_verif) then
	972
[1]	973	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	974	IF(itau.EQ.itaufin) THEN
	975	iav=1
	976	ELSE
	977	iav=0
	978	ENDIF
	979
[1572]	980	! ! Ehouarn: re-compute geopotential for outputs
[1703]	981	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
[1572]	982
[1]	983	IF (ok_dynzon) THEN
	984	#ifdef CPP_IOIPSL
[6]	985	c les traceurs ne sont pas sortis, trop lourd.
	986	c Peut changer eventuellement si besoin.
	987	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
	988	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
[108]	989	& du,dudis,dutop,dufi)
[1]	990	#endif
	991	ENDIF
	992	IF (ok_dyn_ave) THEN
	993	#ifdef CPP_IOIPSL
	994	CALL writedynav(itau,vcov,
	995	& ucov,teta,pk,phi,q,masse,ps,phis)
	996	#endif
	997	ENDIF
	998
	999	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
	1000
	1001	IF(MOD(itau,iecri ).EQ.0) THEN
	1002	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
[5]	1003	! ADAPTATION GCM POUR CP(T)
	1004	call tpot2t(ijp1llm,teta,temp,pk)
	1005	tsurpk = cpp*temp/pk
	1006	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
[1]	1007	unat=0.
	1008	do l=1,llm
	1009	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	1010	vnat(:,l)=vcov(:,l)/cv(:)
	1011	enddo
	1012	#ifdef CPP_IOIPSL
	1013	if (ok_dyn_ins) then
	1014	! write(lunout,*) "leapfrog: call writehist (b)",
	1015	! & itau,iecri
	1016	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	1017	endif ! of if (ok_dyn_ins)
	1018	#endif
	1019	! For some Grads outputs
	1020	if (output_grads_dyn) then
	1021	#include "write_grads_dyn.h"
	1022	endif
	1023
	1024	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
	1025
[1824]	1026	c Determine whether to write to the restart.nc file
	1027	c Decision can't be made in one IF statement as if
	1028	c ecritstart==0 there will be a divide-by-zero error
	1029	lrestart = .false.
	1030	IF (itau.EQ.itaufin) THEN
	1031	lrestart = .true.
	1032	ELSE IF (ecritstart.GT.0) THEN
	1033	IF (MOD(itau,ecritstart).EQ.0) lrestart = .true.
	1034	ENDIF
	1035
	1036	c Write to the restart.nc file
	1037	IF (lrestart) THEN
[1107]	1038	if (planet_type=="mars") then
	1039	CALL dynredem1("restart.nc",REAL(itau)/REAL(day_step),
	1040	& vcov,ucov,teta,q,masse,ps)
[6]	1041	else
[492]	1042	CALL dynredem1("restart.nc",start_time,
[6]	1043	& vcov,ucov,teta,q,masse,ps)
[1107]	1044	endif
[1824]	1045	ENDIF ! of IF (lrestart)
[1]	1046
	1047	forward = .TRUE.
	1048	GO TO 1
	1049
	1050	ENDIF ! of IF (forward)
	1051
	1052	END IF ! of IF(.not.purmats)
	1053
	1054	STOP
	1055	END

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

Download in other formats:

Original Format