Context Navigation

leapfrog.F @ 1685

Last change on this file since 1685 was 1572, checked in by emillour, 9 years ago

All GCMs:
Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation (up to rev r2500 of LMDZ5)

arch:
remove ifort debug option '-check all', replace it with '-check bounds,format,output_conversion,pointers,uninit' (i.e. get it to stop complaining about copying into temporary arrays)

dyn3d_common:
comconst_mod.F90 : add ngroup

dyn3d:
gcm.F90 : minor bug fix (arguments to a call_abort())
leapfrog.F90 : recompute geopotential for bilan_dyn outputs
conf_gcm.F90 : read "ngroup" from run.def
groupe.F , groupeun.F : ngroup no longer a local parameter

dyn3d_par:
conf_gcm.F90 : read "ngroup" from run.def
groupe_p.F , groupeun_p.F : ngroup no longer a local parameter

misc:
regr1_step_av_m.F90 : removed (not used)

phy_common:
mod_phys_lmdz_mpi_transfert.F90 , mod_phys_lmdz_mpi_data.F90 : change is_north_pole and is_south_pole to is_north_pole_dyn and is_south_pole_dyn
mod_phys_lmdz_omp_data.F90 : introduce is_nort_pole_phy and is_south_pole_phy

dynphy_lonlat:
mod_interface_dyn_phys.F90 : use is_north_pole_dyn and is_south_pole_dyn
calfis_p.F : use is_north_pole_dyn and is_south_pole_dyn

phyvenus:
physiq_mod , write_hist*.h : use is_north_pole_phy and is_south_pole_phy to correctly compute mesh area at poles to send to hist*nc files.

phytitan:
physiq_mod , write_hist*.h : use is_north_pole_phy and is_south_pole_phy to correctly compute mesh area at poles to send to hist*nc files.

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

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

Download in other formats:

Original Format