Context Navigation

leapfrog_nogcm.F @ 2919

Last change on this file since 2919 was 2108, checked in by emillour, 6 years ago
Follow-up of r2016 commit of nogcm main program: adapt leapfrog_gcm so that it does not explicitely depend on "phymars" and thus causes compilation issues when using other physics packages. EM
File size: 33.5 KB

Rev	Line
[2106]	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
[2108]	14	USE infotrac, ONLY: nqtot,ok_iso_verif,tname
[2106]	15	USE guide_mod, ONLY : guide_main
	16	USE write_field, ONLY: writefield
	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 sponge_mod, only: callsponge,mode_sponge,sponge
	26	use comuforc_h
	27	USE comvert_mod, ONLY: ap,bp,pressure_exner,presnivs,
	28	& aps,bps,presnivs,pseudoalt,preff,scaleheight
	29	USE comconst_mod, ONLY: daysec,dtvr,dtphys,dtdiss,
[2108]	30	. cpp,ihf,iflag_top_bound,pi,kappa,r
[2106]	31	USE logic_mod, ONLY: iflag_phys,ok_guide,forward,leapf,apphys,
	32	. statcl,conser,purmats,tidal,ok_strato
	33	USE temps_mod, ONLY: jD_ref,jH_ref,itaufin,day_ini,day_ref,
	34	. start_time,dt
	35
	36
	37
	38
	39	IMPLICIT NONE
	40
	41	c ...... Version du 10/01/98 ..........
	42
	43	c avec coordonnees verticales hybrides
	44	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
	45
	46	c=======================================================================
	47	c
	48	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
	49	c -------
	50	c
	51	c Objet:
	52	c ------
	53	c
	54	c GCM LMD nouvelle grille
	55	c
	56	c=======================================================================
	57	c
	58	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
	59	c et possibilite d'appeler une fonction f(y) a derivee tangente
	60	c hyperbolique a la place de la fonction a derivee sinusoidale.
	61
	62	c ... Possibilite de choisir le shema pour l'advection de
	63	c q , en modifiant iadv dans traceur.def (10/02) .
	64	c
	65	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
	66	c Pour Van-Leer iadv=10
	67	c
	68	c-----------------------------------------------------------------------
	69	c Declarations:
	70	c -------------
	71
	72	#include "dimensions.h"
	73	#include "paramet.h"
	74	#include "comdissnew.h"
	75	#include "comgeom.h"
	76	!#include "com_io_dyn.h"
	77	#include "iniprint.h"
	78	#include "academic.h"
	79
	80	! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
	81	! #include "clesphys.h"
	82
	83	REAL,INTENT(IN) :: time_0 ! not used
	84
	85	c dynamical variables:
	86	REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
	87	REAL,INTENT(INOUT) :: vcov(ip1jm,llm) ! meridional covariant wind
	88	REAL,INTENT(INOUT) :: teta(ip1jmp1,llm) ! potential temperature
	89	REAL,INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
	90	REAL,INTENT(INOUT) :: masse(ip1jmp1,llm) ! air mass
	91	REAL,INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
	92	REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
	93
	94	REAL p (ip1jmp1,llmp1 ) ! interlayer pressure
	95	REAL pks(ip1jmp1) ! exner at the surface
	96	REAL pk(ip1jmp1,llm) ! exner at mid-layer
	97	REAL pkf(ip1jmp1,llm) ! filtered exner at mid-layer
	98	REAL phi(ip1jmp1,llm) ! geopotential
	99	REAL w(ip1jmp1,llm) ! vertical velocity
	100	REAL kpd(ip1jmp1) ! TB15 = exp (-z/H)
	101
	102	! ADAPTATION GCM POUR CP(T)
	103	REAL temp(ip1jmp1,llm) ! temperature
	104	REAL tsurpk(ip1jmp1,llm) ! cpp*T/pk
	105
	106	! real zqmin,zqmax
	107
	108
	109	! ED18 nogcm
	110	REAL tau_ps
	111	REAL tau_co2
	112	REAL tau_teta
	113	REAL tetadpmean
	114	REAL tetadp(ip1jmp1,llm)
	115	REAL dpmean
	116	REAL dp(ip1jmp1,llm)
	117	REAL tetamean
	118	real globaverage2d
	119	! LOGICAL firstcall_globaverage2d
	120
	121
	122	c variables dynamiques intermediaire pour le transport
	123	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
	124
	125	c variables dynamiques au pas -1
	126	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
	127	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
	128	REAL massem1(ip1jmp1,llm)
	129
	130	c tendances dynamiques en */s
	131	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
	132	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot)
	133
	134	c tendances de la dissipation en */s
	135	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
	136	REAL dtetadis(ip1jmp1,llm)
	137
	138	c tendances de la couche superieure */s
	139	c REAL dvtop(ip1jm,llm)
	140	REAL dutop(ip1jmp1,llm)
	141	c REAL dtetatop(ip1jmp1,llm)
	142	c REAL dqtop(ip1jmp1,llm,nqtot),dptop(ip1jmp1)
	143
	144	c TITAN : tendances due au forces de marees */s
	145	REAL dvtidal(ip1jm,llm),dutidal(ip1jmp1,llm)
	146
	147	c tendances physiques */s
	148	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
	149	REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
	150
	151	c variables pour le fichier histoire
	152	! REAL dtav ! intervalle de temps elementaire
	153
	154	REAL tppn(iim),tpps(iim),tpn,tps
	155	c
	156	INTEGER itau,itaufinp1,iav
	157	! INTEGER iday ! jour julien
	158	REAL time
	159
	160	REAL SSUM
	161	! REAL finvmaold(ip1jmp1,llm)
	162
	163	cym LOGICAL lafin
	164	LOGICAL :: lafin=.false.
	165	INTEGER ij,iq,l
	166	! INTEGER ik
	167
	168	! real time_step, t_wrt, t_ops
	169
	170	REAL rdaym_ini
	171	! jD_cur: jour julien courant
	172	! jH_cur: heure julienne courante
	173	REAL :: jD_cur, jH_cur
	174	! INTEGER :: an, mois, jour
	175	! REAL :: secondes
	176
	177	LOGICAL first,callinigrads
	178	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
	179	save first
	180	data first/.true./
	181	! real dt_cum
	182	! character*10 infile
	183	! integer zan, tau0, thoriid
	184	! integer nid_ctesGCM
	185	! save nid_ctesGCM
	186	! real degres
	187	! real rlong(iip1), rlatg(jjp1)
	188	! real zx_tmp_2d(iip1,jjp1)
	189	! integer ndex2d(iip1*jjp1)
	190	logical ok_sync
	191	parameter (ok_sync = .true.)
	192	logical physics
	193
	194	data callinigrads/.true./
	195	character*10 string10
	196
	197	! REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
	198	REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
	199
	200	REAL psmean ! pression moyenne
	201	REAL pqco2mean ! moyenne globale ps*qco
	202	REAL p0 ! pression de reference
	203	REAL p00d ! globalaverage(kpd)
	204	REAL qmean_co2,qmean_co2_vert ! mass mean mixing ratio vap co2
	205	REAL pqco2(ip1jmp1) ! average co2 mass index : ps*q_co2
	206	REAL oldps(ip1jmp1) ! saving old pressure ps to calculate qch4
	207
	208	c+jld variables test conservation energie
	209	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
	210	C Tendance de la temp. potentiel d (theta)/ d t due a la
	211	C tansformation d'energie cinetique en energie thermique
	212	C cree par la dissipation
	213	REAL dtetaecdt(ip1jmp1,llm)
	214	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
	215	REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
	216	! REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
	217	CHARACTER*15 ztit
	218	!IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
	219	!IM SAVE ip_ebil_dyn
	220	!IM DATA ip_ebil_dyn/0/
	221	c-jld
	222
	223	! integer :: itau_w ! for write_paramLMDZ_dyn.h
	224
	225	! character*80 dynhist_file, dynhistave_file
	226	character(len=*),parameter :: modname="leapfrog"
	227	character*80 abort_message
	228
	229	logical dissip_conservative
	230	save dissip_conservative
	231	data dissip_conservative/.true./
	232
	233	INTEGER testita
	234	PARAMETER (testita = 9)
	235
	236	logical , parameter :: flag_verif = .false.
	237
	238	! for CP(T)
	239	real :: dtec
	240	real :: ztetaec(ip1jmp1,llm)
	241
	242	! TEMP : diagnostic mass
	243	real :: co2mass(iip1,jjp1)
	244	real :: co2ice_ij(iip1,jjp1)
[2108]	245	integer,save :: igcm_co2=0 ! index of CO2 tracer (if any)
[2106]	246	integer :: i,j,ig
	247	integer, parameter :: ngrid = 2+(jjm-1)*iim
	248	! local mass
	249	real :: mq(ip1jmp1,llm)
	250
	251	if (nday>=0) then
	252	itaufin = nday*day_step
	253	else
	254	! to run a given (-nday) number of dynamical steps
	255	itaufin = -nday
	256	endif
	257	if (less1day) then
	258	c MODIF VENUS: to run less than one day:
	259	itaufin = int(fractday*day_step)
	260	endif
	261	if (ndynstep.gt.0) then
	262	! running a given number of dynamical steps
	263	itaufin=ndynstep
	264	endif
	265	itaufinp1 = itaufin +1
	266
	267
	268
	269
	270	c INITIALISATIONS
	271	dudis(:,:) =0.
	272	dvdis(:,:) =0.
	273	dtetadis(:,:)=0.
	274	dutop(:,:) =0.
	275	c dvtop(:,:) =0.
	276	c dtetatop(:,:)=0.
	277	c dqtop(:,:,:) =0.
	278	c dptop(:) =0.
	279	dufi(:,:) =0.
	280	dvfi(:,:) =0.
	281	dtetafi(:,:)=0.
	282	dqfi(:,:,:) =0.
	283	dpfi(:) =0.
	284
	285	itau = 0
	286	physics=.true.
	287	if (iflag_phys==0.or.iflag_phys==2) physics=.false.
	288
	289	! ED18 nogcm
	290	! firstcall_globaverage2d = .true.
	291
	292	c iday = day_ini+itau/day_step
	293	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	294	c IF(time.GT.1.) THEN
	295	c time = time-1.
	296	c iday = iday+1
	297	c ENDIF
	298
	299
	300	c-----------------------------------------------------------------------
	301	c On initialise la pression et la fonction d'Exner :
	302	c --------------------------------------------------
	303
	304	dq(:,:,:)=0.
	305	CALL pression ( ip1jmp1, ap, bp, ps, p )
	306	if (pressure_exner) then
	307	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
	308	else
	309	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
	310	endif
	311
	312	c------------------
	313	c TEST PK MONOTONE
	314	c------------------
	315	write(,) "Test PK"
	316	do ij=1,ip1jmp1
	317	do l=2,llm
	318	! PRINT*,'pk(ij,l) = ',pk(ij,l)
	319	! PRINT*,'pk(ij,l-1) = ',pk(ij,l-1)
	320	if(pk(ij,l).gt.pk(ij,l-1)) then
	321	c write(,) ij,l,pk(ij,l)
	322	abort_message = 'PK non strictement decroissante'
	323	call abort_gcm(modname,abort_message,1)
	324	c write(,) "ATTENTION, Test PK deconnecte..."
	325	endif
	326	enddo
	327	enddo
	328	write(,) "Fin Test PK"
	329	c stop
	330	c------------------
	331	c Preparing mixing of pressure and tracers in nogcm
	332	c kpd=exp(-z/H) Needed to define a reference pressure :
	333	c P0=pmean/globalaverage(kpd)
	334	c thus P(i) = P0exp(-z(i)/H) = P0kpd(i)
	335	c It is checked that globalaverage2d(Pi)=pmean
	336	DO ij=1,ip1jmp1
	337	kpd(ij) = exp(-phis(ij)/(r*200.))
	338	ENDDO
	339	p00d=globaverage2d(kpd) ! mean pres at ref level
	340	tau_ps = 1. ! constante de rappel for pressure (s)
	341	tau_co2 = 1.E5 !E5 ! constante de rappel for mix ratio qco2 (s)
	342	tau_teta = 1.E7 !constante de rappel for potentiel temperature
	343
	344	! ED18 TEST
	345	! PRINT*,'igcm_co2 = ',igcm_co2
[2108]	346	! Locate tracer "co2" and set igcm_co2:
	347	do iq=1,nqtot
	348	if (tname(iq)=="co2") then
	349	igcm_co2=iq
	350	exit
	351	endif
	352	enddo
[2106]	353
	354	c-----------------------------------------------------------------------
	355	c Debut de l'integration temporelle:
	356	c ----------------------------------
	357
	358	1 CONTINUE ! Matsuno Forward step begins here
	359
	360	c date: (NB: date remains unchanged for Backward step)
	361	c -----
	362
	363	jD_cur = jD_ref + day_ini - day_ref + &
	364	& (itau+1)/day_step
	365	jH_cur = jH_ref + start_time + &
	366	& mod(itau+1,day_step)/float(day_step)
	367	jD_cur = jD_cur + int(jH_cur)
	368	jH_cur = jH_cur - int(jH_cur)
	369
	370	c
	371
	372	! Save fields obtained at previous time step as '...m1'
	373	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
	374	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
	375	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
	376	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
	377	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
	378
	379	forward = .TRUE.
	380	leapf = .FALSE.
	381	dt = dtvr
	382
	383	2 CONTINUE ! Matsuno backward or leapfrog step begins here
	384
	385	c-----------------------------------------------------------------------
	386
	387	c date: (NB: only leapfrog step requires recomputing date)
	388	c -----
	389
	390	IF (leapf) THEN
	391	jD_cur = jD_ref + day_ini - day_ref +
	392	& (itau+1)/day_step
	393	jH_cur = jH_ref + start_time +
	394	& mod(itau+1,day_step)/float(day_step)
	395	jD_cur = jD_cur + int(jH_cur)
	396	jH_cur = jH_cur - int(jH_cur)
	397	ENDIF
	398
	399
	400	c gestion des appels de la physique et des dissipations:
	401	c ------------------------------------------------------
	402	c
	403	c ... P.Le Van ( 6/02/95 ) ....
	404
	405	! ED18: suppression des mentions de la variable apdiss dans le cas
	406	! 'nogcm'
	407
	408	apphys = .FALSE.
	409	statcl = .FALSE.
	410	conser = .FALSE.
	411
	412
	413	IF( purmats ) THEN
	414	! Purely Matsuno time stepping
	415	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
	416
	417
	418	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
	419	s .and. physics ) apphys = .TRUE.
	420	ELSE
	421	! Leapfrog/Matsuno time stepping
	422	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
	423
	424
	425	IF( MOD(itau+1,iphysiq).EQ.0.AND.physics) apphys=.TRUE.
	426	END IF
	427
	428	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
	429	! supress dissipation step
	430
	431
	432
	433	c-----------------------------------------------------------------------
	434	c calcul des tendances dynamiques:
	435	c --------------------------------
	436
	437	! ED18: suppression de l'onglet pour le cas nogcm
	438
	439
	440	dv(:,:) = 0.
	441	du(:,:) = 0.
	442	dteta(:,:) = 0.
	443	dq(:,:,:) = 0.
	444
	445
	446
	447	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
	448	c
	449	c-----------------------------------------------------------------------
	450	c calcul des tendances physiques:
	451	c -------------------------------
	452	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
	453	c
	454	IF( purmats ) THEN
	455	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
	456	ELSE
	457	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
	458	ENDIF
	459	c
	460	c
	461	IF( apphys ) THEN
	462	c
	463	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
	464	c
	465
	466	CALL pression ( ip1jmp1, ap, bp, ps, p )
	467	if (pressure_exner) then
	468	CALL exner_hyb( ip1jmp1, ps, p,pks, pk, pkf )
	469	else
	470	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
	471	endif
	472
	473	! Compute geopotential (physics might need it)
	474	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	475
	476	jD_cur = jD_ref + day_ini - day_ref + &
	477	& (itau+1)/day_step
	478
	479	IF ((planet_type .eq."generic").or.
	480	& (planet_type .eq."mars")) THEN
	481	! AS: we make jD_cur to be pday
	482	jD_cur = int(day_ini + itau/day_step)
	483	ENDIF
	484
	485	! print*,'itau =',itau
	486	! print*,'day_step =',day_step
	487	! print*,'itau/day_step =',itau/day_step
	488
	489
	490	jH_cur = jH_ref + start_time + &
	491	& mod(itau+1,day_step)/float(day_step)
	492	IF ((planet_type .eq."generic").or.
	493	& (planet_type .eq."mars")) THEN
	494	jH_cur = jH_ref + start_time + &
	495	& mod(itau,day_step)/float(day_step)
	496	ENDIF
	497	jD_cur = jD_cur + int(jH_cur)
	498	jH_cur = jH_cur - int(jH_cur)
	499
	500
	501
	502	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
	503	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
	504	! write(lunout,*)'current date = ',an, mois, jour, secondes
	505
	506	c rajout debug
	507	c lafin = .true.
	508
	509
	510	c Interface avec les routines de phylmd (phymars ... )
	511	c -----------------------------------------------------
	512
	513	c+jld
	514
	515	c Diagnostique de conservation de l'Energie : initialisation
	516	IF (ip_ebil_dyn.ge.1 ) THEN
	517	ztit='bil dyn'
	518	! Ehouarn: be careful, diagedyn is Earth-specific!
	519	IF (planet_type.eq."earth") THEN
	520	CALL diagedyn(ztit,2,1,1,dtphys
	521	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	522	ENDIF
	523	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
	524	c-jld
	525	#ifdef CPP_IOIPSL
	526	cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
	527	cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
	528	c IF (first) THEN
	529	c first=.false.
	530	c#include "ini_paramLMDZ_dyn.h"
	531	c ENDIF
	532	c
	533	c#include "write_paramLMDZ_dyn.h"
	534	c
	535	#endif
	536	! #endif of #ifdef CPP_IOIPSL
	537
	538	c call WriteField('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
	539	! print*,'---------LEAPFROG---------------'
	540	! print*,''
	541	! print*,'> AVANT CALFIS'
	542	! print*,''
	543	! print*,'teta(3049,:) = ',teta(3049,:)
	544	! print*,''
	545	! print*,'dteta(3049,:) = ',dteta(3049,:)
	546	! print*,''
	547	! print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
	548	! print*,''
	549
	550
	551	CALL calfis( lafin , jD_cur, jH_cur,
	552	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
	553	$ du,dv,dteta,dq,
	554	$ flxw,
	555	$ dufi,dvfi,dtetafi,dqfi,dpfi )
	556
	557	c call WriteField('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
	558	c call WriteField('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
	559	c call WriteField('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
	560	! print*,'> APRES CALFIS (AVANT ADDFI)'
	561	! print*,''
	562	! print*,'teta(3049,:) = ',teta(3049,:)
	563	! print*,''
	564	! print*,'dteta(3049,:) = ',dteta(3049,:)
	565	! print*,''
	566	! print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
	567	! print*,''
	568
	569
	570	c ajout des tendances physiques
	571	c ------------------------------
	572
	573	CALL addfi( dtphys, leapf, forward ,
	574	$ ucov, vcov, teta , q ,ps ,
	575	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	576
	577	! print*,'> APRES ADDFI'
	578	! print*,''
	579	! print*,'teta(3049,:) = ',teta(3049,:)
	580	! print*,''
	581	! print*,'dteta(3049,:) = ',dteta(3049,:)
	582	! print*,''
	583	! print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
	584	! print*,''
	585
	586	CALL pression (ip1jmp1,ap,bp,ps,p)
	587	CALL massdair(p,masse)
	588
	589	! Mass of tracers in each mess (kg) before Horizontal mixing of pressure
	590	c -------------------------------------------------------------------
	591
	592	DO l=1, llm
	593	DO ij=1,ip1jmp1
	594	mq(ij,l) = masse(ij,l)*q(ij,l,igcm_co2)
	595	ENDDO
	596	ENDDO
	597
	598	c Horizontal mixing of pressure
	599	c ------------------------------
	600	! Rappel newtonien vers psmean
	601	psmean= globaverage2d(ps) ! mean pressure
	602	p0=psmean/p00d
	603	DO ij=1,ip1jmp1
	604	oldps(ij)=ps(ij)
	605	ps(ij)=ps(ij) +(p0*kpd(ij)
	606	& -ps(ij))*(1-exp(-dtphys/tau_ps))
	607	ENDDO
	608
	609	write(72,*) 'psmean = ',psmean ! mean pressure redistributed
	610
	611	! security check: test pressure negative
	612	DO ij=1,ip1jmp1
	613	IF (ps(ij).le.0.) then
	614	!write(,) 'Negative pressure :'
	615	!write(,) 'ps = ',ps(ij),' ig = ',ij
	616	!write(,) 'pmean = ',p0*kpd(ij)
	617	!write(,) 'tau_ps = ',tau_ps
	618	!STOP
	619	ps(ij)=0.0000001*kpd(ij)/p00d
	620	ENDIF
	621	ENDDO
	622	!***********************
	623	! Correction on teta due to surface pressure changes
	624	DO l = 1,llm
	625	DO ij = 1,ip1jmp1
	626	teta(ij,l)= teta(ij,l)(ps(ij)/oldps(ij))*kappa
	627	ENDDO
	628	ENDDO
	629	!***********************
	630
	631
	632
	633	! ! update pressure and update p and pk
	634	! DO ij=1,ip1jmp1
	635	! ps(ij) = ps(ij) + dpfi(ij)*dtphys
	636	! ENDDO
	637	CALL pression (ip1jmp1,ap,bp,ps,p)
	638	CALL massdair(p,masse)
	639	if (pressure_exner) then
	640	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
	641	else
	642	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
	643	endif
	644
	645	! Update tracers after Horizontal mixing of pressure to ! conserve tracer mass
	646	c -------------------------------------------------------------------
	647	DO l=1, llm
	648	DO ij=1,ip1jmp1
	649	q(ij,l,igcm_co2) = mq(ij,l)/ masse(ij,l)
	650	ENDDO
	651	ENDDO
	652
	653
	654
	655	c Horizontal mixing of pressure
	656	c ------------------------------
	657	! Rappel newtonien vers psmean
	658	psmean= globaverage2d(ps) ! mean pressure
	659	! ! increment q_co2 with physical tendancy
	660	! IF (igcm_co2.ne.0) then
	661	! DO l=1, llm
	662	! DO ij=1,ip1jmp1
	663	! q(ij,l,igcm_co2)=q(ij,l,igcm_co2)+
	664	! & dqfi(ij,l,igcm_co2)*dtphys
	665	! ENDDO
	666	! ENDDO
	667	! ENDIF
	668
	669	c Mixing CO2 vertically
	670	c --------------------------
	671	if (igcm_co2.ne.0) then
	672	DO ij=1,ip1jmp1
	673	qmean_co2_vert=0.
	674	DO l=1, llm
	675	qmean_co2_vert= qmean_co2_vert
	676	& + q(ij,l,igcm_co2)*( p(ij,l) - p(ij,l+1))
	677	END DO
	678	qmean_co2_vert= qmean_co2_vert/ps(ij)
	679	DO l=1, llm
	680	q(ij,l,igcm_co2)= qmean_co2_vert
	681	END DO
	682	END DO
	683	end if
	684
	685
	686
	687	c Horizontal mixing of pressure
	688	c ------------------------------
	689	! Rappel newtonien vers psmean
	690	psmean= globaverage2d(ps) ! mean pressure
	691
	692	c Horizontal mixing tracers, and temperature (replace dynamics in nogcm)
	693	c ---------------------------------------------------------------------------------
	694
	695	! Simulate redistribution by dynamics for qco2
	696	if (igcm_co2.ne.0) then
	697
	698	DO ij=1,ip1jmp1
	699	pqco2(ij)= ps(ij) * q(ij,1,igcm_co2)
	700	ENDDO
	701	pqco2mean=globaverage2d(pqco2)
	702
	703	! Rappel newtonien vers qco2_mean
	704	qmean_co2= pqco2mean / psmean
	705
	706	DO ij=1,ip1jmp1
	707	q(ij,1,igcm_co2)=q(ij,1,igcm_co2)+
	708	& (qmean_co2-q(ij,1,igcm_co2))*
	709	& (1.-exp(-dtphys/tau_co2))
	710	ENDDO
	711
	712	DO l=2, llm
	713	DO ij=1,ip1jmp1
	714	q(ij,l,igcm_co2)=q(ij,1,igcm_co2)
	715	END DO
	716	END DO
	717
	718	! TEMPORAIRE (ED)
	719	! PRINT*,'psmean = ',psmean
	720	! PRINT*,'qmean_co2 = ',qmean_co2
	721	! PRINT*,'pqco2mean = ',pqco2mean
	722	! PRINT*,'q(50,1,igcm_co2) = ',q(50,1,igcm_co2)
	723	! PRINT*,'q(50,2,igcm_co2) = ',q(50,2,igcm_co2)
	724	! PRINT*,'q(50,3,igcm_co2) = ',q(50,3,igcm_co2)
	725
	726	endif ! igcm_co2.ne.0
	727
	728
	729	! ********************************************************
	730
	731	c Horizontal mixing of pressure
	732	c ------------------------------
	733	! Rappel newtonien vers psmean
	734	psmean= globaverage2d(ps) ! mean pressure
	735
	736
	737	c Mixing Temperature horizontally
	738	c -------------------------------
	739	! initialize variables that will be averaged
	740	DO l=1,llm
	741	DO ij=1,ip1jmp1
	742	dp(ij,l) = p(ij,l) - p(ij,l+1)
	743	tetadp(ij,l) = teta(ij,l)*dp(ij,l)
	744	ENDDO
	745	ENDDO
	746
	747	DO l=1,llm
	748	tetadpmean = globaverage2d(tetadp(:,l))
	749	dpmean = globaverage2d(dp(:,l))
	750	tetamean = tetadpmean / dpmean
	751	DO ij=1,ip1jmp1
	752	teta(ij,l) = teta(ij,l) + (tetamean - teta(ij,l)) *
	753	& (1 - exp(-dtphys/tau_teta))
	754	ENDDO
	755	ENDDO
	756
	757
	758	ENDIF ! of IF( apphys )
	759
	760
	761	c ********************************************************************
	762	c ********************************************************************
	763	c .... fin de l'integration dynamique et physique pour le pas itau ..
	764	c ********************************************************************
	765	c ********************************************************************
	766
	767
	768	IF ( .NOT.purmats ) THEN
	769	c ........................................................
	770	c .............. schema matsuno + leapfrog ..............
	771	c ........................................................
	772
	773	IF(forward. OR. leapf) THEN
	774	itau= itau + 1
	775	c iday= day_ini+itau/day_step
	776	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	777	c IF(time.GT.1.) THEN
	778	c time = time-1.
	779	c iday = iday+1
	780	c ENDIF
	781	ENDIF
	782
	783
	784	IF( itau. EQ. itaufinp1 ) then
	785	if (flag_verif) then
	786	write(79,*) 'ucov',ucov
	787	write(80,*) 'vcov',vcov
	788	write(81,*) 'teta',teta
	789	write(82,*) 'ps',ps
	790	write(83,*) 'q',q
	791	WRITE(85,*) 'q1 = ',q(:,:,1)
	792	WRITE(86,*) 'q3 = ',q(:,:,3)
	793	endif
	794
	795	abort_message = 'Simulation finished'
	796
	797	call abort_gcm(modname,abort_message,0)
	798	ENDIF
	799	c-----------------------------------------------------------------------
	800	c ecriture du fichier histoire moyenne:
	801	c -------------------------------------
	802
	803	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	804	IF(itau.EQ.itaufin) THEN
	805	iav=1
	806	ELSE
	807	iav=0
	808	ENDIF
	809
	810	! ! Ehouarn: re-compute geopotential for outputs
	811	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	812
	813	IF (ok_dynzon) THEN
	814	#ifdef CPP_IOIPSL
	815	c les traceurs ne sont pas sortis, trop lourd.
	816	c Peut changer eventuellement si besoin.
	817	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
	818	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
	819	& du,dudis,dutop,dufi)
	820	#endif
	821	END IF
	822	IF (ok_dyn_ave) THEN
	823	#ifdef CPP_IOIPSL
	824	CALL writedynav(itau,vcov,
	825	& ucov,teta,pk,phi,q,masse,ps,phis)
	826	#endif
	827	ENDIF
	828
	829	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
	830
	831	if (ok_iso_verif) then
	832	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
	833	endif !if (ok_iso_verif) then
	834
	835	c-----------------------------------------------------------------------
	836	c ecriture de la bande histoire:
	837	c ------------------------------
	838
	839	IF( MOD(itau,iecri).EQ.0) THEN
	840	! Ehouarn: output only during LF or Backward Matsuno
	841	if (leapf.or.(.not.leapf.and.(.not.forward))) then
	842	! ADAPTATION GCM POUR CP(T)
	843	call tpot2t(ijp1llm,teta,temp,pk)
	844	tsurpk = cpp*temp/pk
	845	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
	846	unat=0.
	847	do l=1,llm
	848	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	849	vnat(:,l)=vcov(:,l)/cv(:)
	850	enddo
	851	#ifdef CPP_IOIPSL
	852	if (ok_dyn_ins) then
	853	! write(lunout,*) "leapfrog: call writehist, itau=",itau
	854	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	855	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	856	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
	857	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
	858	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
	859	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
	860	endif ! of if (ok_dyn_ins)
	861	#endif
	862	! For some Grads outputs of fields
	863	if (output_grads_dyn) then
	864	#include "write_grads_dyn.h"
	865	endif
	866	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
	867	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
	868
	869	IF(itau.EQ.itaufin) THEN
	870
	871	if (planet_type=="mars") then
	872	CALL dynredem1("restart.nc",REAL(itau)/REAL(day_step),
	873	& vcov,ucov,teta,q,masse,ps)
	874	else
	875	CALL dynredem1("restart.nc",start_time,
	876	& vcov,ucov,teta,q,masse,ps)
	877	endif
	878	CLOSE(99)
	879	!!! Ehouarn: Why not stop here and now?
	880	ENDIF ! of IF (itau.EQ.itaufin)
	881
	882	c-----------------------------------------------------------------------
	883	c gestion de l'integration temporelle:
	884	c ------------------------------------
	885
	886	IF( MOD(itau,iperiod).EQ.0 ) THEN
	887	GO TO 1
	888	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	889
	890	IF( forward ) THEN
	891	c fin du pas forward et debut du pas backward
	892
	893	forward = .FALSE.
	894	leapf = .FALSE.
	895	GO TO 2
	896
	897	ELSE
	898	c fin du pas backward et debut du premier pas leapfrog
	899
	900	leapf = .TRUE.
	901	dt = 2.*dtvr
	902	GO TO 2
	903	END IF ! of IF (forward)
	904	ELSE
	905
	906	c ...... pas leapfrog .....
	907
	908	leapf = .TRUE.
	909	dt = 2.*dtvr
	910	GO TO 2
	911	END IF ! of IF (MOD(itau,iperiod).EQ.0)
	912	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
	913
	914	ELSE ! of IF (.not.purmats)
	915
	916	if (ok_iso_verif) then
	917	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
	918	endif !if (ok_iso_verif) then
	919
	920	c ........................................................
	921	c .............. schema matsuno ...............
	922	c ........................................................
	923	IF( forward ) THEN
	924
	925	itau = itau + 1
	926	c iday = day_ini+itau/day_step
	927	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	928	c
	929	c IF(time.GT.1.) THEN
	930	c time = time-1.
	931	c iday = iday+1
	932	c ENDIF
	933
	934	forward = .FALSE.
	935	IF( itau. EQ. itaufinp1 ) then
	936	abort_message = 'Simulation finished'
	937	call abort_gcm(modname,abort_message,0)
	938	ENDIF
	939	GO TO 2
	940
	941	ELSE ! of IF(forward) i.e. backward step
	942
	943	if (ok_iso_verif) then
	944	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
	945	endif !if (ok_iso_verif) then
	946
	947	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	948	IF(itau.EQ.itaufin) THEN
	949	iav=1
	950	ELSE
	951	iav=0
	952	ENDIF
	953
	954	! ! Ehouarn: re-compute geopotential for outputs
	955	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	956
	957	IF (ok_dynzon) THEN
	958	#ifdef CPP_IOIPSL
	959	c les traceurs ne sont pas sortis, trop lourd.
	960	c Peut changer eventuellement si besoin.
	961	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
	962	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
	963	& du,dudis,dutop,dufi)
	964	#endif
	965	ENDIF
	966	IF (ok_dyn_ave) THEN
	967	#ifdef CPP_IOIPSL
	968	CALL writedynav(itau,vcov,
	969	& ucov,teta,pk,phi,q,masse,ps,phis)
	970	#endif
	971	ENDIF
	972
	973	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
	974
	975	IF(MOD(itau,iecri ).EQ.0) THEN
	976	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
	977	! ADAPTATION GCM POUR CP(T)
	978	call tpot2t(ijp1llm,teta,temp,pk)
	979	tsurpk = cpp*temp/pk
	980	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
	981	unat=0.
	982	do l=1,llm
	983	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	984	vnat(:,l)=vcov(:,l)/cv(:)
	985	enddo
	986	#ifdef CPP_IOIPSL
	987	if (ok_dyn_ins) then
	988	! write(lunout,*) "leapfrog: call writehist (b)",
	989	! & itau,iecri
	990	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	991	endif ! of if (ok_dyn_ins)
	992	#endif
	993	! For some Grads outputs
	994	if (output_grads_dyn) then
	995	#include "write_grads_dyn.h"
	996	endif
	997
	998	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
	999
	1000	IF(itau.EQ.itaufin) THEN
	1001	if (planet_type=="mars") then
	1002	CALL dynredem1("restart.nc",REAL(itau)/REAL(day_step),
	1003	& vcov,ucov,teta,q,masse,ps)
	1004	else
	1005	CALL dynredem1("restart.nc",start_time,
	1006	& vcov,ucov,teta,q,masse,ps)
	1007	endif
	1008	ENDIF ! of IF(itau.EQ.itaufin)
	1009
	1010	forward = .TRUE.
	1011	GO TO 1
	1012
	1013	ENDIF ! of IF (forward)
	1014
	1015	END IF ! of IF(.not.purmats)
	1016
	1017	STOP
	1018	END
	1019
	1020	! ************************************************************
	1021	! globalaverage VERSION PLuto
	1022	! ************************************************************
	1023	FUNCTION globaverage2d(var)
	1024	! ************************************************************
	1025	IMPLICIT NONE
	1026	#include "dimensions.h"
	1027	#include "paramet.h"
	1028	#include "comgeom2.h"
	1029	!ip1jmp1 called in paramet.h = iip1 x jjp1
	1030	REAL var(iip1,jjp1), globaverage2d
	1031	integer i,j
	1032	REAL airetot
	1033	save airetot
	1034	logical firstcall
	1035	data firstcall/.true./
	1036	save firstcall
	1037
	1038	if (firstcall) then
	1039	firstcall=.false.
	1040	airetot =0.
	1041	DO j=2,jjp1-1 ! lat
	1042	DO i = 1, iim ! lon
	1043	airetot = airetot + aire(i,j)
	1044	END DO
	1045	END DO
	1046	DO i=1,iim
	1047	airetot=airetot + aire(i,1)+aire(i,jjp1)
	1048	ENDDO
	1049	end if
	1050
	1051	globaverage2d = 0.
	1052	DO j=2,jjp1-1
	1053	DO i = 1, iim
	1054	globaverage2d = globaverage2d + var(i,j)*aire(i,j)
	1055	END DO
	1056	END DO
	1057
	1058	DO i=1,iim
	1059	globaverage2d = globaverage2d + var(i,1)*aire(i,1)
	1060	globaverage2d = globaverage2d + var(i,jjp1)*aire(i,jjp1)
	1061	ENDDO
	1062
	1063	globaverage2d = globaverage2d/airetot
	1064	return
	1065	end
	1066

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

Download in other formats:

Original Format