Context Navigation

bilan_dyn_p.F @ 2231

Last change on this file since 2231 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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.1 KB

Rev	Line
[630]	1	!
[1279]	2	! $Id: bilan_dyn_p.F 1907 2013-11-26 13:10:46Z oboucher $
[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
[630]	19	IMPLICIT NONE
	20
	21	#include "dimensions.h"
	22	#include "paramet.h"
	23	#include "comconst.h"
	24	#include "comvert.h"
	25	#include "comgeom2.h"
	26	#include "temps.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
[1885]	59	integer,save :: icum,ncum
	60	!$OMP THREADPRIVATE(icum,ncum)
	61	logical,SAVE :: first=.true.
	62	!$OMP THREADPRIVATE(first)
[630]	63
[1885]	64	real zz,zqy
	65	real,save :: zfactv(jjm,llm)
[630]	66
[1885]	67	integer,parameter :: nQ=7
[630]	68
[1885]	69
[630]	70	cym character*6 nom(nQ)
	71	cym character*6 unites(nQ)
[1885]	72	character(len=6),save :: nom(nQ)
	73	character(len=6),save :: unites(nQ)
[630]	74
[1885]	75	character(len=10) file
[630]	76	integer ifile
	77	parameter (ifile=4)
	78
[1885]	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
[630]	82
[1885]	83	real,SAVE :: time=0.
	84	integer,SAVE :: itau=0.
	85	!$OMP THREADPRIVATE(time,itau)
[630]	86
	87	real ww
	88
	89	c variables dynamiques intermédiaires
[1885]	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)
[630]	97
	98	c champ contenant les scalaires advectés.
[1885]	99	real,save :: Q(iip1,jjp1,llm,nQ)
[630]	100
	101	c champs cumulés
[1885]	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)
[630]	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
[1885]	124	INTEGER,PARAMETER :: iave=1,itot=2,immc=3,itrs=4,istn=5
	125
[630]	126	character*3 ctrs(ntr)
	127	data ctrs/' ','TOT','MMC','TRS','STN'/
	128
[1885]	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)
[630]	132
[1885]	133	real,save :: zv(jjm,llm),psi(jjm,llm+1)
[630]	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
[1885]	147	integer,save :: ndex3d(jjm*llm)
[630]	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
[1885]	158	real,save :: rlong(jjm),rlatg(jjm)
[630]	159	integer :: jjb,jje,jjn,ijb,ije
[1885]	160	type(Request),SAVE :: Req
	161	!$OMP THREADPRIVATE(Req)
[630]	162
[985]	163	! definition du domaine d'ecriture pour le rebuild
[630]	164
[985]	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
[630]	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
[1885]	186	ndex3d=0
[630]	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
[1885]	199	else
	200	write(lunout,*) "bilan_dyn_p: ncum=",ncum
[630]	201	endif
	202
[1885]	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
[630]	212
[1885]	213	!$OMP MASTER
[630]	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
[985]	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)
[630]	262
[985]	263	call histbeg(trim(infile),
[630]	264	. 1, rlong(jjb:jje), jjn, rlatg(jjb:jje),
	265	. 1, 1, 1, jjn,
[985]	266	. tau0, zjulian, dt_cum, thoriid, fileid,
	267	. bilan_dyn_domain_id)
[630]	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
[1885]	339	!$OMP END MASTER
	340	!$OMP BARRIER
[630]	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
[1885]	352	! ucont(:,jjb:jje,:)=0
[630]	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)
[1885]	357	c$OMP BARRIER
[630]	358	call WaitRequest(Req)
[1885]	359	c$OMP BARRIER
[630]	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
[1885]	365	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	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
[1885]	370	!$OMP END DO
[630]	371
[1885]	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
[630]	383
	384	c=====================================================================
	385	c Cumul
	386	c=====================================================================
	387	c
	388	if(icum.EQ.0) then
[985]	389	jjb=jj_begin
	390	jje=jj_end
	391
[1885]	392	!$OMP MASTER
[985]	393	ps_cum(:,jjb:jje)=0.
[1885]	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
[630]	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
[985]	415	jjb=jj_begin
	416	jje=jj_end
	417
[1885]	418	!$OMP MASTER
[985]	419	ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)+ps(:,jjb:jje)
[1885]	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
[985]	431	if (pole_sud) jje=jj_end-1
	432
[1885]	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
[985]	440	jjb=jj_begin
	441	jje=jj_end
	442
[630]	443	do iQ=1,nQ
[1885]	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
[630]	450	enddo
	451
	452	c=====================================================================
	453	c FLUX ET TENDANCES
	454	c=====================================================================
	455
	456	c Flux longitudinal
	457	c -----------------
	458	do iQ=1,nQ
[1885]	459	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	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
[1885]	469	!$OMP END DO NOWAIT
[630]	470	enddo
	471
	472	c flux méridien
	473	c -------------
	474	do iQ=1,nQ
[985]	475	call Register_Hallo(Q(1,1,1,iQ),ip1jmp1,llm,0,1,1,0,Req)
	476	enddo
	477	call SendRequest(Req)
[1885]	478	!$OMP BARRIER
[985]	479	call WaitRequest(Req)
[1885]	480	!$OMP BARRIER
	481
[985]	482	jjb=jj_begin
	483	jje=jj_end
	484	if (pole_sud) jje=jj_end-1
	485
	486	do iQ=1,nQ
[1885]	487	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	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
[1885]	496	!$OMP END DO NOWAIT
[630]	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)
[1885]	507	!$OMP BARRIER
[630]	508	call WaitRequest(Req)
[1885]	509	c$OMP BARRIER
[630]	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)
[1885]	517	!$OMP BARRIER
[630]	518	call WaitRequest(Req)
[1885]	519	c$OMP BARRIER
[630]	520
	521	call convmas_p(flux_u_cum,flux_v_cum,convm)
	522	CALL vitvert_p(convm,w)
[1885]	523	!$OMP BARRIER
[630]	524
[985]	525	jjb=jj_begin
	526	jje=jj_end
	527
[630]	528	do iQ=1,nQ
[1885]	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
[630]	548	enddo
[1885]	549	!$OMP ENDDO NOWAIT
[630]	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
[1885]	562	jjb=jj_begin
	563	jje=jj_end
	564
[630]	565	c Normalisation
	566	do iQ=1,nQ
[1885]	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
[630]	573	enddo
[1885]	574
[1403]	575	zz=1./REAL(ncum)
[630]	576
[1885]	577	!$OMP MASTER
	578	ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)*zz
	579	!$OMP END MASTER
[630]	580
[1885]	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
[985]	590
	591	IF (pole_sud) jje=jj_end-1
[1885]	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
[985]	598
	599	jjb=jj_begin
	600	jje=jj_end
	601
	602
[630]	603	c A retravailler eventuellement
	604	c division de dQ par la masse pour revenir aux bonnes grandeurs
	605	do iQ=1,nQ
[1885]	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
[630]	611	enddo
	612
	613	c=====================================================================
	614	c Transport méridien
	615	c=====================================================================
	616
	617	c cumul zonal des masses des mailles
	618	c ----------------------------------
[985]	619	jjb=jj_begin
	620	jje=jj_end
	621	if (pole_sud) jje=jj_end-1
	622
[1885]	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
[985]	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)
[1885]	636	!$OMP BARRIER
[985]	637	call WaitRequest(Req)
[1885]	638	c$OMP BARRIER
[985]	639
	640	call massbar_p(masse_cum,massebx,masseby)
[630]	641
	642	jjb=jj_begin
	643	jje=jj_end
	644	if (pole_sud) jje=jj_end-1
	645
[1885]	646	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	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
[1885]	656	!$OMP ENDDO
[630]	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
[985]	688	jjb=jj_begin
	689	jje=jj_end
	690	if (pole_sud) jje=jj_end-1
	691
[630]	692	zvQ=0.
	693	psiQ=0.
	694	do iQ=1,nQ
[1885]	695	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	696	do l=1,llm
[1885]	697	zvQtmp(:,l)=0.
[630]	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
[1885]	720	!$OMP ENDDO NOWAIT
[630]	721	c fonction de courant meridienne pour la quantite Q
[1885]	722	!$OMP BARRIER
	723	!$OMP MASTER
[630]	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
[1885]	729	!$OMP END MASTER
	730	!$OMP BARRIER
	731	enddo ! of do iQ=1,nQ
[630]	732
	733	c fonction de courant pour la circulation meridienne moyenne
[1885]	734	!$OMP BARRIER
	735	!$OMP MASTER
[985]	736	psi(jjb:jje,:)=0.
[630]	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
[1885]	743	!$OMP END MASTER
	744	!$OMP BARRIER
[630]	745
	746	c print*,'4OK'
	747	c sorties proprement dites
[1885]	748	!$OMP MASTER
[630]	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
[985]	766	call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm)
[630]	767	s ,jjn*llm,ndex3d)
[985]	768	call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm)
[630]	769	s ,jjn*llm,ndex3d)
[985]	770	psi(jjb:jje,:)=psi(jjb:jje,:)*1.e-9
[630]	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
[985]	781	zamasse(jjb:jje)=0.
[630]	782	do l=1,llm
	783	zamasse(jjb:jje)=zamasse(jjb:jje)+zmasse(jjb:jje,l)
	784	enddo
[985]	785
	786	zavQ(jjb:jje,:,:)=0.
[630]	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
[1885]	799	!$OMP END MASTER
	800	!$OMP BARRIER
[630]	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

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

Download in other formats:

Original Format