/[lmdze]/trunk/Sources/phylmd/Orography/grid_noro_m.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/Orography/grid_noro_m.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/phylmd/Orography/grid_noro_m.f90 revision 68 by guez, Wed Nov 14 16:59:30 2012 UTC trunk/Sources/phylmd/Orography/grid_noro_m.f revision 147 by guez, Wed Jun 17 14:20:14 2015 UTC
# Line 9  contains Line 9  contains
9    
10      ! From dyn3d/grid_noro.F, version 1.1.1.1 2004/05/19 12:53:06      ! From dyn3d/grid_noro.F, version 1.1.1.1 2004/05/19 12:53:06
11    
12      ! Authors: F. Lott, Z. X. Li, A. Harzallah and L. Fairhead      ! Authors: Fran\c{}cois Lott, Laurent Li, A. Harzallah and Laurent
13        ! Fairhead
14    
15        ! Compute the parameters of the sub-grid scale orography scheme as
16        ! described in Lott and Miller (1997) and Lott (1999).
17    
     ! Compute the parameters of the SSO scheme as described in  
     ! Lott and Miller (1997) and Lott (1999).  
18      ! Target points are on a rectangular grid:      ! Target points are on a rectangular grid:
19      ! jjm + 1 latitudes including North and South Poles;      ! jjm + 1 latitudes including North and South Poles;
20      ! iim + 1 longitudes, with periodicity: longitude(iim + 1) = longitude(1)      ! iim + 1 longitudes, with periodicity: longitude(iim + 1) = longitude(1)
# Line 20  contains Line 22  contains
22    
23      ! The parameters a, b, c, d represent the limite of the target      ! The parameters a, b, c, d represent the limite of the target
24      ! gridpoint region. The means over this region are calculated from      ! gridpoint region. The means over this region are calculated from
25      ! USN data, ponderated by a weight proportional to the surface      ! US Navy data, ponderated by a weight proportional to the surface
26      ! occupied by the data inside the model gridpoint area. In most      ! occupied by the data inside the model gridpoint area. In most
27      ! circumstances, this weight is the ratio between the surface of      ! circumstances, this weight is the ratio between the surface of
28      ! the USN gridpoint area and the surface of the model gridpoint      ! the US Navy gridpoint area and the surface of the model gridpoint
29      ! area. See "grid_noto.txt".      ! area. See "grid_noto.txt".
30    
31      use dimens_m, only: iim, jjm      use dimens_m, only: iim, jjm
     use nr_util, only: assert, pi  
32      use mva9_m, only: mva9      use mva9_m, only: mva9
33        use nr_util, only: assert, pi
34    
35        ! 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        REAL, intent(in):: x(:), y(:) ! coordinates of output field
41    
42        ! Correlations of US Navy orography gradients:
43        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    
49        real, intent(out):: zthe(:, :) ! (iim + 1, jjm + 1)
50        ! Orientation of the small axis
51    
52      REAL, intent(in):: xdata(:), ydata(:) ! coordinates of input field      REAL, intent(out):: zpic(:, :) ! (iim + 1, jjm + 1) Maximum altitude
53      REAL, intent(in):: zdata(:, :) ! input field      real, intent(out):: zval(:, :) ! (iim + 1, jjm + 1) Minimum altitude
     REAL, intent(in):: x(:), y(:) ! ccordinates output field  
   
     ! Correlations of USN orography gradients:  
   
     REAL, intent(out):: zphi(:, :)  
     real, intent(out):: zmea(:, :) ! Mean orography  
     real, intent(out):: zstd(:, :) ! Standard deviation  
     REAL zsig(:, :) ! Slope  
     real zgam(:, :) ! Anisotropy  
     real zthe(:, :) ! Orientation of the small axis  
     REAL, intent(out):: zpic(:, :) ! Maximum altitude  
     real, intent(out):: zval(:, :) ! Minimum altitude  
54    
55      real, intent(out):: mask(:, :) ! fraction of land      real, intent(out):: mask(:, :) ! (iim + 1, jjm + 1) fraction of land
56    
57      ! Variables local to the procedure:      ! Variables local to the procedure:
58    
# Line 57  contains Line 64  contains
64      REAL zusn(iusn + 2 * iext, jusn + 2)      REAL zusn(iusn + 2 * iext, jusn + 2)
65    
66      ! Intermediate fields (correlations of orography gradient)      ! Intermediate fields (correlations of orography gradient)
67        REAL, dimension(iim + 1, jjm + 1):: ztz, zxtzx, zytzy, zxtzy, weight
68    
69      REAL ztz(iim + 1, jjm + 1), zxtzx(iim + 1, jjm + 1)      ! Correlations of US Navy orography gradients:
     REAL zytzy(iim + 1, jjm + 1), zxtzy(iim + 1, jjm + 1)  
     REAL weight(iim + 1, jjm + 1)  
   
     ! Correlations of USN orography gradients:  
70      REAL, dimension(iusn + 2 * iext, jusn + 2):: zxtzxusn, zytzyusn, zxtzyusn      REAL, dimension(iusn + 2 * iext, jusn + 2):: zxtzxusn, zytzyusn, zxtzyusn
71    
72      real mask_tmp(size(x), size(y))      real, dimension(iim + 1, jjm + 1):: mask_tmp, num_tot, num_lan
     real num_tot(iim + 1, jjm + 1), num_lan(iim + 1, jjm + 1)  
   
73      REAL a(iim + 1), b(iim + 1), c(jjm + 1), d(jjm + 1)      REAL a(iim + 1), b(iim + 1), c(jjm + 1), d(jjm + 1)
74      real rad, weighx, weighy, xincr, xk, xp, xm, xw, xq, xl      real weighx, weighy, xincr, xk, xp, xm, xw, xq, xl
75      real zbordnor, zdeltax, zbordsud, zdeltay, zbordoue, zlenx, zleny, zmeasud      real zbordnor, zdeltax, zbordsud, zdeltay, zbordoue, zlenx, zleny, zmeasud
76      real zllmpic, zllmmea, zllmgam, zllmthe, zllmstd, zllmsig, zllmval      real zweinor, zweisud, zmeanor, zbordest
     real zpicnor, zminthe, zsigsud, zstdnor, zstdsud, zvalsud, zvalnor  
     real zweinor, zweisud, zsignor, zpicsud, zmeanor, zbordest  
77      integer ii, i, jj, j      integer ii, i, jj, j
78        real, parameter:: rad = 6371229.
79    
80      !-------------------------------      !-------------------------------
81    
# Line 91  contains Line 92  contains
92           size(zsig, 2), size(zgam, 2), size(zthe, 2), size(zpic, 2), &           size(zsig, 2), size(zgam, 2), size(zthe, 2), size(zpic, 2), &
93           size(zval, 2), size(mask, 2)/) == jjm + 1, "grid_noro jjm")           size(zval, 2), size(mask, 2)/) == jjm + 1, "grid_noro jjm")
94    
95      print *, "Paramètres de l'orographie à l'échelle sous-maille"      print *, "Parameters of subgrid-scale orography"
     rad = 6371229.  
96      zdeltay = 2. * pi / real(jusn) * rad      zdeltay = 2. * pi / real(jusn) * rad
97    
98      ! Extension of the USN database to POCEED computations at boundaries:      ! Extension of the US Navy database for computations at boundaries:
99    
100      DO j = 1, jusn      DO j = 1, jusn
101         yusn(j + 1) = ydata(j)         yusn(j + 1) = ydata(j)
# Line 120  contains Line 120  contains
120         zusn(i + iusn / 2 + iext, jusn + 2) = zusn(i, jusn + 1)         zusn(i + iusn / 2 + iext, jusn + 2) = zusn(i, jusn + 1)
121      ENDDO      ENDDO
122    
123      ! COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)      ! Compute limits of model gridpoint area (regular grid)
124    
125      a(1) = x(1) - (x(2) - x(1)) / 2.0      a(1) = x(1) - (x(2) - x(1)) / 2.0
126      b(1) = (x(1) + x(2)) / 2.0      b(1) = (x(1) + x(2)) / 2.0
# Line 150  contains Line 150  contains
150      zpic = - 1E10      zpic = - 1E10
151      zval = 1E10      zval = 1E10
152    
153      ! COMPUTE SLOPES CORRELATIONS ON USN GRID      ! Compute slopes correlations on US Navy grid
154    
155      zytzyusn = 0.      zytzyusn = 0.
156      zxtzxusn = 0.      zxtzxusn = 0.
# Line 166  contains Line 166  contains
166         ENDDO         ENDDO
167      ENDDO      ENDDO
168    
169      ! SUMMATION OVER GRIDPOINT AREA      ! Summation over gridpoint area
170    
171      zleny = pi / real(jusn) * rad      zleny = pi / real(jusn) * rad
172      xincr = pi / 2. / real(jusn)      xincr = pi / 2. / real(jusn)
# Line 179  contains Line 179  contains
179               zdeltax = zdeltay * cos(yusn(j))               zdeltax = zdeltay * cos(yusn(j))
180               zbordnor = (c(jj) - yusn(j) + xincr) * rad               zbordnor = (c(jj) - yusn(j) + xincr) * rad
181               zbordsud = (yusn(j) - d(jj) + xincr) * rad               zbordsud = (yusn(j) - d(jj) + xincr) * rad
182               weighy = AMAX1(0., amin1(zbordnor, zbordsud, zleny))               weighy = MAX(0., min(zbordnor, zbordsud, zleny))
183               IF (weighy /= 0) THEN               IF (weighy /= 0) THEN
184                  DO i = 2, iusn + 2 * iext - 1                  DO i = 2, iusn + 2 * iext - 1
185                     zbordest = (xusn(i) - a(ii) + xincr) * rad * cos(yusn(j))                     zbordest = (xusn(i) - a(ii) + xincr) * rad * cos(yusn(j))
186                     zbordoue = (b(ii) + xincr - xusn(i)) * rad * cos(yusn(j))                     zbordoue = (b(ii) + xincr - xusn(i)) * rad * cos(yusn(j))
187                     weighx = AMAX1(0., amin1(zbordest, zbordoue, zlenx))                     weighx = MAX(0., min(zbordest, zbordoue, zlenx))
188                     IF (weighx /= 0) THEN                     IF (weighx /= 0) THEN
189                        num_tot(ii, jj) = num_tot(ii, jj) + 1.                        num_tot(ii, jj) = num_tot(ii, jj) + 1.
190                        if (zusn(i, j) >= 1.) then                        if (zusn(i, j) >= 1.) then
# Line 202  contains Line 202  contains
202                        ! mean                        ! mean
203                        zmea(ii, jj) = zmea(ii, jj) + zusn(i, j) * weighx * weighy                        zmea(ii, jj) = zmea(ii, jj) + zusn(i, j) * weighx * weighy
204                        ! peacks                        ! peacks
205                        zpic(ii, jj) = amax1(zpic(ii, jj), zusn(i, j))                        zpic(ii, jj) = max(zpic(ii, jj), zusn(i, j))
206                        ! valleys                        ! valleys
207                        zval(ii, jj) = amin1(zval(ii, jj), zusn(i, j))                        zval(ii, jj) = min(zval(ii, jj), zusn(i, j))
208                     ENDIF                     ENDIF
209                  ENDDO                  ENDDO
210               ENDIF               ENDIF
# Line 212  contains Line 212  contains
212         ENDDO         ENDDO
213      ENDDO      ENDDO
214    
215      if (any(weight == 0.)) stop "zero weight in grid_noro"      if (any(weight == 0.)) then
216           print *, "zero weight in grid_noro"
217           stop 1
218        end if
219    
220      ! COMPUTE PARAMETERS NEEDED BY THE LOTT & MILLER (1997) AND      ! Compute parameters needed by the Lott & Miller (1997) and Lott
221      ! LOTT (1999) SSO SCHEME.      ! (1999) subgrid-scale orographic scheme.
222    
     zllmmea = 0.  
     zllmstd = 0.  
     zllmsig = 0.  
     zllmgam = 0.  
     zllmpic = 0.  
     zllmval = 0.  
     zllmthe = 0.  
     zminthe = 0.  
223      DO ii = 1, iim + 1      DO ii = 1, iim + 1
224         DO jj = 1, jjm + 1         DO jj = 1, jjm + 1
225            mask(ii, jj) = num_lan(ii, jj) / num_tot(ii, jj)            mask(ii, jj) = num_lan(ii, jj) / num_tot(ii, jj)
226            ! Mean Orography:            ! Mean orography:
227            zmea (ii, jj) = zmea (ii, jj) / weight(ii, jj)            zmea (ii, jj) = zmea (ii, jj) / weight(ii, jj)
228            zxtzx(ii, jj) = zxtzx(ii, jj) / weight(ii, jj)            zxtzx(ii, jj) = zxtzx(ii, jj) / weight(ii, jj)
229            zytzy(ii, jj) = zytzy(ii, jj) / weight(ii, jj)            zytzy(ii, jj) = zytzy(ii, jj) / weight(ii, jj)
# Line 239  contains Line 234  contains
234         ENDDO         ENDDO
235      ENDDO      ENDDO
236    
237      ! CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:      ! Correct values of horizontal slope near the poles:
238      DO ii = 1, iim + 1      DO ii = 1, iim + 1
239         zxtzx(ii, 1) = zxtzx(ii, 2)         zxtzx(ii, 1) = zxtzx(ii, 2)
240         zxtzx(ii, jjm + 1) = zxtzx(ii, jjm)         zxtzx(ii, jjm + 1) = zxtzx(ii, jjm)
# Line 249  contains Line 244  contains
244         zytzy(ii, jjm + 1) = zytzy(ii, jjm)         zytzy(ii, jjm + 1) = zytzy(ii, jjm)
245      ENDDO      ENDDO
246    
247      ! FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.      ! 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    
252        zphi(:iim, :) = zmea(:iim, :) * mask_tmp(:iim, :)
253        ! (zmea is not yet smoothed)
254    
255      ! FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.      ! Filters to smooth out fields for input into subgrid-scale
256        ! orographic scheme.
257    
258        ! First filter, moving average over 9 points.
259      CALL MVA9(zmea)      CALL MVA9(zmea)
260      CALL MVA9(zstd)      CALL MVA9(zstd)
261      CALL MVA9(zpic)      CALL MVA9(zpic)
# Line 260  contains Line 264  contains
264      CALL MVA9(zxtzy)      CALL MVA9(zxtzy)
265      CALL MVA9(zytzy)      CALL MVA9(zytzy)
266    
     ! Masque prenant en compte maximum de terre. On seuille à 10 % de  
     ! terre car en dessous les paramètres de surface n'ont pas de  
     ! sens.  
     mask_tmp = 0.  
     WHERE (mask >= 0.1) mask_tmp = 1.  
   
267      DO ii = 1, iim      DO ii = 1, iim
268         DO jj = 1, jjm + 1         DO jj = 1, jjm + 1
269            ! Coefficients K, L et M:            ! Coefficients K, L et M:
# Line 275  contains Line 273  contains
273            xp = xk - sqrt(xl**2 + xm**2)            xp = xk - sqrt(xl**2 + xm**2)
274            xq = xk + sqrt(xl**2 + xm**2)            xq = xk + sqrt(xl**2 + xm**2)
275            xw = 1e-8            xw = 1e-8
276            if(xp.le.xw) xp = 0.            if (xp <= xw) xp = 0.
277            if(xq.le.xw) xq = xw            if (xq <= xw) xq = xw
278            if(abs(xm).le.xw) xm = xw * sign(1., xm)            if (abs(xm) <= xw) xm = xw * sign(1., xm)
279            ! modification pour masque de terre fractionnaire            ! modification pour masque de terre fractionnaire
280            ! slope:            ! slope:
281            zsig(ii, jj) = sqrt(xq) * mask_tmp(ii, jj)            zsig(ii, jj) = sqrt(xq) * mask_tmp(ii, jj)
# Line 285  contains Line 283  contains
283            zgam(ii, jj) = xp / xq * mask_tmp(ii, jj)            zgam(ii, jj) = xp / xq * mask_tmp(ii, jj)
284            ! angle theta:            ! angle theta:
285            zthe(ii, jj) = 57.29577951 * atan2(xm, xl) / 2. * mask_tmp(ii, jj)            zthe(ii, jj) = 57.29577951 * atan2(xm, xl) / 2. * mask_tmp(ii, jj)
           zphi(ii, jj) = zmea(ii, jj) * mask_tmp(ii, jj)  
           zmea(ii, jj) = zmea(ii, jj) * mask_tmp(ii, jj)  
           zpic(ii, jj) = zpic(ii, jj) * mask_tmp(ii, jj)  
           zval(ii, jj) = zval(ii, jj) * mask_tmp(ii, jj)  
           zstd(ii, jj) = zstd(ii, jj) * mask_tmp(ii, jj)  
           zllmmea = AMAX1(zmea(ii, jj), zllmmea)  
           zllmstd = AMAX1(zstd(ii, jj), zllmstd)  
           zllmsig = AMAX1(zsig(ii, jj), zllmsig)  
           zllmgam = AMAX1(zgam(ii, jj), zllmgam)  
           zllmthe = AMAX1(zthe(ii, jj), zllmthe)  
           zminthe = amin1(zthe(ii, jj), zminthe)  
           zllmpic = AMAX1(zpic(ii, jj), zllmpic)  
           zllmval = AMAX1(zval(ii, jj), zllmval)  
286         ENDDO         ENDDO
287      ENDDO      ENDDO
288    
289      print *, 'MEAN ORO: ', zllmmea      zmea(:iim, :) = zmea(:iim, :) * mask_tmp(:iim, :)
290      print *, 'ST. DEV.: ', zllmstd      zpic(:iim, :) = zpic(:iim, :) * mask_tmp(:iim, :)
291      print *, 'PENTE: ', zllmsig      zval(:iim, :) = zval(:iim, :) * mask_tmp(:iim, :)
292      print *, 'ANISOTROP: ', zllmgam      zstd(:iim, :) = zstd(:iim, :) * mask_tmp(:iim, :)
293      print *, 'ANGLE: ', zminthe, zllmthe  
294      print *, 'pic: ', zllmpic      print *, 'MEAN ORO: ', MAXVAL(zmea(:iim, :))
295      print *, 'val: ', zllmval      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      ! gamma and theta at 1. and 0. at poles      ! gamma and theta at 1. and 0. at poles
303      zmea(iim + 1, :) = zmea(1, :)      zmea(iim + 1, :) = zmea(1, :)
# Line 319  contains Line 309  contains
309      zgam(iim + 1, :) = zgam(1, :)      zgam(iim + 1, :) = zgam(1, :)
310      zthe(iim + 1, :) = zthe(1, :)      zthe(iim + 1, :) = zthe(1, :)
311    
312      zmeanor = 0.      zweinor = sum(weight(:iim, 1))
313      zmeasud = 0.      zweisud = sum(weight(:iim, jjm + 1))
314      zstdnor = 0.      zmeanor = sum(zmea(:iim, 1) * weight(:iim, 1))
315      zstdsud = 0.      zmeasud = sum(zmea(:iim, jjm + 1) * weight(:iim, jjm + 1))
     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  
316    
317      zmea(:, 1) = zmeanor / zweinor      zmea(:, 1) = zmeanor / zweinor
318      zmea(:, jjm + 1) = zmeasud / zweisud      zmea(:, jjm + 1) = zmeasud / zweisud
# Line 353  contains Line 320  contains
320      zphi(:, 1) = zmeanor / zweinor      zphi(:, 1) = zmeanor / zweinor
321      zphi(:, jjm + 1) = zmeasud / zweisud      zphi(:, jjm + 1) = zmeasud / zweisud
322    
323      zpic(:, 1) = zpicnor / zweinor      zpic(:, 1) = sum(zpic(:iim, 1) * weight(:iim, 1)) / zweinor
324      zpic(:, jjm + 1) = zpicsud / zweisud      zpic(:, jjm + 1) = sum(zpic(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
325             / zweisud
326      zval(:, 1) = zvalnor / zweinor  
327      zval(:, jjm + 1) = zvalsud / zweisud      zval(:, 1) = sum(zval(:iim, 1) * weight(:iim, 1)) / zweinor
328        zval(:, jjm + 1) = sum(zval(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
329      zstd(:, 1) = zstdnor / zweinor           / zweisud
330      zstd(:, jjm + 1) = zstdsud / zweisud  
331        zstd(:, 1) = sum(zstd(:iim, 1) * weight(:iim, 1)) / zweinor
332      zsig(:, 1) = zsignor / zweinor      zstd(:, jjm + 1) = sum(zstd(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
333      zsig(:, jjm + 1) = zsigsud / zweisud           / zweisud
334    
335        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    
339      zgam(:, 1) = 1.      zgam(:, 1) = 1.
340      zgam(:, jjm + 1) = 1.      zgam(:, jjm + 1) = 1.
Legend:
Removed from v.68  
changed lines
  Added in v.147
  ViewVC Help
Powered by ViewVC 1.1.21