source: LMDZ5/trunk/libf/phylmd/grid_noro_m.F90 @ 2306

Last change on this file since 2306 was 2293, checked in by dcugnet, 9 years ago

Initial states creation routines have been reorganized and simplified.
As far as possible, dynamics and physics related routines have been
separated.
Some routines have been converted to fortran 90 and repeated codes sections
have been "factorized".
Array/vector arguments have become implicit in some routines to avoid usage
of "dimensions.h" ; possible for routines with explicit interfaces and if
iim and jjm can be deduced from arguments sizes.

  • dynlonlat_phylonlat/ce0l.F90 calls now phylmd/etat0phys_netcdf.F90 and dyn3d/etat0dyn_netcdf.F90 that replace phylmd/etat0_netcdf.F90. start.nc and startphy.nc creations are now independant.
  • startvar.F90 has been suppressed ; corresponding operations have been simplified and embedded in etat0*_netcdf.F90 routines as internal procedures.
  • Routines converted to fortran 90 and "factorized":
    • dyn3d_common/conf_dat_m.F90 (replaces dyn3d_common/conf_dat2d.F

and dyn3d_common/conf_dat3d.F)

  • dyn3d/dynredem.F90 (replaces dyn3d/dynredem.F)
  • dyn3d/dynetat0.F90 (replaces dyn3d/dynetat0.F)
  • phylmd/grid_noro_m.F90 (replaces dyn3d_common/grid_noro.F)
  • dynlonlat_phylonlat/grid_atob_m.F90 (replaces dyn3d_common/grid_atob.F)
  • dyn3d_common/caldyn0.F90 (replaces dyn3d_common/caldyn0.F)
  • dyn3d_common/covcont.F90 (replaces dyn3d_common/covcont.F)
  • dyn3d_common/pression.F90 (replaces dyn3d_common/pression.F)
  • phylmd/phyredem.F90 and phylmd/limit_netcdf.F90 have been slightly factorized.

TO DO:

  • little fix needed in grid_noro_m.F90 ; untouched yet to ensure results are exactly the same as before. Unsmoothed orography is used to compute "zphi", but smoothed (should be unsmoothed) one is used at poles.
  • add the dyn3dmem versions of dynredem.F90 and dynetat0.F90 (dynredem_loc.F90 and dynetat0_loc.F90, untested yet).
  • test compilation in parallel mode for a single processor.
File size: 14.7 KB
Line 
1MODULE grid_noro_m
2!
3!*******************************************************************************
4
5  PRIVATE
6  PUBLIC :: grid_noro
7
8
9CONTAINS
10
11
12!-------------------------------------------------------------------------------
13!
14SUBROUTINE grid_noro(xd,yd,zd,x,y,zphi,zmea,zstd,zsig,zgam,zthe,zpic,zval,mask)
15!
16!-------------------------------------------------------------------------------
17! Author: F. Lott (see also Z.X. Li, A. Harzallah et L. Fairhead)
18!-------------------------------------------------------------------------------
19! Purpose: Compute the Parameters of the SSO scheme as described in LOTT &MILLER
20!         (1997) and LOTT(1999).
21!-------------------------------------------------------------------------------
22! Comments:
23!  * Target points are on a rectangular grid:
24!      iim+1 latitudes including North and South Poles;
25!      jjm+1 longitudes, with periodicity jjm+1=1.
26!  * At the poles, the fields value is repeated jjm+1 time.
27!  * The parameters a,b,c,d represent the limits of the target gridpoint region.
28!    The means over this region are calculated from USN data, ponderated by a
29!    weight proportional to the surface occupated by the data inside the model
30!    gridpoint area. In most circumstances, this weight is the ratio between the
31!    surfaces of the USN gridpoint area and the model gridpoint area.
32!
33!           (c)
34!        ----d-----
35!        | . . . .|
36!        |        |
37!     (b)a . * . .b(a)
38!        |        |
39!        | . . . .|
40!        ----c-----
41!           (d)
42!  * Hard-coded US Navy dataset dimensions (imdp=2160 ; jmdp=1080) have been
43!    removed (ALLOCATABLE used).
44!  * iext (currently 10% of imdp) represents the margin to ensure output cells
45!    on the edge are contained in input cells.
46!===============================================================================
47  USE assert_eq_m, ONLY: assert_eq
48  IMPLICIT NONE
49!  include "dimensions.h"
50  include "iniprint.h"
51  REAL, PARAMETER :: epsfra = 1.e-5
52!-------------------------------------------------------------------------------
53! Arguments:
54  REAL, INTENT(IN)  :: xd(:), yd(:)  !--- INPUT  COORDINATES     (imdp) (jmdp)
55  REAL, INTENT(IN)  :: zd(:,:)       !--- INPUT  FIELD           (imdp,jmdp)
56  REAL, INTENT(IN)  :: x(:), y(:)    !--- OUTPUT COORDINATES     (imar+1) (jmar)
57  REAL, INTENT(OUT) :: zphi(:,:)     !--- GEOPOTENTIAL           (imar+1,jmar)
58  REAL, INTENT(OUT) :: zmea(:,:)     !--- MEAN OROGRAPHY         (imar+1,jmar)
59  REAL, INTENT(OUT) :: zstd(:,:)     !--- STANDARD DEVIATION     (imar+1,jmar)
60  REAL, INTENT(OUT) :: zsig(:,:)     !--- SLOPE                  (imar+1,jmar)
61  REAL, INTENT(OUT) :: zgam(:,:)     !--- ANISOTROPY             (imar+1,jmar)
62  REAL, INTENT(OUT) :: zthe(:,:)     !--- SMALL AXIS ORIENTATION (imar+1,jmar)
63  REAL, INTENT(OUT) :: zpic(:,:)     !--- MAXIMUM ALTITITUDE     (imar+1,jmar)
64  REAL, INTENT(OUT) :: zval(:,:)     !--- MINIMUM ALTITITUDE     (imar+1,jmar)
65  REAL, INTENT(OUT) :: mask(:,:)     !--- MASK                   (imar+1,jmar)
66!-------------------------------------------------------------------------------
67! Local variables:
68  CHARACTER(LEN=256) :: modname="grid_noro"
69  REAL, ALLOCATABLE :: xusn(:), yusn(:)           ! dim (imdp+2*iext) (jmdp+2)
70  REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,jmdp+2)
71! CORRELATIONS OF OROGRAPHY GRADIENT              ! dim (imar+1,jmar)
72  REAL, ALLOCATABLE :: ztz(:,:), zxtzx(:,:), zytzy(:,:), zxtzy(:,:), weight(:,:)
73! CORRELATIONS OF USN OROGRAPHY GRADIENTS         ! dim (imar+2*iext,jmdp+2)
74  REAL, ALLOCATABLE :: zxtzxusn(:,:), zytzyusn(:,:), zxtzyusn(:,:)
75  REAL, ALLOCATABLE :: mask_tmp(:,:), zmea0(:,:)  ! dim (imar+1,jmar)
76  REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
77  REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imax)
78  REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmax)
79  LOGICAL :: masque_lu
80  INTEGER :: i, ii, imdp, imar, iext
81  INTEGER :: j, jj, jmdp, jmar, nn
82  REAL    :: xpi, zdeltax, zlenx, weighx, xincr,  zmeanor0
83  REAL    :: rad, zdeltay, zleny, weighy, masque, zmeasud0
84  REAL    :: zbordnor, zmeanor, zstdnor, zsignor, zweinor, zpicnor, zvalnor
85  REAL    :: zbordsud, zmeasud, zstdsud, zsigsud, zweisud, zpicsud, zvalsud
86  REAL    :: zbordest, zbordoue, xk, xl, xm, xp, xq, xw
87!-------------------------------------------------------------------------------
88  imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
89  jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
90  imar=assert_eq([SIZE(x),SIZE(zphi,1),SIZE(zmea,1),SIZE(zstd,1),SIZE(zsig,1), &
91                          SIZE(zgam,1),SIZE(zthe,1),SIZE(zpic,1),SIZE(zval,1), &
92                          SIZE(mask,1)],TRIM(modname)//" imar")-1
93  jmar=assert_eq([SIZE(y),SIZE(zphi,2),SIZE(zmea,2),SIZE(zstd,2),SIZE(zsig,2), &
94                          SIZE(zgam,2),SIZE(zthe,2),SIZE(zpic,2),SIZE(zval,2), &
95                          SIZE(mask,2)],TRIM(modname)//" jmar")
96!  IF(imar/=iim)   CALL abort_gcm(TRIM(modname),'imar/=iim'  ,1)
97!  IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
98  iext=imdp/10                                !--- OK up to 36 degrees cell
99  xpi = ACOS(-1.)
100  rad = 6371229.
101  zdeltay=2.*xpi/REAL(jmdp)*rad
102  WRITE(lunout,*)"*** Orography parameters at sub-cell scale ***"
103
104!--- ARE WE USING A READ MASK ?
105  masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
106  WRITE(lunout,*)'Masque lu: ',masque_lu
107
108!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
109  ALLOCATE(xusn(imdp+2*iext))
110  xusn(1     +iext:imdp  +iext)=xd(:)
111  xusn(1          :       iext)=xd(1+imdp-iext:imdp)-2.*xpi
112  xusn(1+imdp+iext:imdp+2*iext)=xd(1          :iext)+2.*xpi
113
114  ALLOCATE(yusn(jmdp+2))
115  yusn(1       )=yd(1)   +(yd(1)   -yd(2))
116  yusn(2:jmdp+1)=yd(:)
117  yusn(  jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
118
119  ALLOCATE(zusn(imdp+2*iext,jmdp+2))
120  zusn(1       +iext:imdp  +iext,2:jmdp+1)=zd  (:                   ,     :)
121  zusn(1            :       iext,2:jmdp+1)=zd  (imdp-iext+1:imdp    ,     :)
122  zusn(1+imdp  +iext:imdp+2*iext,2:jmdp+1)=zd  (1:iext              ,     :)
123  zusn(1            :imdp/2+iext,       1)=zusn(1+imdp/2:imdp  +iext,     2)
124  zusn(1+imdp/2+iext:imdp+2*iext,       1)=zusn(1       :imdp/2+iext,     2)
125  zusn(1            :imdp/2+iext,  jmdp+2)=zusn(1+imdp/2:imdp  +iext,jmdp+1)
126  zusn(1+imdp/2+iext:imdp+2*iext,  jmdp+2)=zusn(1       :imdp/2+iext,jmdp+1)
127
128!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
129  ALLOCATE(a(imar+1),b(imar+1))
130  b(1:imar)=(x(1:imar  )+ x(2:imar+1))/2.0
131  b(imar+1)= x(  imar+1)+(x(  imar+1)-x(imar))/2.0
132  a(1)=x(1)-(x(2)-x(1))/2.0
133  a(2:imar+1)= b(1:imar)
134
135  ALLOCATE(c(jmar),d(jmar))
136  d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
137  d(  jmar  )= y(  jmar  )+(y(  jmar)-y(jmar-1))/2.0
138  c(1)=y(1)-(y(2)-y(1))/2.0
139  c(2:jmar)=d(1:jmar-1)
140
141!--- INITIALIZATIONS:
142  ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
143  ALLOCATE(zxtzx (imar+1,jmar)); zxtzx (:,:)= 0.0
144  ALLOCATE(zytzy (imar+1,jmar)); zytzy (:,:)= 0.0
145  ALLOCATE(zxtzy (imar+1,jmar)); zxtzy (:,:)= 0.0
146  ALLOCATE(ztz   (imar+1,jmar)); ztz   (:,:)= 0.0
147  zmea(:,:)= 0.0
148  zpic(:,:)=-1.E+10
149  zval(:,:)= 1.E+10
150
151!--- COMPUTE SLOPES CORRELATIONS ON USN GRID
152! CORRELATIONS OF USN OROGRAPHY GRADIENTS  ! dim (imdp+2*iext,jmdp+2)
153  ALLOCATE(zytzyusn(imdp+2*iext,jmdp+2)); zytzyusn(:,:)=0.0
154  ALLOCATE(zxtzxusn(imdp+2*iext,jmdp+2)); zxtzxusn(:,:)=0.0
155  ALLOCATE(zxtzyusn(imdp+2*iext,jmdp+2)); zxtzyusn(:,:)=0.0
156  DO j = 2, jmdp+1
157    zdeltax=zdeltay*cos(yusn(j))
158    DO i = 2, imdp+2*iext-1
159      zytzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))**2/zdeltay**2
160      zxtzxusn(i,j)=(zusn(i+1,j)-zusn(i-1,j))**2/zdeltax**2
161      zxtzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))   /zdeltay  &
162     &             *(zusn(i+1,j)-zusn(i-1,j))   /zdeltax
163    END DO
164  END DO
165
166!--- SUMMATION OVER GRIDPOINT AREA
167  zleny=xpi/REAL(jmdp)*rad
168  xincr=xpi/REAL(jmdp)/2.
169  ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
170  ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
171  DO ii = 1, imar+1
172    DO jj = 1, jmar
173      DO j = 2,jmdp+1
174        zlenx  =zleny  *COS(yusn(j))
175        zdeltax=zdeltay*COS(yusn(j))
176        zbordnor=(xincr+c(jj)-yusn(j))*rad
177        zbordsud=(xincr-d(jj)+yusn(j))*rad
178        weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
179        IF(weighy==0.) CYCLE
180        DO i = 2, imdp+2*iext-1
181          zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j))
182          zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j))
183          weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
184          IF(weighx==0.) CYCLE
185          num_tot(ii,jj)=num_tot(ii,jj)+1.0
186          IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
187          weight(ii,jj)=weight(ii,jj)+weighx*weighy
188          zxtzx(ii,jj)=zxtzx(ii,jj)+zxtzxusn(i,j)*weighx*weighy
189          zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)*weighx*weighy
190          zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)*weighx*weighy
191          ztz  (ii,jj)=  ztz(ii,jj)+zusn(i,j)*zusn(i,j)*weighx*weighy
192          zmea (ii,jj)= zmea(ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
193          zpic (ii,jj)=AMAX1(zpic(ii,jj),zusn(i,j))         !--- PEAKS
194          zval (ii,jj)=AMIN1(zval(ii,jj),zusn(i,j))         !--- VALLEYS
195        END DO
196      END DO
197    END DO
198  END DO
199
200!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
201  IF(.NOT.masque_lu) THEN
202    WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
203  END IF
204  nn=COUNT(weight(:,1:jmar-1)==0.0)
205  IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
206  WHERE(weight(:,:)/=0.0)
207    zmea (:,:)=zmea (:,:)/weight(:,:)
208    zxtzx(:,:)=zxtzx(:,:)/weight(:,:)
209    zytzy(:,:)=zytzy(:,:)/weight(:,:)
210    zxtzy(:,:)=zxtzy(:,:)/weight(:,:)
211    ztz  (:,:)=ztz  (:,:)/weight(:,:)
212    zstd (:,:)=ztz  (:,:)-zmea(:,:)**2
213  END WHERE
214  WHERE(zstd(:,:)<0) zstd(:,:)=0.
215  zstd (:,:)=SQRT(zstd(:,:))
216
217!--- CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:
218  zxtzx(:,   1)=zxtzx(:,     2)
219  zxtzx(:,jmar)=zxtzx(:,jmar-1)
220  zxtzy(:,   1)=zxtzy(:,     2)
221  zxtzy(:,jmar)=zxtzy(:,jmar-1)
222  zytzy(:,   1)=zytzy(:,     2)
223  zytzy(:,jmar)=zytzy(:,jmar-1)
224
225!=== FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.
226!--- FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.
227!-------------------------------------------------------------------------------
228  ALLOCATE(zmea0(imar+1,jmar))
229  zmea0(:,:)=zmea(:,:)                           ! GK211005 (CG) UNSMOOTHED TOPO
230  CALL MVA9(zmea);  CALL MVA9(zstd);  CALL MVA9(zpic);  CALL MVA9(zval)
231  CALL MVA9(zxtzx); CALL MVA9(zxtzy); CALL MVA9(zytzy)
232
233!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD. (SURFACE PARAMS MEANINGLESS)
234  ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
235  WHERE(mask>=0.1) mask_tmp = 1.
236  WHERE(weight(:,:)/=0.0)
237!   zphi (:,:)= mask_tmp(:,:)*zmea (:,:) ! GK211005 (CG) not necessarly smoothed
238    zphi (:,:)= mask_tmp(:,:)*zmea0(:,:)
239    zmea0(:,:)= mask_tmp(:,:)*zmea0(:,:)
240    zmea (:,:)= mask_tmp(:,:)*zmea (:,:)
241    zpic (:,:)= mask_tmp(:,:)*zpic (:,:)
242    zval (:,:)= mask_tmp(:,:)*zval (:,:)
243    zstd (:,:)= mask_tmp(:,:)*zstd (:,:)
244  END WHERE
245  DO ii = 1, imar
246    DO jj = 1, jmar
247      IF (weight(ii,jj)/=0.0) THEN
248      !--- Coefficients K, L et M:
249        xk=(zxtzx(ii,jj)+zytzy(ii,jj))/2.
250        xl=(zxtzx(ii,jj)-zytzy(ii,jj))/2.
251        xm=zxtzy(ii,jj)
252        xp=xk-SQRT(xl**2+xm**2)
253        xq=xk+SQRT(xl**2+xm**2)
254        xw=1.e-8
255        IF(xp<=xw) xp=0.
256        IF(xq<=xw) xq=xw
257        IF(ABS(xm)<=xw) xm=xw*SIGN(1.,xm)
258      !--- SLOPE
259        zsig(ii,jj)=SQRT(xq)*mask_tmp(ii,jj)
260      !---ISOTROPY
261        zgam(ii,jj)=xp/xq*mask_tmp(ii,jj)
262      !--- THETA ANGLE
263        zthe(ii,jj)=57.29577951*ATAN2(xm,xl)/2.*mask_tmp(ii,jj)
264      END IF
265    END DO
266  END DO
267  WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
268  WRITE(lunout,*)'  ST. DEV.:' ,MAXVAL(zstd)
269  WRITE(lunout,*)'  PENTE:'    ,MAXVAL(zsig)
270  WRITE(lunout,*)'  ANISOTROP:',MAXVAL(zgam)
271  WRITE(lunout,*)'  ANGLE:'    ,MINVAL(zthe),MAXVAL(zthe)
272  WRITE(lunout,*)'  pic:'      ,MAXVAL(zpic)
273  WRITE(lunout,*)'  val:'      ,MAXVAL(zval)
274     
275!--- Values at poles
276  zmea0(imar+1,:)=zmea0(1,:)
277  zmea (imar+1,:)=zmea (1,:)
278  zphi (imar+1,:)=zphi (1,:)
279  zpic (imar+1,:)=zpic (1,:)
280  zval (imar+1,:)=zval (1,:)
281  zstd (imar+1,:)=zstd (1,:)
282  zsig (imar+1,:)=zsig (1,:)
283  zgam (imar+1,:)=zgam (1,:)
284  zthe (imar+1,:)=zthe (1,:)
285
286  zweinor =SUM(weight(1:imar,   1),DIM=1)
287  zweisud =SUM(weight(1:imar,jmar),DIM=1)
288  zmeanor0=SUM(weight(1:imar,   1)*zmea0(1:imar,   1),DIM=1)
289  zmeasud0=SUM(weight(1:imar,jmar)*zmea0(1:imar,jmar),DIM=1)
290  zmeanor =SUM(weight(1:imar,   1)*zmea (1:imar,   1),DIM=1)
291  zmeasud =SUM(weight(1:imar,jmar)*zmea (1:imar,jmar),DIM=1)
292  zstdnor =SUM(weight(1:imar,   1)*zstd (1:imar,   1),DIM=1)
293  zstdsud =SUM(weight(1:imar,jmar)*zstd (1:imar,jmar),DIM=1)
294  zsignor =SUM(weight(1:imar,   1)*zsig (1:imar,   1),DIM=1)
295  zsigsud =SUM(weight(1:imar,jmar)*zsig (1:imar,jmar),DIM=1)
296  zpicnor =SUM(weight(1:imar,   1)*zpic (1:imar,   1),DIM=1)
297  zpicsud =SUM(weight(1:imar,jmar)*zpic (1:imar,jmar),DIM=1)
298  zvalnor =SUM(weight(1:imar,   1)*zval (1:imar,   1),DIM=1)
299  zvalsud =SUM(weight(1:imar,jmar)*zval (1:imar,jmar),DIM=1)
300
301  zmea(:,1)=zmeanor /zweinor; zmea(:,jmar)=zmeasud /zweisud
302!  zphi(:,1)=zmeanor0/zweinor; zphi(:,jmar)=zmeasud0/zweisud   TO COMMIT
303  zphi(:,1)=zmeanor /zweinor; zphi(:,jmar)=zmeasud /zweisud
304  zpic(:,1)=zpicnor /zweinor; zpic(:,jmar)=zpicsud /zweisud
305  zval(:,1)=zvalnor /zweinor; zval(:,jmar)=zvalsud /zweisud
306  zstd(:,1)=zstdnor /zweinor; zstd(:,jmar)=zstdsud /zweisud
307  zsig(:,1)=zsignor /zweinor; zsig(:,jmar)=zsigsud /zweisud
308  zgam(:,1)=1.;               zgam(:,jmar)=1.
309  zthe(:,1)=0.;               zthe(:,jmar)=0.
310
311END SUBROUTINE grid_noro
312!
313!-------------------------------------------------------------------------------
314
315
316!-------------------------------------------------------------------------------
317!
318SUBROUTINE MVA9(x)
319!
320!-------------------------------------------------------------------------------
321  IMPLICIT NONE
322! MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
323!-------------------------------------------------------------------------------
324! Arguments:
325  REAL, INTENT(INOUT) :: x(:,:)
326!-------------------------------------------------------------------------------
327! Local variables:
328  REAL    :: xf(SIZE(x,DIM=1),SIZE(x,DIM=2)), WEIGHTpb(-1:1,-1:1)
329  INTEGER :: i, j, imar, jmar
330!-------------------------------------------------------------------------------
331  WEIGHTpb=RESHAPE([((1./REAL((1+i**2)*(1+j**2)),i=-1,1),j=-1,1)],SHAPE=[3,3])
332  WEIGHTpb=WEIGHTpb/SUM(WEIGHTpb)
333  imar=SIZE(X,DIM=1); jmar=SIZE(X,DIM=2)
334  DO j=2,jmar-1
335    DO i=2,imar-1
336      xf(i,j)=SUM(x(i-1:i+1,j-1:j+1)*WEIGHTpb(:,:))
337    END DO
338  END DO
339  DO j=2,jmar-1
340    xf(1,j)=SUM(x(imar-1,j-1:j+1)*WEIGHTpb(-1,:))
341    xf(1,j)=xf(1,j)+SUM(x(1:2,j-1:j+1)*WEIGHTpb(0:1,-1:1))
342    xf(imar,j)=xf(1,j)
343  END DO
344  xf(:,   1)=xf(:,     2)
345  xf(:,jmar)=xf(:,jmar-1)
346  x(:,:)=xf(:,:)
347
348END SUBROUTINE MVA9
349!
350!-------------------------------------------------------------------------------
351
352
353END MODULE grid_noro_m
354
Note: See TracBrowser for help on using the repository browser.