Context Navigation

bilan_dyn_loc.F @ 5443

Last change on this file since 5443 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
File size: 23.8 KB

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

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

Download in other formats:

Original Format