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: 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 mva9_m, only: mva9
    use nr_util, only: assert, pi

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

    ! Correlations of US Navy orography gradients:
    REAL, intent(out):: zphi(:, :) ! orography not smoothed
    real, intent(out):: zmea(:, :) ! smoothed 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, dimension(iim + 1, jjm + 1):: ztz, zxtzx, zytzy, zxtzy, weight

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

    real, dimension(iim + 1, jjm + 1):: mask_tmp, num_tot, num_lan, zmea0
    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) subgrid-scale orographic 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

    zmea0 = zmea ! not smoothed

    ! Filters to smooth out fields for input into subgrid-scale
    ! orographic 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 = merge(1., 0., mask >= 0.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) = zmea0(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	95	use mva9_m, only: mva9
33	guez	39	use nr_util, only: assert, pi
34	guez	3
35			REAL, intent(in):: xdata(:), ydata(:) ! coordinates of input field
36			REAL, intent(in):: zdata(:, :) ! input field
37	guez	78	REAL, intent(in):: x(:), y(:) ! coordinates of output field
38	guez	3
39	guez	70	! Correlations of US Navy orography gradients:
40	guez	95	REAL, intent(out):: zphi(:, :) ! orography not smoothed
41			real, intent(out):: zmea(:, :) ! smoothed orography
42	guez	3	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	78	REAL, dimension(iim + 1, jjm + 1):: ztz, zxtzx, zytzy, zxtzy, weight
62	guez	3
63	guez	70	! Correlations of US Navy orography gradients:
64	guez	68	REAL, dimension(iusn + 2 * iext, jusn + 2):: zxtzxusn, zytzyusn, zxtzyusn
65	guez	3
66	guez	78	real, dimension(iim + 1, jjm + 1):: mask_tmp, num_tot, num_lan, zmea0
67	guez	68	REAL a(iim + 1), b(iim + 1), c(jjm + 1), d(jjm + 1)
68	guez	70	real weighx, weighy, xincr, xk, xp, xm, xw, xq, xl
69	guez	3	real zbordnor, zdeltax, zbordsud, zdeltay, zbordoue, zlenx, zleny, zmeasud
70			real zllmpic, zllmmea, zllmgam, zllmthe, zllmstd, zllmsig, zllmval
71			real zpicnor, zminthe, zsigsud, zstdnor, zstdsud, zvalsud, zvalnor
72			real zweinor, zweisud, zsignor, zpicsud, zmeanor, zbordest
73			integer ii, i, jj, j
74	guez	70	real, parameter:: rad = 6371229.
75	guez	3
76			!-------------------------------
77
78			print *, "Call sequence information: grid_noro"
79
80			call assert((/size(xdata), size(zdata, 1)/) == iusn, "grid_noro iusn")
81			call assert((/size(ydata), size(zdata, 2)/) == jusn, "grid_noro jusn")
82
83			call assert((/size(x), size(zphi, 1), size(zmea, 1), size(zstd, 1), &
84			size(zsig, 1), size(zgam, 1), size(zthe, 1), size(zpic, 1), &
85			size(zval, 1), size(mask, 1)/) == iim + 1, "grid_noro iim")
86
87			call assert((/size(y), size(zphi, 2), size(zmea, 2), size(zstd, 2), &
88			size(zsig, 2), size(zgam, 2), size(zthe, 2), size(zpic, 2), &
89			size(zval, 2), size(mask, 2)/) == jjm + 1, "grid_noro jjm")
90
91			print *, "Paramètres de l'orographie à l'échelle sous-maille"
92			zdeltay = 2. * pi / real(jusn) * rad
93
94	guez	70	! Extension of the US Navy database for computations at boundaries:
95	guez	3
96	guez	68	DO j = 1, jusn
97			yusn(j + 1) = ydata(j)
98			DO i = 1, iusn
99			zusn(i + iext, j + 1) = zdata(i, j)
100			xusn(i + iext) = xdata(i)
101	guez	3	ENDDO
102	guez	68	DO i = 1, iext
103			zusn(i, j + 1) = zdata(iusn - iext + i, j)
104			xusn(i) = xdata(iusn - iext + i) - 2. * pi
105			zusn(iusn + iext + i, j + 1) = zdata(i, j)
106			xusn(iusn + iext + i) = xdata(i) + 2. * pi
107	guez	3	ENDDO
108			ENDDO
109
110	guez	68	yusn(1) = ydata(1) + (ydata(1) - ydata(2))
111			yusn(jusn + 2) = ydata(jusn) + (ydata(jusn) - ydata(jusn - 1))
112			DO i = 1, iusn / 2 + iext
113			zusn(i, 1) = zusn(i + iusn / 2, 2)
114			zusn(i + iusn / 2 + iext, 1) = zusn(i, 2)
115			zusn(i, jusn + 2) = zusn(i + iusn / 2, jusn + 1)
116			zusn(i + iusn / 2 + iext, jusn + 2) = zusn(i, jusn + 1)
117	guez	3	ENDDO
118	guez	15
119	guez	78	! Compute limits of model gridpoint area (regular grid)
120	guez	3
121	guez	68	a(1) = x(1) - (x(2) - x(1)) / 2.0
122			b(1) = (x(1) + x(2)) / 2.0
123	guez	3	DO i = 2, iim
124	guez	68	a(i) = b(i - 1)
125			b(i) = (x(i) + x(i + 1)) / 2.0
126	guez	3	ENDDO
127	guez	68	a(iim + 1) = b(iim)
128			b(iim + 1) = x(iim + 1) + (x(iim + 1) - x(iim)) / 2.0
129	guez	3
130	guez	68	c(1) = y(1) - (y(2) - y(1)) / 2.0
131			d(1) = (y(1) + y(2)) / 2.0
132	guez	3	DO j = 2, jjm
133	guez	68	c(j) = d(j - 1)
134			d(j) = (y(j) + y(j + 1)) / 2.0
135	guez	3	ENDDO
136			c(jjm + 1) = d(jjm)
137	guez	68	d(jjm + 1) = y(jjm + 1) + (y(jjm + 1) - y(jjm)) / 2.0
138	guez	3
139			! Initialisations :
140	guez	68	weight = 0.
141			zxtzx = 0.
142			zytzy = 0.
143			zxtzy = 0.
144			ztz = 0.
145			zmea = 0.
146			zpic = - 1E10
147			zval = 1E10
148	guez	3
149	guez	70	! Compute slopes correlations on US Navy grid
150	guez	3
151	guez	68	zytzyusn = 0.
152			zxtzxusn = 0.
153			zxtzyusn = 0.
154	guez	3
155	guez	68	DO j = 2, jusn + 1
156			zdeltax = zdeltay * cos(yusn(j))
157			DO i = 2, iusn + 2 * iext - 1
158			zytzyusn(i, j) = (zusn(i, j + 1) - zusn(i, j - 1))2 / zdeltay2
159			zxtzxusn(i, j) = (zusn(i + 1, j) - zusn(i - 1, j))2 / zdeltax2
160			zxtzyusn(i, j) = (zusn(i, j + 1) - zusn(i, j - 1)) / zdeltay &
161			* (zusn(i + 1, j) - zusn(i - 1, j)) / zdeltax
162	guez	3	ENDDO
163			ENDDO
164
165	guez	78	! Summation over gridpoint area
166	guez	15
167	guez	68	zleny = pi / real(jusn) * rad
168			xincr = pi / 2. / real(jusn)
169			DO ii = 1, iim + 1
170	guez	3	DO jj = 1, jjm + 1
171	guez	68	num_tot(ii, jj) = 0.
172			num_lan(ii, jj) = 0.
173			DO j = 2, jusn + 1
174			zlenx = zleny * cos(yusn(j))
175			zdeltax = zdeltay * cos(yusn(j))
176			zbordnor = (c(jj) - yusn(j) + xincr) * rad
177			zbordsud = (yusn(j) - d(jj) + xincr) * rad
178	guez	70	weighy = MAX(0., min(zbordnor, zbordsud, zleny))
179	guez	3	IF (weighy /= 0) THEN
180	guez	68	DO i = 2, iusn + 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	guez	70	weighx = MAX(0., min(zbordest, zbordoue, zlenx))
184	guez	3	IF (weighx /= 0) THEN
185			num_tot(ii, jj) = num_tot(ii, jj) + 1.
186			if (zusn(i, j) >= 1.) then
187			num_lan(ii, jj) = num_lan(ii, jj) + 1.
188			end if
189			weight(ii, jj) = weight(ii, jj) + weighx * weighy
190	guez	68	zxtzx(ii, jj) = zxtzx(ii, jj) &
191			+ zxtzxusn(i, j) * weighx * weighy
192			zytzy(ii, jj) = zytzy(ii, jj) &
193			+ zytzyusn(i, j) * weighx * weighy
194			zxtzy(ii, jj) = zxtzy(ii, jj) &
195			+ zxtzyusn(i, j) * weighx * weighy
196	guez	3	ztz(ii, jj) = ztz(ii, jj) &
197			+ zusn(i, j) * zusn(i, j) * weighx * weighy
198			! mean
199	guez	68	zmea(ii, jj) = zmea(ii, jj) + zusn(i, j) * weighx * weighy
200	guez	3	! peacks
201	guez	70	zpic(ii, jj) = max(zpic(ii, jj), zusn(i, j))
202	guez	3	! valleys
203	guez	70	zval(ii, jj) = min(zval(ii, jj), zusn(i, j))
204	guez	3	ENDIF
205			ENDDO
206			ENDIF
207			ENDDO
208			ENDDO
209			ENDDO
210
211	guez	70	if (any(weight == 0.)) then
212			print *, "zero weight in grid_noro"
213			stop 1
214			end if
215	guez	3
216	guez	78	! Compute parameters needed by the Lott & Miller (1997) and Lott
217			! (1999) subgrid-scale orographic scheme.
218	guez	3
219	guez	68	zllmmea = 0.
220			zllmstd = 0.
221			zllmsig = 0.
222			zllmgam = 0.
223			zllmpic = 0.
224			zllmval = 0.
225			zllmthe = 0.
226			zminthe = 0.
227			DO ii = 1, iim + 1
228	guez	3	DO jj = 1, jjm + 1
229	guez	68	mask(ii, jj) = num_lan(ii, jj) / num_tot(ii, jj)
230	guez	78	! Mean orography:
231	guez	68	zmea (ii, jj) = zmea (ii, jj) / weight(ii, jj)
232			zxtzx(ii, jj) = zxtzx(ii, jj) / weight(ii, jj)
233			zytzy(ii, jj) = zytzy(ii, jj) / weight(ii, jj)
234			zxtzy(ii, jj) = zxtzy(ii, jj) / weight(ii, jj)
235			ztz(ii, jj) = ztz(ii, jj) / weight(ii, jj)
236			! Standard deviation:
237			zstd(ii, jj) = sqrt(MAX(0., ztz(ii, jj) - zmea(ii, jj)**2))
238	guez	3	ENDDO
239			ENDDO
240
241	guez	78	! Correct values of horizontal slope near the poles:
242	guez	68	DO ii = 1, iim + 1
243			zxtzx(ii, 1) = zxtzx(ii, 2)
244			zxtzx(ii, jjm + 1) = zxtzx(ii, jjm)
245			zxtzy(ii, 1) = zxtzy(ii, 2)
246			zxtzy(ii, jjm + 1) = zxtzy(ii, jjm)
247			zytzy(ii, 1) = zytzy(ii, 2)
248			zytzy(ii, jjm + 1) = zytzy(ii, jjm)
249	guez	3	ENDDO
250
251	guez	78	zmea0 = zmea ! not smoothed
252	guez	3
253	guez	78	! Filters to smooth out fields for input into subgrid-scale
254			! orographic scheme.
255
256			! First filter, moving average over 9 points.
257	guez	47	CALL MVA9(zmea)
258			CALL MVA9(zstd)
259			CALL MVA9(zpic)
260			CALL MVA9(zval)
261			CALL MVA9(zxtzx)
262			CALL MVA9(zxtzy)
263			CALL MVA9(zytzy)
264	guez	3
265	guez	70	! Masque prenant en compte maximum de terre. On met un seuil à 10
266			! % de terre car en dessous les paramètres de surface n'ont pas de
267	guez	68	! sens.
268	guez	78	mask_tmp = merge(1., 0., mask >= 0.1)
269	guez	3
270			DO ii = 1, iim
271			DO jj = 1, jjm + 1
272	guez	68	! Coefficients K, L et M:
273			xk = (zxtzx(ii, jj) + zytzy(ii, jj)) / 2.
274			xl = (zxtzx(ii, jj) - zytzy(ii, jj)) / 2.
275			xm = zxtzy(ii, jj)
276			xp = xk - sqrt(xl2 + xm2)
277			xq = xk + sqrt(xl2 + xm2)
278			xw = 1e-8
279			if(xp.le.xw) xp = 0.
280			if(xq.le.xw) xq = xw
281			if(abs(xm).le.xw) xm = xw * sign(1., xm)
282			! modification pour masque de terre fractionnaire
283			! slope:
284			zsig(ii, jj) = sqrt(xq) * mask_tmp(ii, jj)
285			! isotropy:
286			zgam(ii, jj) = xp / xq * mask_tmp(ii, jj)
287			! angle theta:
288			zthe(ii, jj) = 57.29577951 * atan2(xm, xl) / 2. * mask_tmp(ii, jj)
289	guez	78	zphi(ii, jj) = zmea0(ii, jj) * mask_tmp(ii, jj)
290	guez	68	zmea(ii, jj) = zmea(ii, jj) * mask_tmp(ii, jj)
291			zpic(ii, jj) = zpic(ii, jj) * mask_tmp(ii, jj)
292			zval(ii, jj) = zval(ii, jj) * mask_tmp(ii, jj)
293			zstd(ii, jj) = zstd(ii, jj) * mask_tmp(ii, jj)
294	guez	70	zllmmea = MAX(zmea(ii, jj), zllmmea)
295			zllmstd = MAX(zstd(ii, jj), zllmstd)
296			zllmsig = MAX(zsig(ii, jj), zllmsig)
297			zllmgam = MAX(zgam(ii, jj), zllmgam)
298			zllmthe = MAX(zthe(ii, jj), zllmthe)
299			zminthe = min(zthe(ii, jj), zminthe)
300			zllmpic = MAX(zpic(ii, jj), zllmpic)
301			zllmval = MAX(zval(ii, jj), zllmval)
302	guez	3	ENDDO
303			ENDDO
304	guez	68
305	guez	15	print *, 'MEAN ORO: ', zllmmea
306			print *, 'ST. DEV.: ', zllmstd
307			print *, 'PENTE: ', zllmsig
308			print *, 'ANISOTROP: ', zllmgam
309			print *, 'ANGLE: ', zminthe, zllmthe
310			print *, 'pic: ', zllmpic
311			print *, 'val: ', zllmval
312	guez	3
313	guez	68	! gamma and theta at 1. and 0. at poles
314			zmea(iim + 1, :) = zmea(1, :)
315			zphi(iim + 1, :) = zphi(1, :)
316			zpic(iim + 1, :) = zpic(1, :)
317			zval(iim + 1, :) = zval(1, :)
318			zstd(iim + 1, :) = zstd(1, :)
319			zsig(iim + 1, :) = zsig(1, :)
320			zgam(iim + 1, :) = zgam(1, :)
321			zthe(iim + 1, :) = zthe(1, :)
322	guez	3
323	guez	68	zmeanor = 0.
324			zmeasud = 0.
325			zstdnor = 0.
326			zstdsud = 0.
327			zsignor = 0.
328			zsigsud = 0.
329			zweinor = 0.
330			zweisud = 0.
331			zpicnor = 0.
332			zpicsud = 0.
333			zvalnor = 0.
334			zvalsud = 0.
335	guez	3
336	guez	68	DO ii = 1, iim
337			zweinor = zweinor + weight(ii, 1)
338			zweisud = zweisud + weight(ii, jjm + 1)
339			zmeanor = zmeanor + zmea(ii, 1) * weight(ii, 1)
340			zmeasud = zmeasud + zmea(ii, jjm + 1) * weight(ii, jjm + 1)
341			zstdnor = zstdnor + zstd(ii, 1) * weight(ii, 1)
342			zstdsud = zstdsud + zstd(ii, jjm + 1) * weight(ii, jjm + 1)
343			zsignor = zsignor + zsig(ii, 1) * weight(ii, 1)
344			zsigsud = zsigsud + zsig(ii, jjm + 1) * weight(ii, jjm + 1)
345			zpicnor = zpicnor + zpic(ii, 1) * weight(ii, 1)
346			zpicsud = zpicsud + zpic(ii, jjm + 1) * weight(ii, jjm + 1)
347			zvalnor = zvalnor + zval(ii, 1) * weight(ii, 1)
348			zvalsud = zvalsud + zval(ii, jjm + 1) * weight(ii, jjm + 1)
349	guez	3	ENDDO
350
351	guez	68	zmea(:, 1) = zmeanor / zweinor
352			zmea(:, jjm + 1) = zmeasud / zweisud
353	guez	3
354	guez	68	zphi(:, 1) = zmeanor / zweinor
355			zphi(:, jjm + 1) = zmeasud / zweisud
356	guez	3
357	guez	68	zpic(:, 1) = zpicnor / zweinor
358			zpic(:, jjm + 1) = zpicsud / zweisud
359	guez	3
360	guez	68	zval(:, 1) = zvalnor / zweinor
361			zval(:, jjm + 1) = zvalsud / zweisud
362	guez	3
363	guez	68	zstd(:, 1) = zstdnor / zweinor
364			zstd(:, jjm + 1) = zstdsud / zweisud
365	guez	3
366	guez	68	zsig(:, 1) = zsignor / zweinor
367			zsig(:, jjm + 1) = zsigsud / zweisud
368	guez	3
369	guez	68	zgam(:, 1) = 1.
370			zgam(:, jjm + 1) = 1.
371	guez	3
372	guez	68	zthe(:, 1) = 0.
373			zthe(:, jjm + 1) = 0.
374	guez	3
375			END SUBROUTINE grid_noro
376
377			end module grid_noro_m