Context Navigation

bilan_dyn_loc.F @ 1721

Last change on this file since 1721 was 1673, checked in by Laurent Fairhead, 12 years ago

Fin du phasage de la dynamique parallele localisee (petite memoire) avec le tronc LMDZ5 r1671
Il reste quelques routines a verifier (en particulier ce qui touche a l'etude des cas academiques)
et la validation a effectuer

End of the phasing of the localised (low memory) parallel dynamics package with the
LMDZ5 trunk (r1671)
Some routines still need some checking (in particular the academic cases) and some
validation is still required

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

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

Download in other formats:

Original Format