Context Navigation

bilan_dyn.F @ 9

Last change on this file since 9 was 6, checked in by slebonnois, 15 years ago

cf commit_v6.log :

manipulation traceurs
homogeneisation .def
bilan_dyn
etats initiaux start.nc
appels specifiques pour physique

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

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

Download in other formats:

Original Format