Context Navigation

leapfrog_nogcm.F @ 3303

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

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

Download in other formats:

Original Format