Context Navigation

bilan_dyn_p.F @ 3567

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

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

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

Download in other formats:

Original Format