Context Navigation

bilan_dyn_p.F @ 3856

Last change on this file since 3856 was 2601, checked in by Ehouarn Millour, 8 years ago
Cleanup in the dynamics: turn temps.h into module temps_mod.F90 EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 22.2 KB

Rev	Line
[630]	1	!
[1279]	2	! $Id: bilan_dyn_p.F 2601 2016-07-24 09:51:55Z fhourdin $
[630]	3	!
	4	SUBROUTINE bilan_dyn_p (ntrac,dt_app,dt_cum,
	5	s ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac)
	6
	7	c AFAIRE
	8	c Prevoir en champ nq+1 le diagnostique de l'energie
	9	c en faisant Qzon=Cv T + L * ...
	10	c vQ..A=Cp T + L * ...
	11
[1279]	12	#ifdef CPP_IOIPSL
[630]	13	USE IOIPSL
[1279]	14	#endif
[1823]	15	USE parallel_lmdz
[630]	16	USE mod_hallo
	17	use misc_mod
[1885]	18	use write_field_p
[2597]	19	USE comconst_mod, ONLY: cpp, pi
[2600]	20	USE comvert_mod, ONLY: presnivs
[2601]	21	USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn
[2600]	22
[630]	23	IMPLICIT NONE
	24
	25	#include "dimensions.h"
	26	#include "paramet.h"
	27	#include "comgeom2.h"
	28	#include "iniprint.h"
	29
	30	c====================================================================
	31	c
	32	c Sous-programme consacre à des diagnostics dynamiques de base
	33	c
	34	c
	35	c De facon generale, les moyennes des scalaires Q sont ponderees par
	36	c la masse.
	37	c
	38	c Les flux de masse sont eux simplement moyennes.
	39	c
	40	c====================================================================
	41
	42	c Arguments :
	43	c ===========
	44
	45	integer ntrac
	46	real dt_app,dt_cum
	47	real ps(iip1,jjp1)
	48	real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
	49	real flux_u(iip1,jjp1,llm)
	50	real flux_v(iip1,jjm,llm)
	51	real teta(iip1,jjp1,llm)
	52	real phi(iip1,jjp1,llm)
	53	real ucov(iip1,jjp1,llm)
	54	real vcov(iip1,jjm,llm)
	55	real trac(iip1,jjp1,llm,ntrac)
	56
	57	c Local :
	58	c =======
	59
[1885]	60	integer,save :: icum,ncum
	61	!$OMP THREADPRIVATE(icum,ncum)
	62	logical,SAVE :: first=.true.
	63	!$OMP THREADPRIVATE(first)
[630]	64
[1885]	65	real zz,zqy
	66	real,save :: zfactv(jjm,llm)
[630]	67
[1885]	68	integer,parameter :: nQ=7
[630]	69
[1885]	70
[630]	71	cym character*6 nom(nQ)
	72	cym character*6 unites(nQ)
[1885]	73	character(len=6),save :: nom(nQ)
	74	character(len=6),save :: unites(nQ)
[630]	75
[1885]	76	character(len=10) file
[630]	77	integer ifile
	78	parameter (ifile=4)
	79
[1885]	80	integer,PARAMETER :: itemp=1,igeop=2,iecin=3,iang=4,iu=5
	81	INTEGER,PARAMETER :: iovap=6,iun=7
	82	integer,PARAMETER :: i_sortie=1
[630]	83
[1885]	84	real,SAVE :: time=0.
	85	integer,SAVE :: itau=0.
	86	!$OMP THREADPRIVATE(time,itau)
[630]	87
	88	real ww
	89
	90	c variables dynamiques intermédiaires
[1885]	91	REAL,save :: vcont(iip1,jjm,llm),ucont(iip1,jjp1,llm)
	92	REAL,save :: ang(iip1,jjp1,llm),unat(iip1,jjp1,llm)
	93	REAL,save :: massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
	94	REAL,save :: vorpot(iip1,jjm,llm)
	95	REAL,save :: w(iip1,jjp1,llm),ecin(iip1,jjp1,llm)
	96	REAL,save ::convm(iip1,jjp1,llm)
	97	REAL,save :: bern(iip1,jjp1,llm)
[630]	98
	99	c champ contenant les scalaires advectés.
[1885]	100	real,save :: Q(iip1,jjp1,llm,nQ)
[630]	101
	102	c champs cumulés
[1885]	103	real,save :: ps_cum(iip1,jjp1)
	104	real,save :: masse_cum(iip1,jjp1,llm)
	105	real,save :: flux_u_cum(iip1,jjp1,llm)
	106	real,save :: flux_v_cum(iip1,jjm,llm)
	107	real,save :: Q_cum(iip1,jjp1,llm,nQ)
	108	real,save :: flux_uQ_cum(iip1,jjp1,llm,nQ)
	109	real,save :: flux_vQ_cum(iip1,jjm,llm,nQ)
	110	real,save :: flux_wQ_cum(iip1,jjp1,llm,nQ)
	111	real,save :: dQ(iip1,jjp1,llm,nQ)
[630]	112
	113
	114	c champs de tansport en moyenne zonale
	115	integer ntr,itr
	116	parameter (ntr=5)
	117
	118	cym character*10 znom(ntr,nQ)
	119	cym character*20 znoml(ntr,nQ)
	120	cym character*10 zunites(ntr,nQ)
	121	character*10,save :: znom(ntr,nQ)
	122	character*20,save :: znoml(ntr,nQ)
	123	character*10,save :: zunites(ntr,nQ)
	124
[1885]	125	INTEGER,PARAMETER :: iave=1,itot=2,immc=3,itrs=4,istn=5
	126
[630]	127	character*3 ctrs(ntr)
	128	data ctrs/' ','TOT','MMC','TRS','STN'/
	129
[1885]	130	real,save :: zvQ(jjm,llm,ntr,nQ),zvQtmp(jjm,llm)
	131	real,save :: zavQ(jjm,ntr,nQ),psiQ(jjm,llm+1,nQ)
	132	real,save :: zmasse(jjm,llm),zamasse(jjm)
[630]	133
[1885]	134	real,save :: zv(jjm,llm),psi(jjm,llm+1)
[630]	135
	136	integer i,j,l,iQ
	137
	138
	139	c Initialisation du fichier contenant les moyennes zonales.
	140	c ---------------------------------------------------------
	141
	142	character*10 infile
	143
	144	integer fileid
	145	integer thoriid, zvertiid
	146	save fileid
	147
[1885]	148	integer,save :: ndex3d(jjm*llm)
[630]	149
	150	C Variables locales
	151	C
	152	integer tau0
	153	real zjulian
	154	character*3 str
	155	character*10 ctrac
	156	integer ii,jj
	157	integer zan, dayref
	158	C
[1885]	159	real,save :: rlong(jjm),rlatg(jjm)
[630]	160	integer :: jjb,jje,jjn,ijb,ije
[1885]	161	type(Request),SAVE :: Req
	162	!$OMP THREADPRIVATE(Req)
[630]	163
[985]	164	! definition du domaine d'ecriture pour le rebuild
[630]	165
[985]	166	INTEGER,DIMENSION(1) :: ddid
	167	INTEGER,DIMENSION(1) :: dsg
	168	INTEGER,DIMENSION(1) :: dsl
	169	INTEGER,DIMENSION(1) :: dpf
	170	INTEGER,DIMENSION(1) :: dpl
	171	INTEGER,DIMENSION(1) :: dhs
	172	INTEGER,DIMENSION(1) :: dhe
	173
	174	INTEGER :: bilan_dyn_domain_id
	175
	176
[630]	177	c=====================================================================
	178	c Initialisation
	179	c=====================================================================
	180	if (adjust) return
	181
	182	time=time+dt_app
	183	itau=itau+1
	184
	185	if (first) then
	186
[1885]	187	ndex3d=0
[630]	188
	189	icum=0
	190	c initialisation des fichiers
	191	first=.false.
	192	c ncum est la frequence de stokage en pas de temps
	193	ncum=dt_cum/dt_app
	194	if (abs(ncumdt_app-dt_cum).gt.1.e-5dt_app) then
	195	WRITE(lunout,*)
	196	. 'Pb : le pas de cumule doit etre multiple du pas'
	197	WRITE(lunout,*)'dt_app=',dt_app
	198	WRITE(lunout,*)'dt_cum=',dt_cum
	199	stop
[1885]	200	else
	201	write(lunout,*) "bilan_dyn_p: ncum=",ncum
[630]	202	endif
	203
[1885]	204	! if (i_sortie.eq.1) then
[2600]	205	! file='dynzon'
[1885]	206	! if (mpi_rank==0) then
[2600]	207	! call inigrads(ifile,1
[1885]	208	! s ,0.,180./pi,0.,0.,jjm,rlatv,-90.,90.,180./pi
	209	! s ,llm,presnivs,1.
	210	! s ,dt_cum,file,'dyn_zon ')
	211	! endif
	212	! endif
[630]	213
[1885]	214	!$OMP MASTER
[630]	215	nom(itemp)='T'
	216	nom(igeop)='gz'
	217	nom(iecin)='K'
	218	nom(iang)='ang'
	219	nom(iu)='u'
	220	nom(iovap)='ovap'
	221	nom(iun)='un'
	222
	223	unites(itemp)='K'
	224	unites(igeop)='m2/s2'
	225	unites(iecin)='m2/s2'
	226	unites(iang)='ang'
	227	unites(iu)='m/s'
	228	unites(iovap)='kg/kg'
	229	unites(iun)='un'
	230
	231
	232	c Initialisation du fichier contenant les moyennes zonales.
	233	c ---------------------------------------------------------
	234
	235	infile='dynzon'
	236
	237	zan = annee_ref
	238	dayref = day_ref
	239	CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
	240	tau0 = itau_dyn
	241
	242	rlong=0.
	243	rlatg=rlatv*180./pi
	244
	245	jjb=jj_begin
	246	jje=jj_end
	247	jjn=jj_nb
[985]	248	IF (pole_sud) THEN
	249	jjn=jj_nb-1
	250	jje=jj_end-1
	251	ENDIF
	252
	253	ddid=(/ 2 /)
	254	dsg=(/ jjm /)
	255	dsl=(/ jjn /)
	256	dpf=(/ jjb /)
	257	dpl=(/ jje /)
	258	dhs=(/ 0 /)
	259	dhe=(/ 0 /)
	260
	261	call flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe,
	262	. 'box',bilan_dyn_domain_id)
[630]	263
[985]	264	call histbeg(trim(infile),
[630]	265	. 1, rlong(jjb:jje), jjn, rlatg(jjb:jje),
	266	. 1, 1, 1, jjn,
[985]	267	. tau0, zjulian, dt_cum, thoriid, fileid,
	268	. bilan_dyn_domain_id)
[630]	269
	270	C
	271	C Appel a histvert pour la grille verticale
	272	C
	273	call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
	274	. llm, presnivs, zvertiid)
	275	C
	276	C Appels a histdef pour la definition des variables a sauvegarder
	277	do iQ=1,nQ
	278	do itr=1,ntr
	279	if(itr.eq.1) then
	280	znom(itr,iQ)=nom(iQ)
	281	znoml(itr,iQ)=nom(iQ)
	282	zunites(itr,iQ)=unites(iQ)
	283	else
	284	znom(itr,iQ)=ctrs(itr)//'v'//nom(iQ)
	285	znoml(itr,iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)
	286	zunites(itr,iQ)='m/s * '//unites(iQ)
	287	endif
	288	enddo
	289	enddo
	290
	291	c Declarations des champs avec dimension verticale
	292	c print*,'1HISTDEF'
	293	do iQ=1,nQ
	294	do itr=1,ntr
	295	IF (prt_level > 5)
	296	. WRITE(lunout,*)'var ',itr,iQ
	297	. ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
	298	call histdef(fileid,znom(itr,iQ),znoml(itr,iQ),
	299	. zunites(itr,iQ),1,jjn,thoriid,llm,1,llm,zvertiid,
	300	. 32,'ave(X)',dt_cum,dt_cum)
	301	enddo
	302	c Declarations pour les fonctions de courant
	303	c print*,'2HISTDEF'
	304	call histdef(fileid,'psi'//nom(iQ)
	305	. ,'stream fn. '//znoml(itot,iQ),
	306	. zunites(itot,iQ),1,jjn,thoriid,llm,1,llm,zvertiid,
	307	. 32,'ave(X)',dt_cum,dt_cum)
	308	enddo
	309
	310
	311	c Declarations pour les champs de transport d'air
	312	c print*,'3HISTDEF'
	313	call histdef(fileid, 'masse', 'masse',
	314	. 'kg', 1, jjn, thoriid, llm, 1, llm, zvertiid,
	315	. 32, 'ave(X)', dt_cum, dt_cum)
	316	call histdef(fileid, 'v', 'v',
	317	. 'm/s', 1, jjn, thoriid, llm, 1, llm, zvertiid,
	318	. 32, 'ave(X)', dt_cum, dt_cum)
	319	c Declarations pour les fonctions de courant
	320	c print*,'4HISTDEF'
	321	call histdef(fileid,'psi','stream fn. MMC ','mega t/s',
	322	. 1,jjn,thoriid,llm,1,llm,zvertiid,
	323	. 32,'ave(X)',dt_cum,dt_cum)
	324
	325
	326	c Declaration des champs 1D de transport en latitude
	327	c print*,'5HISTDEF'
	328	do iQ=1,nQ
	329	do itr=2,ntr
	330	call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ),
	331	. zunites(itr,iQ),1,jjn,thoriid,1,1,1,-99,
	332	. 32,'ave(X)',dt_cum,dt_cum)
	333	enddo
	334	enddo
	335
	336
	337	c print*,'8HISTDEF'
	338	CALL histend(fileid)
	339
[1885]	340	!$OMP END MASTER
	341	!$OMP BARRIER
[630]	342	endif
	343
	344
	345	c=====================================================================
	346	c Calcul des champs dynamiques
	347	c ----------------------------
	348
	349	jjb=jj_begin
	350	jje=jj_end
	351
	352	c énergie cinétique
[1885]	353	! ucont(:,jjb:jje,:)=0
[630]	354
	355	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Req)
	356	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Req)
	357	call SendRequest(Req)
[1885]	358	c$OMP BARRIER
[630]	359	call WaitRequest(Req)
[1885]	360	c$OMP BARRIER
[630]	361
	362	CALL covcont_p(llm,ucov,vcov,ucont,vcont)
	363	CALL enercin_p(vcov,ucov,vcont,ucont,ecin)
	364
	365	c moment cinétique
[1885]	366	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	367	do l=1,llm
	368	ang(:,jjb:jje,l)=ucov(:,jjb:jje,l)+constang(:,jjb:jje)
	369	unat(:,jjb:jje,l)=ucont(:,jjb:jje,l)*cu(:,jjb:jje)
	370	enddo
[1885]	371	!$OMP END DO
[630]	372
[1885]	373	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	374	DO l=1,llm
	375	Q(:,jjb:jje,l,itemp)=teta(:,jjb:jje,l)*pk(:,jjb:jje,l)/cpp
	376	Q(:,jjb:jje,l,igeop)=phi(:,jjb:jje,l)
	377	Q(:,jjb:jje,l,iecin)=ecin(:,jjb:jje,l)
	378	Q(:,jjb:jje,l,iang)=ang(:,jjb:jje,l)
	379	Q(:,jjb:jje,l,iu)=unat(:,jjb:jje,l)
	380	Q(:,jjb:jje,l,iovap)=trac(:,jjb:jje,l,1)
	381	Q(:,jjb:jje,l,iun)=1.
	382	ENDDO
	383	!$OMP END DO NOWAIT
[630]	384
	385	c=====================================================================
	386	c Cumul
	387	c=====================================================================
	388	c
	389	if(icum.EQ.0) then
[985]	390	jjb=jj_begin
	391	jje=jj_end
	392
[1885]	393	!$OMP MASTER
[985]	394	ps_cum(:,jjb:jje)=0.
[1885]	395	!$OMP END MASTER
	396	!$OMP BARRIER
	397
	398	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	399	DO l=1,llm
	400	masse_cum(:,jjb:jje,l)=0.
	401	flux_u_cum(:,jjb:jje,l)=0.
	402	Q_cum(:,jjb:jje,l,:)=0.
	403	flux_uQ_cum(:,jjb:jje,l,:)=0.
	404	if (pole_sud) jje=jj_end-1
	405	flux_v_cum(:,jjb:jje,l)=0.
	406	flux_vQ_cum(:,jjb:jje,l,:)=0.
	407	ENDDO
	408	!$OMP END DO NOWAIT
[630]	409	endif
	410
	411	IF (prt_level > 5)
	412	. WRITE(lunout,*)'dans bilan_dyn ',icum,'->',icum+1
	413	icum=icum+1
	414
	415	c accumulation des flux de masse horizontaux
[985]	416	jjb=jj_begin
	417	jje=jj_end
	418
[1885]	419	!$OMP MASTER
[985]	420	ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)+ps(:,jjb:jje)
[1885]	421	!$OMP END MASTER
	422	!$OMP BARRIER
	423
	424	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	425	DO l=1,llm
	426	masse_cum(:,jjb:jje,l)=masse_cum(:,jjb:jje,l)+masse(:,jjb:jje,l)
	427	flux_u_cum(:,jjb:jje,l)=flux_u_cum(:,jjb:jje,l)
	428	. +flux_u(:,jjb:jje,l)
	429	ENDDO
	430	!$OMP END DO NOWAIT
	431
[985]	432	if (pole_sud) jje=jj_end-1
	433
[1885]	434	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	435	DO l=1,llm
	436	flux_v_cum(:,jjb:jje,l)=flux_v_cum(:,jjb:jje,l)
	437	. +flux_v(:,jjb:jje,l)
	438	ENDDO
	439	!$OMP END DO NOWAIT
	440
[985]	441	jjb=jj_begin
	442	jje=jj_end
	443
[630]	444	do iQ=1,nQ
[1885]	445	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	446	DO l=1,llm
	447	Q_cum(:,jjb:jje,l,iQ)=Q_cum(:,jjb:jje,l,iQ)
	448	. +Q(:,jjb:jje,l,iQ)*masse(:,jjb:jje,l)
	449	ENDDO
	450	!$OMP END DO NOWAIT
[630]	451	enddo
	452
	453	c=====================================================================
	454	c FLUX ET TENDANCES
	455	c=====================================================================
	456
	457	c Flux longitudinal
	458	c -----------------
	459	do iQ=1,nQ
[1885]	460	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	461	do l=1,llm
	462	do j=jjb,jje
	463	do i=1,iim
	464	flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ)
	465	s +flux_u(i,j,l)0.5(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
	466	enddo
	467	flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
	468	enddo
	469	enddo
[1885]	470	!$OMP END DO NOWAIT
[630]	471	enddo
	472
	473	c flux méridien
	474	c -------------
	475	do iQ=1,nQ
[985]	476	call Register_Hallo(Q(1,1,1,iQ),ip1jmp1,llm,0,1,1,0,Req)
	477	enddo
	478	call SendRequest(Req)
[1885]	479	!$OMP BARRIER
[985]	480	call WaitRequest(Req)
[1885]	481	!$OMP BARRIER
	482
[985]	483	jjb=jj_begin
	484	jje=jj_end
	485	if (pole_sud) jje=jj_end-1
	486
	487	do iQ=1,nQ
[1885]	488	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	489	do l=1,llm
	490	do j=jjb,jje
	491	do i=1,iip1
	492	flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ)
	493	s +flux_v(i,j,l)0.5(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
	494	enddo
	495	enddo
	496	enddo
[1885]	497	!$OMP END DO NOWAIT
[630]	498	enddo
	499
	500
	501	c tendances
	502	c ---------
	503
	504	c convergence horizontale
	505	call Register_Hallo(flux_uQ_cum,ip1jmp1,llm,2,2,2,2,Req)
	506	call Register_Hallo(flux_vQ_cum,ip1jm,llm,2,2,2,2,Req)
	507	call SendRequest(Req)
[1885]	508	!$OMP BARRIER
[630]	509	call WaitRequest(Req)
[1885]	510	c$OMP BARRIER
[630]	511
	512	call convflu_p(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
	513
	514	c calcul de la vitesse verticale
	515	call Register_Hallo(flux_u_cum,ip1jmp1,llm,2,2,2,2,Req)
	516	call Register_Hallo(flux_v_cum,ip1jm,llm,2,2,2,2,Req)
	517	call SendRequest(Req)
[1885]	518	!$OMP BARRIER
[630]	519	call WaitRequest(Req)
[1885]	520	c$OMP BARRIER
[630]	521
	522	call convmas_p(flux_u_cum,flux_v_cum,convm)
	523	CALL vitvert_p(convm,w)
[1885]	524	!$OMP BARRIER
[630]	525
[985]	526	jjb=jj_begin
	527	jje=jj_end
	528
[630]	529	do iQ=1,nQ
[1885]	530	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	531	do l=1,llm
	532	IF (l<llm) THEN
	533	do j=jjb,jje
	534	do i=1,iip1
	535	ww=-0.5w(i,j,l+1)(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
	536	dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww
	537	dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
	538	enddo
	539	enddo
	540	ENDIF
	541	IF (l>2) THEN
	542	do j=jjb,jje
	543	do i=1,iip1
	544	ww=-0.5w(i,j,l)(Q(i,j,l-1,iQ)+Q(i,j,l,iQ))
	545	dQ(i,j,l,iQ)=dQ(i,j,l,iQ)+ww
	546	enddo
	547	enddo
	548	ENDIF
[630]	549	enddo
[1885]	550	!$OMP ENDDO NOWAIT
[630]	551	enddo
	552	IF (prt_level > 5)
	553	. WRITE(lunout,*)'Apres les calculs fait a chaque pas'
	554	c=====================================================================
	555	c PAS DE TEMPS D'ECRITURE
	556	c=====================================================================
	557	if (icum.eq.ncum) then
	558	c=====================================================================
	559
	560	IF (prt_level > 5)
	561	. WRITE(lunout,*)'Pas d ecriture'
	562
[1885]	563	jjb=jj_begin
	564	jje=jj_end
	565
[630]	566	c Normalisation
	567	do iQ=1,nQ
[1885]	568	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	569	do l=1,llm
	570	Q_cum(:,jjb:jje,l,iQ)=Q_cum(:,jjb:jje,l,iQ)
[2600]	571	. /masse_cum(:,jjb:jje,l)
[1885]	572	enddo
	573	!$OMP ENDDO NOWAIT
[630]	574	enddo
[1885]	575
[1403]	576	zz=1./REAL(ncum)
[630]	577
[1885]	578	!$OMP MASTER
	579	ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)*zz
	580	!$OMP END MASTER
[630]	581
[1885]	582	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	583	DO l=1,llm
	584	masse_cum(:,jjb:jje,l)=masse_cum(:,jjb:jje,l)*zz
	585	flux_u_cum(:,jjb:jje,l)=flux_u_cum(:,jjb:jje,l)*zz
	586	flux_uQ_cum(:,jjb:jje,l,:)=flux_uQ_cum(:,jjb:jje,l,:)*zz
	587	dQ(:,jjb:jje,l,:)=dQ(:,jjb:jje,l,:)*zz
	588	ENDDO
	589	!$OMP ENDDO NOWAIT
	590
[985]	591
	592	IF (pole_sud) jje=jj_end-1
[1885]	593	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	594	DO l=1,llm
	595	flux_v_cum(:,jjb:jje,l)=flux_v_cum(:,jjb:jje,l)*zz
	596	flux_vQ_cum(:,jjb:jje,l,:)=flux_vQ_cum(:,jjb:jje,l,:)*zz
	597	ENDDO
	598	!$OMP ENDDO
[985]	599
	600	jjb=jj_begin
	601	jje=jj_end
	602
	603
[630]	604	c A retravailler eventuellement
	605	c division de dQ par la masse pour revenir aux bonnes grandeurs
	606	do iQ=1,nQ
[1885]	607	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	608	DO l=1,llm
	609	dQ(:,jjb:jje,l,iQ)=dQ(:,jjb:jje,l,iQ)/masse_cum(:,jjb:jje,l)
	610	ENDDO
	611	!$OMP ENDDO NOWAIT
[630]	612	enddo
	613
	614	c=====================================================================
	615	c Transport méridien
	616	c=====================================================================
	617
	618	c cumul zonal des masses des mailles
	619	c ----------------------------------
[985]	620	jjb=jj_begin
	621	jje=jj_end
	622	if (pole_sud) jje=jj_end-1
	623
[1885]	624	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	625	DO l=1,llm
	626	zv(jjb:jje,l)=0.
	627	zmasse(jjb:jje,l)=0.
	628	ENDDO
	629	!$OMP ENDDO NOWAIT
[985]	630
	631	call Register_Hallo(masse_cum,ip1jmp1,llm,1,1,1,1,Req)
	632	do iQ=1,nQ
	633	call Register_Hallo(Q_cum(1,1,1,iQ),ip1jmp1,llm,0,1,1,0,Req)
	634	enddo
	635
	636	call SendRequest(Req)
[1885]	637	!$OMP BARRIER
[985]	638	call WaitRequest(Req)
[1885]	639	c$OMP BARRIER
[985]	640
	641	call massbar_p(masse_cum,massebx,masseby)
[630]	642
	643	jjb=jj_begin
	644	jje=jj_end
	645	if (pole_sud) jje=jj_end-1
	646
[1885]	647	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	648	do l=1,llm
	649	do j=jjb,jje
	650	do i=1,iim
	651	zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
	652	zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
	653	enddo
	654	zfactv(j,l)=cv(1,j)/zmasse(j,l)
	655	enddo
	656	enddo
[1885]	657	!$OMP ENDDO
[630]	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
	681	c
	682	c --------------------------------------------------------------
	683
	684
	685	c ----------------------------------------
	686	c Transport dans le plan latitude-altitude
	687	c ----------------------------------------
	688
[985]	689	jjb=jj_begin
	690	jje=jj_end
	691	if (pole_sud) jje=jj_end-1
	692
[630]	693	zvQ=0.
	694	psiQ=0.
	695	do iQ=1,nQ
[1885]	696	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	697	do l=1,llm
[1885]	698	zvQtmp(:,l)=0.
[630]	699	do j=jjb,jje
	700	c print*,'j,l,iQ=',j,l,iQ
	701	c Calcul des moyennes zonales du transort total et de zvQtmp
	702	do i=1,iim
	703	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
	704	s +flux_vQ_cum(i,j,l,iQ)
	705	zqy= 0.5(Q_cum(i,j,l,iQ)masse_cum(i,j,l)+
	706	s Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
	707	zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
	708	s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
	709	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
	710	enddo
	711	c print*,'aOK'
	712	c Decomposition
	713	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
	714	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
	715	zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
	716	zvQ(j,l,immc,iQ)=zv(j,l)zvQ(j,l,iave,iQ)zfactv(j,l)
	717	zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
	718	zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
	719	enddo
	720	enddo
[1885]	721	!$OMP ENDDO NOWAIT
[630]	722	c fonction de courant meridienne pour la quantite Q
[1885]	723	!$OMP BARRIER
	724	!$OMP MASTER
[630]	725	do l=llm,1,-1
	726	do j=jjb,jje
	727	psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
	728	enddo
	729	enddo
[1885]	730	!$OMP END MASTER
	731	!$OMP BARRIER
	732	enddo ! of do iQ=1,nQ
[630]	733
	734	c fonction de courant pour la circulation meridienne moyenne
[1885]	735	!$OMP BARRIER
	736	!$OMP MASTER
[985]	737	psi(jjb:jje,:)=0.
[630]	738	do l=llm,1,-1
	739	do j=jjb,jje
	740	psi(j,l)=psi(j,l+1)+zv(j,l)
	741	zv(j,l)=zv(j,l)*zfactv(j,l)
	742	enddo
	743	enddo
[1885]	744	!$OMP END MASTER
	745	!$OMP BARRIER
[630]	746
	747	c print*,'4OK'
	748	c sorties proprement dites
[1885]	749	!$OMP MASTER
[630]	750	if (i_sortie.eq.1) then
	751	jjb=jj_begin
	752	jje=jj_end
	753	jjn=jj_nb
	754	if (pole_sud) jje=jj_end-1
	755	if (pole_sud) jjn=jj_nb-1
	756
	757	do iQ=1,nQ
	758	do itr=1,ntr
	759	call histwrite(fileid,znom(itr,iQ),itau,
	760	s zvQ(jjb:jje,:,itr,iQ)
	761	s ,jjn*llm,ndex3d)
	762	enddo
	763	call histwrite(fileid,'psi'//nom(iQ),
	764	s itau,psiQ(jjb:jje,1:llm,iQ)
	765	s ,jjn*llm,ndex3d)
	766	enddo
[985]	767	call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm)
[630]	768	s ,jjn*llm,ndex3d)
[985]	769	call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm)
[630]	770	s ,jjn*llm,ndex3d)
[985]	771	psi(jjb:jje,:)=psi(jjb:jje,:)*1.e-9
[630]	772	call histwrite(fileid,'psi',itau,psi(jjb:jje,1:llm),
	773	s jjn*llm,ndex3d)
	774
	775	endif
	776
	777
	778	c -----------------
	779	c Moyenne verticale
	780	c -----------------
	781
[985]	782	zamasse(jjb:jje)=0.
[630]	783	do l=1,llm
	784	zamasse(jjb:jje)=zamasse(jjb:jje)+zmasse(jjb:jje,l)
	785	enddo
[985]	786
	787	zavQ(jjb:jje,:,:)=0.
[630]	788	do iQ=1,nQ
	789	do itr=2,ntr
	790	do l=1,llm
	791	zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)
	792	s +zvQ(jjb:jje,l,itr,iQ)
	793	s *zmasse(jjb:jje,l)
	794	enddo
	795	zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)/zamasse(jjb:jje)
	796	call histwrite(fileid,'a'//znom(itr,iQ),itau,
	797	s zavQ(jjb:jje,itr,iQ),jjn*llm,ndex3d)
	798	enddo
	799	enddo
[1885]	800	!$OMP END MASTER
	801	!$OMP BARRIER
[630]	802	c on doit pouvoir tracer systematiquement la fonction de courant.
	803
	804	c=====================================================================
	805	c/////////////////////////////////////////////////////////////////////
	806	icum=0 !///////////////////////////////////////
	807	endif ! icum.eq.ncum !///////////////////////////////////////
	808	c/////////////////////////////////////////////////////////////////////
	809	c=====================================================================
	810	return
	811	end

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

Download in other formats:

Original Format