Context Navigation

bilan_dyn_loc.F @ 2173

Last change on this file since 2173 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

File size: 23.7 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
	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
[1848]	159	type(Request),SAVE :: Req
	160	!$OMP THREADPRIVATE(Req)
[1632]	161
	162	! definition du domaine d'ecriture pour le rebuild
	163
	164	INTEGER,DIMENSION(1) :: ddid
	165	INTEGER,DIMENSION(1) :: dsg
	166	INTEGER,DIMENSION(1) :: dsl
	167	INTEGER,DIMENSION(1) :: dpf
	168	INTEGER,DIMENSION(1) :: dpl
	169	INTEGER,DIMENSION(1) :: dhs
	170	INTEGER,DIMENSION(1) :: dhe
	171
	172	INTEGER :: bilan_dyn_domain_id
	173
	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
	218	ALLOCATE(rlong(jjb_v:jje_v))
	219	ALLOCATE(rlatg(jjb_v:jje_v))
	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),
	287	. 1, rlong(jjb:jje), jjn, rlatg(jjb:jje),
	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
	516	!$OMP END DO NOWAIT
	517	enddo
	518
	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
	543	jjb=jj_begin
	544	jje=jj_end
	545
	546	! do iQ=1,nQ
	547	! do l=1,llm-1
	548	! do j=jjb,jje
	549	! do i=1,iip1
	550	! ww=-0.5w(i,j,l+1)(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
	551	! dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww
	552	! dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
	553	! enddo
	554	! enddo
	555	! enddo
	556	! enddo
	557
	558	do iQ=1,nQ
	559	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	560	do l=1,llm
	561	IF (l<llm) THEN
	562	do j=jjb,jje
	563	do i=1,iip1
	564	ww=-0.5w(i,j,l+1)(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
	565	dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww
	566	dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
	567	enddo
	568	enddo
	569	ENDIF
	570	IF (l>2) THEN
	571	do j=jjb,jje
	572	do i=1,iip1
	573	ww=-0.5w(i,j,l)(Q(i,j,l-1,iQ)+Q(i,j,l,iQ))
	574	dQ(i,j,l,iQ)=dQ(i,j,l,iQ)+ww
	575	enddo
	576	enddo
	577	ENDIF
	578	enddo
	579	!$OMP ENDDO NOWAIT
	580	enddo
	581	IF (prt_level > 5)
	582	. WRITE(lunout,*)'Apres les calculs fait a chaque pas'
	583	c=====================================================================
	584	c PAS DE TEMPS D'ECRITURE
	585	c=====================================================================
	586	if (icum.eq.ncum) then
	587	c=====================================================================
	588
	589	IF (prt_level > 5)
	590	. WRITE(lunout,*)'Pas d ecriture'
	591
	592	jjb=jj_begin
	593	jje=jj_end
	594
	595	c Normalisation
	596	do iQ=1,nQ
	597	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	598	do l=1,llm
	599	Q_cum(:,jjb:jje,l,iQ)=Q_cum(:,jjb:jje,l,iQ)
	600	. /masse_cum(:,jjb:jje,l)
	601	enddo
	602	!$OMP ENDDO NOWAIT
	603	enddo
	604
	605	zz=1./REAL(ncum)
	606
	607	!$OMP MASTER
	608	ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)*zz
	609	!$OMP END MASTER
	610
	611	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	612	DO l=1,llm
	613	masse_cum(:,jjb:jje,l)=masse_cum(:,jjb:jje,l)*zz
	614	flux_u_cum(:,jjb:jje,l)=flux_u_cum(:,jjb:jje,l)*zz
	615	flux_uQ_cum(:,jjb:jje,l,:)=flux_uQ_cum(:,jjb:jje,l,:)*zz
	616	dQ(:,jjb:jje,l,:)=dQ(:,jjb:jje,l,:)*zz
	617	ENDDO
	618	!$OMP ENDDO NOWAIT
	619
	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
	628
	629	jjb=jj_begin
	630	jje=jj_end
	631
	632
	633	c A retravailler eventuellement
	634	c division de dQ par la masse pour revenir aux bonnes grandeurs
	635	do iQ=1,nQ
	636	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	637	DO l=1,llm
	638	dQ(:,jjb:jje,l,iQ)=dQ(:,jjb:jje,l,iQ)/masse_cum(:,jjb:jje,l)
	639	ENDDO
	640	!$OMP ENDDO NOWAIT
	641	enddo
	642
	643	c=====================================================================
	644	c Transport méridien
	645	c=====================================================================
	646
	647	c cumul zonal des masses des mailles
	648	c ----------------------------------
	649	jjb=jj_begin
	650	jje=jj_end
	651	if (pole_sud) jje=jj_end-1
	652
	653	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	654	DO l=1,llm
	655	zv(jjb:jje,l)=0.
	656	zmasse(jjb:jje,l)=0.
	657	ENDDO
	658	!$OMP ENDDO NOWAIT
	659
	660	call Register_Hallo_u(masse_cum,llm,1,1,1,1,Req)
	661	do iQ=1,nQ
[1874]	662	call Register_Hallo_u(Q_cum(1,jjb_u,1,iQ),llm,0,1,1,0,Req)
[1632]	663	enddo
	664
	665	call SendRequest(Req)
	666	!$OMP BARRIER
	667	call WaitRequest(Req)
	668
	669	call massbar_loc(masse_cum,massebx,masseby)
	670
	671	jjb=jj_begin
	672	jje=jj_end
	673	if (pole_sud) jje=jj_end-1
	674
	675	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	676	do l=1,llm
	677	do j=jjb,jje
	678	do i=1,iim
	679	zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
	680	zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
	681	enddo
	682	zfactv(j,l)=cv(1,j)/zmasse(j,l)
	683	enddo
	684	enddo
	685	!$OMP ENDDO NOWAIT
	686
	687	c print*,'3OK'
	688	c --------------------------------------------------------------
	689	c calcul de la moyenne zonale du transport :
	690	c ------------------------------------------
	691	c
	692	c --
	693	c TOT : la circulation totale [ vq ]
	694	c
	695	c - -
	696	c MMC : mean meridional circulation [ v ] [ q ]
	697	c
	698	c ---- -- - -
	699	c TRS : transitoires [ v'q'] = [ vq ] - [ v q ]
	700	c
	701	c - * - * - - - -
	702	c STT : stationaires [ v q ] = [ v q ] - [ v ] [ q ]
	703	c
	704	c - -
	705	c on utilise aussi l'intermediaire TMP : [ v q ]
	706	c
	707	c la variable zfactv transforme un transport meridien cumule
	708	c en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
	709	c
	710	c --------------------------------------------------------------
	711
	712
	713	c ----------------------------------------
	714	c Transport dans le plan latitude-altitude
	715	c ----------------------------------------
	716
	717	jjb=jj_begin
	718	jje=jj_end
	719	if (pole_sud) jje=jj_end-1
	720
	721	zvQ=0.
	722	psiQ=0.
	723	do iQ=1,nQ
	724	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	725	do l=1,llm
	726	zvQtmp(:,l)=0.
	727	do j=jjb,jje
	728	c print*,'j,l,iQ=',j,l,iQ
	729	c Calcul des moyennes zonales du transort total et de zvQtmp
	730	do i=1,iim
	731	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
	732	s +flux_vQ_cum(i,j,l,iQ)
	733	zqy= 0.5(Q_cum(i,j,l,iQ)masse_cum(i,j,l)+
	734	s Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
	735	zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
	736	s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
	737	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
	738	enddo
	739	c print*,'aOK'
	740	c Decomposition
	741	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
	742	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
	743	zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
	744	zvQ(j,l,immc,iQ)=zv(j,l)zvQ(j,l,iave,iQ)zfactv(j,l)
	745	zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
	746	zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
	747	enddo
	748	enddo
	749	!$OMP ENDDO NOWAIT
	750	c fonction de courant meridienne pour la quantite Q
	751	!$OMP BARRIER
	752	!$OMP MASTER
	753	do l=llm,1,-1
	754	do j=jjb,jje
	755	psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
	756	enddo
	757	enddo
	758	!$OMP END MASTER
	759	!$OMP BARRIER
	760	enddo
	761
	762	c fonction de courant pour la circulation meridienne moyenne
	763	!$OMP BARRIER
	764	!$OMP MASTER
	765	psi(jjb:jje,:)=0.
	766	do l=llm,1,-1
	767	do j=jjb,jje
	768	psi(j,l)=psi(j,l+1)+zv(j,l)
	769	zv(j,l)=zv(j,l)*zfactv(j,l)
	770	enddo
	771	enddo
	772	!$OMP END MASTER
	773	!$OMP BARRIER
	774
	775	c print*,'4OK'
	776	c sorties proprement dites
	777	!$OMP MASTER
	778	if (i_sortie.eq.1) then
	779	jjb=jj_begin
	780	jje=jj_end
	781	jjn=jj_nb
	782	if (pole_sud) jje=jj_end-1
	783	if (pole_sud) jjn=jj_nb-1
	784	do iQ=1,nQ
	785	do itr=1,ntr
	786	call histwrite(fileid,znom(itr,iQ),itau,
	787	s zvQ(jjb:jje,:,itr,iQ)
	788	s ,jjn*llm,ndex3d)
	789	enddo
	790	call histwrite(fileid,'psi'//nom(iQ),
	791	s itau,psiQ(jjb:jje,1:llm,iQ)
	792	s ,jjn*llm,ndex3d)
	793	enddo
	794
	795	call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm)
	796	s ,jjn*llm,ndex3d)
	797	call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm)
	798	s ,jjn*llm,ndex3d)
	799	psi(jjb:jje,:)=psi(jjb:jje,:)*1.e-9
	800	call histwrite(fileid,'psi',itau,psi(jjb:jje,1:llm),
	801	s jjn*llm,ndex3d)
	802
	803	endif
	804
	805
	806	c -----------------
	807	c Moyenne verticale
	808	c -----------------
	809
	810	zamasse(jjb:jje)=0.
	811	do l=1,llm
	812	zamasse(jjb:jje)=zamasse(jjb:jje)+zmasse(jjb:jje,l)
	813	enddo
	814
	815	zavQ(jjb:jje,:,:)=0.
	816	do iQ=1,nQ
	817	do itr=2,ntr
	818	do l=1,llm
	819	zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)
	820	s +zvQ(jjb:jje,l,itr,iQ)
	821	s *zmasse(jjb:jje,l)
	822	enddo
	823	zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)/zamasse(jjb:jje)
	824	call histwrite(fileid,'a'//znom(itr,iQ),itau,
	825	s zavQ(jjb:jje,itr,iQ),jjn*llm,ndex3d)
	826	enddo
	827	enddo
	828	!$OMP END MASTER
	829	c on doit pouvoir tracer systematiquement la fonction de courant.
	830
	831	c=====================================================================
	832	c/////////////////////////////////////////////////////////////////////
	833	icum=0 !///////////////////////////////////////
	834	endif ! icum.eq.ncum !///////////////////////////////////////
	835	c/////////////////////////////////////////////////////////////////////
	836	c=====================================================================
	837
	838	return
	839	end

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

Download in other formats:

Original Format