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\c{}cois 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 dimensions, only: iim, jjm
    use mva9_m, only: mva9
    use nr_util, only: assert, pi

    ! Coordinates of input field:
    REAL, intent(in):: xdata(:) ! (iusn)
    REAL, intent(in):: ydata(:) ! (jusn)

    REAL, intent(in):: zdata(:, :) ! (iusn, jusn) input field, in m
    REAL, intent(in):: x(:), y(:) ! coordinates of output field

    ! Correlations of US Navy orography gradients:

    REAL, intent(out):: zphi(:, :) ! (iim + 1, jjm + 1)
    ! geoptential height of orography, not smoothed, in m
    
    real, intent(out):: zmea(:, :) ! (iim + 1, jjm + 1) smoothed orography
    real, intent(out):: zstd(:, :) ! (iim + 1, jjm + 1) Standard deviation
    REAL, intent(out):: zsig(:, :) ! (iim + 1, jjm + 1) Slope
    real, intent(out):: zgam(:, :) ! (iim + 1, jjm + 1) Anisotropy

    real, intent(out):: zthe(:, :) ! (iim + 1, jjm + 1)
    ! Orientation of the small axis

    REAL, intent(out):: zpic(:, :) ! (iim + 1, jjm + 1) Maximum altitude
    real, intent(out):: zval(:, :) ! (iim + 1, jjm + 1) Minimum altitude

    real, intent(out):: mask(:, :) ! (iim + 1, jjm + 1) fraction of land

    ! Local:

    ! 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) ! in m

    ! 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
    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 zweinor, zweisud, 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 *, "Parameters of subgrid-scale orography" 
    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.

    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

    ! Masque prenant en compte maximum de terre. On met un seuil \`a 10
    ! % de terre car en dessous les param\`etres de surface n'ont pas de
    ! sens.
    mask_tmp = merge(1., 0., mask >= 0.1)

    zphi(:iim, :) = zmea(:iim, :) * mask_tmp(:iim, :)
    ! (zmea is not yet 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)

    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 <= xw) xp = 0.
          if (xq <= xw) xq = xw
          if (abs(xm) <= 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)
       ENDDO
    ENDDO

    zmea(:iim, :) = zmea(:iim, :) * mask_tmp(:iim, :)
    zpic(:iim, :) = zpic(:iim, :) * mask_tmp(:iim, :)
    zval(:iim, :) = zval(:iim, :) * mask_tmp(:iim, :)
    zstd(:iim, :) = zstd(:iim, :) * mask_tmp(:iim, :)

    print *, 'MEAN ORO: ', MAXVAL(zmea(:iim, :))
    print *, 'ST. DEV.: ', MAXVAL(zstd(:iim, :))
    print *, 'PENTE: ', MAXVAL(zsig(:iim, :))
    print *, 'ANISOTROP: ', MAXVAL(zgam(:iim, :))
    print *, 'ANGLE: ', minval(zthe(:iim, :)), MAXVAL(zthe(:iim, :))
    print *, 'pic: ', MAXVAL(zpic(:iim, :))
    print *, 'val: ', MAXVAL(zval(:iim, :))

    ! 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, :)

    zweinor = sum(weight(:iim, 1))
    zweisud = sum(weight(:iim, jjm + 1))
    zmeanor = sum(zmea(:iim, 1) * weight(:iim, 1))
    zmeasud = sum(zmea(:iim, jjm + 1) * weight(:iim, jjm + 1))

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

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

    zpic(:, 1) = sum(zpic(:iim, 1) * weight(:iim, 1)) / zweinor
    zpic(:, jjm + 1) = sum(zpic(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
         / zweisud

    zval(:, 1) = sum(zval(:iim, 1) * weight(:iim, 1)) / zweinor
    zval(:, jjm + 1) = sum(zval(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
         / zweisud

    zstd(:, 1) = sum(zstd(:iim, 1) * weight(:iim, 1)) / zweinor
    zstd(:, jjm + 1) = sum(zstd(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
         / zweisud

    zsig(:, 1) = sum(zsig(:iim, 1) * weight(:iim, 1)) / zweinor
    zsig(:, jjm + 1) = sum(zsig(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
         / 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	147	! Authors: Fran\c{}cois Lott, Laurent Li, A. Harzallah and Laurent
13	guez	70	! 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	guez	265	use dimensions, only: iim, jjm
32	guez	95	use mva9_m, only: mva9
33	guez	39	use nr_util, only: assert, pi
34	guez	3
35	guez	147	! Coordinates of input field:
36			REAL, intent(in):: xdata(:) ! (iusn)
37			REAL, intent(in):: ydata(:) ! (jusn)
38
39	guez	227	REAL, intent(in):: zdata(:, :) ! (iusn, jusn) input field, in m
40	guez	78	REAL, intent(in):: x(:), y(:) ! coordinates of output field
41	guez	3
42	guez	70	! Correlations of US Navy orography gradients:
43	guez	227
44			REAL, intent(out):: zphi(:, :) ! (iim + 1, jjm + 1)
45			! geoptential height of orography, not smoothed, in m
46
47	guez	147	real, intent(out):: zmea(:, :) ! (iim + 1, jjm + 1) smoothed orography
48			real, intent(out):: zstd(:, :) ! (iim + 1, jjm + 1) Standard deviation
49			REAL, intent(out):: zsig(:, :) ! (iim + 1, jjm + 1) Slope
50			real, intent(out):: zgam(:, :) ! (iim + 1, jjm + 1) Anisotropy
51	guez	3
52	guez	147	real, intent(out):: zthe(:, :) ! (iim + 1, jjm + 1)
53			! Orientation of the small axis
54	guez	3
55	guez	147	REAL, intent(out):: zpic(:, :) ! (iim + 1, jjm + 1) Maximum altitude
56			real, intent(out):: zval(:, :) ! (iim + 1, jjm + 1) Minimum altitude
57
58			real, intent(out):: mask(:, :) ! (iim + 1, jjm + 1) fraction of land
59
60	guez	227	! Local:
61	guez	3
62			! In this version it is assumed that the input data come from
63			! the US Navy dataset:
64	guez	68	integer, parameter:: iusn = 2160, jusn = 1080
65			integer, parameter:: iext = 216
66			REAL xusn(iusn + 2 * iext), yusn(jusn + 2)
67	guez	227	REAL zusn(iusn + 2 * iext, jusn + 2) ! in m
68	guez	3
69	guez	68	! Intermediate fields (correlations of orography gradient)
70	guez	78	REAL, dimension(iim + 1, jjm + 1):: ztz, zxtzx, zytzy, zxtzy, weight
71	guez	3
72	guez	70	! Correlations of US Navy orography gradients:
73	guez	68	REAL, dimension(iusn + 2 * iext, jusn + 2):: zxtzxusn, zytzyusn, zxtzyusn
74	guez	3
75	guez	147	real, dimension(iim + 1, jjm + 1):: mask_tmp, num_tot, num_lan
76	guez	68	REAL a(iim + 1), b(iim + 1), c(jjm + 1), d(jjm + 1)
77	guez	70	real weighx, weighy, xincr, xk, xp, xm, xw, xq, xl
78	guez	3	real zbordnor, zdeltax, zbordsud, zdeltay, zbordoue, zlenx, zleny, zmeasud
79	guez	147	real zweinor, zweisud, zmeanor, zbordest
80	guez	3	integer ii, i, jj, j
81	guez	70	real, parameter:: rad = 6371229.
82	guez	3
83			!-------------------------------
84
85			print *, "Call sequence information: grid_noro"
86
87			call assert((/size(xdata), size(zdata, 1)/) == iusn, "grid_noro iusn")
88			call assert((/size(ydata), size(zdata, 2)/) == jusn, "grid_noro jusn")
89
90			call assert((/size(x), size(zphi, 1), size(zmea, 1), size(zstd, 1), &
91			size(zsig, 1), size(zgam, 1), size(zthe, 1), size(zpic, 1), &
92			size(zval, 1), size(mask, 1)/) == iim + 1, "grid_noro iim")
93
94			call assert((/size(y), size(zphi, 2), size(zmea, 2), size(zstd, 2), &
95			size(zsig, 2), size(zgam, 2), size(zthe, 2), size(zpic, 2), &
96			size(zval, 2), size(mask, 2)/) == jjm + 1, "grid_noro jjm")
97
98	guez	147	print *, "Parameters of subgrid-scale orography"
99	guez	3	zdeltay = 2. * pi / real(jusn) * rad
100
101	guez	70	! Extension of the US Navy database for computations at boundaries:
102	guez	3
103	guez	68	DO j = 1, jusn
104			yusn(j + 1) = ydata(j)
105			DO i = 1, iusn
106			zusn(i + iext, j + 1) = zdata(i, j)
107			xusn(i + iext) = xdata(i)
108	guez	3	ENDDO
109	guez	68	DO i = 1, iext
110			zusn(i, j + 1) = zdata(iusn - iext + i, j)
111			xusn(i) = xdata(iusn - iext + i) - 2. * pi
112			zusn(iusn + iext + i, j + 1) = zdata(i, j)
113			xusn(iusn + iext + i) = xdata(i) + 2. * pi
114	guez	3	ENDDO
115			ENDDO
116
117	guez	68	yusn(1) = ydata(1) + (ydata(1) - ydata(2))
118			yusn(jusn + 2) = ydata(jusn) + (ydata(jusn) - ydata(jusn - 1))
119			DO i = 1, iusn / 2 + iext
120			zusn(i, 1) = zusn(i + iusn / 2, 2)
121			zusn(i + iusn / 2 + iext, 1) = zusn(i, 2)
122			zusn(i, jusn + 2) = zusn(i + iusn / 2, jusn + 1)
123			zusn(i + iusn / 2 + iext, jusn + 2) = zusn(i, jusn + 1)
124	guez	3	ENDDO
125	guez	15
126	guez	78	! Compute limits of model gridpoint area (regular grid)
127	guez	3
128	guez	68	a(1) = x(1) - (x(2) - x(1)) / 2.0
129			b(1) = (x(1) + x(2)) / 2.0
130	guez	3	DO i = 2, iim
131	guez	68	a(i) = b(i - 1)
132			b(i) = (x(i) + x(i + 1)) / 2.0
133	guez	3	ENDDO
134	guez	68	a(iim + 1) = b(iim)
135			b(iim + 1) = x(iim + 1) + (x(iim + 1) - x(iim)) / 2.0
136	guez	3
137	guez	68	c(1) = y(1) - (y(2) - y(1)) / 2.0
138			d(1) = (y(1) + y(2)) / 2.0
139	guez	3	DO j = 2, jjm
140	guez	68	c(j) = d(j - 1)
141			d(j) = (y(j) + y(j + 1)) / 2.0
142	guez	3	ENDDO
143			c(jjm + 1) = d(jjm)
144	guez	68	d(jjm + 1) = y(jjm + 1) + (y(jjm + 1) - y(jjm)) / 2.0
145	guez	3
146			! Initialisations :
147	guez	68	weight = 0.
148			zxtzx = 0.
149			zytzy = 0.
150			zxtzy = 0.
151			ztz = 0.
152			zmea = 0.
153			zpic = - 1E10
154			zval = 1E10
155	guez	3
156	guez	70	! Compute slopes correlations on US Navy grid
157	guez	3
158	guez	68	zytzyusn = 0.
159			zxtzxusn = 0.
160			zxtzyusn = 0.
161	guez	3
162	guez	68	DO j = 2, jusn + 1
163			zdeltax = zdeltay * cos(yusn(j))
164			DO i = 2, iusn + 2 * iext - 1
165			zytzyusn(i, j) = (zusn(i, j + 1) - zusn(i, j - 1))2 / zdeltay2
166			zxtzxusn(i, j) = (zusn(i + 1, j) - zusn(i - 1, j))2 / zdeltax2
167			zxtzyusn(i, j) = (zusn(i, j + 1) - zusn(i, j - 1)) / zdeltay &
168			* (zusn(i + 1, j) - zusn(i - 1, j)) / zdeltax
169	guez	3	ENDDO
170			ENDDO
171
172	guez	78	! Summation over gridpoint area
173	guez	15
174	guez	68	zleny = pi / real(jusn) * rad
175			xincr = pi / 2. / real(jusn)
176			DO ii = 1, iim + 1
177	guez	3	DO jj = 1, jjm + 1
178	guez	68	num_tot(ii, jj) = 0.
179			num_lan(ii, jj) = 0.
180			DO j = 2, jusn + 1
181			zlenx = zleny * cos(yusn(j))
182			zdeltax = zdeltay * cos(yusn(j))
183			zbordnor = (c(jj) - yusn(j) + xincr) * rad
184			zbordsud = (yusn(j) - d(jj) + xincr) * rad
185	guez	70	weighy = MAX(0., min(zbordnor, zbordsud, zleny))
186	guez	3	IF (weighy /= 0) THEN
187	guez	68	DO i = 2, iusn + 2 * iext - 1
188			zbordest = (xusn(i) - a(ii) + xincr) * rad * cos(yusn(j))
189			zbordoue = (b(ii) + xincr - xusn(i)) * rad * cos(yusn(j))
190	guez	70	weighx = MAX(0., min(zbordest, zbordoue, zlenx))
191	guez	3	IF (weighx /= 0) THEN
192			num_tot(ii, jj) = num_tot(ii, jj) + 1.
193			if (zusn(i, j) >= 1.) then
194			num_lan(ii, jj) = num_lan(ii, jj) + 1.
195			end if
196			weight(ii, jj) = weight(ii, jj) + weighx * weighy
197	guez	68	zxtzx(ii, jj) = zxtzx(ii, jj) &
198			+ zxtzxusn(i, j) * weighx * weighy
199			zytzy(ii, jj) = zytzy(ii, jj) &
200			+ zytzyusn(i, j) * weighx * weighy
201			zxtzy(ii, jj) = zxtzy(ii, jj) &
202			+ zxtzyusn(i, j) * weighx * weighy
203	guez	3	ztz(ii, jj) = ztz(ii, jj) &
204			+ zusn(i, j) * zusn(i, j) * weighx * weighy
205			! mean
206	guez	68	zmea(ii, jj) = zmea(ii, jj) + zusn(i, j) * weighx * weighy
207	guez	3	! peacks
208	guez	70	zpic(ii, jj) = max(zpic(ii, jj), zusn(i, j))
209	guez	3	! valleys
210	guez	70	zval(ii, jj) = min(zval(ii, jj), zusn(i, j))
211	guez	3	ENDIF
212			ENDDO
213			ENDIF
214			ENDDO
215			ENDDO
216			ENDDO
217
218	guez	70	if (any(weight == 0.)) then
219			print *, "zero weight in grid_noro"
220			stop 1
221			end if
222	guez	3
223	guez	78	! Compute parameters needed by the Lott & Miller (1997) and Lott
224			! (1999) subgrid-scale orographic scheme.
225	guez	3
226	guez	68	DO ii = 1, iim + 1
227	guez	3	DO jj = 1, jjm + 1
228	guez	68	mask(ii, jj) = num_lan(ii, jj) / num_tot(ii, jj)
229	guez	78	! Mean orography:
230	guez	68	zmea (ii, jj) = zmea (ii, jj) / weight(ii, jj)
231			zxtzx(ii, jj) = zxtzx(ii, jj) / weight(ii, jj)
232			zytzy(ii, jj) = zytzy(ii, jj) / weight(ii, jj)
233			zxtzy(ii, jj) = zxtzy(ii, jj) / weight(ii, jj)
234			ztz(ii, jj) = ztz(ii, jj) / weight(ii, jj)
235			! Standard deviation:
236			zstd(ii, jj) = sqrt(MAX(0., ztz(ii, jj) - zmea(ii, jj)**2))
237	guez	3	ENDDO
238			ENDDO
239
240	guez	78	! Correct values of horizontal slope near the poles:
241	guez	68	DO ii = 1, iim + 1
242			zxtzx(ii, 1) = zxtzx(ii, 2)
243			zxtzx(ii, jjm + 1) = zxtzx(ii, jjm)
244			zxtzy(ii, 1) = zxtzy(ii, 2)
245			zxtzy(ii, jjm + 1) = zxtzy(ii, jjm)
246			zytzy(ii, 1) = zytzy(ii, 2)
247			zytzy(ii, jjm + 1) = zytzy(ii, jjm)
248	guez	3	ENDDO
249
250	guez	147	! Masque prenant en compte maximum de terre. On met un seuil \`a 10
251			! % de terre car en dessous les param\`etres de surface n'ont pas de
252			! sens.
253			mask_tmp = merge(1., 0., mask >= 0.1)
254	guez	3
255	guez	147	zphi(:iim, :) = zmea(:iim, :) * mask_tmp(:iim, :)
256			! (zmea is not yet smoothed)
257
258	guez	78	! Filters to smooth out fields for input into subgrid-scale
259			! orographic scheme.
260
261			! First filter, moving average over 9 points.
262	guez	47	CALL MVA9(zmea)
263			CALL MVA9(zstd)
264			CALL MVA9(zpic)
265			CALL MVA9(zval)
266			CALL MVA9(zxtzx)
267			CALL MVA9(zxtzy)
268			CALL MVA9(zytzy)
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	guez	147	if (xp <= xw) xp = 0.
280			if (xq <= xw) xq = xw
281			if (abs(xm) <= xw) xm = xw * sign(1., xm)
282	guez	68	! 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	3	ENDDO
290			ENDDO
291	guez	68
292	guez	147	zmea(:iim, :) = zmea(:iim, :) * mask_tmp(:iim, :)
293			zpic(:iim, :) = zpic(:iim, :) * mask_tmp(:iim, :)
294			zval(:iim, :) = zval(:iim, :) * mask_tmp(:iim, :)
295			zstd(:iim, :) = zstd(:iim, :) * mask_tmp(:iim, :)
296	guez	3
297	guez	147	print *, 'MEAN ORO: ', MAXVAL(zmea(:iim, :))
298			print *, 'ST. DEV.: ', MAXVAL(zstd(:iim, :))
299			print *, 'PENTE: ', MAXVAL(zsig(:iim, :))
300			print *, 'ANISOTROP: ', MAXVAL(zgam(:iim, :))
301			print *, 'ANGLE: ', minval(zthe(:iim, :)), MAXVAL(zthe(:iim, :))
302			print *, 'pic: ', MAXVAL(zpic(:iim, :))
303			print *, 'val: ', MAXVAL(zval(:iim, :))
304
305	guez	68	! gamma and theta at 1. and 0. at poles
306			zmea(iim + 1, :) = zmea(1, :)
307			zphi(iim + 1, :) = zphi(1, :)
308			zpic(iim + 1, :) = zpic(1, :)
309			zval(iim + 1, :) = zval(1, :)
310			zstd(iim + 1, :) = zstd(1, :)
311			zsig(iim + 1, :) = zsig(1, :)
312			zgam(iim + 1, :) = zgam(1, :)
313			zthe(iim + 1, :) = zthe(1, :)
314	guez	3
315	guez	147	zweinor = sum(weight(:iim, 1))
316			zweisud = sum(weight(:iim, jjm + 1))
317			zmeanor = sum(zmea(:iim, 1) * weight(:iim, 1))
318			zmeasud = sum(zmea(:iim, jjm + 1) * weight(:iim, jjm + 1))
319	guez	3
320	guez	68	zmea(:, 1) = zmeanor / zweinor
321			zmea(:, jjm + 1) = zmeasud / zweisud
322	guez	3
323	guez	68	zphi(:, 1) = zmeanor / zweinor
324			zphi(:, jjm + 1) = zmeasud / zweisud
325	guez	3
326	guez	147	zpic(:, 1) = sum(zpic(:iim, 1) * weight(:iim, 1)) / zweinor
327			zpic(:, jjm + 1) = sum(zpic(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
328			/ zweisud
329	guez	3
330	guez	147	zval(:, 1) = sum(zval(:iim, 1) * weight(:iim, 1)) / zweinor
331			zval(:, jjm + 1) = sum(zval(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
332			/ zweisud
333	guez	3
334	guez	147	zstd(:, 1) = sum(zstd(:iim, 1) * weight(:iim, 1)) / zweinor
335			zstd(:, jjm + 1) = sum(zstd(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
336			/ zweisud
337	guez	3
338	guez	147	zsig(:, 1) = sum(zsig(:iim, 1) * weight(:iim, 1)) / zweinor
339			zsig(:, jjm + 1) = sum(zsig(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
340			/ zweisud
341	guez	3
342	guez	68	zgam(:, 1) = 1.
343			zgam(:, jjm + 1) = 1.
344	guez	3
345	guez	68	zthe(:, 1) = 0.
346			zthe(:, jjm + 1) = 0.
347	guez	3
348			END SUBROUTINE grid_noro
349
350			end module grid_noro_m