Context Navigation

pentes_ini.F @ 1506

Last change on this file since 1506 was 764, checked in by Laurent Fairhead, 17 years ago
Merge entre la version V3_conv et le HEAD YM, JG, LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 13.3 KB

Rev	Line
[630]	1	!
	2	! $Header$
	3	!
	4	SUBROUTINE pentes_ini (q,w,masse,pbaru,pbarv,mode)
	5	IMPLICIT NONE
	6
	7	c=======================================================================
	8	c Adaptation LMDZ: A.Armengaud (LGGE)
	9	c ----------------
	10	c
	11	c ********************************************************************
	12	c Transport des traceurs par la methode des pentes
	13	c ********************************************************************
	14	c Reference possible : Russel. G.L., Lerner J.A.:
	15	c A new Finite-Differencing Scheme for Traceur Transport
	16	c Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81
	17	c ********************************************************************
	18	c q,w,masse,pbaru et pbarv
	19	c sont des arguments d'entree pour le s-pg ....
	20	c
	21	c=======================================================================
	22
	23
	24	#include "dimensions.h"
	25	#include "paramet.h"
	26	#include "comconst.h"
	27	#include "comvert.h"
	28	#include "comgeom2.h"
	29
	30	c Arguments:
	31	c ----------
	32	integer mode
	33	REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
	34	REAL q( iip1,jjp1,llm,0:3)
	35	REAL w( ip1jmp1,llm )
	36	REAL masse( iip1,jjp1,llm)
	37	c Local:
	38	c ------
	39	LOGICAL limit
	40	REAL sm ( iip1,jjp1, llm )
	41	REAL s0( iip1,jjp1,llm ), sx( iip1,jjp1,llm )
	42	REAL sy( iip1,jjp1,llm ), sz( iip1,jjp1,llm )
	43	real masn,mass,zz
	44	INTEGER i,j,l,iq
	45
	46	c modif Fred 24 03 96
	47
	48	real sinlon(iip1),sinlondlon(iip1)
	49	real coslon(iip1),coslondlon(iip1)
	50	save sinlon,coslon,sinlondlon,coslondlon
	51	real dyn1,dyn2,dys1,dys2
	52	real qpn,qps,dqzpn,dqzps
	53	real smn,sms,s0n,s0s,sxn(iip1),sxs(iip1)
	54	real qmin,zq,pente_max
	55	c
	56	REAL SSUM
	57	integer ismax,ismin,lati,latf
[764]	58	EXTERNAL SSUM, ismin,ismax
[630]	59	logical first
	60	save first
	61	c fin modif
	62
	63
	64	c modif Fred 24 03 96
	65	data first/.true./
	66
	67	limit = .TRUE.
	68	pente_max=2
	69	c if (mode.eq.1.or.mode.eq.3) then
	70	c if (mode.eq.1) then
	71	if (mode.ge.1) then
	72	lati=2
	73	latf=jjm
	74	else
	75	lati=1
	76	latf=jjp1
	77	endif
	78
	79	qmin=0.4995
	80	qmin=0.
	81	if(first) then
	82	print*,'SCHEMA AMONT NOUVEAU'
	83	first=.false.
	84	do i=2,iip1
	85	coslon(i)=cos(rlonv(i))
	86	sinlon(i)=sin(rlonv(i))
	87	coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
	88	sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
	89	print*,coslondlon(i),sinlondlon(i)
	90	enddo
	91	coslon(1)=coslon(iip1)
	92	coslondlon(1)=coslondlon(iip1)
	93	sinlon(1)=sinlon(iip1)
	94	sinlondlon(1)=sinlondlon(iip1)
	95	print*,'sum sinlondlon ',ssum(iim,sinlondlon,1)/sinlondlon(1)
	96	print*,'sum coslondlon ',ssum(iim,coslondlon,1)/coslondlon(1)
	97	DO l = 1,llm
	98	DO j = 1,jjp1
	99	DO i = 1,iip1
	100	q ( i,j,l,1 )=0.
	101	q ( i,j,l,2 )=0.
	102	q ( i,j,l,3 )=0.
	103	ENDDO
	104	ENDDO
	105	ENDDO
	106
	107	endif
	108	c Fin modif Fred
	109
	110	c *** q contient les qqtes de traceur avant l'advection
	111
	112	c *** Affectation des tableaux S a partir de Q
	113	c *** Rem : utilisation de SCOPY ulterieurement
	114
	115	DO l = 1,llm
	116	DO j = 1,jjp1
	117	DO i = 1,iip1
	118	s0( i,j,llm+1-l ) = q ( i,j,l,0 )
	119	sx( i,j,llm+1-l ) = q ( i,j,l,1 )
	120	sy( i,j,llm+1-l ) = q ( i,j,l,2 )
	121	sz( i,j,llm+1-l ) = q ( i,j,l,3 )
	122	ENDDO
	123	ENDDO
	124	ENDDO
	125
	126	c PRINT*,'----- S0 just before conversion -------'
	127	c PRINT*,'S0(16,12,1)=',s0(16,12,1)
	128	c PRINT*,'Q(16,12,1,4)=',q(16,12,1,4)
	129
	130	c *** On calcule la masse d'air en kg
	131
	132	DO l = 1,llm
	133	DO j = 1,jjp1
	134	DO i = 1,iip1
	135	sm ( i,j,llm+1-l)=masse( i,j,l )
	136	ENDDO
	137	ENDDO
	138	ENDDO
	139
	140	c *** On converti les champs S en atome (resp. kg)
	141	c *** Les routines d'advection traitent les champs
	142	c *** a advecter si ces derniers sont en atome (resp. kg)
	143	c *** A optimiser !!!
	144
	145	DO l = 1,llm
	146	DO j = 1,jjp1
	147	DO i = 1,iip1
	148	s0(i,j,l) = s0(i,j,l) * sm ( i,j,l )
	149	sx(i,j,l) = sx(i,j,l) * sm ( i,j,l )
	150	sy(i,j,l) = sy(i,j,l) * sm ( i,j,l )
	151	sz(i,j,l) = sz(i,j,l) * sm ( i,j,l )
	152	ENDDO
	153	ENDDO
	154	ENDDO
	155
	156	c ss0 = 0.
	157	c DO l = 1,llm
	158	c DO j = 1,jjp1
	159	c DO i = 1,iim
	160	c ss0 = ss0 + s0 ( i,j,l )
	161	c ENDDO
	162	c ENDDO
	163	c ENDDO
	164	c PRINT*, 'valeur tot s0 avant advection=',ss0
	165
	166	c *** Appel des subroutines d'advection en X, en Y et en Z
	167	c *** Advection avec "time-splitting"
	168
	169	c-----------------------------------------------------------
	170	c PRINT*,'----- S0 just before ADVX -------'
	171	c PRINT*,'S0(16,12,1)=',s0(16,12,1)
	172
	173	c-----------------------------------------------------------
	174	c do l=1,llm
	175	c do j=1,jjp1
	176	c do i=1,iip1
	177	c zq=s0(i,j,l)/sm(i,j,l)
	178	c if(zq.lt.qmin)
	179	c , print*,'avant advx1, s0(',i,',',j,',',l,')=',zq
	180	c enddo
	181	c enddo
	182	c enddo
	183	CCC
	184	if(mode.eq.2) then
	185	do l=1,llm
	186	s0s=0.
	187	s0n=0.
	188	dyn1=0.
	189	dys1=0.
	190	dyn2=0.
	191	dys2=0.
	192	smn=0.
	193	sms=0.
	194	do i=1,iim
	195	smn=smn+sm(i,1,l)
	196	sms=sms+sm(i,jjp1,l)
	197	s0n=s0n+s0(i,1,l)
	198	s0s=s0s+s0(i,jjp1,l)
	199	zz=sy(i,1,l)/sm(i,1,l)
	200	dyn1=dyn1+sinlondlon(i)*zz
	201	dyn2=dyn2+coslondlon(i)*zz
	202	zz=sy(i,jjp1,l)/sm(i,jjp1,l)
	203	dys1=dys1+sinlondlon(i)*zz
	204	dys2=dys2+coslondlon(i)*zz
	205	enddo
	206	do i=1,iim
	207	sy(i,1,l)=dyn1sinlon(i)+dyn2coslon(i)
	208	sy(i,jjp1,l)=dys1sinlon(i)+dys2coslon(i)
	209	enddo
	210	do i=1,iim
	211	s0(i,1,l)=s0n/smn+sy(i,1,l)
	212	s0(i,jjp1,l)=s0s/sms-sy(i,jjp1,l)
	213	enddo
	214
	215	s0(iip1,1,l)=s0(1,1,l)
	216	s0(iip1,jjp1,l)=s0(1,jjp1,l)
	217
	218	do i=1,iim
	219	sxn(i)=s0(i+1,1,l)-s0(i,1,l)
	220	sxs(i)=s0(i+1,jjp1,l)-s0(i,jjp1,l)
	221	c on rerentre les masses
	222	enddo
	223	do i=1,iim
	224	sy(i,1,l)=sy(i,1,l)*sm(i,1,l)
	225	sy(i,jjp1,l)=sy(i,jjp1,l)*sm(i,jjp1,l)
	226	s0(i,1,l)=s0(i,1,l)*sm(i,1,l)
	227	s0(i,jjp1,l)=s0(i,jjp1,l)*sm(i,jjp1,l)
	228	enddo
	229	sxn(iip1)=sxn(1)
	230	sxs(iip1)=sxs(1)
	231	do i=1,iim
	232	sx(i+1,1,l)=0.25(sxn(i)+sxn(i+1))sm(i+1,1,l)
	233	sx(i+1,jjp1,l)=0.25(sxs(i)+sxs(i+1))sm(i+1,jjp1,l)
	234	enddo
	235	s0(iip1,1,l)=s0(1,1,l)
	236	s0(iip1,jjp1,l)=s0(1,jjp1,l)
	237	sy(iip1,1,l)=sy(1,1,l)
	238	sy(iip1,jjp1,l)=sy(1,jjp1,l)
	239	sx(1,1,l)=sx(iip1,1,l)
	240	sx(1,jjp1,l)=sx(iip1,jjp1,l)
	241	enddo
	242	endif
	243
	244	if (mode.eq.4) then
	245	do l=1,llm
	246	do i=1,iip1
	247	sx(i,1,l)=0.
	248	sx(i,jjp1,l)=0.
	249	sy(i,1,l)=0.
	250	sy(i,jjp1,l)=0.
	251	enddo
	252	enddo
	253	endif
	254	call limx(s0,sx,sm,pente_max)
	255	c call minmaxq(zq,1.e33,-1.e33,'avant advx ')
	256	call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
	257	c call minmaxq(zq,1.e33,-1.e33,'avant advy ')
	258	if (mode.eq.4) then
	259	do l=1,llm
	260	do i=1,iip1
	261	sx(i,1,l)=0.
	262	sx(i,jjp1,l)=0.
	263	sy(i,1,l)=0.
	264	sy(i,jjp1,l)=0.
	265	enddo
	266	enddo
	267	endif
	268	call limy(s0,sy,sm,pente_max)
	269	call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )
	270	c call minmaxq(zq,1.e33,-1.e33,'avant advz ')
	271	do j=1,jjp1
	272	do i=1,iip1
	273	sz(i,j,1)=0.
	274	sz(i,j,llm)=0.
	275	enddo
	276	enddo
	277	call limz(s0,sz,sm,pente_max)
	278	call advz( limit,dtvr,w,sm,s0,sx,sy,sz )
	279	if (mode.eq.4) then
	280	do l=1,llm
	281	do i=1,iip1
	282	sx(i,1,l)=0.
	283	sx(i,jjp1,l)=0.
	284	sy(i,1,l)=0.
	285	sy(i,jjp1,l)=0.
	286	enddo
	287	enddo
	288	endif
	289	call limy(s0,sy,sm,pente_max)
	290	call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )
	291	do l=1,llm
	292	do j=1,jjp1
	293	sm(iip1,j,l)=sm(1,j,l)
	294	s0(iip1,j,l)=s0(1,j,l)
	295	sx(iip1,j,l)=sx(1,j,l)
	296	sy(iip1,j,l)=sy(1,j,l)
	297	sz(iip1,j,l)=sz(1,j,l)
	298	enddo
	299	enddo
	300
	301
	302	c call minmaxq(zq,1.e33,-1.e33,'avant advx ')
	303	if (mode.eq.4) then
	304	do l=1,llm
	305	do i=1,iip1
	306	sx(i,1,l)=0.
	307	sx(i,jjp1,l)=0.
	308	sy(i,1,l)=0.
	309	sy(i,jjp1,l)=0.
	310	enddo
	311	enddo
	312	endif
	313	call limx(s0,sx,sm,pente_max)
	314	call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
	315	c call minmaxq(zq,1.e33,-1.e33,'apres advx ')
	316	c do l=1,llm
	317	c do j=1,jjp1
	318	c do i=1,iip1
	319	c zq=s0(i,j,l)/sm(i,j,l)
	320	c if(zq.lt.qmin)
	321	c , print*,'apres advx2, s0(',i,',',j,',',l,')=',zq
	322	c enddo
	323	c enddo
	324	c enddo
	325	c *** On repasse les S dans la variable q directement 14/10/94
	326	c On revient a des rapports de melange en divisant par la masse
	327
	328	c En dehors des poles:
	329
	330	DO l = 1,llm
	331	DO j = 1,jjp1
	332	DO i = 1,iim
	333	q(i,j,llm+1-l,0)=s0(i,j,l)/sm(i,j,l)
	334	q(i,j,llm+1-l,1)=sx(i,j,l)/sm(i,j,l)
	335	q(i,j,llm+1-l,2)=sy(i,j,l)/sm(i,j,l)
	336	q(i,j,llm+1-l,3)=sz(i,j,l)/sm(i,j,l)
	337	ENDDO
	338	ENDDO
	339	ENDDO
	340
	341	c Traitements specifiques au pole
	342
	343	if(mode.ge.1) then
	344	DO l=1,llm
	345	c filtrages aux poles
	346	masn=ssum(iim,sm(1,1,l),1)
	347	mass=ssum(iim,sm(1,jjp1,l),1)
	348	qpn=ssum(iim,s0(1,1,l),1)/masn
	349	qps=ssum(iim,s0(1,jjp1,l),1)/mass
	350	dqzpn=ssum(iim,sz(1,1,l),1)/masn
	351	dqzps=ssum(iim,sz(1,jjp1,l),1)/mass
	352	do i=1,iip1
	353	q( i,1,llm+1-l,3)=dqzpn
	354	q( i,jjp1,llm+1-l,3)=dqzps
	355	q( i,1,llm+1-l,0)=qpn
	356	q( i,jjp1,llm+1-l,0)=qps
	357	enddo
	358	if(mode.eq.3) then
	359	dyn1=0.
	360	dys1=0.
	361	dyn2=0.
	362	dys2=0.
	363	do i=1,iim
	364	dyn1=dyn1+sinlondlon(i)*sy(i,1,l)/sm(i,1,l)
	365	dyn2=dyn2+coslondlon(i)*sy(i,1,l)/sm(i,1,l)
	366	dys1=dys1+sinlondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)
	367	dys2=dys2+coslondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)
	368	enddo
	369	do i=1,iim
	370	q(i,1,llm+1-l,2)=
	371	s (sinlon(i)dyn1+coslon(i)dyn2)
	372	q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)
	373	q(i,jjp1,llm+1-l,2)=
	374	s (sinlon(i)dys1+coslon(i)dys2)
	375	q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
	376	s -q(i,jjp1,llm+1-l,2)
	377	enddo
	378	endif
	379	if(mode.eq.1) then
	380	c on filtre les valeurs au bord de la "grande maille pole"
	381	dyn1=0.
	382	dys1=0.
	383	dyn2=0.
	384	dys2=0.
	385	do i=1,iim
	386	zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)
	387	dyn1=dyn1+sinlondlon(i)*zz
	388	dyn2=dyn2+coslondlon(i)*zz
	389	zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)
	390	dys1=dys1+sinlondlon(i)*zz
	391	dys2=dys2+coslondlon(i)*zz
	392	enddo
	393	do i=1,iim
	394	q(i,1,llm+1-l,2)=
	395	s (sinlon(i)dyn1+coslon(i)dyn2)/2.
	396	q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)
	397	q(i,jjp1,llm+1-l,2)=
	398	s (sinlon(i)dys1+coslon(i)dys2)/2.
	399	q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
	400	s -q(i,jjp1,llm+1-l,2)
	401	enddo
	402	q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)
	403	q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)
	404
	405	do i=1,iim
	406	sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)
	407	sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)
	408	enddo
	409	sxn(iip1)=sxn(1)
	410	sxs(iip1)=sxs(1)
	411	do i=1,iim
	412	q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))
	413	q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))
	414	enddo
	415	q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)
	416	q(1,jjp1,llm+1-l,1)=q(iip1,jjp1,llm+1-l,1)
	417
	418	endif
	419
	420	ENDDO
	421	endif
	422
	423	c bouclage en longitude
	424	do iq=0,3
	425	do l=1,llm
	426	do j=1,jjp1
	427	q(iip1,j,l,iq)=q(1,j,l,iq)
	428	enddo
	429	enddo
	430	enddo
	431
	432	c PRINT*, ' SORTIE DE PENTES --- ca peut glisser ....'
	433
	434	DO l = 1,llm
	435	DO j = 1,jjp1
	436	DO i = 1,iip1
	437	IF (q(i,j,l,0).lt.0.) THEN
	438	c PRINT*,'------------ BIP-----------'
	439	c PRINT*,'Q0(',i,j,l,')=',q(i,j,l,0)
	440	c PRINT*,'QX(',i,j,l,')=',q(i,j,l,1)
	441	c PRINT*,'QY(',i,j,l,')=',q(i,j,l,2)
	442	c PRINT*,'QZ(',i,j,l,')=',q(i,j,l,3)
	443	c PRINT*,' PBL EN SORTIE DE PENTES'
	444	q(i,j,l,0)=0.
	445	c STOP
	446	ENDIF
	447	ENDDO
	448	ENDDO
	449	ENDDO
	450
	451	c PRINT*, '-------------------------------------------'
	452
	453	do l=1,llm
	454	do j=1,jjp1
	455	do i=1,iip1
	456	if(q(i,j,l,0).lt.qmin)
	457	, print*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0)
	458	enddo
	459	enddo
	460	enddo
	461	RETURN
	462	END
	463
	464
	465
	466
	467
	468
	469
	470
	471
	472
	473
	474

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

Download in other formats:

Original Format