Context Navigation

leapfrog.F @ 1540

Last change on this file since 1540 was 1508, checked in by emillour, 10 years ago

Common dynamics:
Updates in the dynamics (seq and ) to keep up with updates
in LMDZ5 (up to LMDZ5 trunk, rev 2325):
IMPORTANT: Modifications for isotopes are only done in dyn3d, not in dyn3dpar

as in LMDZ5 these modifications were done in dyn3dmem.
Related LMDZ5 revisions are r2270 and r2281

in dynlonlat_phylonlat:
add module "grid_atob_m.F90" (a regridding utility so far only used by phylmd/ce0l.F90, used to be dyn3d_common/grid_atob.F)

in misc:
follow up updates on wxios.F (add missing_val module variable)

in dyn3d_common:
pression.F => pression.F90
misc_mod.F90: moved from misc to dyn3d_common
added new iso_verif_dyn.F
covcont.F => covcont.F90
infotrac.F90 : add handling of isotopes (reading of corresponding traceur.def for planets not implemented)
dynetat0.F => dynetat0.F90 with some code factorization
dynredem.F => dynredem.F90 with some code factorization
added dynredem_mod.F90: routines used by dynredem
iniacademic.F90 : added isotopes-related initialization for Earth case

in dyn3d:
added check_isotopes.F
modified (isotopes) advtrac.F90, caladvtrac.F
guide_mod.F90: ported updates
leapfrog.F : (isotopes) updates (NB: call integrd with nqtot tracers)
qminimium.F : adaptations for isotopes (copied over, except that #include comvert.h is not needed).
vlsplt.F: adaptations for isotopes (copied over, except than #include logic.h, comvert.h not needed, and replace "include comconst.h" with use comconst_mod, ONLY: pi)
vlspltqs.F : same as vlsplt.F, but also keeping added modification for CP(T)

in dyn3dpar:
leapfrog_p.F: remove unecessary #ifdef CPP_EARTH cpp flag. and call integrd_p with nqtot tracers (only important for Earth)
dynredem_p.F => dynredem_p.F90 and some code factorization
and no isotopes-relates changes in dyn3dpar (since these changes have been made in LMDZ5 dyn3dmem).

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

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

Download in other formats:

Original Format