Context Navigation

bilan_dyn.F @ 91

Last change on this file since 91 was 37, checked in by emillour, 14 years ago

Remise en route chantier compilation -- Ehouarn

Modifs et corrections pour pouvoir compiler le gcm (en séentiel, avec

makelmdz_fcm pour l'instant):

ajout de fichiers 'arch' pour linux-64 (pour Bellonzi, avec ioipsl et en r8)
modification de makelmdz_fcm, ajout de la cléPP_PHYS si on compile avec une physique
correction de quelques typos/bugs réléàa compilation:
infotrac.F90 : supression des appels àlnblnk' (remplacépar len_trim)
bilan_dyn.F : déaration des variables znom3,znom3l,zunites3, planet_type
cpdet.F : "use control_mod, ONLY: planet_type" mis aux bons endroits

(idem sur cpdet.F dans dyn3dpar)

leapfrog.F : declaration de ztetaec(), dtec, cpdet , itau_w, duspg()
diagedyn.F : correction typo; attention dans diagedyn.F il y a du

include "../phylmd/YOMCST.h"
Ca va poser problè qd on change de physique....

Avec ces modifs, la compilation marche sans physque, avec ou sans ioipsl et avec la physique terrestre phylmd.

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
	917	do j=1,jjp1
	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