Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

leapfrog_nogcm_pluto.F @ 3493

Last change on this file since 3493 was 3229, checked in by tbertrand, 9 months ago
Pluto PCM: Adding Nogcm.F : 2D surface model with 1 atm layer
File size: 32.7 KB

Rev	Line
[3229]	1	!
	2	! $Id: leapfrog.F 1446 2010-10-22 09:27:25Z emillour $
	3	!
	4	c
	5	c
	6	SUBROUTINE leapfrog_nogcm_pluto(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 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,
	30	. cpp,ihf,iflag_top_bound,pi,kappa,r
	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_n2
	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 pqn2mean ! moyenne globale ps*qco
	202	REAL p0 ! pression de reference
	203	REAL p00d ! globalaverage(kpd)
	204	REAL qmean_n2,qmean_n2_vert ! mass mean mixing ratio vap n2
	205	REAL pqn2(ip1jmp1) ! average n2 mass index : ps*q_n2
	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 :: n2mass(iip1,jjp1)
	244	real :: n2ice_ij(iip1,jjp1)
	245	integer,save :: igcm_n2=0 ! index of CO2 tracer (if any)
	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_n2 = 1.E5 !E5 ! constante de rappel for mix ratio qn2 (s)
	342	tau_teta = 1.E7 !constante de rappel for potentiel temperature
	343
	344	PRINT*,'igcm_n2 = ',igcm_n2
	345	do iq=1,nqtot
	346	if (tname(iq)=="n2") then
	347	igcm_n2=iq
	348	exit
	349	endif
	350	enddo
	351
	352	c-----------------------------------------------------------------------
	353	c Debut de l'integration temporelle:
	354	c ----------------------------------
	355
	356	1 CONTINUE ! Matsuno Forward step begins here
	357
	358	c date: (NB: date remains unchanged for Backward step)
	359	c -----
	360
	361	jD_cur = jD_ref + day_ini - day_ref + &
	362	& (itau+1)/day_step
	363	jH_cur = jH_ref + start_time + &
	364	& mod(itau+1,day_step)/float(day_step)
	365	jD_cur = jD_cur + int(jH_cur)
	366	jH_cur = jH_cur - int(jH_cur)
	367
	368	c
	369
	370	! Save fields obtained at previous time step as '...m1'
	371	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
	372	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
	373	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
	374	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
	375	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
	376
	377	forward = .TRUE.
	378	leapf = .FALSE.
	379	dt = dtvr
	380
	381	2 CONTINUE ! Matsuno backward or leapfrog step begins here
	382
	383	c-----------------------------------------------------------------------
	384
	385	c date: (NB: only leapfrog step requires recomputing date)
	386	c -----
	387
	388	IF (leapf) THEN
	389	jD_cur = jD_ref + day_ini - day_ref +
	390	& (itau+1)/day_step
	391	jH_cur = jH_ref + start_time +
	392	& mod(itau+1,day_step)/float(day_step)
	393	jD_cur = jD_cur + int(jH_cur)
	394	jH_cur = jH_cur - int(jH_cur)
	395	ENDIF
	396
	397
	398	c gestion des appels de la physique et des dissipations:
	399	c ------------------------------------------------------
	400	c
	401	c ... P.Le Van ( 6/02/95 ) ....
	402
	403	! ED18: suppression des mentions de la variable apdiss dans le cas
	404	! 'nogcm'
	405
	406	apphys = .FALSE.
	407	statcl = .FALSE.
	408	conser = .FALSE.
	409
	410
	411	IF( purmats ) THEN
	412	! Purely Matsuno time stepping
	413	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
	414
	415
	416	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
	417	s .and. physics ) apphys = .TRUE.
	418	ELSE
	419	! Leapfrog/Matsuno time stepping
	420	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
	421
	422
	423	IF( MOD(itau+1,iphysiq).EQ.0.AND.physics) apphys=.TRUE.
	424	END IF
	425
	426	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
	427	! supress dissipation step
	428
	429
	430
	431	c-----------------------------------------------------------------------
	432	c calcul des tendances dynamiques:
	433	c --------------------------------
	434
	435	! ED18: suppression de l'onglet pour le cas nogcm
	436
	437
	438	dv(:,:) = 0.
	439	du(:,:) = 0.
	440	dteta(:,:) = 0.
	441	dq(:,:,:) = 0.
	442
	443
	444
	445	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
	446	c
	447	c-----------------------------------------------------------------------
	448	c calcul des tendances physiques:
	449	c -------------------------------
	450	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
	451	c
	452	IF( purmats ) THEN
	453	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
	454	ELSE
	455	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
	456	ENDIF
	457	c
	458	c
	459	IF( apphys ) THEN
	460	c
	461	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
	462	c
	463
	464	CALL pression ( ip1jmp1, ap, bp, ps, p )
	465	if (pressure_exner) then
	466	CALL exner_hyb( ip1jmp1, ps, p,pks, pk, pkf )
	467	else
	468	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
	469	endif
	470
	471	! Compute geopotential (physics might need it)
	472	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	473
	474	jD_cur = jD_ref + day_ini - day_ref + &
	475	& (itau+1)/day_step
	476
	477	! AS: we make jD_cur to be pday
	478	jD_cur = int(day_ini + itau/day_step)
	479
	480	! print*,'itau =',itau
	481	! print*,'day_step =',day_step
	482	! print*,'itau/day_step =',itau/day_step
	483
	484
	485	jH_cur = jH_ref + start_time + &
	486	& mod(itau+1,day_step)/float(day_step)
	487	jH_cur = jH_ref + start_time + &
	488	& mod(itau,day_step)/float(day_step)
	489	jD_cur = jD_cur + int(jH_cur)
	490	jH_cur = jH_cur - int(jH_cur)
	491
	492
	493
	494	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
	495	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
	496	! write(lunout,*)'current date = ',an, mois, jour, secondes
	497
	498	c rajout debug
	499	c lafin = .true.
	500
	501
	502	c Interface avec les routines de phylmd (phymars ... )
	503	c -----------------------------------------------------
	504
	505	c+jld
	506
	507	c Diagnostique de conservation de l'Energie : initialisation
	508	IF (ip_ebil_dyn.ge.1 ) THEN
	509	ztit='bil dyn'
	510	! Ehouarn: be careful, diagedyn is Earth-specific!
	511	IF (planet_type.eq."earth") THEN
	512	CALL diagedyn(ztit,2,1,1,dtphys
	513	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
	514	ENDIF
	515	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
	516	c-jld
	517	#ifdef CPP_IOIPSL
	518	cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
	519	cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
	520	c IF (first) THEN
	521	c first=.false.
	522	c#include "ini_paramLMDZ_dyn.h"
	523	c ENDIF
	524	c
	525	c#include "write_paramLMDZ_dyn.h"
	526	c
	527	#endif
	528	! #endif of #ifdef CPP_IOIPSL
	529
	530	c call WriteField('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
	531	! print*,'---------LEAPFROG---------------'
	532	! print*,''
	533	! print*,'> AVANT CALFIS'
	534	! print*,''
	535	! print*,'teta(3049,:) = ',teta(3049,:)
	536	! print*,''
	537	! print*,'dteta(3049,:) = ',dteta(3049,:)
	538	! print*,''
	539	! print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
	540	! print*,''
	541
	542
	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,
	547	$ dufi,dvfi,dtetafi,dqfi,dpfi )
	548
	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	! print*,'> APRES CALFIS (AVANT ADDFI)'
	553	! print*,''
	554	! print*,'teta(3049,:) = ',teta(3049,:)
	555	! print*,''
	556	! print*,'dteta(3049,:) = ',dteta(3049,:)
	557	! print*,''
	558	! print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
	559	! print*,''
	560
	561
	562	c ajout des tendances physiques
	563	c ------------------------------
	564
	565	CALL addfi( dtphys, leapf, forward ,
	566	$ ucov, vcov, teta , q ,ps ,
	567	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
	568
	569	! print*,'> APRES ADDFI'
	570	! print*,''
	571	! print*,'teta(3049,:) = ',teta(3049,:)
	572	! print*,''
	573	! print*,'dteta(3049,:) = ',dteta(3049,:)
	574	! print*,''
	575	! print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
	576	! print*,''
	577
	578	CALL pression (ip1jmp1,ap,bp,ps,p)
	579	CALL massdair(p,masse)
	580
	581	! Mass of tracers in each mess (kg) before Horizontal mixing of pressure
	582	c -------------------------------------------------------------------
	583
	584	DO l=1, llm
	585	DO ij=1,ip1jmp1
	586	mq(ij,l) = masse(ij,l)*q(ij,l,igcm_n2)
	587	ENDDO
	588	ENDDO
	589
	590	c Horizontal mixing of pressure
	591	c ------------------------------
	592	! Rappel newtonien vers psmean
	593	psmean= globaverage2d(ps) ! mean pressure
	594	p0=psmean/p00d
	595	DO ij=1,ip1jmp1
	596	oldps(ij)=ps(ij)
	597	ps(ij)=ps(ij) +(p0*kpd(ij)
	598	& -ps(ij))*(1-exp(-dtphys/tau_ps))
	599	ENDDO
	600
	601	write(72,*) 'psmean = ',psmean ! mean pressure redistributed
	602
	603	! security check: test pressure negative
	604	DO ij=1,ip1jmp1
	605	IF (ps(ij).le.0.) then
	606	!write(,) 'Negative pressure :'
	607	!write(,) 'ps = ',ps(ij),' ig = ',ij
	608	!write(,) 'pmean = ',p0*kpd(ij)
	609	!write(,) 'tau_ps = ',tau_ps
	610	!STOP
	611	ps(ij)=0.0000001*kpd(ij)/p00d
	612	ENDIF
	613	ENDDO
	614	!***********************
	615	! Correction on teta due to surface pressure changes
	616	DO l = 1,llm
	617	DO ij = 1,ip1jmp1
	618	teta(ij,l)= teta(ij,l)(ps(ij)/oldps(ij))*kappa
	619	ENDDO
	620	ENDDO
	621	!***********************
	622
	623
	624
	625	! ! update pressure and update p and pk
	626	! DO ij=1,ip1jmp1
	627	! ps(ij) = ps(ij) + dpfi(ij)*dtphys
	628	! ENDDO
	629	CALL pression (ip1jmp1,ap,bp,ps,p)
	630	CALL massdair(p,masse)
	631	if (pressure_exner) then
	632	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
	633	else
	634	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
	635	endif
	636
	637	! Update tracers after Horizontal mixing of pressure to ! conserve tracer mass
	638	c -------------------------------------------------------------------
	639	DO l=1, llm
	640	DO ij=1,ip1jmp1
	641	q(ij,l,igcm_n2) = mq(ij,l)/ masse(ij,l)
	642	ENDDO
	643	ENDDO
	644
	645
	646
	647	c Horizontal mixing of pressure
	648	c ------------------------------
	649	! Rappel newtonien vers psmean
	650	psmean= globaverage2d(ps) ! mean pressure
	651	! ! increment q_n2 with physical tendancy
	652	! IF (igcm_n2.ne.0) then
	653	! DO l=1, llm
	654	! DO ij=1,ip1jmp1
	655	! q(ij,l,igcm_n2)=q(ij,l,igcm_n2)+
	656	! & dqfi(ij,l,igcm_n2)*dtphys
	657	! ENDDO
	658	! ENDDO
	659	! ENDIF
	660
	661	c Mixing CO2 vertically
	662	c --------------------------
	663	if (igcm_n2.ne.0) then
	664	DO ij=1,ip1jmp1
	665	qmean_n2_vert=0.
	666	DO l=1, llm
	667	qmean_n2_vert= qmean_n2_vert
	668	& + q(ij,l,igcm_n2)*( p(ij,l) - p(ij,l+1))
	669	END DO
	670	qmean_n2_vert= qmean_n2_vert/ps(ij)
	671	DO l=1, llm
	672	q(ij,l,igcm_n2)= qmean_n2_vert
	673	END DO
	674	END DO
	675	end if
	676
	677
	678
	679	c Horizontal mixing of pressure
	680	c ------------------------------
	681	! Rappel newtonien vers psmean
	682	psmean= globaverage2d(ps) ! mean pressure
	683
	684	c Horizontal mixing tracers, and temperature (replace dynamics in nogcm)
	685	c ---------------------------------------------------------------------------------
	686
	687	! Simulate redistribution by dynamics for qn2
	688	if (igcm_n2.ne.0) then
	689
	690	DO ij=1,ip1jmp1
	691	pqn2(ij)= ps(ij) * q(ij,1,igcm_n2)
	692	ENDDO
	693	pqn2mean=globaverage2d(pqn2)
	694
	695	! Rappel newtonien vers qn2_mean
	696	qmean_n2= pqn2mean / psmean
	697
	698	DO ij=1,ip1jmp1
	699	q(ij,1,igcm_n2)=q(ij,1,igcm_n2)+
	700	& (qmean_n2-q(ij,1,igcm_n2))*
	701	& (1.-exp(-dtphys/tau_n2))
	702	ENDDO
	703
	704	DO l=2, llm
	705	DO ij=1,ip1jmp1
	706	q(ij,l,igcm_n2)=q(ij,1,igcm_n2)
	707	END DO
	708	END DO
	709
	710	! TEMPORAIRE (ED)
	711	! PRINT*,'psmean = ',psmean
	712	! PRINT*,'qmean_n2 = ',qmean_n2
	713	! PRINT*,'pqn2mean = ',pqn2mean
	714	! PRINT*,'q(50,1,igcm_n2) = ',q(50,1,igcm_n2)
	715	! PRINT*,'q(50,2,igcm_n2) = ',q(50,2,igcm_n2)
	716	! PRINT*,'q(50,3,igcm_n2) = ',q(50,3,igcm_n2)
	717
	718	endif ! igcm_n2.ne.0
	719
	720
	721	! ********************************************************
	722
	723	c Horizontal mixing of pressure
	724	c ------------------------------
	725	! Rappel newtonien vers psmean
	726	psmean= globaverage2d(ps) ! mean pressure
	727
	728
	729	c Mixing Temperature horizontally
	730	c -------------------------------
	731	! initialize variables that will be averaged
	732	DO l=1,llm
	733	DO ij=1,ip1jmp1
	734	dp(ij,l) = p(ij,l) - p(ij,l+1)
	735	tetadp(ij,l) = teta(ij,l)*dp(ij,l)
	736	ENDDO
	737	ENDDO
	738
	739	DO l=1,llm
	740	tetadpmean = globaverage2d(tetadp(:,l))
	741	dpmean = globaverage2d(dp(:,l))
	742	tetamean = tetadpmean / dpmean
	743	DO ij=1,ip1jmp1
	744	teta(ij,l) = teta(ij,l) + (tetamean - teta(ij,l)) *
	745	& (1 - exp(-dtphys/tau_teta))
	746	ENDDO
	747	ENDDO
	748
	749
	750	ENDIF ! of IF( apphys )
	751
	752
	753	c ********************************************************************
	754	c ********************************************************************
	755	c .... fin de l'integration dynamique et physique pour le pas itau ..
	756	c ********************************************************************
	757	c ********************************************************************
	758
	759
	760	IF ( .NOT.purmats ) THEN
	761	c ........................................................
	762	c .............. schema matsuno + leapfrog ..............
	763	c ........................................................
	764
	765	IF(forward. OR. leapf) THEN
	766	itau= itau + 1
	767	c iday= day_ini+itau/day_step
	768	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	769	c IF(time.GT.1.) THEN
	770	c time = time-1.
	771	c iday = iday+1
	772	c ENDIF
	773	ENDIF
	774
	775
	776	IF( itau. EQ. itaufinp1 ) then
	777	if (flag_verif) then
	778	write(79,*) 'ucov',ucov
	779	write(80,*) 'vcov',vcov
	780	write(81,*) 'teta',teta
	781	write(82,*) 'ps',ps
	782	write(83,*) 'q',q
	783	WRITE(85,*) 'q1 = ',q(:,:,1)
	784	WRITE(86,*) 'q3 = ',q(:,:,3)
	785	endif
	786
	787	abort_message = 'Simulation finished'
	788
	789	call abort_gcm(modname,abort_message,0)
	790	ENDIF
	791	c-----------------------------------------------------------------------
	792	c ecriture du fichier histoire moyenne:
	793	c -------------------------------------
	794
	795	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	796	IF(itau.EQ.itaufin) THEN
	797	iav=1
	798	ELSE
	799	iav=0
	800	ENDIF
	801
	802	! ! Ehouarn: re-compute geopotential for outputs
	803	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	804
	805	IF (ok_dynzon) THEN
	806	#ifdef CPP_IOIPSL
	807	c les traceurs ne sont pas sortis, trop lourd.
	808	c Peut changer eventuellement si besoin.
	809	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
	810	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
	811	& du,dudis,dutop,dufi)
	812	#endif
	813	END IF
	814	IF (ok_dyn_ave) THEN
	815	#ifdef CPP_IOIPSL
	816	CALL writedynav(itau,vcov,
	817	& ucov,teta,pk,phi,q,masse,ps,phis)
	818	#endif
	819	ENDIF
	820
	821	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
	822
	823	if (ok_iso_verif) then
	824	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
	825	endif !if (ok_iso_verif) then
	826
	827	c-----------------------------------------------------------------------
	828	c ecriture de la bande histoire:
	829	c ------------------------------
	830
	831	IF( MOD(itau,iecri).EQ.0) THEN
	832	! Ehouarn: output only during LF or Backward Matsuno
	833	if (leapf.or.(.not.leapf.and.(.not.forward))) then
	834	! ADAPTATION GCM POUR CP(T)
	835	call tpot2t(ijp1llm,teta,temp,pk)
	836	tsurpk = cpp*temp/pk
	837	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
	838	unat=0.
	839	do l=1,llm
	840	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	841	vnat(:,l)=vcov(:,l)/cv(:)
	842	enddo
	843	#ifdef CPP_IOIPSL
	844	if (ok_dyn_ins) then
	845	! write(lunout,*) "leapfrog: call writehist, itau=",itau
	846	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	847	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
	848	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
	849	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
	850	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
	851	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
	852	endif ! of if (ok_dyn_ins)
	853	#endif
	854	! For some Grads outputs of fields
	855	if (output_grads_dyn) then
	856	#include "write_grads_dyn.h"
	857	endif
	858	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
	859	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
	860
	861	IF(itau.EQ.itaufin) THEN
	862
	863	CALL dynredem1("restart.nc",start_time,
	864	& vcov,ucov,teta,q,masse,ps)
	865	CLOSE(99)
	866	!!! Ehouarn: Why not stop here and now?
	867	ENDIF ! of IF (itau.EQ.itaufin)
	868
	869	c-----------------------------------------------------------------------
	870	c gestion de l'integration temporelle:
	871	c ------------------------------------
	872
	873	IF( MOD(itau,iperiod).EQ.0 ) THEN
	874	GO TO 1
	875	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	876
	877	IF( forward ) THEN
	878	c fin du pas forward et debut du pas backward
	879
	880	forward = .FALSE.
	881	leapf = .FALSE.
	882	GO TO 2
	883
	884	ELSE
	885	c fin du pas backward et debut du premier pas leapfrog
	886
	887	leapf = .TRUE.
	888	dt = 2.*dtvr
	889	GO TO 2
	890	END IF ! of IF (forward)
	891	ELSE
	892
	893	c ...... pas leapfrog .....
	894
	895	leapf = .TRUE.
	896	dt = 2.*dtvr
	897	GO TO 2
	898	END IF ! of IF (MOD(itau,iperiod).EQ.0)
	899	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
	900
	901	ELSE ! of IF (.not.purmats)
	902
	903	if (ok_iso_verif) then
	904	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
	905	endif !if (ok_iso_verif) then
	906
	907	c ........................................................
	908	c .............. schema matsuno ...............
	909	c ........................................................
	910	IF( forward ) THEN
	911
	912	itau = itau + 1
	913	c iday = day_ini+itau/day_step
	914	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
	915	c
	916	c IF(time.GT.1.) THEN
	917	c time = time-1.
	918	c iday = iday+1
	919	c ENDIF
	920
	921	forward = .FALSE.
	922	IF( itau. EQ. itaufinp1 ) then
	923	abort_message = 'Simulation finished'
	924	call abort_gcm(modname,abort_message,0)
	925	ENDIF
	926	GO TO 2
	927
	928	ELSE ! of IF(forward) i.e. backward step
	929
	930	if (ok_iso_verif) then
	931	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
	932	endif !if (ok_iso_verif) then
	933
	934	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	935	IF(itau.EQ.itaufin) THEN
	936	iav=1
	937	ELSE
	938	iav=0
	939	ENDIF
	940
	941	! ! Ehouarn: re-compute geopotential for outputs
	942	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
	943
	944	IF (ok_dynzon) THEN
	945	#ifdef CPP_IOIPSL
	946	c les traceurs ne sont pas sortis, trop lourd.
	947	c Peut changer eventuellement si besoin.
	948	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
	949	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
	950	& du,dudis,dutop,dufi)
	951	#endif
	952	ENDIF
	953	IF (ok_dyn_ave) THEN
	954	#ifdef CPP_IOIPSL
	955	CALL writedynav(itau,vcov,
	956	& ucov,teta,pk,phi,q,masse,ps,phis)
	957	#endif
	958	ENDIF
	959
	960	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
	961
	962	IF(MOD(itau,iecri ).EQ.0) THEN
	963	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
	964	! ADAPTATION GCM POUR CP(T)
	965	call tpot2t(ijp1llm,teta,temp,pk)
	966	tsurpk = cpp*temp/pk
	967	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
	968	unat=0.
	969	do l=1,llm
	970	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
	971	vnat(:,l)=vcov(:,l)/cv(:)
	972	enddo
	973	#ifdef CPP_IOIPSL
	974	if (ok_dyn_ins) then
	975	! write(lunout,*) "leapfrog: call writehist (b)",
	976	! & itau,iecri
	977	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
	978	endif ! of if (ok_dyn_ins)
	979	#endif
	980	! For some Grads outputs
	981	if (output_grads_dyn) then
	982	#include "write_grads_dyn.h"
	983	endif
	984
	985	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
	986
	987	IF(itau.EQ.itaufin) THEN
	988	CALL dynredem1("restart.nc",start_time,
	989	& vcov,ucov,teta,q,masse,ps)
	990	ENDIF ! of IF(itau.EQ.itaufin)
	991
	992	forward = .TRUE.
	993	GO TO 1
	994
	995	ENDIF ! of IF (forward)
	996
	997	END IF ! of IF(.not.purmats)
	998
	999	STOP
	1000	END
	1001
	1002	! ************************************************************
	1003	! globalaverage VERSION PLuto
	1004	! ************************************************************
	1005	FUNCTION globaverage2d(var)
	1006	! ************************************************************
	1007	IMPLICIT NONE
	1008	#include "dimensions.h"
	1009	#include "paramet.h"
	1010	#include "comgeom2.h"
	1011	!ip1jmp1 called in paramet.h = iip1 x jjp1
	1012	REAL var(iip1,jjp1), globaverage2d
	1013	integer i,j
	1014	REAL airetot
	1015	save airetot
	1016	logical firstcall
	1017	data firstcall/.true./
	1018	save firstcall
	1019
	1020	if (firstcall) then
	1021	firstcall=.false.
	1022	airetot =0.
	1023	DO j=2,jjp1-1 ! lat
	1024	DO i = 1, iim ! lon
	1025	airetot = airetot + aire(i,j)
	1026	END DO
	1027	END DO
	1028	DO i=1,iim
	1029	airetot=airetot + aire(i,1)+aire(i,jjp1)
	1030	ENDDO
	1031	end if
	1032
	1033	globaverage2d = 0.
	1034	DO j=2,jjp1-1
	1035	DO i = 1, iim
	1036	globaverage2d = globaverage2d + var(i,j)*aire(i,j)
	1037	END DO
	1038	END DO
	1039
	1040	DO i=1,iim
	1041	globaverage2d = globaverage2d + var(i,1)*aire(i,1)
	1042	globaverage2d = globaverage2d + var(i,jjp1)*aire(i,jjp1)
	1043	ENDDO
	1044
	1045	globaverage2d = globaverage2d/airetot
	1046	return
	1047	end
	1048

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

Download in other formats:

Original Format