phylmd/Orography/grid_noro_m.f

module grid_noro_m

  implicit none

contains

  SUBROUTINE grid_noro(xdata, ydata, zdata, x, y, zphi, zmea, zstd, zsig, &
       zgam, zthe, zpic, zval, mask)

    ! From dyn3d/grid_noro.F, version 1.1.1.1 2004/05/19 12:53:06

    ! Authors: F. Lott, Z. X. Li, A. Harzallah and L. Fairhead

    ! Compute the parameters of the SSO scheme as described in
    ! Lott and Miller (1997) and Lott (1999).
    ! Target points are on a rectangular grid:
    ! jjm + 1 latitudes including North and South Poles;
    ! iim + 1 longitudes, with periodicity: longitude(iim + 1) = longitude(1)
    ! At the poles the field value is repeated iim + 1 times.

    ! The parameters a, b, c, d represent the limite of the target
    ! gridpoint region. The means over this region are calculated from
    ! USN data, ponderated by a weight proportional to the surface
    ! occupied by the data inside the model gridpoint area. In most
    ! circumstances, this weight is the ratio between the surface of
    ! the USN gridpoint area and the surface of the model gridpoint
    ! area. See "grid_noto.txt".

    use dimens_m, only: iim, jjm
    use nr_util, only: assert, pi
    use mva9_m, only: mva9

    REAL, intent(in):: xdata(:), ydata(:) ! coordinates of input field
    REAL, intent(in):: zdata(:, :) ! input field
    REAL, intent(in):: x(:), y(:) ! ccordinates output field

    ! Correlations of USN orography gradients:

    REAL, intent(out):: zphi(:, :)
    real, intent(out):: zmea(:, :) ! Mean orography
    real, intent(out):: zstd(:, :) ! Standard deviation
    REAL zsig(:, :) ! Slope
    real zgam(:, :) ! Anisotropy
    real zthe(:, :) ! Orientation of the small axis
    REAL, intent(out):: zpic(:, :) ! Maximum altitude
    real, intent(out):: zval(:, :) ! Minimum altitude

    real, intent(out):: mask(:, :) ! fraction of land

    ! Variables local to the procedure:

    ! In this version it is assumed that the input data come from
    ! the US Navy dataset:
    integer, parameter:: iusn = 2160, jusn = 1080
    integer, parameter:: iext = 216
    REAL xusn(iusn + 2 * iext), yusn(jusn + 2)
    REAL zusn(iusn + 2 * iext, jusn + 2)

    ! Intermediate fields (correlations of orography gradient)

    REAL ztz(iim + 1, jjm + 1), zxtzx(iim + 1, jjm + 1)
    REAL zytzy(iim + 1, jjm + 1), zxtzy(iim + 1, jjm + 1)
    REAL weight(iim + 1, jjm + 1)

    ! Correlations of USN orography gradients:
    REAL, dimension(iusn + 2 * iext, jusn + 2):: zxtzxusn, zytzyusn, zxtzyusn

    real mask_tmp(size(x), size(y))
    real num_tot(iim + 1, jjm + 1), num_lan(iim + 1, jjm + 1)

    REAL a(iim + 1), b(iim + 1), c(jjm + 1), d(jjm + 1)
    real rad, weighx, weighy, xincr, xk, xp, xm, xw, xq, xl
    real zbordnor, zdeltax, zbordsud, zdeltay, zbordoue, zlenx, zleny, zmeasud
    real zllmpic, zllmmea, zllmgam, zllmthe, zllmstd, zllmsig, zllmval
    real zpicnor, zminthe, zsigsud, zstdnor, zstdsud, zvalsud, zvalnor
    real zweinor, zweisud, zsignor, zpicsud, zmeanor, zbordest
    integer ii, i, jj, j

    !-------------------------------

    print *, "Call sequence information: grid_noro"

    call assert((/size(xdata), size(zdata, 1)/) == iusn, "grid_noro iusn")
    call assert((/size(ydata), size(zdata, 2)/) == jusn, "grid_noro jusn")

    call assert((/size(x), size(zphi, 1), size(zmea, 1), size(zstd, 1), &
         size(zsig, 1), size(zgam, 1), size(zthe, 1), size(zpic, 1), &
         size(zval, 1), size(mask, 1)/) == iim + 1, "grid_noro iim")

    call assert((/size(y), size(zphi, 2), size(zmea, 2), size(zstd, 2), &
         size(zsig, 2), size(zgam, 2), size(zthe, 2), size(zpic, 2), &
         size(zval, 2), size(mask, 2)/) == jjm + 1, "grid_noro jjm")

    print *, "Paramètres de l'orographie à l'échelle sous-maille" 
    rad = 6371229.
    zdeltay = 2. * pi / real(jusn) * rad

    ! Extension of the USN database to POCEED computations at boundaries:

    DO j = 1, jusn
       yusn(j + 1) = ydata(j)
       DO i = 1, iusn
          zusn(i + iext, j + 1) = zdata(i, j)
          xusn(i + iext) = xdata(i)
       ENDDO
       DO i = 1, iext
          zusn(i, j + 1) = zdata(iusn - iext + i, j)
          xusn(i) = xdata(iusn - iext + i) - 2. * pi
          zusn(iusn + iext + i, j + 1) = zdata(i, j)
          xusn(iusn + iext + i) = xdata(i) + 2. * pi
       ENDDO
    ENDDO

    yusn(1) = ydata(1) + (ydata(1) - ydata(2))
    yusn(jusn + 2) = ydata(jusn) + (ydata(jusn) - ydata(jusn - 1))
    DO i = 1, iusn / 2 + iext
       zusn(i, 1) = zusn(i + iusn / 2, 2)
       zusn(i + iusn / 2 + iext, 1) = zusn(i, 2)
       zusn(i, jusn + 2) = zusn(i + iusn / 2, jusn + 1)
       zusn(i + iusn / 2 + iext, jusn + 2) = zusn(i, jusn + 1)
    ENDDO

    ! COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)

    a(1) = x(1) - (x(2) - x(1)) / 2.0
    b(1) = (x(1) + x(2)) / 2.0
    DO i = 2, iim
       a(i) = b(i - 1)
       b(i) = (x(i) + x(i + 1)) / 2.0
    ENDDO
    a(iim + 1) = b(iim)
    b(iim + 1) = x(iim + 1) + (x(iim + 1) - x(iim)) / 2.0

    c(1) = y(1) - (y(2) - y(1)) / 2.0
    d(1) = (y(1) + y(2)) / 2.0
    DO j = 2, jjm
       c(j) = d(j - 1)
       d(j) = (y(j) + y(j + 1)) / 2.0
    ENDDO
    c(jjm + 1) = d(jjm)
    d(jjm + 1) = y(jjm + 1) + (y(jjm + 1) - y(jjm)) / 2.0

    ! Initialisations :
    weight = 0.
    zxtzx = 0.
    zytzy = 0.
    zxtzy = 0.
    ztz = 0.
    zmea = 0.
    zpic = - 1E10
    zval = 1E10

    ! COMPUTE SLOPES CORRELATIONS ON USN GRID

    zytzyusn = 0.
    zxtzxusn = 0.
    zxtzyusn = 0.

    DO j = 2, jusn + 1 
       zdeltax = zdeltay * cos(yusn(j))
       DO i = 2, iusn + 2 * iext - 1
          zytzyusn(i, j) = (zusn(i, j + 1) - zusn(i, j - 1))**2 / zdeltay**2
          zxtzxusn(i, j) = (zusn(i + 1, j) - zusn(i - 1, j))**2 / zdeltax**2
          zxtzyusn(i, j) = (zusn(i, j + 1) - zusn(i, j - 1)) / zdeltay &
               * (zusn(i + 1, j) - zusn(i - 1, j)) / zdeltax
       ENDDO
    ENDDO

    ! SUMMATION OVER GRIDPOINT AREA

    zleny = pi / real(jusn) * rad
    xincr = pi / 2. / real(jusn)
    DO ii = 1, iim + 1
       DO jj = 1, jjm + 1
          num_tot(ii, jj) = 0.
          num_lan(ii, jj) = 0.
          DO j = 2, jusn + 1 
             zlenx = zleny * cos(yusn(j))
             zdeltax = zdeltay * cos(yusn(j))
             zbordnor = (c(jj) - yusn(j) + xincr) * rad
             zbordsud = (yusn(j) - d(jj) + xincr) * rad
             weighy = AMAX1(0., amin1(zbordnor, zbordsud, zleny))
             IF (weighy /= 0) THEN
                DO i = 2, iusn + 2 * iext - 1
                   zbordest = (xusn(i) - a(ii) + xincr) * rad * cos(yusn(j))
                   zbordoue = (b(ii) + xincr - xusn(i)) * rad * cos(yusn(j))
                   weighx = AMAX1(0., amin1(zbordest, zbordoue, zlenx))
                   IF (weighx /= 0) THEN
                      num_tot(ii, jj) = num_tot(ii, jj) + 1.
                      if (zusn(i, j) >= 1.) then
                         num_lan(ii, jj) = num_lan(ii, jj) + 1.
                      end if
                      weight(ii, jj) = weight(ii, jj) + weighx * weighy
                      zxtzx(ii, jj) = zxtzx(ii, jj) &
                           + zxtzxusn(i, j) * weighx * weighy
                      zytzy(ii, jj) = zytzy(ii, jj) &
                           + zytzyusn(i, j) * weighx * weighy
                      zxtzy(ii, jj) = zxtzy(ii, jj) &
                           + zxtzyusn(i, j) * weighx * weighy
                      ztz(ii, jj) = ztz(ii, jj) &
                           + zusn(i, j) * zusn(i, j) * weighx * weighy
                      ! mean
                      zmea(ii, jj) = zmea(ii, jj) + zusn(i, j) * weighx * weighy
                      ! peacks
                      zpic(ii, jj) = amax1(zpic(ii, jj), zusn(i, j))
                      ! valleys
                      zval(ii, jj) = amin1(zval(ii, jj), zusn(i, j))
                   ENDIF
                ENDDO
             ENDIF
          ENDDO
       ENDDO
    ENDDO

    if (any(weight == 0.)) stop "zero weight in grid_noro"

    ! COMPUTE PARAMETERS NEEDED BY THE LOTT & MILLER (1997) AND
    ! LOTT (1999) SSO SCHEME.

    zllmmea = 0.
    zllmstd = 0.
    zllmsig = 0.
    zllmgam = 0.
    zllmpic = 0.
    zllmval = 0.
    zllmthe = 0.
    zminthe = 0.
    DO ii = 1, iim + 1
       DO jj = 1, jjm + 1
          mask(ii, jj) = num_lan(ii, jj) / num_tot(ii, jj)
          ! Mean Orography:
          zmea (ii, jj) = zmea (ii, jj) / weight(ii, jj)
          zxtzx(ii, jj) = zxtzx(ii, jj) / weight(ii, jj)
          zytzy(ii, jj) = zytzy(ii, jj) / weight(ii, jj)
          zxtzy(ii, jj) = zxtzy(ii, jj) / weight(ii, jj)
          ztz(ii, jj) = ztz(ii, jj) / weight(ii, jj)
          ! Standard deviation:
          zstd(ii, jj) = sqrt(MAX(0., ztz(ii, jj) - zmea(ii, jj)**2))
       ENDDO
    ENDDO

    ! CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:
    DO ii = 1, iim + 1
       zxtzx(ii, 1) = zxtzx(ii, 2)
       zxtzx(ii, jjm + 1) = zxtzx(ii, jjm)
       zxtzy(ii, 1) = zxtzy(ii, 2)
       zxtzy(ii, jjm + 1) = zxtzy(ii, jjm)
       zytzy(ii, 1) = zytzy(ii, 2)
       zytzy(ii, jjm + 1) = zytzy(ii, jjm)
    ENDDO

    ! FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.

    ! FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.
    CALL MVA9(zmea)
    CALL MVA9(zstd)
    CALL MVA9(zpic)
    CALL MVA9(zval)
    CALL MVA9(zxtzx)
    CALL MVA9(zxtzy) 
    CALL MVA9(zytzy)

    ! Masque prenant en compte maximum de terre. On seuille à 10 % de
    ! terre car en dessous les paramètres de surface n'ont pas de
    ! sens.
    mask_tmp = 0.
    WHERE (mask >= 0.1) mask_tmp = 1.

    DO ii = 1, iim
       DO jj = 1, jjm + 1
          ! Coefficients K, L et M:
          xk = (zxtzx(ii, jj) + zytzy(ii, jj)) / 2.
          xl = (zxtzx(ii, jj) - zytzy(ii, jj)) / 2.
          xm = zxtzy(ii, jj)
          xp = xk - sqrt(xl**2 + xm**2)
          xq = xk + sqrt(xl**2 + xm**2)
          xw = 1e-8
          if(xp.le.xw) xp = 0.
          if(xq.le.xw) xq = xw
          if(abs(xm).le.xw) xm = xw * sign(1., xm)
          ! modification pour masque de terre fractionnaire
          ! slope: 
          zsig(ii, jj) = sqrt(xq) * mask_tmp(ii, jj)
          ! isotropy:
          zgam(ii, jj) = xp / xq * mask_tmp(ii, jj)
          ! angle theta:
          zthe(ii, jj) = 57.29577951 * atan2(xm, xl) / 2. * mask_tmp(ii, jj)
          zphi(ii, jj) = zmea(ii, jj) * mask_tmp(ii, jj)
          zmea(ii, jj) = zmea(ii, jj) * mask_tmp(ii, jj)
          zpic(ii, jj) = zpic(ii, jj) * mask_tmp(ii, jj)
          zval(ii, jj) = zval(ii, jj) * mask_tmp(ii, jj)
          zstd(ii, jj) = zstd(ii, jj) * mask_tmp(ii, jj)
          zllmmea = AMAX1(zmea(ii, jj), zllmmea)
          zllmstd = AMAX1(zstd(ii, jj), zllmstd)
          zllmsig = AMAX1(zsig(ii, jj), zllmsig)
          zllmgam = AMAX1(zgam(ii, jj), zllmgam)
          zllmthe = AMAX1(zthe(ii, jj), zllmthe)
          zminthe = amin1(zthe(ii, jj), zminthe)
          zllmpic = AMAX1(zpic(ii, jj), zllmpic)
          zllmval = AMAX1(zval(ii, jj), zllmval)
       ENDDO
    ENDDO

    print *, 'MEAN ORO: ', zllmmea
    print *, 'ST. DEV.: ', zllmstd
    print *, 'PENTE: ', zllmsig
    print *, 'ANISOTROP: ', zllmgam
    print *, 'ANGLE: ', zminthe, zllmthe
    print *, 'pic: ', zllmpic
    print *, 'val: ', zllmval

    ! gamma and theta at 1. and 0. at poles
    zmea(iim + 1, :) = zmea(1, :)
    zphi(iim + 1, :) = zphi(1, :)
    zpic(iim + 1, :) = zpic(1, :)
    zval(iim + 1, :) = zval(1, :)
    zstd(iim + 1, :) = zstd(1, :)
    zsig(iim + 1, :) = zsig(1, :)
    zgam(iim + 1, :) = zgam(1, :)
    zthe(iim + 1, :) = zthe(1, :)

    zmeanor = 0.
    zmeasud = 0.
    zstdnor = 0.
    zstdsud = 0.
    zsignor = 0.
    zsigsud = 0.
    zweinor = 0.
    zweisud = 0.
    zpicnor = 0.
    zpicsud = 0. 
    zvalnor = 0.
    zvalsud = 0. 

    DO ii = 1, iim
       zweinor = zweinor + weight(ii, 1)
       zweisud = zweisud + weight(ii, jjm + 1)
       zmeanor = zmeanor + zmea(ii, 1) * weight(ii, 1)
       zmeasud = zmeasud + zmea(ii, jjm + 1) * weight(ii, jjm + 1)
       zstdnor = zstdnor + zstd(ii, 1) * weight(ii, 1)
       zstdsud = zstdsud + zstd(ii, jjm + 1) * weight(ii, jjm + 1)
       zsignor = zsignor + zsig(ii, 1) * weight(ii, 1)
       zsigsud = zsigsud + zsig(ii, jjm + 1) * weight(ii, jjm + 1)
       zpicnor = zpicnor + zpic(ii, 1) * weight(ii, 1)
       zpicsud = zpicsud + zpic(ii, jjm + 1) * weight(ii, jjm + 1)
       zvalnor = zvalnor + zval(ii, 1) * weight(ii, 1)
       zvalsud = zvalsud + zval(ii, jjm + 1) * weight(ii, jjm + 1)
    ENDDO

    zmea(:, 1) = zmeanor / zweinor
    zmea(:, jjm + 1) = zmeasud / zweisud

    zphi(:, 1) = zmeanor / zweinor
    zphi(:, jjm + 1) = zmeasud / zweisud

    zpic(:, 1) = zpicnor / zweinor
    zpic(:, jjm + 1) = zpicsud / zweisud

    zval(:, 1) = zvalnor / zweinor
    zval(:, jjm + 1) = zvalsud / zweisud

    zstd(:, 1) = zstdnor / zweinor
    zstd(:, jjm + 1) = zstdsud / zweisud

    zsig(:, 1) = zsignor / zweinor
    zsig(:, jjm + 1) = zsigsud / zweisud

    zgam(:, 1) = 1.
    zgam(:, jjm + 1) = 1.

    zthe(:, 1) = 0.
    zthe(:, jjm + 1) = 0.

  END SUBROUTINE grid_noro

end module grid_noro_m
1	guez	3	module grid_noro_m
2
3			implicit none
4
5			contains
6
7			SUBROUTINE grid_noro(xdata, ydata, zdata, x, y, zphi, zmea, zstd, zsig, &
8			zgam, zthe, zpic, zval, mask)
9
10			! From dyn3d/grid_noro.F, version 1.1.1.1 2004/05/19 12:53:06
11
12			! Authors: F. Lott, Z. X. Li, A. Harzallah and L. Fairhead
13
14			! Compute the parameters of the SSO scheme as described in
15			! Lott and Miller (1997) and Lott (1999).
16			! Target points are on a rectangular grid:
17	guez	68	! jjm + 1 latitudes including North and South Poles;
18			! iim + 1 longitudes, with periodicity: longitude(iim + 1) = longitude(1)
19			! At the poles the field value is repeated iim + 1 times.
20	guez	3
21			! The parameters a, b, c, d represent the limite of the target
22	guez	68	! gridpoint region. The means over this region are calculated from
23			! USN data, ponderated by a weight proportional to the surface
24			! occupied by the data inside the model gridpoint area. In most
25			! circumstances, this weight is the ratio between the surface of
26			! the USN gridpoint area and the surface of the model gridpoint
27			! area. See "grid_noto.txt".
28	guez	3
29			use dimens_m, only: iim, jjm
30	guez	39	use nr_util, only: assert, pi
31	guez	47	use mva9_m, only: mva9
32	guez	3
33			REAL, intent(in):: xdata(:), ydata(:) ! coordinates of input field
34			REAL, intent(in):: zdata(:, :) ! input field
35			REAL, intent(in):: x(:), y(:) ! ccordinates output field
36
37			! Correlations of USN orography gradients:
38
39	guez	68	REAL, intent(out):: zphi(:, :)
40	guez	3	real, intent(out):: zmea(:, :) ! Mean orography
41			real, intent(out):: zstd(:, :) ! Standard deviation
42			REAL zsig(:, :) ! Slope
43			real zgam(:, :) ! Anisotropy
44			real zthe(:, :) ! Orientation of the small axis
45			REAL, intent(out):: zpic(:, :) ! Maximum altitude
46			real, intent(out):: zval(:, :) ! Minimum altitude
47
48	guez	15	real, intent(out):: mask(:, :) ! fraction of land
49	guez	3
50			! Variables local to the procedure:
51
52			! In this version it is assumed that the input data come from
53			! the US Navy dataset:
54	guez	68	integer, parameter:: iusn = 2160, jusn = 1080
55			integer, parameter:: iext = 216
56			REAL xusn(iusn + 2 * iext), yusn(jusn + 2)
57			REAL zusn(iusn + 2 * iext, jusn + 2)
58	guez	3
59	guez	68	! Intermediate fields (correlations of orography gradient)
60	guez	3
61	guez	68	REAL ztz(iim + 1, jjm + 1), zxtzx(iim + 1, jjm + 1)
62			REAL zytzy(iim + 1, jjm + 1), zxtzy(iim + 1, jjm + 1)
63			REAL weight(iim + 1, jjm + 1)
64	guez	3
65			! Correlations of USN orography gradients:
66	guez	68	REAL, dimension(iusn + 2 * iext, jusn + 2):: zxtzxusn, zytzyusn, zxtzyusn
67	guez	3
68			real mask_tmp(size(x), size(y))
69	guez	68	real num_tot(iim + 1, jjm + 1), num_lan(iim + 1, jjm + 1)
70	guez	3
71	guez	68	REAL a(iim + 1), b(iim + 1), c(jjm + 1), d(jjm + 1)
72	guez	3	real rad, weighx, weighy, xincr, xk, xp, xm, xw, xq, xl
73			real zbordnor, zdeltax, zbordsud, zdeltay, zbordoue, zlenx, zleny, zmeasud
74			real zllmpic, zllmmea, zllmgam, zllmthe, zllmstd, zllmsig, zllmval
75			real zpicnor, zminthe, zsigsud, zstdnor, zstdsud, zvalsud, zvalnor
76			real zweinor, zweisud, zsignor, zpicsud, zmeanor, zbordest
77			integer ii, i, jj, j
78
79			!-------------------------------
80
81			print *, "Call sequence information: grid_noro"
82
83			call assert((/size(xdata), size(zdata, 1)/) == iusn, "grid_noro iusn")
84			call assert((/size(ydata), size(zdata, 2)/) == jusn, "grid_noro jusn")
85
86			call assert((/size(x), size(zphi, 1), size(zmea, 1), size(zstd, 1), &
87			size(zsig, 1), size(zgam, 1), size(zthe, 1), size(zpic, 1), &
88			size(zval, 1), size(mask, 1)/) == iim + 1, "grid_noro iim")
89
90			call assert((/size(y), size(zphi, 2), size(zmea, 2), size(zstd, 2), &
91			size(zsig, 2), size(zgam, 2), size(zthe, 2), size(zpic, 2), &
92			size(zval, 2), size(mask, 2)/) == jjm + 1, "grid_noro jjm")
93
94			print *, "Paramètres de l'orographie à l'échelle sous-maille"
95			rad = 6371229.
96			zdeltay = 2. * pi / real(jusn) * rad
97
98			! Extension of the USN database to POCEED computations at boundaries:
99
100	guez	68	DO j = 1, jusn
101			yusn(j + 1) = ydata(j)
102			DO i = 1, iusn
103			zusn(i + iext, j + 1) = zdata(i, j)
104			xusn(i + iext) = xdata(i)
105	guez	3	ENDDO
106	guez	68	DO i = 1, iext
107			zusn(i, j + 1) = zdata(iusn - iext + i, j)
108			xusn(i) = xdata(iusn - iext + i) - 2. * pi
109			zusn(iusn + iext + i, j + 1) = zdata(i, j)
110			xusn(iusn + iext + i) = xdata(i) + 2. * pi
111	guez	3	ENDDO
112			ENDDO
113
114	guez	68	yusn(1) = ydata(1) + (ydata(1) - ydata(2))
115			yusn(jusn + 2) = ydata(jusn) + (ydata(jusn) - ydata(jusn - 1))
116			DO i = 1, iusn / 2 + iext
117			zusn(i, 1) = zusn(i + iusn / 2, 2)
118			zusn(i + iusn / 2 + iext, 1) = zusn(i, 2)
119			zusn(i, jusn + 2) = zusn(i + iusn / 2, jusn + 1)
120			zusn(i + iusn / 2 + iext, jusn + 2) = zusn(i, jusn + 1)
121	guez	3	ENDDO
122	guez	15
123	guez	68	! COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
124	guez	3
125	guez	68	a(1) = x(1) - (x(2) - x(1)) / 2.0
126			b(1) = (x(1) + x(2)) / 2.0
127	guez	3	DO i = 2, iim
128	guez	68	a(i) = b(i - 1)
129			b(i) = (x(i) + x(i + 1)) / 2.0
130	guez	3	ENDDO
131	guez	68	a(iim + 1) = b(iim)
132			b(iim + 1) = x(iim + 1) + (x(iim + 1) - x(iim)) / 2.0
133	guez	3
134	guez	68	c(1) = y(1) - (y(2) - y(1)) / 2.0
135			d(1) = (y(1) + y(2)) / 2.0
136	guez	3	DO j = 2, jjm
137	guez	68	c(j) = d(j - 1)
138			d(j) = (y(j) + y(j + 1)) / 2.0
139	guez	3	ENDDO
140			c(jjm + 1) = d(jjm)
141	guez	68	d(jjm + 1) = y(jjm + 1) + (y(jjm + 1) - y(jjm)) / 2.0
142	guez	3
143			! Initialisations :
144	guez	68	weight = 0.
145			zxtzx = 0.
146			zytzy = 0.
147			zxtzy = 0.
148			ztz = 0.
149			zmea = 0.
150			zpic = - 1E10
151			zval = 1E10
152	guez	3
153	guez	68	! COMPUTE SLOPES CORRELATIONS ON USN GRID
154	guez	3
155	guez	68	zytzyusn = 0.
156			zxtzxusn = 0.
157			zxtzyusn = 0.
158	guez	3
159	guez	68	DO j = 2, jusn + 1
160			zdeltax = zdeltay * cos(yusn(j))
161			DO i = 2, iusn + 2 * iext - 1
162			zytzyusn(i, j) = (zusn(i, j + 1) - zusn(i, j - 1))2 / zdeltay2
163			zxtzxusn(i, j) = (zusn(i + 1, j) - zusn(i - 1, j))2 / zdeltax2
164			zxtzyusn(i, j) = (zusn(i, j + 1) - zusn(i, j - 1)) / zdeltay &
165			* (zusn(i + 1, j) - zusn(i - 1, j)) / zdeltax
166	guez	3	ENDDO
167			ENDDO
168
169	guez	68	! SUMMATION OVER GRIDPOINT AREA
170	guez	15
171	guez	68	zleny = pi / real(jusn) * rad
172			xincr = pi / 2. / real(jusn)
173			DO ii = 1, iim + 1
174	guez	3	DO jj = 1, jjm + 1
175	guez	68	num_tot(ii, jj) = 0.
176			num_lan(ii, jj) = 0.
177			DO j = 2, jusn + 1
178			zlenx = zleny * cos(yusn(j))
179			zdeltax = zdeltay * cos(yusn(j))
180			zbordnor = (c(jj) - yusn(j) + xincr) * rad
181			zbordsud = (yusn(j) - d(jj) + xincr) * rad
182			weighy = AMAX1(0., amin1(zbordnor, zbordsud, zleny))
183	guez	3	IF (weighy /= 0) THEN
184	guez	68	DO i = 2, iusn + 2 * iext - 1
185			zbordest = (xusn(i) - a(ii) + xincr) * rad * cos(yusn(j))
186			zbordoue = (b(ii) + xincr - xusn(i)) * rad * cos(yusn(j))
187			weighx = AMAX1(0., amin1(zbordest, zbordoue, zlenx))
188	guez	3	IF (weighx /= 0) THEN
189			num_tot(ii, jj) = num_tot(ii, jj) + 1.
190			if (zusn(i, j) >= 1.) then
191			num_lan(ii, jj) = num_lan(ii, jj) + 1.
192			end if
193			weight(ii, jj) = weight(ii, jj) + weighx * weighy
194	guez	68	zxtzx(ii, jj) = zxtzx(ii, jj) &
195			+ zxtzxusn(i, j) * weighx * weighy
196			zytzy(ii, jj) = zytzy(ii, jj) &
197			+ zytzyusn(i, j) * weighx * weighy
198			zxtzy(ii, jj) = zxtzy(ii, jj) &
199			+ zxtzyusn(i, j) * weighx * weighy
200	guez	3	ztz(ii, jj) = ztz(ii, jj) &
201			+ zusn(i, j) * zusn(i, j) * weighx * weighy
202			! mean
203	guez	68	zmea(ii, jj) = zmea(ii, jj) + zusn(i, j) * weighx * weighy
204	guez	3	! peacks
205	guez	68	zpic(ii, jj) = amax1(zpic(ii, jj), zusn(i, j))
206	guez	3	! valleys
207	guez	68	zval(ii, jj) = amin1(zval(ii, jj), zusn(i, j))
208	guez	3	ENDIF
209			ENDDO
210			ENDIF
211			ENDDO
212			ENDDO
213			ENDDO
214
215			if (any(weight == 0.)) stop "zero weight in grid_noro"
216
217	guez	68	! COMPUTE PARAMETERS NEEDED BY THE LOTT & MILLER (1997) AND
218			! LOTT (1999) SSO SCHEME.
219	guez	3
220	guez	68	zllmmea = 0.
221			zllmstd = 0.
222			zllmsig = 0.
223			zllmgam = 0.
224			zllmpic = 0.
225			zllmval = 0.
226			zllmthe = 0.
227			zminthe = 0.
228			DO ii = 1, iim + 1
229	guez	3	DO jj = 1, jjm + 1
230	guez	68	mask(ii, jj) = num_lan(ii, jj) / num_tot(ii, jj)
231			! Mean Orography:
232			zmea (ii, jj) = zmea (ii, jj) / weight(ii, jj)
233			zxtzx(ii, jj) = zxtzx(ii, jj) / weight(ii, jj)
234			zytzy(ii, jj) = zytzy(ii, jj) / weight(ii, jj)
235			zxtzy(ii, jj) = zxtzy(ii, jj) / weight(ii, jj)
236			ztz(ii, jj) = ztz(ii, jj) / weight(ii, jj)
237			! Standard deviation:
238			zstd(ii, jj) = sqrt(MAX(0., ztz(ii, jj) - zmea(ii, jj)**2))
239	guez	3	ENDDO
240			ENDDO
241
242			! CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:
243	guez	68	DO ii = 1, iim + 1
244			zxtzx(ii, 1) = zxtzx(ii, 2)
245			zxtzx(ii, jjm + 1) = zxtzx(ii, jjm)
246			zxtzy(ii, 1) = zxtzy(ii, 2)
247			zxtzy(ii, jjm + 1) = zxtzy(ii, jjm)
248			zytzy(ii, 1) = zytzy(ii, 2)
249			zytzy(ii, jjm + 1) = zytzy(ii, jjm)
250	guez	3	ENDDO
251
252	guez	68	! FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.
253	guez	3
254	guez	68	! FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.
255	guez	47	CALL MVA9(zmea)
256			CALL MVA9(zstd)
257			CALL MVA9(zpic)
258			CALL MVA9(zval)
259			CALL MVA9(zxtzx)
260			CALL MVA9(zxtzy)
261			CALL MVA9(zytzy)
262	guez	3
263	guez	68	! Masque prenant en compte maximum de terre. On seuille à 10 % de
264			! terre car en dessous les paramètres de surface n'ont pas de
265			! sens.
266			mask_tmp = 0.
267	guez	3	WHERE (mask >= 0.1) mask_tmp = 1.
268
269			DO ii = 1, iim
270			DO jj = 1, jjm + 1
271	guez	68	! Coefficients K, L et M:
272			xk = (zxtzx(ii, jj) + zytzy(ii, jj)) / 2.
273			xl = (zxtzx(ii, jj) - zytzy(ii, jj)) / 2.
274			xm = zxtzy(ii, jj)
275			xp = xk - sqrt(xl2 + xm2)
276			xq = xk + sqrt(xl2 + xm2)
277			xw = 1e-8
278			if(xp.le.xw) xp = 0.
279			if(xq.le.xw) xq = xw
280			if(abs(xm).le.xw) xm = xw * sign(1., xm)
281			! modification pour masque de terre fractionnaire
282			! slope:
283			zsig(ii, jj) = sqrt(xq) * mask_tmp(ii, jj)
284			! isotropy:
285			zgam(ii, jj) = xp / xq * mask_tmp(ii, jj)
286			! angle theta:
287			zthe(ii, jj) = 57.29577951 * atan2(xm, xl) / 2. * mask_tmp(ii, jj)
288			zphi(ii, jj) = zmea(ii, jj) * mask_tmp(ii, jj)
289			zmea(ii, jj) = zmea(ii, jj) * mask_tmp(ii, jj)
290			zpic(ii, jj) = zpic(ii, jj) * mask_tmp(ii, jj)
291			zval(ii, jj) = zval(ii, jj) * mask_tmp(ii, jj)
292			zstd(ii, jj) = zstd(ii, jj) * mask_tmp(ii, jj)
293			zllmmea = AMAX1(zmea(ii, jj), zllmmea)
294			zllmstd = AMAX1(zstd(ii, jj), zllmstd)
295			zllmsig = AMAX1(zsig(ii, jj), zllmsig)
296			zllmgam = AMAX1(zgam(ii, jj), zllmgam)
297			zllmthe = AMAX1(zthe(ii, jj), zllmthe)
298			zminthe = amin1(zthe(ii, jj), zminthe)
299			zllmpic = AMAX1(zpic(ii, jj), zllmpic)
300			zllmval = AMAX1(zval(ii, jj), zllmval)
301	guez	3	ENDDO
302			ENDDO
303	guez	68
304	guez	15	print *, 'MEAN ORO: ', zllmmea
305			print *, 'ST. DEV.: ', zllmstd
306			print *, 'PENTE: ', zllmsig
307			print *, 'ANISOTROP: ', zllmgam
308			print *, 'ANGLE: ', zminthe, zllmthe
309			print *, 'pic: ', zllmpic
310			print *, 'val: ', zllmval
311	guez	3
312	guez	68	! gamma and theta at 1. and 0. at poles
313			zmea(iim + 1, :) = zmea(1, :)
314			zphi(iim + 1, :) = zphi(1, :)
315			zpic(iim + 1, :) = zpic(1, :)
316			zval(iim + 1, :) = zval(1, :)
317			zstd(iim + 1, :) = zstd(1, :)
318			zsig(iim + 1, :) = zsig(1, :)
319			zgam(iim + 1, :) = zgam(1, :)
320			zthe(iim + 1, :) = zthe(1, :)
321	guez	3
322	guez	68	zmeanor = 0.
323			zmeasud = 0.
324			zstdnor = 0.
325			zstdsud = 0.
326			zsignor = 0.
327			zsigsud = 0.
328			zweinor = 0.
329			zweisud = 0.
330			zpicnor = 0.
331			zpicsud = 0.
332			zvalnor = 0.
333			zvalsud = 0.
334	guez	3
335	guez	68	DO ii = 1, iim
336			zweinor = zweinor + weight(ii, 1)
337			zweisud = zweisud + weight(ii, jjm + 1)
338			zmeanor = zmeanor + zmea(ii, 1) * weight(ii, 1)
339			zmeasud = zmeasud + zmea(ii, jjm + 1) * weight(ii, jjm + 1)
340			zstdnor = zstdnor + zstd(ii, 1) * weight(ii, 1)
341			zstdsud = zstdsud + zstd(ii, jjm + 1) * weight(ii, jjm + 1)
342			zsignor = zsignor + zsig(ii, 1) * weight(ii, 1)
343			zsigsud = zsigsud + zsig(ii, jjm + 1) * weight(ii, jjm + 1)
344			zpicnor = zpicnor + zpic(ii, 1) * weight(ii, 1)
345			zpicsud = zpicsud + zpic(ii, jjm + 1) * weight(ii, jjm + 1)
346			zvalnor = zvalnor + zval(ii, 1) * weight(ii, 1)
347			zvalsud = zvalsud + zval(ii, jjm + 1) * weight(ii, jjm + 1)
348	guez	3	ENDDO
349
350	guez	68	zmea(:, 1) = zmeanor / zweinor
351			zmea(:, jjm + 1) = zmeasud / zweisud
352	guez	3
353	guez	68	zphi(:, 1) = zmeanor / zweinor
354			zphi(:, jjm + 1) = zmeasud / zweisud
355	guez	3
356	guez	68	zpic(:, 1) = zpicnor / zweinor
357			zpic(:, jjm + 1) = zpicsud / zweisud
358	guez	3
359	guez	68	zval(:, 1) = zvalnor / zweinor
360			zval(:, jjm + 1) = zvalsud / zweisud
361	guez	3
362	guez	68	zstd(:, 1) = zstdnor / zweinor
363			zstd(:, jjm + 1) = zstdsud / zweisud
364	guez	3
365	guez	68	zsig(:, 1) = zsignor / zweinor
366			zsig(:, jjm + 1) = zsigsud / zweisud
367	guez	3
368	guez	68	zgam(:, 1) = 1.
369			zgam(:, jjm + 1) = 1.
370	guez	3
371	guez	68	zthe(:, 1) = 0.
372			zthe(:, jjm + 1) = 0.
373	guez	3
374			END SUBROUTINE grid_noro
375
376			end module grid_noro_m