Context Navigation

bilan_dyn.F @ 125

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

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

Download in other formats:

Original Format