Context Navigation

grid_noro.F @ 822

Last change on this file since 822 was 404, checked in by lmdzadmin, 22 years ago
seuil a 10% de terre dans parametrisation orographie PB IM/LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 15.2 KB

Rev	Line
[370]	1	c
	2	c $Header$
	3	c
[2]	4	SUBROUTINE grid_noro(imdep, jmdep, xdata, ydata, zdata,
	5	. imar, jmar, x, y,
	6	. zphi,zmea,zstd,zsig,zgam,zthe,
	7	. zpic,zval,mask)
	8	c=======================================================================
	9	c (F. Lott) (voir aussi z.x. Li, A. Harzallah et L. Fairhead)
	10	c
	11	c Compute the Parameters of the SSO scheme as described in
	12	c LOTT & MILLER (1997) and LOTT(1999).
	13	c Target points are on a rectangular grid:
	14	c iim+1 latitudes including North and South Poles;
	15	c jjm+1 longitudes, with periodicity jjm+1=1.
	16	c aux poles. At the poles the fields value is repeated
	17	c jjm+1 time.
	18	c The parameters a,b,c,d represent the limite of the target
	19	c gridpoint region. The means over this region are calculated
	20	c from USN data, ponderated by a weight proportional to the
	21	c surface occupated by the data inside the model gridpoint area.
	22	c In most circumstances, this weight is the ratio between the
	23	c surface of the USN gridpoint area and the surface of the
	24	c model gridpoint area.
	25	c
	26	c (c)
	27	c ----d-----
	28	c \| . . . .\|
	29	c \| \|
	30	c (b)a . * . .b(a)
	31	c \| \|
	32	c \| . . . .\|
	33	c ----c-----
	34	c (d)
	35	C=======================================================================
	36	c INPUT:
	37	c imdep, jmdep: dimensions X and Y input field
	38	c xdata, ydata: coordinates X and Y input field
	39	c zdata: Input field
	40	c In this version it is assumed that the entry data come from
	41	c the USNavy dataset: imdep=iusn=2160, jmdep=jusn=1080.
	42	c OUTPUT:
	43	c imar, jmar: dimensions X and Y Output field
	44	c x, y: ccordinates X and Y Output field.
	45	c zmea: Mean orographie
	46	c zstd: Standard deviation
	47	c zsig: Slope
	48	c zgam: Anisotropy
	49	c zthe: Orientation of the small axis
	50	c zpic: Maximum altitude
	51	c zval: Minimum altitude
	52	C=======================================================================
	53	IMPLICIT INTEGER (I,J)
	54	IMPLICIT REAL(X,Z)
	55
[228]	56	parameter(iusn=2160,jusn=1080,iext=216, epsfra = 1.e-5)
[2]	57	#include "dimensions.h"
	58	REAL xusn(iusn+2*iext),yusn(jusn+2)
	59	REAL zusn(iusn+2*iext,jusn+2)
	60
	61	INTEGER imdep, jmdep
	62	REAL xdata(imdep),ydata(jmdep)
	63	REAL zdata(imdep,jmdep)
	64	c
	65	INTEGER imar, jmar
	66
	67	C INTERMEDIATE FIELDS (CORRELATIONS OF OROGRAPHY GRADIENT)
	68
	69	REAL ztz(iim+1,jjm+1),zxtzx(iim+1,jjm+1)
	70	REAL zytzy(iim+1,jjm+1),zxtzy(iim+1,jjm+1)
	71	REAL weight(iim+1,jjm+1)
	72
	73	C CORRELATIONS OF USN OROGRAPHY GRADIENTS
	74
	75	REAL zxtzxusn(iusn+2iext,jusn+2),zytzyusn(iusn+2iext,jusn+2)
	76	REAL zxtzyusn(iusn+2*iext,jusn+2)
	77	REAL x(imar+1),y(jmar),zphi(imar+1,jmar)
	78	REAL zmea(imar+1,jmar),zstd(imar+1,jmar)
	79	REAL zsig(imar+1,jmar),zgam(imar+1,jmar),zthe(imar+1,jmar)
	80	REAL zpic(imar+1,jmar),zval(imar+1,jmar)
[99]	81	c$$$ PB integer mask(imar+1,jmar)
[228]	82	real mask(imar+1,jmar), mask_tmp(imar+1,jmar)
[2]	83	real num_tot(2200,1100),num_lan(2200,1100)
	84	c
	85	REAL a(2200),b(2200),c(1100),d(1100)
[370]	86	logical masque_lu
[2]	87	c
	88	print *,' parametres de l orographie a l echelle sous maille'
	89	xpi=acos(-1.)
	90	rad = 6 371 229.
	91	zdeltay=2.xpi/float(jusn)rad
	92	c
[370]	93	c utilise-t'on un masque lu?
	94	c
	95	masque_lu = .true.
	96	if (maxval(mask) == -99999 .and. minval(mask) == -99999) then
	97	masque_lu= .false.
	98	masque = 0.0
	99	endif
	100	write(,)'Masque lu', masque_lu
	101	c
[2]	102	c quelques tests de dimensions:
	103	c
	104	c
	105	if(iim.ne.imar) STOP 'Problem dim. x'
	106	if(jjm.ne.jmar-1) STOP 'Problem dim. y'
	107	IF (imar.GT.2200 .OR. jmar.GT.1100) THEN
	108	PRINT*, 'imar or jmar too big', imar, jmar
	109	CALL ABORT
	110	ENDIF
	111
	112	IF(imdep.ne.iusn.or.jmdep.ne.jusn)then
	113	print *,' imdep or jmdep bad dimensions:',imdep,jmdep
	114	call abort
	115	ENDIF
	116
	117	IF(imar+1.ne.iim+1.or.jmar.ne.jjm+1)THEN
	118	print *,' imar or jmar bad dimensions:',imar,jmar
	119	call abort
	120	ENDIF
	121
	122
	123	c print *,'xdata:',xdata
	124	c print *,'ydata:',ydata
	125	c print *,'x:',x
	126	c print *,'y:',y
	127	c
	128	C EXTENSION OF THE USN DATABASE TO POCEED COMPUTATIONS AT
	129	C BOUNDARIES:
	130	c
	131	DO j=1,jusn
	132	yusn(j+1)=ydata(j)
	133	DO i=1,iusn
	134	zusn(i+iext,j+1)=zdata(i,j)
	135	xusn(i+iext)=xdata(i)
	136	ENDDO
	137	DO i=1,iext
	138	zusn(i,j+1)=zdata(iusn-iext+i,j)
	139	xusn(i)=xdata(iusn-iext+i)-2.*xpi
	140	zusn(iusn+iext+i,j+1)=zdata(i,j)
	141	xusn(iusn+iext+i)=xdata(i)+2.*xpi
	142	ENDDO
	143	ENDDO
	144
	145	yusn(1)=ydata(1)+(ydata(1)-ydata(2))
	146	yusn(jusn+2)=ydata(jusn)+(ydata(jusn)-ydata(jusn-1))
	147	DO i=1,iusn/2+iext
	148	zusn(i,1)=zusn(i+iusn/2,2)
	149	zusn(i+iusn/2+iext,1)=zusn(i,2)
	150	zusn(i,jusn+2)=zusn(i+iusn/2,jusn+1)
	151	zusn(i+iusn/2+iext,jusn+2)=zusn(i,jusn+1)
	152	ENDDO
	153	c
	154	c COMPUTE LIMITS OF MODEL GRIDPOINT AREA
	155	C ( REGULAR GRID)
	156	c
	157	a(1) = x(1) - (x(2)-x(1))/2.0
	158	b(1) = (x(1)+x(2))/2.0
	159	DO i = 2, imar
	160	a(i) = b(i-1)
	161	b(i) = (x(i)+x(i+1))/2.0
	162	ENDDO
	163	a(imar+1) = b(imar)
	164	b(imar+1) = x(imar+1) + (x(imar+1)-x(imar))/2.0
	165
	166	c(1) = y(1) - (y(2)-y(1))/2.0
	167	d(1) = (y(1)+y(2))/2.0
	168	DO j = 2, jmar-1
	169	c(j) = d(j-1)
	170	d(j) = (y(j)+y(j+1))/2.0
	171	ENDDO
	172	c(jmar) = d(jmar-1)
	173	d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0
	174	c
	175	c initialisations:
	176	c
	177	DO i = 1, imar+1
	178	DO j = 1, jmar
	179	weight(i,j) = 0.0
	180	zxtzx(i,j) = 0.0
	181	zytzy(i,j) = 0.0
	182	zxtzy(i,j) = 0.0
	183	ztz(i,j) = 0.0
	184	zmea(i,j) = 0.0
	185	zpic(i,j) =-1.E+10
	186	zval(i,j) = 1.E+10
	187	ENDDO
	188	ENDDO
	189	c
	190	c COMPUTE SLOPES CORRELATIONS ON USN GRID
	191	c
	192	DO j = 1,jusn+2
	193	DO i = 1, iusn+2*iext
	194	zytzyusn(i,j)=0.0
	195	zxtzxusn(i,j)=0.0
	196	zxtzyusn(i,j)=0.0
	197	ENDDO
	198	ENDDO
	199
	200
	201	DO j = 2,jusn+1
	202	zdeltax=zdeltay*cos(yusn(j))
	203	DO i = 2, iusn+2*iext-1
	204	zytzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))2/zdeltay2
	205	zxtzxusn(i,j)=(zusn(i+1,j)-zusn(i-1,j))2/zdeltax2
	206	zxtzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))/zdeltay
	207	* *(zusn(i+1,j)-zusn(i-1,j))/zdeltax
	208	ENDDO
	209	ENDDO
	210	c
	211	c SUMMATION OVER GRIDPOINT AREA
	212	c
	213	zleny=xpi/float(jusn)*rad
	214	xincr=xpi/2./float(jusn)
	215	DO ii = 1, imar+1
	216	DO jj = 1, jmar
	217	num_tot(ii,jj)=0.
	218	num_lan(ii,jj)=0.
	219	c PRINT *,' iteration ii jj:',ii,jj
	220	DO j = 2,jusn+1
	221	c DO j = 3,jusn
	222	zlenx=zleny*cos(yusn(j))
	223	zdeltax=zdeltay*cos(yusn(j))
	224	zbordnor=(c(jj)-yusn(j)+xincr)*rad
	225	zbordsud=(yusn(j)-d(jj)+xincr)*rad
	226	weighy=AMAX1(0.,
	227	* amin1(zbordnor,zbordsud,zleny))
	228	IF(weighy.ne.0)THEN
	229	DO i = 2, iusn+2*iext-1
	230	zbordest=(xusn(i)-a(ii)+xincr)radcos(yusn(j))
	231	zbordoue=(b(ii)+xincr-xusn(i))radcos(yusn(j))
	232	weighx=AMAX1(0.,
	233	* amin1(zbordest,zbordoue,zlenx))
	234	IF(weighx.ne.0)THEN
	235	num_tot(ii,jj)=num_tot(ii,jj)+1.0
	236	if(zusn(i,j).ge.1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
	237	weight(ii,jj)=weight(ii,jj)+weighx*weighy
	238	zxtzx(ii,jj)=zxtzx(ii,jj)+zxtzxusn(i,j)weighxweighy
	239	zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)weighxweighy
	240	zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)weighxweighy
	241	ztz(ii,jj) =ztz(ii,jj) +zusn(i,j)zusn(i,j)weighx*weighy
	242	c mean
	243	zmea(ii,jj) =zmea(ii,jj)+zusn(i,j)weighxweighy
	244	c peacks
	245	zpic(ii,jj)=amax1(zpic(ii,jj),zusn(i,j))
	246	c valleys
	247	zval(ii,jj)=amin1(zval(ii,jj),zusn(i,j))
	248	ENDIF
	249	ENDDO
	250	ENDIF
	251	ENDDO
	252	ENDDO
	253	ENDDO
	254	c
	255	c COMPUTE PARAMETERS NEEDED BY THE LOTT & MILLER (1997) AND
	256	C LOTT (1999) SSO SCHEME.
	257	c
	258	zllmmea=0.
	259	zllmstd=0.
	260	zllmsig=0.
	261	zllmgam=0.
	262	zllmpic=0.
	263	zllmval=0.
	264	zllmthe=0.
	265	zminthe=0.
	266	c print 100,' '
	267	c100 format(1X,A1,'II JJ',4X,'H',8X,'SD',8X,'SI',3X,'GA',3X,'TH')
	268	DO ii = 1, imar+1
	269	DO jj = 1, jmar
	270	IF (weight(ii,jj) .NE. 0.0) THEN
	271	c Mask
[99]	272	c$$$ if(num_lan(ii,jj)/num_tot(ii,jj).ge.0.5)then
	273	c$$$ mask(ii,jj)=1
	274	c$$$ else
	275	c$$$ mask(ii,jj)=0
	276	c$$$ ENDIF
[370]	277	if (.not. masque_lu) then
	278	mask(ii,jj) = num_lan(ii,jj)/num_tot(ii,jj)
	279	endif
[2]	280	c Mean Orography:
	281	zmea (ii,jj)=zmea (ii,jj)/weight(ii,jj)
	282	zxtzx(ii,jj)=zxtzx(ii,jj)/weight(ii,jj)
	283	zytzy(ii,jj)=zytzy(ii,jj)/weight(ii,jj)
	284	zxtzy(ii,jj)=zxtzy(ii,jj)/weight(ii,jj)
	285	ztz(ii,jj) =ztz(ii,jj)/weight(ii,jj)
	286	c Standard deviation:
	287	zstd(ii,jj)=sqrt(AMAX1(0.,ztz(ii,jj)-zmea(ii,jj)**2))
	288	ELSE
	289	PRINT*, 'probleme,ii,jj=', ii,jj
	290	ENDIF
	291	ENDDO
	292	ENDDO
	293
	294	C CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:
	295
	296	DO ii = 1, imar+1
	297	zxtzx(ii,1)=zxtzx(ii,2)
	298	zxtzx(ii,jmar)=zxtzx(ii,jmar-1)
	299	zxtzy(ii,1)=zxtzy(ii,2)
	300	zxtzy(ii,jmar)=zxtzy(ii,jmar-1)
	301	zytzy(ii,1)=zytzy(ii,2)
	302	zytzy(ii,jmar)=zytzy(ii,jmar-1)
	303	ENDDO
	304
	305	C FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.
	306
	307	C FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.
	308
	309	CALL MVA9(zmea,iim+1,jjm+1)
	310	CALL MVA9(zstd,iim+1,jjm+1)
	311	CALL MVA9(zpic,iim+1,jjm+1)
	312	CALL MVA9(zval,iim+1,jjm+1)
	313	CALL MVA9(zxtzx,iim+1,jjm+1)
[99]	314	CALL MVA9(zxtzy,iim+1,jjm+1)
[369]	315	CALL MVA9(zytzy,iim+1,jjm+1)
[228]	316	C$$$ Masque prenant en compte maximum de terre
[404]	317	C$$$ On seuil a 10% de terre de terre car en dessous les parametres de surface n'on
	318	C$$$ pas de sens (PB)
[228]	319	mask_tmp= 0.0
[404]	320	WHERE(mask .GE. 0.1) mask_tmp = 1.
[2]	321
	322	DO ii = 1, imar
	323	DO jj = 1, jmar
	324	IF (weight(ii,jj) .NE. 0.0) THEN
	325	c Coefficients K, L et M:
	326	xk=(zxtzx(ii,jj)+zytzy(ii,jj))/2.
	327	xl=(zxtzx(ii,jj)-zytzy(ii,jj))/2.
	328	xm=zxtzy(ii,jj)
	329	xp=xk-sqrt(xl2+xm2)
	330	xq=xk+sqrt(xl2+xm2)
	331	xw=1.e-8
	332	if(xp.le.xw) xp=0.
	333	if(xq.le.xw) xq=xw
	334	if(abs(xm).le.xw) xm=xw*sign(1.,xm)
[228]	335	c slope:
[99]	336	c$$$ zsig(ii,jj)=sqrt(xq)*mask(ii,jj)
	337	c$$$c isotropy:
	338	c$$$ zgam(ii,jj)=xp/xq*mask(ii,jj)
	339	c$$$c angle theta:
	340	c$$$ zthe(ii,jj)=57.29577951atan2(xm,xl)/2.mask(ii,jj)
	341	c$$$ zphi(ii,jj)=zmea(ii,jj)*mask(ii,jj)
	342	c$$$ zmea(ii,jj)=zmea(ii,jj)*mask(ii,jj)
	343	c$$$ zpic(ii,jj)=zpic(ii,jj)*mask(ii,jj)
	344	c$$$ zval(ii,jj)=zval(ii,jj)*mask(ii,jj)
	345	c$$$ zstd(ii,jj)=zstd(ii,jj)*mask(ii,jj)
	346	C$$* PB modif pour maque de terre fractionnaire
[2]	347	c slope:
[228]	348	zsig(ii,jj)=sqrt(xq)*mask_tmp(ii,jj)
[2]	349	c isotropy:
[228]	350	zgam(ii,jj)=xp/xq*mask_tmp(ii,jj)
[2]	351	c angle theta:
[228]	352	zthe(ii,jj)=57.29577951atan2(xm,xl)/2.mask_tmp(ii,jj)
	353	zphi(ii,jj)=zmea(ii,jj)*mask_tmp(ii,jj)
	354	zmea(ii,jj)=zmea(ii,jj)*mask_tmp(ii,jj)
	355	zpic(ii,jj)=zpic(ii,jj)*mask_tmp(ii,jj)
	356	zval(ii,jj)=zval(ii,jj)*mask_tmp(ii,jj)
	357	zstd(ii,jj)=zstd(ii,jj)*mask_tmp(ii,jj)
[2]	358	c print 101,ii,jj,
	359	c * zmea(ii,jj),zstd(ii,jj),zsig(ii,jj),zgam(ii,jj),
	360	c * zthe(ii,jj)
	361	c101 format(1x,2(1x,i2),2(1x,f7.1),1x,f7.4,2x,f4.2,1x,f5.1)
	362	ELSE
	363	c PRINT*, 'probleme,ii,jj=', ii,jj
	364	ENDIF
	365	zllmmea=AMAX1(zmea(ii,jj),zllmmea)
	366	zllmstd=AMAX1(zstd(ii,jj),zllmstd)
	367	zllmsig=AMAX1(zsig(ii,jj),zllmsig)
	368	zllmgam=AMAX1(zgam(ii,jj),zllmgam)
	369	zllmthe=AMAX1(zthe(ii,jj),zllmthe)
	370	zminthe=amin1(zthe(ii,jj),zminthe)
	371	zllmpic=AMAX1(zpic(ii,jj),zllmpic)
	372	zllmval=AMAX1(zval(ii,jj),zllmval)
	373	ENDDO
	374	ENDDO
	375	print *,' MEAN ORO:',zllmmea
	376	print *,' ST. DEV.:',zllmstd
	377	print *,' PENTE:',zllmsig
	378	print *,' ANISOTROP:',zllmgam
	379	print *,' ANGLE:',zminthe,zllmthe
	380	print *,' pic:',zllmpic
	381	print *,' val:',zllmval
	382
	383	C
	384	c gamma and theta a 1. and 0. at poles
	385	c
	386	DO jj=1,jmar
	387	zmea(imar+1,jj)=zmea(1,jj)
[7]	388	zphi(imar+1,jj)=zphi(1,jj)
[2]	389	zpic(imar+1,jj)=zpic(1,jj)
	390	zval(imar+1,jj)=zval(1,jj)
	391	zstd(imar+1,jj)=zstd(1,jj)
	392	zsig(imar+1,jj)=zsig(1,jj)
	393	zgam(imar+1,jj)=zgam(1,jj)
	394	zthe(imar+1,jj)=zthe(1,jj)
	395	ENDDO
	396
	397
	398	zmeanor=0.0
	399	zmeasud=0.0
	400	zstdnor=0.0
	401	zstdsud=0.0
	402	zsignor=0.0
	403	zsigsud=0.0
	404	zweinor=0.0
	405	zweisud=0.0
	406	zpicnor=0.0
	407	zpicsud=0.0
	408	zvalnor=0.0
	409	zvalsud=0.0
	410
	411	DO ii=1,imar
	412	zweinor=zweinor+ weight(ii, 1)
	413	zweisud=zweisud+ weight(ii,jmar)
	414	zmeanor=zmeanor+zmea(ii, 1)*weight(ii, 1)
	415	zmeasud=zmeasud+zmea(ii,jmar)*weight(ii,jmar)
	416	zstdnor=zstdnor+zstd(ii, 1)*weight(ii, 1)
	417	zstdsud=zstdsud+zstd(ii,jmar)*weight(ii,jmar)
	418	zsignor=zsignor+zsig(ii, 1)*weight(ii, 1)
	419	zsigsud=zsigsud+zsig(ii,jmar)*weight(ii,jmar)
	420	zpicnor=zpicnor+zpic(ii, 1)*weight(ii, 1)
	421	zpicsud=zpicsud+zpic(ii,jmar)*weight(ii,jmar)
	422	zvalnor=zvalnor+zval(ii, 1)*weight(ii, 1)
	423	zvalsud=zvalsud+zval(ii,jmar)*weight(ii,jmar)
	424	ENDDO
	425
	426	DO ii=1,imar+1
	427	zmea(ii, 1)=zmeanor/zweinor
	428	zmea(ii,jmar)=zmeasud/zweisud
	429	zphi(ii, 1)=zmeanor/zweinor
	430	zphi(ii,jmar)=zmeasud/zweisud
	431	zpic(ii, 1)=zpicnor/zweinor
	432	zpic(ii,jmar)=zpicsud/zweisud
	433	zval(ii, 1)=zvalnor/zweinor
	434	zval(ii,jmar)=zvalsud/zweisud
	435	zstd(ii, 1)=zstdnor/zweinor
	436	zstd(ii,jmar)=zstdsud/zweisud
	437	zsig(ii, 1)=zsignor/zweinor
	438	zsig(ii,jmar)=zsigsud/zweisud
	439	zgam(ii, 1)=1.
	440	zgam(ii,jmar)=1.
	441	zthe(ii, 1)=0.
	442	zthe(ii,jmar)=0.
	443	ENDDO
	444
	445	RETURN
	446	END
	447
	448	SUBROUTINE MVA9(X,IMAR,JMAR)
	449
	450	C MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
	451
	452	PARAMETER (ISMo=300,JSMo=200)
	453	REAL X(IMAR,JMAR),XF(ISMo,JSMo)
[99]	454	real WEIGHTpb(-1:1,-1:1)
[2]	455
	456	if(imar.gt.ismo) stop'surdimensionner ismo dans mva9 (grid_noro)'
	457	if(jmar.gt.jsmo) stop'surdimensionner jsmo dans mva9 (grid_noro)'
	458
	459	SUM=0.
	460	DO IS=-1,1
[99]	461	DO JS=-1,1
	462	WEIGHTpb(IS,JS)=1./FLOAT((1+IS*2)(1+JS**2))
	463	SUM=SUM+WEIGHTpb(IS,JS)
	464	ENDDO
[2]	465	ENDDO
[99]	466
[325]	467	c WRITE(,) 'MVA9 ', IMAR, JMAR
	468	c WRITE(,) 'MVA9 ', WEIGHTpb
	469	c WRITE(,) 'MVA9 SUM ', SUM
[2]	470	DO IS=-1,1
[99]	471	DO JS=-1,1
	472	WEIGHTpb(IS,JS)=WEIGHTpb(IS,JS)/SUM
	473	ENDDO
[2]	474	ENDDO
	475
	476	DO J=2,JMAR-1
[99]	477	DO I=2,IMAR-1
	478	XF(I,J)=0.
	479	DO IS=-1,1
	480	DO JS=-1,1
	481	XF(I,J)=XF(I,J)+X(I+IS,J+JS)*WEIGHTpb(IS,JS)
	482	ENDDO
	483	ENDDO
	484	ENDDO
[2]	485	ENDDO
	486
	487	DO J=2,JMAR-1
[99]	488	XF(1,J)=0.
	489	IS=IMAR-1
	490	DO JS=-1,1
	491	XF(1,J)=XF(1,J)+X(IS,J+JS)*WEIGHTpb(-1,JS)
	492	ENDDO
	493	DO IS=0,1
	494	DO JS=-1,1
	495	XF(1,J)=XF(1,J)+X(1+IS,J+JS)*WEIGHTpb(IS,JS)
	496	ENDDO
	497	ENDDO
	498	XF(IMAR,J)=XF(1,J)
[2]	499	ENDDO
	500
	501	DO I=1,IMAR
[99]	502	XF(I,1)=XF(I,2)
	503	XF(I,JMAR)=XF(I,JMAR-1)
[2]	504	ENDDO
	505
	506	DO I=1,IMAR
[99]	507	DO J=1,JMAR
	508	X(I,J)=XF(I,J)
	509	ENDDO
[2]	510	ENDDO
	511
	512	RETURN
	513	END
	514
	515
	516

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

Download in other formats:

Original Format