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\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 dimens_m, 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
    REAL, intent(in):: x(:), y(:) ! coordinates of output field

    ! Correlations of US Navy orography gradients:
    REAL, intent(out):: zphi(:, :) ! (iim + 1, jjm + 1) orography not smoothed
    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

    ! 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
    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			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	guez	147	! Coordinates of input field:
36			REAL, intent(in):: xdata(:) ! (iusn)
37			REAL, intent(in):: ydata(:) ! (jusn)
38
39			REAL, intent(in):: zdata(:, :) ! (iusn, jusn) input field
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	147	REAL, intent(out):: zphi(:, :) ! (iim + 1, jjm + 1) orography not smoothed
44			real, intent(out):: zmea(:, :) ! (iim + 1, jjm + 1) smoothed orography
45			real, intent(out):: zstd(:, :) ! (iim + 1, jjm + 1) Standard deviation
46			REAL, intent(out):: zsig(:, :) ! (iim + 1, jjm + 1) Slope
47			real, intent(out):: zgam(:, :) ! (iim + 1, jjm + 1) Anisotropy
48	guez	3
49	guez	147	real, intent(out):: zthe(:, :) ! (iim + 1, jjm + 1)
50			! Orientation of the small axis
51	guez	3
52	guez	147	REAL, intent(out):: zpic(:, :) ! (iim + 1, jjm + 1) Maximum altitude
53			real, intent(out):: zval(:, :) ! (iim + 1, jjm + 1) Minimum altitude
54
55			real, intent(out):: mask(:, :) ! (iim + 1, jjm + 1) fraction of land
56
57	guez	3	! Variables local to the procedure:
58
59			! In this version it is assumed that the input data come from
60			! the US Navy dataset:
61	guez	68	integer, parameter:: iusn = 2160, jusn = 1080
62			integer, parameter:: iext = 216
63			REAL xusn(iusn + 2 * iext), yusn(jusn + 2)
64			REAL zusn(iusn + 2 * iext, jusn + 2)
65	guez	3
66	guez	68	! Intermediate fields (correlations of orography gradient)
67	guez	78	REAL, dimension(iim + 1, jjm + 1):: ztz, zxtzx, zytzy, zxtzy, weight
68	guez	3
69	guez	70	! Correlations of US Navy orography gradients:
70	guez	68	REAL, dimension(iusn + 2 * iext, jusn + 2):: zxtzxusn, zytzyusn, zxtzyusn
71	guez	3
72	guez	147	real, dimension(iim + 1, jjm + 1):: mask_tmp, num_tot, num_lan
73	guez	68	REAL a(iim + 1), b(iim + 1), c(jjm + 1), d(jjm + 1)
74	guez	70	real weighx, weighy, xincr, xk, xp, xm, xw, xq, xl
75	guez	3	real zbordnor, zdeltax, zbordsud, zdeltay, zbordoue, zlenx, zleny, zmeasud
76	guez	147	real zweinor, zweisud, zmeanor, zbordest
77	guez	3	integer ii, i, jj, j
78	guez	70	real, parameter:: rad = 6371229.
79	guez	3
80			!-------------------------------
81
82			print *, "Call sequence information: grid_noro"
83
84			call assert((/size(xdata), size(zdata, 1)/) == iusn, "grid_noro iusn")
85			call assert((/size(ydata), size(zdata, 2)/) == jusn, "grid_noro jusn")
86
87			call assert((/size(x), size(zphi, 1), size(zmea, 1), size(zstd, 1), &
88			size(zsig, 1), size(zgam, 1), size(zthe, 1), size(zpic, 1), &
89			size(zval, 1), size(mask, 1)/) == iim + 1, "grid_noro iim")
90
91			call assert((/size(y), size(zphi, 2), size(zmea, 2), size(zstd, 2), &
92			size(zsig, 2), size(zgam, 2), size(zthe, 2), size(zpic, 2), &
93			size(zval, 2), size(mask, 2)/) == jjm + 1, "grid_noro jjm")
94
95	guez	147	print *, "Parameters of subgrid-scale orography"
96	guez	3	zdeltay = 2. * pi / real(jusn) * rad
97
98	guez	70	! Extension of the US Navy database for computations at boundaries:
99	guez	3
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	78	! 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	70	! Compute slopes correlations on US Navy 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	78	! 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	guez	70	weighy = MAX(0., min(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	guez	70	weighx = MAX(0., min(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	70	zpic(ii, jj) = max(zpic(ii, jj), zusn(i, j))
206	guez	3	! valleys
207	guez	70	zval(ii, jj) = min(zval(ii, jj), zusn(i, j))
208	guez	3	ENDIF
209			ENDDO
210			ENDIF
211			ENDDO
212			ENDDO
213			ENDDO
214
215	guez	70	if (any(weight == 0.)) then
216			print *, "zero weight in grid_noro"
217			stop 1
218			end if
219	guez	3
220	guez	78	! Compute parameters needed by the Lott & Miller (1997) and Lott
221			! (1999) subgrid-scale orographic scheme.
222	guez	3
223	guez	68	DO ii = 1, iim + 1
224	guez	3	DO jj = 1, jjm + 1
225	guez	68	mask(ii, jj) = num_lan(ii, jj) / num_tot(ii, jj)
226	guez	78	! Mean orography:
227	guez	68	zmea (ii, jj) = zmea (ii, jj) / weight(ii, jj)
228			zxtzx(ii, jj) = zxtzx(ii, jj) / weight(ii, jj)
229			zytzy(ii, jj) = zytzy(ii, jj) / weight(ii, jj)
230			zxtzy(ii, jj) = zxtzy(ii, jj) / weight(ii, jj)
231			ztz(ii, jj) = ztz(ii, jj) / weight(ii, jj)
232			! Standard deviation:
233			zstd(ii, jj) = sqrt(MAX(0., ztz(ii, jj) - zmea(ii, jj)**2))
234	guez	3	ENDDO
235			ENDDO
236
237	guez	78	! Correct values of horizontal slope near the poles:
238	guez	68	DO ii = 1, iim + 1
239			zxtzx(ii, 1) = zxtzx(ii, 2)
240			zxtzx(ii, jjm + 1) = zxtzx(ii, jjm)
241			zxtzy(ii, 1) = zxtzy(ii, 2)
242			zxtzy(ii, jjm + 1) = zxtzy(ii, jjm)
243			zytzy(ii, 1) = zytzy(ii, 2)
244			zytzy(ii, jjm + 1) = zytzy(ii, jjm)
245	guez	3	ENDDO
246
247	guez	147	! Masque prenant en compte maximum de terre. On met un seuil \`a 10
248			! % de terre car en dessous les param\`etres de surface n'ont pas de
249			! sens.
250			mask_tmp = merge(1., 0., mask >= 0.1)
251	guez	3
252	guez	147	zphi(:iim, :) = zmea(:iim, :) * mask_tmp(:iim, :)
253			! (zmea is not yet smoothed)
254
255	guez	78	! Filters to smooth out fields for input into subgrid-scale
256			! orographic scheme.
257
258			! 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			DO ii = 1, iim
268			DO jj = 1, jjm + 1
269	guez	68	! Coefficients K, L et M:
270			xk = (zxtzx(ii, jj) + zytzy(ii, jj)) / 2.
271			xl = (zxtzx(ii, jj) - zytzy(ii, jj)) / 2.
272			xm = zxtzy(ii, jj)
273			xp = xk - sqrt(xl2 + xm2)
274			xq = xk + sqrt(xl2 + xm2)
275			xw = 1e-8
276	guez	147	if (xp <= xw) xp = 0.
277			if (xq <= xw) xq = xw
278			if (abs(xm) <= xw) xm = xw * sign(1., xm)
279	guez	68	! modification pour masque de terre fractionnaire
280			! slope:
281			zsig(ii, jj) = sqrt(xq) * mask_tmp(ii, jj)
282			! isotropy:
283			zgam(ii, jj) = xp / xq * mask_tmp(ii, jj)
284			! angle theta:
285			zthe(ii, jj) = 57.29577951 * atan2(xm, xl) / 2. * mask_tmp(ii, jj)
286	guez	3	ENDDO
287			ENDDO
288	guez	68
289	guez	147	zmea(:iim, :) = zmea(:iim, :) * mask_tmp(:iim, :)
290			zpic(:iim, :) = zpic(:iim, :) * mask_tmp(:iim, :)
291			zval(:iim, :) = zval(:iim, :) * mask_tmp(:iim, :)
292			zstd(:iim, :) = zstd(:iim, :) * mask_tmp(:iim, :)
293	guez	3
294	guez	147	print *, 'MEAN ORO: ', MAXVAL(zmea(:iim, :))
295			print *, 'ST. DEV.: ', MAXVAL(zstd(:iim, :))
296			print *, 'PENTE: ', MAXVAL(zsig(:iim, :))
297			print *, 'ANISOTROP: ', MAXVAL(zgam(:iim, :))
298			print *, 'ANGLE: ', minval(zthe(:iim, :)), MAXVAL(zthe(:iim, :))
299			print *, 'pic: ', MAXVAL(zpic(:iim, :))
300			print *, 'val: ', MAXVAL(zval(:iim, :))
301
302	guez	68	! gamma and theta at 1. and 0. at poles
303			zmea(iim + 1, :) = zmea(1, :)
304			zphi(iim + 1, :) = zphi(1, :)
305			zpic(iim + 1, :) = zpic(1, :)
306			zval(iim + 1, :) = zval(1, :)
307			zstd(iim + 1, :) = zstd(1, :)
308			zsig(iim + 1, :) = zsig(1, :)
309			zgam(iim + 1, :) = zgam(1, :)
310			zthe(iim + 1, :) = zthe(1, :)
311	guez	3
312	guez	147	zweinor = sum(weight(:iim, 1))
313			zweisud = sum(weight(:iim, jjm + 1))
314			zmeanor = sum(zmea(:iim, 1) * weight(:iim, 1))
315			zmeasud = sum(zmea(:iim, jjm + 1) * weight(:iim, jjm + 1))
316	guez	3
317	guez	68	zmea(:, 1) = zmeanor / zweinor
318			zmea(:, jjm + 1) = zmeasud / zweisud
319	guez	3
320	guez	68	zphi(:, 1) = zmeanor / zweinor
321			zphi(:, jjm + 1) = zmeasud / zweisud
322	guez	3
323	guez	147	zpic(:, 1) = sum(zpic(:iim, 1) * weight(:iim, 1)) / zweinor
324			zpic(:, jjm + 1) = sum(zpic(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
325			/ zweisud
326	guez	3
327	guez	147	zval(:, 1) = sum(zval(:iim, 1) * weight(:iim, 1)) / zweinor
328			zval(:, jjm + 1) = sum(zval(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
329			/ zweisud
330	guez	3
331	guez	147	zstd(:, 1) = sum(zstd(:iim, 1) * weight(:iim, 1)) / zweinor
332			zstd(:, jjm + 1) = sum(zstd(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
333			/ zweisud
334	guez	3
335	guez	147	zsig(:, 1) = sum(zsig(:iim, 1) * weight(:iim, 1)) / zweinor
336			zsig(:, jjm + 1) = sum(zsig(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
337			/ zweisud
338	guez	3
339	guez	68	zgam(:, 1) = 1.
340			zgam(:, jjm + 1) = 1.
341	guez	3
342	guez	68	zthe(:, 1) = 0.
343			zthe(:, jjm + 1) = 0.
344	guez	3
345			END SUBROUTINE grid_noro
346
347			end module grid_noro_m