Context Navigation

leapfrog.F @ 3020

Last change on this file since 3020 was 2544, checked in by romain.vande, 3 years ago

MARS GCM:

update of r2507: correction for time output in the case of multiple restart files
RV

File size: 35.7 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,
[2507]	33	. start_time,dt,hour_ini,day_end
[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)
[2507]	280	c write(,) "ATTENTION, Test PK deconnecte..."
[97]	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
[2380]	438	flxw(:,:) = 0.D+0
[1024]	439	#else
[1]	440	CALL caldyn
[5]	441	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,tsurpk,phis ,
[1]	442	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
	443
[1024]	444	! Simple zonal wind nudging for generic planetary model
	445	! AS 09/2013
	446	! ---------------------------------------------------
	447	if (planet_type.eq."generic") then
	448	if (ok_guide) then
	449	du(:,:) = du(:,:) + ((uforc(:,:)-ucov(:,:)) / facwind)
	450	endif
	451	endif
[1]	452
	453	c-----------------------------------------------------------------------
	454	c calcul des tendances advection des traceurs (dont l'humidite)
	455	c -------------------------------------------------------------
	456
[1508]	457	if (ok_iso_verif) then
	458	call check_isotopes_seq(q,ip1jmp1,
	459	& 'leapfrog 686: avant caladvtrac')
	460	endif !if (ok_iso_verif) then
	461
[1189]	462	IF( forward. OR . leapf ) THEN
	463	! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
[1]	464	CALL caladvtrac(q,pbaru,pbarv,
	465	* p, masse, dq, teta,
	466	. flxw, pk)
	467
	468	IF (offline) THEN
	469	Cmaf stokage du flux de masse pour traceurs OFF-LINE
	470
	471	#ifdef CPP_IOIPSL
	472	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
	473	. dtvr, itau)
	474	#endif
	475
	476
	477	ENDIF ! of IF (offline)
	478	c
	479	ENDIF ! of IF( forward. OR . leapf )
	480
	481
	482	c-----------------------------------------------------------------------
	483	c integrations dynamique et traceurs:
	484	c ----------------------------------
	485
[1508]	486	if (ok_iso_verif) then
	487	write(,) 'leapfrog 720'
	488	call check_isotopes_seq(q,ip1jmp1,'leapfrog 756')
	489	endif !if (ok_iso_verif) then
	490
	491	CALL integrd ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 ,
[776]	492	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
	493	! $ finvmaold )
[1]	494
[1508]	495	if (ok_iso_verif) then
	496	write(,) 'leapfrog 724'
	497	call check_isotopes_seq(q,ip1jmp1,'leapfrog 762')
	498	endif !if (ok_iso_verif) then
	499
[500]	500	IF ((planet_type.eq."titan").and.(tidal)) then
[495]	501	c-----------------------------------------------------------------------
[2507]	502	c Marees gravitationnelles causees par Saturne
[495]	503	c B. Charnay (28/10/2010)
	504	c ----------------------------------------------------------
	505	CALL tidal_forces(rdaym_ini, dutidal, dvtidal)
	506	ucov=ucov+dutidal*dt
	507	vcov=vcov+dvtidal*dt
	508	ENDIF
[1]	509
[1024]	510	! NODYN precompiling flag
	511	#endif
	512
[1]	513	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
	514	c
	515	c-----------------------------------------------------------------------
	516	c calcul des tendances physiques:
	517	c -------------------------------
	518	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
	519	c
	520	IF( purmats ) THEN
	521	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
	522	ELSE
	523	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
	524	ENDIF
	525	c
	526	c
	527	IF( apphys ) THEN
	528	c
	529	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
	530	c
	531
	532	CALL pression ( ip1jmp1, ap, bp, ps, p )
[776]	533	if (pressure_exner) then
[1302]	534	CALL exner_hyb( ip1jmp1, ps, p,pks, pk, pkf )
[776]	535	else
[1302]	536	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
[109]	537	endif
[1]	538
[1302]	539	! Compute geopotential (physics might need it)
[2195]	540	! GEOP CORRECTION
	541	! ADAPTATION GCM POUR CP(T)
	542	! CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	543	call tpot2t(ijp1llm,teta,temp,pk)
	544	tsurpk = cpp*temp/pk
	545	CALL geopot ( ip1jmp1, tsurpk, pk, pks, phis, phi )
[1302]	546
[492]	547	jD_cur = jD_ref + day_ini - day_ref + &
[1549]	548	& (itau+1)/day_step
[841]	549
[1107]	550	IF ((planet_type .eq."generic").or.
	551	& (planet_type .eq."mars")) THEN
[841]	552	! AS: we make jD_cur to be pday
	553	jD_cur = int(day_ini + itau/day_step)
	554	ENDIF
	555
[1824]	556	IF (planet_type .eq. "mars") THEN
	557	jH_cur = jH_ref + hour_ini + &
	558	& mod(itau,day_step)/float(day_step)
	559	ELSE IF (planet_type .eq. "generic") THEN
[1549]	560	jH_cur = jH_ref + start_time + &
[1824]	561	& mod(itau,day_step)/float(day_step)
	562	ELSE
	563	jH_cur = jH_ref + start_time + &
	564	& mod(itau+1,day_step)/float(day_step)
[1549]	565	ENDIF
[492]	566	jD_cur = jD_cur + int(jH_cur)
	567	jH_cur = jH_cur - int(jH_cur)
[1]	568	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
	569	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
	570	! write(lunout,*)'current date = ',an, mois, jour, secondes
	571
	572	c rajout debug
	573	c lafin = .true.
	574
	575
[6]	576	c Interface avec les routines de phylmd (phymars ... )
[1]	577	c -----------------------------------------------------
	578
	579	c+jld
	580
[2507]	581	c Diagnostique de conservation de l'energie : initialisation
[1]	582	IF (ip_ebil_dyn.ge.1 ) THEN
	583	ztit='bil dyn'
[1508]	584	! Ehouarn: be careful, diagedyn is Earth-specific!
[1]	585	IF (planet_type.eq."earth") THEN
	586	CALL diagedyn(ztit,2,1,1,dtphys
	587	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	588	ENDIF
	589	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
	590	c-jld
	591	#ifdef CPP_IOIPSL
[841]	592	cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
	593	cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
	594	c IF (first) THEN
	595	c first=.false.
	596	c#include "ini_paramLMDZ_dyn.h"
	597	c ENDIF
[1]	598	c
[841]	599	c#include "write_paramLMDZ_dyn.h"
[1]	600	c
	601	#endif
	602	! #endif of #ifdef CPP_IOIPSL
[6]	603
[1056]	604	c call WriteField('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
	605
[1]	606	CALL calfis( lafin , jD_cur, jH_cur,
	607	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
	608	$ du,dv,dteta,dq,
	609	$ flxw,
[97]	610	$ dufi,dvfi,dtetafi,dqfi,dpfi )
[1]	611
[1056]	612	c call WriteField('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
	613	c call WriteField('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
	614	c call WriteField('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
	615
[6]	616	c ajout des tendances physiques:
	617	c ------------------------------
	618	CALL addfi( dtphys, leapf, forward ,
	619	$ ucov, vcov, teta , q ,ps ,
	620	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
[1189]	621	! since addfi updates ps(), also update p(), masse() and pk()
	622	CALL pression (ip1jmp1,ap,bp,ps,p)
	623	CALL massdair(p,masse)
	624	if (pressure_exner) then
[1302]	625	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
[1189]	626	else
[1302]	627	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
[1189]	628	endif
	629
[6]	630	c Couche superieure :
	631	c -------------------
[1012]	632	IF (iflag_top_bound > 0) THEN
[1010]	633	CALL top_bound(vcov,ucov,teta,masse,dtphys,dutop)
	634	dutop(:,:)=dutop(:,:)/dtphys ! convert to a tendency in (m/s)/s
[108]	635	ENDIF
	636
[2507]	637	c Diagnostique de conservation de l'energie : difference
[1]	638	IF (ip_ebil_dyn.ge.1 ) THEN
	639	ztit='bil phys'
	640	IF (planet_type.eq."earth") THEN
	641	CALL diagedyn(ztit,2,1,1,dtphys
	642	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	643	ENDIF
	644	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
	645
	646	ENDIF ! of IF( apphys )
	647
	648	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
	649	! Academic case : Simple friction and Newtonan relaxation
	650	! -------------------------------------------------------
	651	DO l=1,llm
	652	DO ij=1,ip1jmp1
	653	teta(ij,l)=teta(ij,l)-dtvr*
	654	& (teta(ij,l)-tetarappel(ij,l))(knewt_g+knewt_t(l)clat4(ij))
	655	ENDDO
	656	ENDDO ! of DO l=1,llm
[53]	657
[270]	658	if (planet_type.eq."giant") then
	659	! add an intrinsic heat flux at the base of the atmosphere
	660	teta(:,1)=teta(:,1)+dtvraire(:)ihf/cpp/masse(:,1)
	661	endif
[53]	662
[270]	663	call friction(ucov,vcov,dtvr)
	664
[53]	665
[1]	666	! Sponge layer (if any)
	667	IF (ok_strato) THEN
[1010]	668	CALL top_bound(vcov,ucov,teta,masse,dtvr,dutop)
	669	dutop(:,:)=dutop(:,:)/dtvr ! convert to a tendency in (m/s)/s
[1]	670	ENDIF ! of IF (ok_strato)
	671	ENDIF ! of IF (iflag_phys.EQ.2)
	672
	673
	674	c-jld
	675
	676	CALL pression ( ip1jmp1, ap, bp, ps, p )
[776]	677	if (pressure_exner) then
[1302]	678	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
[776]	679	else
[1302]	680	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
[109]	681	endif
[1189]	682	CALL massdair(p,masse)
[1]	683
[1508]	684	if (ok_iso_verif) then
	685	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196')
	686	endif !if (ok_iso_verif) then
	687
[1]	688	c-----------------------------------------------------------------------
	689	c dissipation horizontale et verticale des petites echelles:
	690	c ----------------------------------------------------------
	691
	692	IF(apdiss) THEN
	693
[1017]	694	! sponge layer
	695	if (callsponge) then
	696	CALL sponge(ucov,vcov,teta,ps,dtdiss,mode_sponge)
	697	endif
[1]	698
	699	c calcul de l'energie cinetique avant dissipation
	700	call covcont(llm,ucov,vcov,ucont,vcont)
	701	call enercin(vcov,ucov,vcont,ucont,ecin0)
	702
	703	c dissipation
[5]	704	! ADAPTATION GCM POUR CP(T)
	705	call tpot2t(ijp1llm,teta,temp,pk)
	706
[1]	707	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
	708	ucov=ucov+dudis
	709	vcov=vcov+dvdis
[6]	710	dudis=dudis/dtdiss ! passage en (m/s)/s
	711	dvdis=dvdis/dtdiss ! passage en (m/s)/s
[1]	712
	713	c------------------------------------------------------------------------
	714	if (dissip_conservative) then
	715	C On rajoute la tendance due a la transform. Ec -> E therm. cree
	716	C lors de la dissipation
	717	call covcont(llm,ucov,vcov,ucont,vcont)
	718	call enercin(vcov,ucov,vcont,ucont,ecin)
[5]	719	! ADAPTATION GCM POUR CP(T)
	720	do ij=1,ip1jmp1
	721	do l=1,llm
	722	dtec = (ecin0(ij,l)-ecin(ij,l))/cpdet(temp(ij,l))
	723	temp(ij,l) = temp(ij,l) + dtec
	724	enddo
	725	enddo
	726	call t2tpot(ijp1llm,temp,ztetaec,pk)
	727	dtetaecdt=ztetaec-teta
[1]	728	dtetadis=dtetadis+dtetaecdt
	729	endif
	730	teta=teta+dtetadis
[6]	731	dtetadis=dtetadis/dtdiss ! passage en K/s
[1]	732	c------------------------------------------------------------------------
	733
	734
	735	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
	736	c ... Calcul de la valeur moyenne, unique de h aux poles .....
	737	c
	738
	739	DO l = 1, llm
	740	DO ij = 1,iim
	741	tppn(ij) = aire( ij ) * teta( ij ,l)
	742	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	743	ENDDO
	744	tpn = SSUM(iim,tppn,1)/apoln
	745	tps = SSUM(iim,tpps,1)/apols
	746
	747	DO ij = 1, iip1
	748	teta( ij ,l) = tpn
	749	teta(ij+ip1jm,l) = tps
	750	ENDDO
	751	ENDDO
	752
[776]	753	if (1 == 0) then
	754	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
	755	!!! 2) should probably not be here anyway
	756	!!! but are kept for those who would want to revert to previous behaviour
	757	DO ij = 1,iim
	758	tppn(ij) = aire( ij ) * ps ( ij )
	759	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
	760	ENDDO
	761	tpn = SSUM(iim,tppn,1)/apoln
	762	tps = SSUM(iim,tpps,1)/apols
[1]	763
[776]	764	DO ij = 1, iip1
	765	ps( ij ) = tpn
	766	ps(ij+ip1jm) = tps
	767	ENDDO
	768	endif ! of if (1 == 0)
[1]	769
	770	END IF ! of IF(apdiss)
	771
	772	c ajout debug
	773	c IF( lafin ) then
	774	c abort_message = 'Simulation finished'
	775	c call abort_gcm(modname,abort_message,0)
	776	c ENDIF
	777
	778	c ********************************************************************
	779	c ********************************************************************
	780	c .... fin de l'integration dynamique et physique pour le pas itau ..
	781	c ********************************************************************
	782	c ********************************************************************
	783
	784	c preparation du pas d'integration suivant ......
	785
[1508]	786	if (ok_iso_verif) then
	787	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509')
	788	endif !if (ok_iso_verif) then
	789
[1]	790	IF ( .NOT.purmats ) THEN
	791	c ........................................................
	792	c .............. schema matsuno + leapfrog ..............
	793	c ........................................................
	794
	795	IF(forward. OR. leapf) THEN
	796	itau= itau + 1
	797	c iday= day_ini+itau/day_step
	798	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	799	c IF(time.GT.1.) THEN
	800	c time = time-1.
	801	c iday = iday+1
	802	c ENDIF
	803	ENDIF
	804
	805
	806	IF( itau. EQ. itaufinp1 ) then
	807	if (flag_verif) then
	808	write(79,*) 'ucov',ucov
	809	write(80,*) 'vcov',vcov
	810	write(81,*) 'teta',teta
	811	write(82,*) 'ps',ps
	812	write(83,*) 'q',q
	813	WRITE(85,*) 'q1 = ',q(:,:,1)
	814	WRITE(86,*) 'q3 = ',q(:,:,3)
	815	endif
	816
	817	abort_message = 'Simulation finished'
	818
	819	call abort_gcm(modname,abort_message,0)
	820	ENDIF
	821	c-----------------------------------------------------------------------
	822	c ecriture du fichier histoire moyenne:
	823	c -------------------------------------
	824
	825	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	826	IF(itau.EQ.itaufin) THEN
	827	iav=1
	828	ELSE
	829	iav=0
	830	ENDIF
	831
[1572]	832	! ! Ehouarn: re-compute geopotential for outputs
[2195]	833	! GEOP CORRECTION
	834	! ADAPTATION GCM POUR CP(T)
	835	! CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	836	call tpot2t(ijp1llm,teta,temp,pk)
	837	tsurpk = cpp*temp/pk
	838	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
[1572]	839
[1]	840	IF (ok_dynzon) THEN
	841	#ifdef CPP_IOIPSL
[6]	842	c les traceurs ne sont pas sortis, trop lourd.
	843	c Peut changer eventuellement si besoin.
	844	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
	845	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
[108]	846	& du,dudis,dutop,dufi)
[1]	847	#endif
	848	END IF
	849	IF (ok_dyn_ave) THEN
	850	#ifdef CPP_IOIPSL
	851	CALL writedynav(itau,vcov,
	852	& ucov,teta,pk,phi,q,masse,ps,phis)
	853	#endif
	854	ENDIF
	855
	856	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
	857
[1508]	858	if (ok_iso_verif) then
	859	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
	860	endif !if (ok_iso_verif) then
	861
[1]	862	c-----------------------------------------------------------------------
	863	c ecriture de la bande histoire:
	864	c ------------------------------
	865
	866	IF( MOD(itau,iecri).EQ.0) THEN
	867	! Ehouarn: output only during LF or Backward Matsuno
	868	if (leapf.or.(.not.leapf.and.(.not.forward))) then
[5]	869	! ADAPTATION GCM POUR CP(T)
	870	call tpot2t(ijp1llm,teta,temp,pk)
	871	tsurpk = cpp*temp/pk
	872	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
[1]	873	unat=0.
	874	do l=1,llm
	875	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	876	vnat(:,l)=vcov(:,l)/cv(:)
	877	enddo
	878	#ifdef CPP_IOIPSL
	879	if (ok_dyn_ins) then
	880	! write(lunout,*) "leapfrog: call writehist, itau=",itau
	881	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	882	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	883	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
	884	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
	885	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
	886	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
	887	endif ! of if (ok_dyn_ins)
	888	#endif
	889	! For some Grads outputs of fields
	890	if (output_grads_dyn) then
	891	#include "write_grads_dyn.h"
	892	endif
	893	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
	894	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
	895
[1824]	896	c Determine whether to write to the restart.nc file
	897	c Decision can't be made in one IF statement as if
	898	c ecritstart==0 there will be a divide-by-zero error
	899	lrestart = .false.
	900	IF (itau.EQ.itaufin) THEN
	901	lrestart = .true.
	902	ELSE IF (ecritstart.GT.0) THEN
	903	IF (MOD(itau,ecritstart).EQ.0) lrestart = .true.
	904	ENDIF
[1]	905
[1824]	906	c Write to the restart.nc file
	907	IF (lrestart) THEN
[1107]	908	if (planet_type=="mars") then
[2507]	909	if(ecritstart.GT.0) then
	910	CALL dynredem1("restart.nc",
[2544]	911	& REAL(itau)/REAL(day_step)+
	912	& time_0-floor(time_0),
[2507]	913	& vcov,ucov,teta,q,masse,ps)
	914	else
	915	CALL dynredem1("restart.nc",
[2544]	916	& REAL(itau)/REAL(day_step)-(day_end-day_ini)
	917	& +time_0-floor(time_0),
[2507]	918	& vcov,ucov,teta,q,masse,ps)
	919	endif
[6]	920	else
[492]	921	CALL dynredem1("restart.nc",start_time,
[1]	922	& vcov,ucov,teta,q,masse,ps)
[1107]	923	endif
[1]	924	CLOSE(99)
[776]	925	!!! Ehouarn: Why not stop here and now?
[1824]	926	ENDIF ! of IF (lrestart)
[1]	927
	928	c-----------------------------------------------------------------------
	929	c gestion de l'integration temporelle:
	930	c ------------------------------------
	931
	932	IF( MOD(itau,iperiod).EQ.0 ) THEN
	933	GO TO 1
	934	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	935
	936	IF( forward ) THEN
	937	c fin du pas forward et debut du pas backward
	938
	939	forward = .FALSE.
	940	leapf = .FALSE.
	941	GO TO 2
	942
	943	ELSE
	944	c fin du pas backward et debut du premier pas leapfrog
	945
	946	leapf = .TRUE.
	947	dt = 2.*dtvr
	948	GO TO 2
	949	END IF ! of IF (forward)
	950	ELSE
	951
	952	c ...... pas leapfrog .....
	953
	954	leapf = .TRUE.
	955	dt = 2.*dtvr
	956	GO TO 2
	957	END IF ! of IF (MOD(itau,iperiod).EQ.0)
	958	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
	959
	960	ELSE ! of IF (.not.purmats)
	961
[1508]	962	if (ok_iso_verif) then
	963	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
	964	endif !if (ok_iso_verif) then
	965
[1]	966	c ........................................................
	967	c .............. schema matsuno ...............
	968	c ........................................................
	969	IF( forward ) THEN
	970
	971	itau = itau + 1
	972	c iday = day_ini+itau/day_step
	973	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	974	c
	975	c IF(time.GT.1.) THEN
	976	c time = time-1.
	977	c iday = iday+1
	978	c ENDIF
	979
	980	forward = .FALSE.
	981	IF( itau. EQ. itaufinp1 ) then
	982	abort_message = 'Simulation finished'
	983	call abort_gcm(modname,abort_message,0)
	984	ENDIF
	985	GO TO 2
	986
	987	ELSE ! of IF(forward) i.e. backward step
	988
[1508]	989	if (ok_iso_verif) then
	990	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
	991	endif !if (ok_iso_verif) then
	992
[1]	993	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	994	IF(itau.EQ.itaufin) THEN
	995	iav=1
	996	ELSE
	997	iav=0
	998	ENDIF
	999
[1572]	1000	! ! Ehouarn: re-compute geopotential for outputs
[2195]	1001	! GEOP CORRECTION
	1002	! ADAPTATION GCM POUR CP(T)
	1003	! CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	1004	call tpot2t(ijp1llm,teta,temp,pk)
	1005	tsurpk = cpp*temp/pk
	1006	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
[1572]	1007
[1]	1008	IF (ok_dynzon) THEN
	1009	#ifdef CPP_IOIPSL
[6]	1010	c les traceurs ne sont pas sortis, trop lourd.
	1011	c Peut changer eventuellement si besoin.
	1012	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
	1013	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
[108]	1014	& du,dudis,dutop,dufi)
[1]	1015	#endif
	1016	ENDIF
	1017	IF (ok_dyn_ave) THEN
	1018	#ifdef CPP_IOIPSL
	1019	CALL writedynav(itau,vcov,
	1020	& ucov,teta,pk,phi,q,masse,ps,phis)
	1021	#endif
	1022	ENDIF
	1023
	1024	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
	1025
	1026	IF(MOD(itau,iecri ).EQ.0) THEN
	1027	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
[5]	1028	! ADAPTATION GCM POUR CP(T)
	1029	call tpot2t(ijp1llm,teta,temp,pk)
	1030	tsurpk = cpp*temp/pk
	1031	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
[1]	1032	unat=0.
	1033	do l=1,llm
	1034	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	1035	vnat(:,l)=vcov(:,l)/cv(:)
	1036	enddo
	1037	#ifdef CPP_IOIPSL
	1038	if (ok_dyn_ins) then
	1039	! write(lunout,*) "leapfrog: call writehist (b)",
	1040	! & itau,iecri
	1041	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	1042	endif ! of if (ok_dyn_ins)
	1043	#endif
	1044	! For some Grads outputs
	1045	if (output_grads_dyn) then
	1046	#include "write_grads_dyn.h"
	1047	endif
	1048
	1049	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
	1050
[1824]	1051	c Determine whether to write to the restart.nc file
	1052	c Decision can't be made in one IF statement as if
	1053	c ecritstart==0 there will be a divide-by-zero error
	1054	lrestart = .false.
	1055	IF (itau.EQ.itaufin) THEN
	1056	lrestart = .true.
	1057	ELSE IF (ecritstart.GT.0) THEN
	1058	IF (MOD(itau,ecritstart).EQ.0) lrestart = .true.
	1059	ENDIF
	1060
	1061	c Write to the restart.nc file
	1062	IF (lrestart) THEN
[1107]	1063	if (planet_type=="mars") then
[2507]	1064	if(ecritstart.GT.0) then
	1065	CALL dynredem1("restart.nc",
	1066	& REAL(itau)/REAL(day_step),
[1107]	1067	& vcov,ucov,teta,q,masse,ps)
[2507]	1068	else
	1069	CALL dynredem1("restart.nc",
	1070	& REAL(itau)/REAL(day_step)-(day_end-day_ini),
	1071	& vcov,ucov,teta,q,masse,ps)
	1072	endif
[6]	1073	else
[492]	1074	CALL dynredem1("restart.nc",start_time,
[6]	1075	& vcov,ucov,teta,q,masse,ps)
[1107]	1076	endif
[1824]	1077	ENDIF ! of IF (lrestart)
[1]	1078
	1079	forward = .TRUE.
	1080	GO TO 1
	1081
	1082	ENDIF ! of IF (forward)
	1083
	1084	END IF ! of IF(.not.purmats)
	1085
	1086	STOP
	1087	END

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

Download in other formats:

Original Format