Context Navigation

bilan_dyn.F @ 3537

Last change on this file since 3537 was 1422, checked in by milmd, 10 years ago
In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).
File size: 27.7 KB

Rev	Line
[1]	1	!
	2	! $Id: bilan_dyn.F 1403 2010-07-01 09:02:53Z fairhead $
	3	!
[6]	4	SUBROUTINE bilan_dyn (dt_app,dt_cum,
	5	s ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,
	6	s ducovdyn,ducovdis,ducovspg,ducovphy)
	7	c si besoin des traceurs:
	8	c SUBROUTINE bilan_dyn (ntrac,dt_app,dt_cum,
	9	c s ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac,
	10	c s ducovdyn,ducovdis,ducovspg,ducovphy)
[1]	11
	12	c AFAIRE
	13	c Prevoir en champ nq+1 le diagnostique de l'energie
	14	c en faisant Qzon=Cv T + L * ...
	15	c vQ..A=Cp T + L * ...
	16
	17	#ifdef CPP_IOIPSL
	18	USE IOIPSL
	19	#endif
	20
[37]	21	USE control_mod, ONLY: planet_type
[1017]	22	USE cpdet_mod, only: tpot2t
[1422]	23	USE comvert_mod, ONLY: ap,bp,presnivs
	24	USE comconst_mod, ONLY: rad,omeg,pi
	25	USE temps_mod, ONLY: annee_ref,day_ref
[37]	26
[1]	27	IMPLICIT NONE
	28
	29	#include "dimensions.h"
	30	#include "paramet.h"
	31	#include "comgeom2.h"
	32	#include "iniprint.h"
	33
	34	c====================================================================
	35	c
	36	c Sous-programme consacre à des diagnostics dynamiques de base
	37	c
	38	c
	39	c De facon generale, les moyennes des scalaires Q sont ponderees par
	40	c la masse.
	41	c
	42	c Les flux de masse sont eux simplement moyennes.
	43	c
	44	c====================================================================
	45
	46	c Arguments :
	47	c ===========
	48
[6]	49	c integer ntrac
[1]	50	real dt_app,dt_cum
	51	real ps(iip1,jjp1)
	52	real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
	53	real flux_u(iip1,jjp1,llm)
	54	real flux_v(iip1,jjm,llm)
	55	real teta(iip1,jjp1,llm)
	56	real phi(iip1,jjp1,llm)
	57	real ucov(iip1,jjp1,llm)
	58	real vcov(iip1,jjm,llm)
[6]	59	c real trac(iip1,jjp1,llm,ntrac)
	60	c Tendances en m/s2 :
	61	real ducovdyn(iip1,jjp1,llm)
	62	real ducovdis(iip1,jjp1,llm)
	63	real ducovspg(iip1,jjp1,llm)
	64	real ducovphy(iip1,jjp1,llm)
[1]	65
	66	c Local :
	67	c =======
	68
	69	integer icum,ncum
[6]	70	save icum,ncum
[1]	71
[6]	72	integer i,j,l,iQ,num
	73	real zz,zqy,zfactv(jjm,llm),zfactw(jjm,llm)
	74	character*2 strd2
	75	real ww
[1]	76
[6]	77	logical first
	78	save first
	79	data first/.true./
[1]	80
	81	integer i_sortie
	82	save i_sortie
[6]	83	data i_sortie/1/
[1]	84
	85	real time
	86	integer itau
	87	save time,itau
	88	data time,itau/0.,0/
	89
[6]	90	! facteur = -1. pour Venus
	91	real fact_geovenus
	92	save fact_geovenus
[1]	93
	94	c variables dynamiques intermédiaires
[6]	95	c -----------------------------------
[1]	96	REAL vcont(iip1,jjm,llm),ucont(iip1,jjp1,llm)
	97	REAL ang(iip1,jjp1,llm),unat(iip1,jjp1,llm)
	98	REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
	99	REAL vorpot(iip1,jjm,llm)
	100	REAL w(iip1,jjp1,llm),ecin(iip1,jjp1,llm),convm(iip1,jjp1,llm)
[6]	101	real temp(iip1,jjp1,llm)
	102	real dudyn(iip1,jjp1,llm)
	103	real dudis(iip1,jjp1,llm)
	104	real duspg(iip1,jjp1,llm)
	105	real duphy(iip1,jjp1,llm)
[1]	106
[6]	107	c CHAMPS SCALAIRES Q ADVECTES
	108	c ----------------------------
	109	integer nQ
	110	c avec tous les composes, ca fait trop.... Je les enleve
	111	c parameter (nQ=6+nqmx)
	112	parameter (nQ=6)
	113
	114	character*6,save :: nom(nQ)
	115	character*6,save :: unites(nQ)
	116
	117	integer itemp,igeop,iecin,iang,iu,iun
	118	save itemp,igeop,iecin,iang,iu,iun
	119	data itemp,igeop,iecin,iang,iu,iun/1,2,3,4,5,6/
	120
[1]	121	c champ contenant les scalaires advectés.
	122	real Q(iip1,jjp1,llm,nQ)
	123
	124	c champs cumulés
	125	real ps_cum(iip1,jjp1)
	126	real masse_cum(iip1,jjp1,llm)
	127	real flux_u_cum(iip1,jjp1,llm)
	128	real flux_v_cum(iip1,jjm,llm)
[6]	129	real flux_w_cum(iip1,jjp1,llm)
[1]	130	real Q_cum(iip1,jjp1,llm,nQ)
	131	real flux_uQ_cum(iip1,jjp1,llm,nQ)
	132	real flux_vQ_cum(iip1,jjm,llm,nQ)
	133	real flux_wQ_cum(iip1,jjp1,llm,nQ)
	134	real dQ(iip1,jjp1,llm,nQ)
	135
	136	save ps_cum,masse_cum,flux_u_cum,flux_v_cum
	137	save Q_cum,flux_uQ_cum,flux_vQ_cum
[6]	138	save flux_w_cum,flux_wQ_cum
[1]	139
[6]	140	c champs de transport en moyenne zonale
[1]	141	integer ntr,itr
	142	parameter (ntr=5)
	143
	144	character*10,save :: znom(ntr,nQ)
	145	character*20,save :: znoml(ntr,nQ)
	146	character*10,save :: zunites(ntr,nQ)
[6]	147	character*10,save :: znom2(ntr,nQ)
	148	character*20,save :: znom2l(ntr,nQ)
	149	character*10,save :: zunites2(ntr,nQ)
[37]	150	character*10,save :: znom3(nQ)
	151	character*20,save :: znom3l(nQ)
	152	character*10,save :: zunites3(nQ)
[1]	153
	154	integer iave,itot,immc,itrs,istn
	155	data iave,itot,immc,itrs,istn/1,2,3,4,5/
	156	character*3 ctrs(ntr)
	157	data ctrs/' ','TOT','MMC','TRS','STN'/
	158
	159	real zvQ(jjm,llm,ntr,nQ),zvQtmp(jjm,llm)
[6]	160	real zwQ(jjm,llm,ntr,nQ),zwQtmp(jjm,llm)
[1]	161	real zavQ(jjm,ntr,nQ),psiQ(jjm,llm+1,nQ)
[6]	162	real zawQ(jjm,llm,ntr,nQ)
	163	real zdQ(jjm,llm,nQ)
	164	real zmasse(jjm,llm),zavmasse(jjm),zawmasse(llm)
	165	real psbar(jjm)
[1]	166
[6]	167	real zv(jjm,llm),zw(jjp1,llm),psi(jjm,llm+1)
[1]	168
[6]	169	c TENDANCES POUR MOMENT CINETIQUE
	170	c -------------------------------
[1]	171
[6]	172	integer ntdc,itdc
	173	parameter (ntdc=4)
[1]	174
[6]	175	integer itdcdyn,itdcdis,itdcspg,itdcphy
	176	data itdcdyn,itdcdis,itdcspg,itdcphy/1,2,3,4/
	177
	178	character*6,save :: nomtdc(ntdc)
	179
	180	c champ contenant les tendances du moment cinetique.
	181	real tdc(iip1,jjp1,llm,ntdc)
	182	real ztdc(jjm,llm,ntdc) ! moyenne zonale
	183
	184	c champs cumulés
	185	real tdc_cum(iip1,jjp1,llm,ntdc)
	186	save tdc_cum
	187
	188	c integrations completes
	189	real mtot,mctot,dmctot(ntdc)
	190
[1]	191	c Initialisation du fichier contenant les moyennes zonales.
	192	c ---------------------------------------------------------
	193
	194	character*10 infile
	195
	196	integer fileid
	197	integer thoriid, zvertiid
	198	save fileid
	199
	200	integer ndex3d(jjm*llm)
[6]	201	real ztmp3d(jjm,llm)
[1]	202
	203	C Variables locales
	204	C
	205	integer tau0
	206	real zjulian
	207	integer zan, dayref
	208	C
	209	real rlong(jjm),rlatg(jjm)
	210
	211
	212
	213	c=====================================================================
	214	c Initialisation
	215	c=====================================================================
	216
	217	ndex3d=0
	218
	219	if (first) then
	220
[6]	221	if (planet_type.eq."venus") then
	222	fact_geovenus = -1.
	223	else
	224	fact_geovenus = 1.
	225	endif
[1]	226
	227	icum=0
	228	c initialisation des fichiers
	229	first=.false.
	230	c ncum est la frequence de stokage en pas de temps
	231	ncum=dt_cum/dt_app
[6]	232	if (abs(ncumdt_app-dt_cum).gt.1.e-2dt_app) then
	233	if (abs((ncum+1)dt_app-dt_cum).lt.1.e-2dt_app) then
	234	ncum = ncum+1
	235	elseif (abs((ncum-1)dt_app-dt_cum).lt.1.e-2dt_app) then
	236	ncum = ncum-1
	237	else
[1]	238	WRITE(lunout,*)
	239	. 'Pb : le pas de cumule doit etre multiple du pas'
	240	WRITE(lunout,*)'dt_app=',dt_app
	241	WRITE(lunout,*)'dt_cum=',dt_cum
[6]	242	WRITE(lunout,)'ncumdt_app=',ncum*dt_app
	243	WRITE(lunout,*)'ncum=',ncum
[1]	244	stop
[6]	245	endif
[1]	246	endif
	247
[6]	248	c VARIABLES ADVECTEES:
[1]	249
[6]	250	nom(itemp)='temp'
[1]	251	nom(igeop)='gz'
[6]	252	nom(iecin)='ecin'
[1]	253	nom(iang)='ang'
	254	nom(iu)='u'
	255	nom(iun)='un'
	256
	257	unites(itemp)='K'
	258	unites(igeop)='m2/s2'
	259	unites(iecin)='m2/s2'
	260	unites(iang)='ang'
	261	unites(iu)='m/s'
	262	unites(iun)='un'
	263
[6]	264	c avec tous les composes, ca fait trop.... Je les enleve
	265	c do num=1,ntrac
	266	c write(strd2,'(i2.2)') num
	267	c nom(6+num)='trac'//strd2
	268	c unites(6+num)='kg/kg'
	269	c enddo
[1]	270
[6]	271	c TENDANCES MOMENT CIN:
	272
	273	nomtdc(itdcdyn) ='dmcdyn'
	274	nomtdc(itdcdis) ='dmcdis'
	275	nomtdc(itdcspg) ='dmcspg'
	276	nomtdc(itdcphy) ='dmcphy'
	277
[1]	278	c Initialisation du fichier contenant les moyennes zonales.
	279	c ---------------------------------------------------------
	280
	281	infile='dynzon'
	282
	283	zan = annee_ref
	284	dayref = day_ref
	285	CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
[6]	286	c tau0 = itau_dyn
	287	tau0 = 0
	288	itau = tau0
[1]	289
	290	rlong=0.
[6]	291	rlatg=rlatv180./pifact_geovenus
[1]	292
	293	call histbeg(infile, 1, rlong, jjm, rlatg,
	294	. 1, 1, 1, jjm,
	295	. tau0, zjulian, dt_cum, thoriid, fileid)
	296
	297	C
	298	C Appel a histvert pour la grille verticale
	299	C
	300	call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
	301	. llm, presnivs, zvertiid)
	302	C
	303	C Appels a histdef pour la definition des variables a sauvegarder
[6]	304
[1]	305	do iQ=1,nQ
	306	do itr=1,ntr
	307	if(itr.eq.1) then
[6]	308	znom(itr,iQ) =nom(iQ)
	309	znoml(itr,iQ) =nom(iQ)
	310	zunites(itr,iQ) =unites(iQ)
[1]	311	else
[6]	312	znom(itr,iQ) =ctrs(itr)//'v'//nom(iQ)
	313	znoml(itr,iQ) ='transport : v * '//nom(iQ)//' '//ctrs(itr)
	314	zunites(itr,iQ) ='m/s * '//unites(iQ)
	315	znom2(itr,iQ) =ctrs(itr)//'w'//nom(iQ)
	316	znom2l(itr,iQ) ='transport: w * '//nom(iQ)//' '//ctrs(itr)
	317	zunites2(itr,iQ)='Pa/s * '//unites(iQ)
[1]	318	endif
	319	enddo
[6]	320	znom3(iQ)='d'//nom(iQ)
	321	znom3l(iQ)='convergence: '//nom(iQ)
	322	zunites3(iQ)=unites(iQ)//' /s'
	323	c print*,'DEBUG:',znom3(iQ),znom3l(iQ),zunites3(iQ)
[1]	324	enddo
	325
	326	c Declarations des champs avec dimension verticale
[6]	327
	328	if (1.eq.0) then ! on les sort, ou pas...
	329
	330	c do iQ=1,nQ
	331	c !!!! JE NE SORS ICI QUE temp et ang POUR CAUSE DE PLACE !
	332	do iQ=1,4,3
[1]	333	do itr=1,ntr
	334	IF (prt_level > 5)
	335	. WRITE(lunout,*)'var ',itr,iQ
	336	. ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
	337	call histdef(fileid,znom(itr,iQ),znoml(itr,iQ),
	338	. zunites(itr,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
	339	. 32,'ave(X)',dt_cum,dt_cum)
	340	enddo
[6]	341	c transport vertical:
	342	do itr=2,ntr
	343	IF (prt_level > 5)
	344	. WRITE(lunout,*)'var ',itr,iQ
	345	. ,znom2(itr,iQ),znom2l(itr,iQ),zunites2(itr,iQ)
	346	call histdef(fileid,znom2(itr,iQ),znom2l(itr,iQ),
	347	. zunites2(itr,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
	348	. 32,'ins(X)',dt_cum,dt_cum)
	349	enddo
	350
	351	c Declarations pour convergences
	352	IF (prt_level > 5)
	353	. WRITE(lunout,*)'var ',iQ
	354	. ,znom3(iQ),znom3l(iQ),zunites3(iQ)
	355	call histdef(fileid,znom3(iQ),znom3l(iQ),
	356	. zunites3(iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
	357	. 32,'ins(X)',dt_cum,dt_cum)
	358
[1]	359	c Declarations pour les fonctions de courant
[6]	360	c Non sorties ici...
	361	c call histdef(fileid,'psi'//nom(iQ)
	362	c . ,'stream fn. '//znoml(itot,iQ),
	363	c . zunites(itot,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
	364	c . 32,'ave(X)',dt_cum,dt_cum)
[1]	365
[6]	366	enddo ! iQ
[1]	367
[6]	368	endif ! 1=1 sortie ou non...
	369
[1]	370	c Declarations pour les champs de transport d'air
	371	call histdef(fileid, 'masse', 'masse',
	372	. 'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid,
	373	. 32, 'ave(X)', dt_cum, dt_cum)
	374	call histdef(fileid, 'v', 'v',
	375	. 'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
	376	. 32, 'ave(X)', dt_cum, dt_cum)
[6]	377	call histdef(fileid, 'w', 'w',
	378	. 'Pa/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
	379	. 32, 'ins(X)', dt_cum, dt_cum)
	380
	381	c Declarations pour la fonction de courant
[1]	382	call histdef(fileid,'psi','stream fn. MMC ','mega t/s',
	383	. 1,jjm,thoriid,llm,1,llm,zvertiid,
	384	. 32,'ave(X)',dt_cum,dt_cum)
	385
	386
[6]	387	c Declarations pour les tendances de moment cinetique
	388	do itdc=1,ntdc
	389	call histdef(fileid, nomtdc(itdc), nomtdc(itdc),
	390	. 'ang/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
	391	. 32, 'ins(X)', dt_cum, dt_cum)
	392	enddo
	393
[1]	394	c Declaration des champs 1D de transport en latitude
[6]	395	c do iQ=1,nQ
	396	c !!!! JE NE SORS ICI QUE temp et ang POUR CAUSE DE PLACE !
	397	do iQ=1,4,3
[1]	398	do itr=2,ntr
	399	call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ),
	400	. zunites(itr,iQ),1,jjm,thoriid,1,1,1,-99,
	401	. 32,'ave(X)',dt_cum,dt_cum)
[6]	402	c JE VIRE LE VERTICAL POUR L'INSTANT
	403	c call histdef(fileid,'a'//znom2(itr,iQ),znom2l(itr,iQ),
	404	c . zunites2(itr,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
	405	c . 32,'ins(X)',dt_cum,dt_cum)
[1]	406	enddo
	407	enddo
	408
	409	CALL histend(fileid)
	410
	411
[6]	412	endif ! first
[1]	413
	414
	415	c=====================================================================
	416	c Calcul des champs dynamiques
	417	c ----------------------------
	418
	419	c énergie cinétique
	420	ucont(:,:,:)=0
	421	CALL covcont(llm,ucov,vcov,ucont,vcont)
	422	CALL enercin(vcov,ucov,vcont,ucont,ecin)
	423
[6]	424	c moment cinétique et tendances
	425	dudyn = 0.
	426	dudis = 0.
	427	duspg = 0.
	428	duphy = 0.
[1]	429	do l=1,llm
	430	unat(:,:,l)=ucont(:,:,l)*cu(:,:)
[6]	431	dudyn(:,2:jjm,l) = ducovdyn(:,2:jjm,l)/cu(:,2:jjm)
	432	dudis(:,2:jjm,l) = ducovdis(:,2:jjm,l)/cu(:,2:jjm)
	433	duspg(:,2:jjm,l) = ducovspg(:,2:jjm,l)/cu(:,2:jjm)
	434	duphy(:,2:jjm,l) = ducovphy(:,2:jjm,l)/cu(:,2:jjm)
	435	do j=1,jjp1
	436	ang(:,j,l)= radcos(rlatu(j))
	437	. ( unat(:,j,l) + radcos(rlatu(j))omeg )
	438	tdc(:,j,l,1) = radcos(rlatu(j))dudyn(:,j,l)
	439	tdc(:,j,l,2) = radcos(rlatu(j))dudis(:,j,l)
	440	tdc(:,j,l,3) = radcos(rlatu(j))duspg(:,j,l)
	441	tdc(:,j,l,4) = radcos(rlatu(j))duphy(:,j,l)
	442	enddo
[1]	443	enddo
[6]	444	c Normalisation:
	445	ang = ang / (2./3. radrad*omeg)
	446	do itdc=1,ntdc
	447	tdc(:,:,:,itdc)=tdc(:,:,:,itdc) / (2./3. radrad*omeg)
	448	enddo
[1]	449
[6]	450	! ADAPTATION GCM POUR CP(T)
	451	call tpot2t(ip1jmp1*llm,teta,temp,pk)
	452	Q(:,:,:,itemp) = temp(:,:,:)
	453	Q(:,:,:,igeop) =phi(:,:,:)
	454	Q(:,:,:,iecin) =ecin(:,:,:)
	455	Q(:,:,:,iang) =ang(:,:,:)
	456	Q(:,:,:,iu) =unat(:,:,:)
	457	Q(:,:,:,iun) =1.
	458	c avec tous les composes, ca fait trop.... Je les enleve
	459	c do num=1,ntrac
	460	c Q(:,:,:,6+num)=trac(:,:,:,num)
	461	c enddo
[1]	462
[6]	463	c calcul du flux de masse vertical (+ vers le bas)
	464	call convmas(flux_u,flux_v,convm)
	465	CALL vitvert(convm,w)
[1]	466
	467	c=====================================================================
	468	c Cumul
	469	c=====================================================================
	470	c
	471	if(icum.EQ.0) then
[6]	472	ps_cum = 0.
	473	masse_cum = 0.
	474	flux_u_cum = 0.
	475	flux_v_cum = 0.
	476	flux_w_cum = 0.
	477	Q_cum = 0.
	478	flux_vQ_cum = 0.
	479	flux_uQ_cum = 0.
	480	flux_wQ_cum = 0.
	481	tdc_cum = 0.
[1]	482	endif
	483
	484	IF (prt_level > 5)
	485	. WRITE(lunout,*)'dans bilan_dyn ',icum,'->',icum+1
	486	icum=icum+1
	487
	488	c accumulation des flux de masse horizontaux
[6]	489	ps_cum = ps_cum + ps
	490	masse_cum = masse_cum + masse
	491	flux_u_cum = flux_u_cum + flux_u
	492	flux_v_cum = flux_v_cum + flux_v
	493	flux_w_cum = flux_w_cum + w
[1]	494	do iQ=1,nQ
[6]	495	Q_cum(:,:,:,iQ) = Q_cum(:,:,:,iQ) + Q(:,:,:,iQ)*masse(:,:,:)
[1]	496	enddo
[6]	497	do itdc=1,ntdc
	498	tdc_cum(:,:,:,itdc) =
	499	. tdc_cum(:,:,:,itdc) + tdc(:,:,:,itdc)*masse(:,:,:)
	500	enddo
[1]	501
	502	c=====================================================================
	503	c FLUX ET TENDANCES
	504	c=====================================================================
	505
	506	c Flux longitudinal
	507	c -----------------
	508	do iQ=1,nQ
	509	do l=1,llm
	510	do j=1,jjp1
	511	do i=1,iim
	512	flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ)
	513	s +flux_u(i,j,l)0.5(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
	514	enddo
	515	flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
	516	enddo
	517	enddo
	518	enddo
	519
	520	c flux méridien
	521	c -------------
	522	do iQ=1,nQ
	523	do l=1,llm
	524	do j=1,jjm
	525	do i=1,iip1
	526	flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ)
	527	s +flux_v(i,j,l)0.5(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
	528	enddo
	529	enddo
	530	enddo
	531	enddo
	532
[6]	533	c Flux vertical
	534	c -------------
	535	do iQ=1,nQ
	536	do l=2,llm
	537	do j=1,jjp1
	538	do i=1,iip1
	539	flux_wQ_cum(i,j,l,iQ)=flux_wQ_cum(i,j,l,iQ)
	540	s +w(i,j,l)0.5(Q(i,j,l-1,iQ)+Q(i,j,l,iQ))
	541	enddo
	542	enddo
	543	enddo
	544	flux_wQ_cum(:,:,1,iQ)=0.0e0
	545	enddo
[1]	546
	547	c tendances
	548	c ---------
	549
	550	c convergence horizontale
	551	call convflu(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
	552
	553	c calcul de la vitesse verticale
	554	call convmas(flux_u_cum,flux_v_cum,convm)
	555	CALL vitvert(convm,w)
	556
[6]	557	c ajustement tendances (vitesse verticale)
[1]	558	do iQ=1,nQ
	559	do l=1,llm-1
	560	do j=1,jjp1
	561	do i=1,iip1
	562	ww=-0.5w(i,j,l+1)(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
	563	dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww
	564	dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
	565	enddo
	566	enddo
	567	enddo
	568	enddo
	569	IF (prt_level > 5)
	570	. WRITE(lunout,*)'Apres les calculs fait a chaque pas'
[6]	571
[1]	572	c=====================================================================
	573	c PAS DE TEMPS D'ECRITURE
	574	c=====================================================================
	575	if (icum.eq.ncum) then
	576	c=====================================================================
	577
[6]	578	time=time+dt_cum
	579	itau=itau+1
	580
[1]	581	IF (prt_level > 5)
	582	. WRITE(lunout,*)'Pas d ecriture'
	583
	584	c Normalisation
	585	do iQ=1,nQ
[6]	586	Q_cum(:,:,:,iQ) = Q_cum(:,:,:,iQ)/masse_cum(:,:,:)
	587	dQ(:,:,:,iQ) = dQ(:,:,:,iQ) /masse_cum(:,:,:)
[1]	588	enddo
[6]	589	do itdc=1,ntdc
	590	tdc_cum(:,:,:,itdc) = tdc_cum(:,:,:,itdc)/masse_cum(:,:,:)
	591	enddo
	592
[1]	593	zz=1./REAL(ncum)
[6]	594	ps_cum = ps_cum *zz
	595	masse_cum = masse_cum *zz
	596	flux_u_cum = flux_u_cum *zz
	597	flux_v_cum = flux_v_cum *zz
	598	flux_w_cum = flux_w_cum *zz
	599	flux_uQ_cum = flux_uQ_cum *zz
	600	flux_vQ_cum = flux_vQ_cum *zz
	601	flux_wQ_cum = flux_wQ_cum *zz
[1]	602
[6]	603	c Integration complete
	604	mtot = 0.
	605	mctot = 0.
	606	dmctot = 0.
	607	do l=1,llm
	608	do j=2,jjm
	609	do i=1,iim
	610	mtot = mtot + masse_cum(i,j,l)
	611	mctot = mctot + Q_cum(i,j,l,iang)*masse_cum(i,j,l)
	612	enddo
	613	enddo
	614	enddo
	615	mctot = mctot/mtot
	616	do itdc=1,ntdc
	617	do l=1,llm
	618	do j=2,jjm
	619	do i=1,iim
	620	dmctot(itdc) = dmctot(itdc)
	621	. + tdc_cum(i,j,l,itdc)*masse_cum(i,j,l)/mtot
	622	enddo
	623	enddo
	624	enddo
	625	enddo
[1]	626
	627	c=====================================================================
	628	c Transport méridien
	629	c=====================================================================
	630
	631	c cumul zonal des masses des mailles
	632	c ----------------------------------
	633	zv=0.
[6]	634	zw=0.
[1]	635	zmasse=0.
	636	call massbar(masse_cum,massebx,masseby)
[6]	637
	638	c moy zonale de la ps cumulee
	639	do j=1,jjm
	640	psbar(j)=0.
	641	do i=1,iim
	642	psbar(j)=psbar(j)+ps_cum(i,j)/iim
	643	enddo
	644	enddo
	645
[1]	646	do l=1,llm
	647	do j=1,jjm
	648	do i=1,iim
	649	zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
	650	zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
[6]	651	zw(j,l)=zw(j,l)+flux_w_cum(i,j,l)
[1]	652	enddo
	653	zfactv(j,l)=cv(1,j)/zmasse(j,l)
[6]	654	zfactw(j,l)=((ap(l)-ap(l+1))+psbar(j)*(bp(l)-bp(l+1)))
	655	. /zmasse(j,l)
[1]	656	enddo
[6]	657	do i=1,iim
	658	zw(jjp1,l)=zw(jjp1,l)+flux_w_cum(i,jjp1,l)
	659	enddo
[1]	660	enddo
	661
	662	c print*,'3OK'
	663	c --------------------------------------------------------------
	664	c calcul de la moyenne zonale du transport :
	665	c ------------------------------------------
	666	c
	667	c --
	668	c TOT : la circulation totale [ vq ]
	669	c
	670	c - -
	671	c MMC : mean meridional circulation [ v ] [ q ]
	672	c
	673	c ---- -- - -
	674	c TRS : transitoires [ v'q'] = [ vq ] - [ v q ]
	675	c
	676	c - * - * - - - -
	677	c STT : stationaires [ v q ] = [ v q ] - [ v ] [ q ]
	678	c
	679	c - -
	680	c on utilise aussi l'intermediaire TMP : [ v q ]
	681	c
	682	c la variable zfactv transforme un transport meridien cumule
	683	c en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
[6]	684	c la variable zfactw transforme un transport vertical cumule
	685	c en kg/s * unte-du-champ-transporte en Pa/s * unite-du-champ-transporte
[1]	686	c
	687	c --------------------------------------------------------------
	688
	689
	690	c ----------------------------------------
	691	c Transport dans le plan latitude-altitude
	692	c ----------------------------------------
	693
	694	zvQ=0.
[6]	695	zwQ=0.
	696	zdQ=0.
[1]	697	psiQ=0.
[6]	698
[1]	699	do iQ=1,nQ
[6]	700
	701	c transport meridien
[1]	702	zvQtmp=0.
	703	do l=1,llm
	704	do j=1,jjm
	705	c print*,'j,l,iQ=',j,l,iQ
[6]	706	c Calcul des moyennes zonales du transport total et de zvQtmp
[1]	707	do i=1,iim
	708	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
	709	s +flux_vQ_cum(i,j,l,iQ)
	710	zqy= 0.5(Q_cum(i,j,l,iQ)masse_cum(i,j,l)+
	711	s Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
	712	zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
	713	s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
	714	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
	715	enddo
	716	c print*,'aOK'
	717	c Decomposition
	718	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
	719	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
	720	zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
	721	zvQ(j,l,immc,iQ)=zv(j,l)zvQ(j,l,iave,iQ)zfactv(j,l)
	722	zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
	723	zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
	724	enddo
	725	enddo
	726	c fonction de courant meridienne pour la quantite Q
	727	do l=llm,1,-1
	728	do j=1,jjm
[6]	729	psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)/zfactv(j,l)
[1]	730	enddo
	731	enddo
[37]	732	!! enddo
[1]	733
[6]	734	c transport vertical
	735	zwQtmp=0.
	736	do l=1,llm
	737	do j=1,jjm
	738	c print*,'j,l,iQ=',j,l,iQ
	739	c Calcul des moyennes zonales du transport vertical total et de zwQtmp
	740	do i=1,iim
	741	zwQ(j,l,itot,iQ)=zwQ(j,l,itot,iQ)
	742	s +flux_wQ_cum(i,j,l,iQ)
	743	zqy= 0.5(Q_cum(i,j,l,iQ)masse_cum(i,j,l)+
	744	s Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
	745	zwQtmp(j,l)=zwQtmp(j,l)+flux_w_cum(i,j,l)*zqy
	746	s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
	747	zwQ(j,l,iave,iQ)=zwQ(j,l,iave,iQ)+zqy
	748	enddo
	749	c Decomposition
	750	zwQ(j,l,iave,iQ)=zwQ(j,l,iave,iQ)/zmasse(j,l)
	751	zwQ(j,l,itot,iQ)=zwQ(j,l,itot,iQ)*zfactw(j,l)
	752	zwQtmp(j,l)=zwQtmp(j,l)*zfactw(j,l)
	753	zwQ(j,l,immc,iQ)=zw(j,l)zwQ(j,l,iave,iQ)zfactw(j,l)
	754	zwQ(j,l,itrs,iQ)=zwQ(j,l,itot,iQ)-zwQtmp(j,l)
	755	zwQ(j,l,istn,iQ)=zwQtmp(j,l)-zwQ(j,l,immc,iQ)
	756	enddo
	757	enddo
	758
	759	c convergence
	760	c Calcul moyenne zonale de la convergence totale
	761	do l=1,llm
	762	do j=1,jjm
	763	c print*,'j,l,iQ=',j,l,iQ
	764	do i=1,iim
	765	zdQ(j,l,iQ)=zdQ(j,l,iQ) +
	766	. ( dQ(i,j,l,iQ) * masse_cum(i,j,l)
	767	. + dQ(i,j+1,l,iQ) * masse_cum(i,j+1,l))
	768	. / ( masse_cum(i,j,l)+masse_cum(i,j+1,l))
	769	enddo
	770	enddo
	771	enddo
[37]	772	enddo ! of do iQ=1,nQ
[6]	773
[1]	774	c fonction de courant pour la circulation meridienne moyenne
	775	psi=0.
	776	do l=llm,1,-1
	777	do j=1,jjm
[6]	778	psi(j,l)= psi(j,l+1)+zv(j,l)
	779	zv(j,l) = zv(j,l)*zfactv(j,l)
	780	zw(j,l) = 0.5(zw(j,l)+zw(j+1,l))zfactw(j,l)
[1]	781	enddo
	782	enddo
	783
[6]	784	c Calcul moyenne zonale des tendances moment cin.
	785	ztdc=0.
	786	do itdc=1,ntdc
	787	do l=1,llm
	788	do j=1,jjm
	789	do i=1,iim
	790	ztdc(j,l,itdc)=ztdc(j,l,itdc) +
	791	. ( tdc_cum(i,j,l,itdc) * masse_cum(i,j,l)
	792	. + tdc_cum(i,j+1,l,itdc) * masse_cum(i,j+1,l))
	793	. / ( masse_cum(i,j,l)+masse_cum(i,j+1,l))
	794	enddo
	795	enddo
	796	enddo
	797	enddo
	798
[1]	799	c print*,'4OK'
[6]	800
	801	c--------------------------------------
	802	c--------------------------------------
[1]	803	c sorties proprement dites
[6]	804	c--------------------------------------
	805	c--------------------------------------
	806
[1]	807	if (i_sortie.eq.1) then
[6]	808
	809	c sortie des integrations completes dans le listing
	810	write(*,'(A12,5(1PE11.4,X))') "BILANMCDYN ",mctot,dmctot
	811
	812	c sorties dans fichier dynzon
	813
	814	if (1.eq.0) then ! on les sort, ou pas...
	815
	816	c avec tous les composes, ca fait trop.... Je les enleve
	817	c do iQ=1,nQ
	818	c do iQ=1,6
	819	c !!!! JE NE SORS ICI QUE temp et ang POUR CAUSE DE PLACE !
	820	do iQ=1,4,3
	821
	822	ztmp3d(:,:)= zvQ(:,:,1,iQ) ! valeur moyenne
	823	call histwrite(fileid,znom(1,iQ),itau,ztmp3d
[1]	824	s ,jjm*llm,ndex3d)
[6]	825	do itr=2,ntr
	826	ztmp3d(:,:)= zvQ(:,:,itr,iQ)*fact_geovenus ! transport horizontal
	827	call histwrite(fileid,znom(itr,iQ),itau,ztmp3d
[1]	828	s ,jjm*llm,ndex3d)
[6]	829	enddo
	830
	831	do itr=2,ntr
	832	ztmp3d(:,:)=zwQ(:,:,itr,iQ)
	833	call histwrite(fileid,znom2(itr,iQ),itau,ztmp3d
	834	s ,jjm*llm,ndex3d)
	835	enddo
	836
	837	ztmp3d(:,:)= zdQ(:,:,iQ)
	838	call histwrite(fileid,znom3(iQ),itau,ztmp3d
	839	s ,jjm*llm,ndex3d)
	840
	841	c ztmp3d(:,:)= psiQ(:,1:llm,iQ)*fact_geovenus
	842	c call histwrite(fileid,'psi'//nom(iQ),itau,ztmp3d
	843	c s ,jjm*llm,ndex3d)
[1]	844	enddo
	845
[6]	846	endif ! 1=1 sortie ou non...
	847
	848	ztmp3d=zmasse
	849	call histwrite(fileid,'masse',itau,ztmp3d
[1]	850	s ,jjm*llm,ndex3d)
[6]	851
	852	ztmp3d= zv*fact_geovenus
	853	call histwrite(fileid,'v',itau,ztmp3d
[1]	854	s ,jjm*llm,ndex3d)
[6]	855	ztmp3d(:,:)=zw(1:jjm,:)
	856	call histwrite(fileid,'w',itau,ztmp3d
	857	s ,jjm*llm,ndex3d)
	858	ztmp3d= psi(:,1:llm)1.e-9fact_geovenus
	859	call histwrite(fileid,'psi',itau,ztmp3d,jjm*llm,ndex3d)
[1]	860
[6]	861	do itdc=1,ntdc
	862	ztmp3d(:,:)= ztdc(:,:,itdc)
	863	call histwrite(fileid,nomtdc(itdc),itau,ztmp3d
	864	s ,jjm*llm,ndex3d)
	865	enddo
[1]	866
[6]	867	endif ! i_sortie
[1]	868
[6]	869
[1]	870	c -----------------
	871	c Moyenne verticale
	872	c -----------------
	873
[6]	874	zavmasse=0.
[1]	875	do l=1,llm
[6]	876	zavmasse(:)=zavmasse(:)+zmasse(:,l)
[1]	877	enddo
	878	zavQ=0.
[6]	879
	880	c avec tous les composes, ca fait trop.... Je les enleve
	881	c do iQ=1,nQ
	882	c do iQ=1,6
	883	c !!!! JE NE SORS ICI QUE temp et ang POUR CAUSE DE PLACE !
	884	do iQ=1,4,3
[1]	885	do itr=2,ntr
	886	do l=1,llm
	887	zavQ(:,itr,iQ)=zavQ(:,itr,iQ)+zvQ(:,l,itr,iQ)*zmasse(:,l)
	888	enddo
[6]	889	zavQ(:,itr,iQ)=zavQ(:,itr,iQ)/zavmasse(:)
	890	if (i_sortie.eq.1) then
	891	ztmp3d=0.0
	892	ztmp3d(:,1)= zavQ(:,itr,iQ)*fact_geovenus
	893	call histwrite(fileid,'a'//znom(itr,iQ),itau,ztmp3d
	894	. ,jjm*llm,ndex3d)
	895	endif
[1]	896	enddo
	897	enddo
	898
[6]	899	c ------------------
	900	c Moyenne meridienne
	901	c ------------------
[1]	902
[6]	903	zawmasse=0.
	904	do j=1,jjm
	905	do l=1,llm
	906	zawmasse(l)=zawmasse(l)+zmasse(j,l)
	907	enddo
	908	enddo
	909	zawQ=0.
	910
	911	c avec tous les composes, ca fait trop.... Je les enleve
	912	c do iQ=1,nQ
	913	c do iQ=1,6
	914	c !!!! JE NE SORS ICI QUE temp et ang POUR CAUSE DE PLACE !
	915	do iQ=1,4,3
	916	do itr=2,ntr
	917	do l=1,llm
[97]	918	do j=1,jjm
[6]	919	zawQ(1,l,itr,iQ)=zawQ(1,l,itr,iQ)+zwQ(j,l,itr,iQ)*zmasse(j,l)
	920	enddo
	921	zawQ(1,l,itr,iQ)=zawQ(1,l,itr,iQ)/zawmasse(l)
	922	enddo
	923	if (i_sortie.eq.1) then
	924	ztmp3d=0.0
	925	do l=1,llm
	926	ztmp3d(1,l)=zawQ(1,l,itr,iQ)
	927	enddo
	928	c JE VIRE LE VERTICAL POUR L'INSTANT
	929	c call histwrite(fileid,'a'//znom2(itr,iQ),itau,ztmp3d
	930	c . ,jjm*llm,ndex3d)
	931	endif
	932	enddo
	933	enddo
	934
	935	call histsync(fileid)
	936
[1]	937	c=====================================================================
	938	c/////////////////////////////////////////////////////////////////////
	939	icum=0 !///////////////////////////////////////
	940	endif ! icum.eq.ncum !///////////////////////////////////////
	941	c/////////////////////////////////////////////////////////////////////
	942	c=====================================================================
	943
	944	return
	945	end

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

Download in other formats:

Original Format