Context Navigation

leapfrog.F @ 1443

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

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

Download in other formats:

Original Format