Context Navigation

grid_noro.F @ 937

Last change on this file since 937 was 773, checked in by Laurent Fairhead, 17 years ago
Suite du merge entre la version et la HEAD: quelques modifications de Yann sur le LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 15.5 KB

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

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

Download in other formats:

Original Format