Context Navigation

leapfrog_nogcm.F @ 3607

Last change on this file since 3607 was 3539, checked in by tbertrand, 2 months ago
LMDZ.PLUTO: Update and Validation of the Volatile Transport Model (VTM, or "nogcm") Merging missing elements with Pluto.old TB
File size: 26.3 KB

Rev	Line
[3410]	1	!
	2	! $Id: leapfrog.F 1446 2010-10-22 09:27:25Z emillour $
	3	!
	4	c
	5	c
	6	SUBROUTINE leapfrog_nogcm(ucov,vcov,teta,ps,masse,phis,q,
	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
	14	USE infotrac, ONLY: nqtot,ok_iso_verif,tname
	15	USE write_field, ONLY: writefield
[3421]	16	USE callkeys_mod, only: tau_n2, tau_ch4, tau_co
[3410]	17	USE control_mod, ONLY: planet_type,nday,day_step,iperiod,iphysiq,
	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
	22	use exner_hyb_m, only: exner_hyb
	23	use exner_milieu_m, only: exner_milieu
	24	use cpdet_mod, only: cpdet,tpot2t,t2tpot
	25	use comuforc_h
	26	USE comvert_mod, ONLY: ap,bp,pressure_exner,presnivs,
	27	& aps,bps,presnivs,pseudoalt,preff,scaleheight
	28	USE comconst_mod, ONLY: daysec,dtvr,dtphys,dtdiss,
	29	. cpp,ihf,iflag_top_bound,pi,kappa,r
	30	USE logic_mod, ONLY: iflag_phys,ok_guide,forward,leapf,apphys,
	31	. statcl,conser,purmats,tidal,ok_strato
	32	USE temps_mod, ONLY: jD_ref,jH_ref,itaufin,day_ini,day_ref,
	33	. start_time,dt
	34
	35
	36
	37
	38	IMPLICIT NONE
	39
[3539]	40	c ...... Version du 30/11/24 ..........
[3410]	41
	42	c=======================================================================
	43	c
[3539]	44	c Auteur: T. Bertrand, F. Forget, A. Falco
[3410]	45	c -------
	46	c
	47	c Objet:
	48	c ------
	49	c
[3539]	50	c GCM LMD sans dynamique, pour modele saisonnier de surface
[3410]	51	c
	52	c=======================================================================
	53	c
	54	c-----------------------------------------------------------------------
	55	c Declarations:
	56	c -------------
	57
	58	#include "dimensions.h"
	59	#include "paramet.h"
	60	#include "comdissnew.h"
	61	#include "comgeom.h"
	62	!#include "com_io_dyn.h"
	63	#include "iniprint.h"
	64	#include "academic.h"
	65
	66	REAL,INTENT(IN) :: time_0 ! not used
	67
	68	c dynamical variables:
	69	REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
	70	REAL,INTENT(INOUT) :: vcov(ip1jm,llm) ! meridional covariant wind
	71	REAL,INTENT(INOUT) :: teta(ip1jmp1,llm) ! potential temperature
	72	REAL,INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
	73	REAL,INTENT(INOUT) :: masse(ip1jmp1,llm) ! air mass
	74	REAL,INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
	75	REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
	76
	77	REAL p (ip1jmp1,llmp1 ) ! interlayer pressure
	78	REAL pks(ip1jmp1) ! exner at the surface
	79	REAL pk(ip1jmp1,llm) ! exner at mid-layer
	80	REAL pkf(ip1jmp1,llm) ! filtered exner at mid-layer
	81	REAL phi(ip1jmp1,llm) ! geopotential
	82	REAL w(ip1jmp1,llm) ! vertical velocity
	83	REAL kpd(ip1jmp1) ! TB15 = exp (-z/H)
	84
	85	! ADAPTATION GCM POUR CP(T)
	86	REAL temp(ip1jmp1,llm) ! temperature
	87	REAL tsurpk(ip1jmp1,llm) ! cpp*T/pk
	88
[3539]	89	! nogcm
[3410]	90	REAL tau_ps
[3421]	91	REAL tau_x
	92	REAL tau_def
[3410]	93	REAL tau_teta
	94	REAL tetadpmean
	95	REAL tetadp(ip1jmp1,llm)
	96	REAL dpmean
	97	REAL dp(ip1jmp1,llm)
	98	REAL tetamean
	99	real globaverage2d
	100
	101	c variables dynamiques intermediaire pour le transport
	102	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
	103
	104	c variables dynamiques au pas -1
	105	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
	106	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
	107	REAL massem1(ip1jmp1,llm)
	108
	109	c tendances dynamiques en */s
	110	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
	111	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot)
	112
	113	c tendances physiques */s
	114	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
	115	REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
	116
	117	c variables pour le fichier histoire
	118	! REAL dtav ! intervalle de temps elementaire
	119
	120	REAL tppn(iim),tpps(iim),tpn,tps
	121	c
	122	INTEGER itau,itaufinp1,iav
	123	! INTEGER iday ! jour julien
	124	REAL time
	125
	126	REAL SSUM
	127
	128	cym LOGICAL lafin
	129	LOGICAL :: lafin=.false.
	130	INTEGER ij,iq,l
	131
	132	REAL rdaym_ini
	133	! jD_cur: jour julien courant
	134	! jH_cur: heure julienne courante
	135	REAL :: jD_cur, jH_cur
	136	! INTEGER :: an, mois, jour
	137	! REAL :: secondes
	138
	139	LOGICAL first,callinigrads
	140	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
	141	save first
	142	data first/.true./
	143	logical ok_sync
	144	parameter (ok_sync = .true.)
	145	logical physics
	146
	147	data callinigrads/.true./
	148	character*10 string10
	149
	150	! REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
	151	REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
	152
	153	REAL psmean ! pression moyenne
[3421]	154	REAL pqxmean ! moyenne globale ps*qco
[3410]	155	REAL p0 ! pression de reference
	156	REAL p00d ! globalaverage(kpd)
[3421]	157	REAL qmean_x,qmean_x_vert ! mass mean mixing ratio vap n2
	158	REAL pqx(ip1jmp1) ! average n2 mass index : ps*q_n2
[3410]	159	REAL oldps(ip1jmp1) ! saving old pressure ps to calculate qch4
	160
	161	c+jld variables test conservation energie
	162	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
	163	C Tendance de la temp. potentiel d (theta)/ d t due a la
	164	C tansformation d'energie cinetique en energie thermique
	165	C cree par la dissipation
	166	REAL dtetaecdt(ip1jmp1,llm)
	167	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
	168	REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
	169	CHARACTER*15 ztit
	170
	171	character(len=*),parameter :: modname="leapfrog"
	172	character*80 abort_message
	173
	174	logical dissip_conservative
	175	save dissip_conservative
	176	data dissip_conservative/.true./
	177
	178	INTEGER testita
	179	PARAMETER (testita = 9)
	180
	181	logical , parameter :: flag_verif = .false.
	182
	183	! for CP(T)
	184	real :: dtec
	185	real :: ztetaec(ip1jmp1,llm)
	186
	187	! TEMP : diagnostic mass
[3411]	188	real :: n2mass(iip1,jjp1)
	189	real :: n2ice_ij(iip1,jjp1)
[3421]	190	integer,save :: igcm_n2=0 ! index of N2 tracer (if any)
	191	integer,save :: igcm_ch4=0 ! index of CH4 tracer (if any)
	192	integer,save :: igcm_co=0 ! index of CO tracer (if any)
[3410]	193	integer :: i,j,ig
	194	integer, parameter :: ngrid = 2+(jjm-1)*iim
	195	! local mass
	196	real :: mq(ip1jmp1,llm)
	197
	198	if (nday>=0) then
	199	itaufin = nday*day_step
	200	else
	201	! to run a given (-nday) number of dynamical steps
	202	itaufin = -nday
	203	endif
	204	if (less1day) then
	205	c MODIF VENUS: to run less than one day:
	206	itaufin = int(fractday*day_step)
	207	endif
	208	if (ndynstep.gt.0) then
	209	! running a given number of dynamical steps
	210	itaufin=ndynstep
	211	endif
	212	itaufinp1 = itaufin +1
	213
	214
	215
	216
	217	c INITIALISATIONS
	218	dufi(:,:) =0.
	219	dvfi(:,:) =0.
	220	dtetafi(:,:)=0.
	221	dqfi(:,:,:) =0.
	222	dpfi(:) =0.
	223
	224	itau = 0
	225	physics=.true.
	226	if (iflag_phys==0.or.iflag_phys==2) physics=.false.
	227
	228	! ED18 nogcm
	229	! firstcall_globaverage2d = .true.
	230
	231	c iday = day_ini+itau/day_step
	232	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	233	c IF(time.GT.1.) THEN
	234	c time = time-1.
	235	c iday = iday+1
	236	c ENDIF
	237
	238
	239	c-----------------------------------------------------------------------
	240	c On initialise la pression et la fonction d'Exner :
	241	c --------------------------------------------------
	242
	243	dq(:,:,:)=0.
	244	CALL pression ( ip1jmp1, ap, bp, ps, p )
	245	if (pressure_exner) then
	246	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
	247	else
	248	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
	249	endif
	250
	251	c------------------
	252	c TEST PK MONOTONE
	253	c------------------
	254	write(,) "Test PK"
	255	do ij=1,ip1jmp1
	256	do l=2,llm
	257	! PRINT*,'pk(ij,l) = ',pk(ij,l)
	258	! PRINT*,'pk(ij,l-1) = ',pk(ij,l-1)
	259	if(pk(ij,l).gt.pk(ij,l-1)) then
	260	c write(,) ij,l,pk(ij,l)
	261	abort_message = 'PK non strictement decroissante'
	262	call abort_gcm(modname,abort_message,1)
	263	c write(,) "ATTENTION, Test PK deconnecte..."
	264	endif
	265	enddo
	266	enddo
	267	write(,) "Fin Test PK"
	268	c stop
[3539]	269
[3410]	270	c------------------
	271	c Preparing mixing of pressure and tracers in nogcm
	272	c kpd=exp(-z/H) Needed to define a reference pressure :
	273	c P0=pmean/globalaverage(kpd)
	274	c thus P(i) = P0exp(-z(i)/H) = P0kpd(i)
	275	c It is checked that globalaverage2d(Pi)=pmean
	276	DO ij=1,ip1jmp1
[3539]	277	kpd(ij) = exp(-phis(ij)/(r*38.))
[3410]	278	ENDDO
	279	p00d=globaverage2d(kpd) ! mean pres at ref level
[3539]	280	tau_ps = tau_n2 ! constante de rappel for pressure (s)
	281	tau_def = 1. ! default constante de rappel for mix ratio qX (s)
[3410]	282	tau_teta = 1.E7 !constante de rappel for potentiel temperature
	283
	284	do iq=1,nqtot
[3411]	285	if (tname(iq)=="n2") then
	286	igcm_n2=iq
[3421]	287	! exit
	288	else if (tname(iq)=="ch4_gas") then
	289	igcm_ch4=iq
	290	else if (tname(iq)=="co_gas") then
	291	igcm_co=iq
[3410]	292	endif
	293	enddo
[3421]	294	PRINT*,'igcm_n2 = ',igcm_n2
	295	PRINT*,'igcm_ch4 = ',igcm_ch4
	296	PRINT*,'igcm_co = ',igcm_co
[3410]	297
	298	c-----------------------------------------------------------------------
	299	c Debut de l'integration temporelle:
	300	c ----------------------------------
	301
	302	1 CONTINUE ! Matsuno Forward step begins here
	303
	304	c date: (NB: date remains unchanged for Backward step)
	305	c -----
	306
[3539]	307	jD_cur = jD_ref + day_ini - day_ref +
[3410]	308	& (itau+1)/day_step
[3539]	309	jH_cur = jH_ref + start_time +
[3410]	310	& mod(itau+1,day_step)/float(day_step)
	311	jD_cur = jD_cur + int(jH_cur)
	312	jH_cur = jH_cur - int(jH_cur)
	313
	314	! Save fields obtained at previous time step as '...m1'
	315	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
	316	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
	317	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
	318	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
	319	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
	320
	321	forward = .TRUE.
	322	leapf = .FALSE.
	323	dt = dtvr
	324
	325	2 CONTINUE ! Matsuno backward or leapfrog step begins here
	326
	327	c-----------------------------------------------------------------------
	328
	329	c date: (NB: only leapfrog step requires recomputing date)
	330	c -----
	331
	332	IF (leapf) THEN
	333	jD_cur = jD_ref + day_ini - day_ref +
	334	& (itau+1)/day_step
	335	jH_cur = jH_ref + start_time +
	336	& mod(itau+1,day_step)/float(day_step)
	337	jD_cur = jD_cur + int(jH_cur)
	338	jH_cur = jH_cur - int(jH_cur)
	339	ENDIF
	340
	341
	342	c gestion des appels de la physique et des dissipations:
	343	c ------------------------------------------------------
	344
	345	apphys = .FALSE.
	346	statcl = .FALSE.
	347	conser = .FALSE.
	348
	349	IF( purmats ) THEN
	350	! Purely Matsuno time stepping
	351	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
	352
	353	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
	354	s .and. physics ) apphys = .TRUE.
	355	ELSE
[3539]	356	! Leapfrog/Matsuno time stepping
	357	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
	358	IF( MOD(itau+1,iphysiq).EQ.0.AND.physics) apphys=.TRUE.
[3410]	359	END IF
	360
	361	c-----------------------------------------------------------------------
	362	c calcul des tendances dynamiques:
	363	c --------------------------------
	364
[3539]	365	dv(:,:) = 0.
	366	du(:,:) = 0.
	367	dteta(:,:) = 0.
	368	dq(:,:,:) = 0.
[3410]	369
	370	c-----------------------------------------------------------------------
	371	c calcul des tendances physiques:
	372	c -------------------------------
	373	c
	374	IF( purmats ) THEN
	375	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
	376	ELSE
	377	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
	378	ENDIF
[3539]	379
[3410]	380	IF( apphys ) THEN
	381
	382	CALL pression ( ip1jmp1, ap, bp, ps, p )
	383	if (pressure_exner) then
	384	CALL exner_hyb( ip1jmp1, ps, p,pks, pk, pkf )
	385	else
	386	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
	387	endif
	388
	389	! Compute geopotential (physics might need it)
	390	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	391
[3539]	392	jD_cur = jD_ref + day_ini - day_ref +
[3410]	393	& (itau+1)/day_step
	394
[3539]	395	! AS: we make jD_cur to be pday
	396	jD_cur = int(day_ini + itau/day_step)
[3410]	397
[3539]	398	jH_cur = jH_ref + start_time +
[3410]	399	& mod(itau+1,day_step)/float(day_step)
[3539]	400	jH_cur = jH_ref + start_time +
[3410]	401	& mod(itau,day_step)/float(day_step)
[3539]	402	jD_cur = jD_cur + int(jH_cur)
	403	jH_cur = jH_cur - int(jH_cur)
[3410]	404
	405	c Interface avec les routines de phylmd (phymars ... )
	406	c -----------------------------------------------------
	407
	408	CALL calfis( lafin , jD_cur, jH_cur,
	409	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
	410	$ du,dv,dteta,dq,
	411	$ flxw,
	412	$ dufi,dvfi,dtetafi,dqfi,dpfi )
	413
	414	c ajout des tendances physiques
	415	c ------------------------------
	416
[3539]	417	CALL addfi( dtphys, leapf, forward ,
[3410]	418	$ ucov, vcov, teta , q ,ps ,
	419	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	420
[3539]	421	CALL pression (ip1jmp1,ap,bp,ps,p)
	422	CALL massdair(p,masse)
[3410]	423
[3539]	424	! 1) Mass of tracers in each mess (kg) before Horizontal mixing of pressure
[3410]	425
	426	DO l=1, llm
	427	DO ij=1,ip1jmp1
[3411]	428	mq(ij,l) = masse(ij,l)*q(ij,l,igcm_n2)
[3410]	429	ENDDO
	430	ENDDO
	431
[3539]	432	! 2) Horizontal mixing of pressure
[3410]	433	! Rappel newtonien vers psmean
[3539]	434	psmean= globaverage2d(ps) ! mean pressure
	435	p0=psmean/p00d
	436	DO ij=1,ip1jmp1
	437	oldps(ij)=ps(ij)
	438	ps(ij)=ps(ij) +(p0*kpd(ij)
[3410]	439	& -ps(ij))*(1-exp(-dtphys/tau_ps))
[3539]	440	ENDDO
[3410]	441
[3539]	442	write(72,*) 'psmean = ',psmean ! mean pressure redistributed
[3410]	443
[3539]	444	! 3) Security check: test pressure negative
	445	DO ij=1,ip1jmp1
[3410]	446	IF (ps(ij).le.0.) then
	447	!write(,) 'Negative pressure :'
	448	!write(,) 'ps = ',ps(ij),' ig = ',ij
	449	!write(,) 'pmean = ',p0*kpd(ij)
	450	!write(,) 'tau_ps = ',tau_ps
	451	!STOP
[3539]	452	ps(ij)=1.e-6*kpd(ij)/p00d
[3410]	453	ENDIF
[3539]	454	ENDDO
[3410]	455
[3539]	456	! ! 4) Correction on teta due to surface pressure changes
	457	! DO l = 1,llm
	458	! DO ij = 1,ip1jmp1
	459	! teta(ij,l)= teta(ij,l)(ps(ij)/oldps(ij))*kappa
	460	! ENDDO
	461	! ENDDO
[3410]	462
[3539]	463	! 5) Update pressure p and pk
	464	CALL pression (ip1jmp1,ap,bp,ps,p)
	465	CALL massdair(p,masse)
	466	if (pressure_exner) then
[3410]	467	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
[3539]	468	else
[3410]	469	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
[3539]	470	endif
[3410]	471
[3539]	472	! 6) Update tracers after Horizontal mixing of pressure to ! conserve tracer mass
[3410]	473	DO l=1, llm
	474	DO ij=1,ip1jmp1
[3411]	475	q(ij,l,igcm_n2) = mq(ij,l)/ masse(ij,l)
[3410]	476	ENDDO
	477	ENDDO
	478
[3539]	479	! 7) Horizontal mixing tracers, and temperature (replace dynamics in nogcm)
	480	psmean= globaverage2d(ps) ! mean pressure
[3410]	481
[3421]	482	! Simulate redistribution by dynamics for qX
	483	DO iq=1,nqtot
	484	if ((iq.eq.igcm_n2).or.(iq.eq.igcm_ch4).or.
	485	& (iq.eq.igcm_co)) then
[3410]	486
	487	DO ij=1,ip1jmp1
[3421]	488	pqx(ij)= ps(ij) * q(ij,1,iq)
[3410]	489	ENDDO
[3421]	490	pqxmean=globaverage2d(pqx)
[3410]	491
[3539]	492	! Rappel newtonien vers qx_mean
[3421]	493	qmean_x= pqxmean / psmean
	494
	495	tau_x = tau_def
	496	if (iq.eq.igcm_n2) then
	497	tau_x = tau_n2
	498	else if (iq.eq.igcm_ch4) then
	499	tau_x = tau_ch4
	500	else if (iq.eq.igcm_co) then
	501	tau_x = tau_co
	502	end if
	503	PRINT*,' tau_x ',iq,tau_x
[3410]	504
	505	DO ij=1,ip1jmp1
[3421]	506	q(ij,1,iq)=q(ij,1,iq)+
	507	& (qmean_x-q(ij,1,iq))*
	508	& (1.-exp(-dtphys/tau_x))
[3410]	509	ENDDO
	510
	511	DO l=2, llm
	512	DO ij=1,ip1jmp1
[3421]	513	q(ij,l,iq)=q(ij,1,iq)
[3410]	514	END DO
	515	END DO
	516
[3539]	517	! Diagnostics
	518	! PRINT*,'psmean = ',psmean
	519	! PRINT*,'qmean_x = ',qmean_x
	520	! PRINT*,'pqxmean = ',pqxmean
[3421]	521	! PRINT*,' q(50,1) = ',iq,q(50,1,iq)
	522	! PRINT*,' q(50,2) = ',iq,q(50,2,iq)
	523	! PRINT*,' q(50,3) = ',iq,q(50,3,iq)
[3410]	524
[3411]	525	endif ! igcm_n2.ne.0
[3539]	526	ENDDO
[3410]	527
[3539]	528	! 8) Mixing Temperature horizontally
	529	! initialize variables that will be averaged
	530	! DO l=1,llm
	531	! DO ij=1,ip1jmp1
	532	! dp(ij,l) = p(ij,l) - p(ij,l+1)
	533	! tetadp(ij,l) = teta(ij,l)*dp(ij,l)
	534	! ENDDO
	535	! ENDDO
[3410]	536
[3539]	537	! DO l=1,llm
	538	! tetadpmean = globaverage2d(tetadp(:,l))
	539	! dpmean = globaverage2d(dp(:,l))
	540	! tetamean = tetadpmean / dpmean
	541	! DO ij=1,ip1jmp1
	542	! teta(ij,l) = teta(ij,l) + (tetamean - teta(ij,l)) *
	543	! & (1 - exp(-dtphys/tau_teta))
	544	! ENDDO
	545	! ENDDO
[3410]	546
	547	ENDIF ! of IF( apphys )
	548
	549
	550	c ********************************************************************
	551	c ********************************************************************
	552	c .... fin de l'integration dynamique et physique pour le pas itau ..
	553	c ********************************************************************
	554	c ********************************************************************
	555
	556	IF ( .NOT.purmats ) THEN
	557	c ........................................................
	558	c .............. schema matsuno + leapfrog ..............
	559	c ........................................................
	560
	561	IF(forward. OR. leapf) THEN
	562	itau= itau + 1
	563	c iday= day_ini+itau/day_step
	564	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	565	c IF(time.GT.1.) THEN
	566	c time = time-1.
	567	c iday = iday+1
	568	c ENDIF
	569	ENDIF
	570
	571	IF( itau. EQ. itaufinp1 ) then
	572	if (flag_verif) then
	573	write(79,*) 'ucov',ucov
	574	write(80,*) 'vcov',vcov
	575	write(81,*) 'teta',teta
	576	write(82,*) 'ps',ps
	577	write(83,*) 'q',q
	578	WRITE(85,*) 'q1 = ',q(:,:,1)
	579	WRITE(86,*) 'q3 = ',q(:,:,3)
	580	endif
	581
	582	abort_message = 'Simulation finished'
	583
	584	call abort_gcm(modname,abort_message,0)
	585	ENDIF
[3539]	586
[3410]	587	c-----------------------------------------------------------------------
	588	c ecriture du fichier histoire moyenne:
	589	c -------------------------------------
	590
	591	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	592	IF(itau.EQ.itaufin) THEN
	593	iav=1
	594	ELSE
	595	iav=0
	596	ENDIF
	597
	598	! ! Ehouarn: re-compute geopotential for outputs
	599	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	600
	601	IF (ok_dyn_ave) THEN
	602	#ifdef CPP_IOIPSL
	603	CALL writedynav(itau,vcov,
	604	& ucov,teta,pk,phi,q,masse,ps,phis)
	605	#endif
	606	ENDIF
	607
	608	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
	609
	610	if (ok_iso_verif) then
	611	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
	612	endif !if (ok_iso_verif) then
	613
	614	c-----------------------------------------------------------------------
	615	c ecriture de la bande histoire:
	616	c ------------------------------
	617
	618	IF( MOD(itau,iecri).EQ.0) THEN
	619	! Ehouarn: output only during LF or Backward Matsuno
	620	if (leapf.or.(.not.leapf.and.(.not.forward))) then
	621	! ADAPTATION GCM POUR CP(T)
	622	call tpot2t(ijp1llm,teta,temp,pk)
	623	tsurpk = cpp*temp/pk
	624	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
	625	unat=0.
	626	do l=1,llm
	627	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	628	vnat(:,l)=vcov(:,l)/cv(:)
	629	enddo
	630	#ifdef CPP_IOIPSL
	631	if (ok_dyn_ins) then
	632	! write(lunout,*) "leapfrog: call writehist, itau=",itau
	633	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	634	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	635	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
	636	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
	637	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
	638	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
	639	endif ! of if (ok_dyn_ins)
	640	#endif
	641	! For some Grads outputs of fields
	642	if (output_grads_dyn) then
	643	#include "write_grads_dyn.h"
	644	endif
	645	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
	646	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
	647
	648	IF(itau.EQ.itaufin) THEN
	649
[3411]	650	CALL dynredem1("restart.nc",start_time,
[3410]	651	& vcov,ucov,teta,q,masse,ps)
	652	CLOSE(99)
	653	!!! Ehouarn: Why not stop here and now?
	654	ENDIF ! of IF (itau.EQ.itaufin)
	655
	656	c-----------------------------------------------------------------------
	657	c gestion de l'integration temporelle:
	658	c ------------------------------------
	659
	660	IF( MOD(itau,iperiod).EQ.0 ) THEN
	661	GO TO 1
	662	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	663
	664	IF( forward ) THEN
	665	c fin du pas forward et debut du pas backward
	666
	667	forward = .FALSE.
	668	leapf = .FALSE.
	669	GO TO 2
	670
	671	ELSE
	672	c fin du pas backward et debut du premier pas leapfrog
	673
	674	leapf = .TRUE.
	675	dt = 2.*dtvr
	676	GO TO 2
	677	END IF ! of IF (forward)
	678	ELSE
	679
	680	c ...... pas leapfrog .....
	681
	682	leapf = .TRUE.
	683	dt = 2.*dtvr
	684	GO TO 2
	685	END IF ! of IF (MOD(itau,iperiod).EQ.0)
	686	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
	687
	688	ELSE ! of IF (.not.purmats)
	689
	690	if (ok_iso_verif) then
	691	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
	692	endif !if (ok_iso_verif) then
	693
	694	c ........................................................
	695	c .............. schema matsuno ...............
	696	c ........................................................
	697	IF( forward ) THEN
	698
	699	itau = itau + 1
	700	c iday = day_ini+itau/day_step
	701	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	702	c
	703	c IF(time.GT.1.) THEN
	704	c time = time-1.
	705	c iday = iday+1
	706	c ENDIF
	707
	708	forward = .FALSE.
	709	IF( itau. EQ. itaufinp1 ) then
	710	abort_message = 'Simulation finished'
	711	call abort_gcm(modname,abort_message,0)
	712	ENDIF
	713	GO TO 2
	714
	715	ELSE ! of IF(forward) i.e. backward step
	716
	717	if (ok_iso_verif) then
	718	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
	719	endif !if (ok_iso_verif) then
	720
	721	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	722	IF(itau.EQ.itaufin) THEN
	723	iav=1
	724	ELSE
	725	iav=0
	726	ENDIF
	727
	728	! ! Ehouarn: re-compute geopotential for outputs
	729	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	730
	731	IF (ok_dyn_ave) THEN
	732	#ifdef CPP_IOIPSL
	733	CALL writedynav(itau,vcov,
	734	& ucov,teta,pk,phi,q,masse,ps,phis)
	735	#endif
	736	ENDIF
	737
	738	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
	739
	740	IF(MOD(itau,iecri ).EQ.0) THEN
	741	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
	742	! ADAPTATION GCM POUR CP(T)
	743	call tpot2t(ijp1llm,teta,temp,pk)
	744	tsurpk = cpp*temp/pk
	745	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
	746	unat=0.
	747	do l=1,llm
	748	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	749	vnat(:,l)=vcov(:,l)/cv(:)
	750	enddo
	751	#ifdef CPP_IOIPSL
	752	if (ok_dyn_ins) then
	753	! write(lunout,*) "leapfrog: call writehist (b)",
	754	! & itau,iecri
	755	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	756	endif ! of if (ok_dyn_ins)
	757	#endif
	758	! For some Grads outputs
	759	if (output_grads_dyn) then
	760	#include "write_grads_dyn.h"
	761	endif
	762
	763	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
	764
	765	IF(itau.EQ.itaufin) THEN
	766	CALL dynredem1("restart.nc",start_time,
	767	& vcov,ucov,teta,q,masse,ps)
	768	ENDIF ! of IF(itau.EQ.itaufin)
	769
	770	forward = .TRUE.
	771	GO TO 1
	772
	773	ENDIF ! of IF (forward)
	774
	775	END IF ! of IF(.not.purmats)
	776
	777	STOP
	778	END
	779
	780	! ************************************************************
	781	! globalaverage VERSION PLuto
	782	! ************************************************************
	783	FUNCTION globaverage2d(var)
	784	! ************************************************************
	785	IMPLICIT NONE
	786	#include "dimensions.h"
	787	#include "paramet.h"
	788	#include "comgeom2.h"
	789	!ip1jmp1 called in paramet.h = iip1 x jjp1
	790	REAL var(iip1,jjp1), globaverage2d
	791	integer i,j
	792	REAL airetot
	793	save airetot
	794	logical firstcall
	795	data firstcall/.true./
	796	save firstcall
	797
	798	if (firstcall) then
	799	firstcall=.false.
	800	airetot =0.
	801	DO j=2,jjp1-1 ! lat
	802	DO i = 1, iim ! lon
	803	airetot = airetot + aire(i,j)
	804	END DO
	805	END DO
	806	DO i=1,iim
	807	airetot=airetot + aire(i,1)+aire(i,jjp1)
	808	ENDDO
	809	end if
	810
	811	globaverage2d = 0.
	812	DO j=2,jjp1-1
	813	DO i = 1, iim
	814	globaverage2d = globaverage2d + var(i,j)*aire(i,j)
	815	END DO
	816	END DO
	817
	818	DO i=1,iim
	819	globaverage2d = globaverage2d + var(i,1)*aire(i,1)
	820	globaverage2d = globaverage2d + var(i,jjp1)*aire(i,jjp1)
	821	ENDDO
	822
	823	globaverage2d = globaverage2d/airetot
	824	return
	825	end
	826

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

Download in other formats:

Original Format