source: LMDZ5/trunk/libf/dyn3d_common/grid_noro.F @ 2210

Last change on this file since 2210 was 2199, checked in by Ehouarn Millour, 10 years ago

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