phylmd/Orography/grid_noro.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: François Lott, Laurent Li, A. Harzallah and Laurent
    ! Fairhead

    ! Compute the parameters of the sub-grid scale orography 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
    ! US Navy 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 US Navy 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 US Navy orography gradients:
    REAL, intent(out):: zphi(:, :)
    real, intent(out):: zmea(:, :) ! Mean orography
    real, intent(out):: zstd(:, :) ! Standard deviation
    REAL, intent(out):: zsig(:, :) ! Slope
    real, intent(out):: zgam(:, :) ! Anisotropy
    real, intent(out):: 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 US Navy 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 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
    real, parameter:: rad = 6371229.

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

    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" 
    zdeltay = 2. * pi / real(jusn) * rad

    ! Extension of the US Navy database for 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 US Navy 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 = MAX(0., min(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 = MAX(0., min(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) = max(zpic(ii, jj), zusn(i, j))
                      ! valleys
                      zval(ii, jj) = min(zval(ii, jj), zusn(i, j))
                   ENDIF
                ENDDO
             ENDIF
          ENDDO
       ENDDO
    ENDDO

    if (any(weight == 0.)) then
       print *, "zero weight in grid_noro"
       stop 1
    end if

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