/[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/phylmd/Orography/grid_noro_m.f90 revision 76 by guez, Fri Nov 15 18:45:49 2013 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: François Lott, Laurent Li, A. Harzallah and Laurent      ! Authors: Fran\c{}cois Lott, Laurent Li, A. Harzallah and Laurent
13      ! Fairhead      ! Fairhead
14    
15      ! Compute the parameters of the sub-grid scale orography scheme as      ! Compute the parameters of the sub-grid scale orography scheme as
# Line 29  contains Line 29  contains
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):: xdata(:), ydata(:) ! coordinates of input field      REAL, intent(in):: zdata(:, :) ! (iusn, jusn) input field
40      REAL, intent(in):: zdata(:, :) ! input field      REAL, intent(in):: x(:), y(:) ! coordinates of output field
     REAL, intent(in):: x(:), y(:) ! ccordinates output field  
41    
42      ! Correlations of US Navy orography gradients:      ! Correlations of US Navy orography gradients:
43      REAL, intent(out):: zphi(:, :)      REAL, intent(out):: zphi(:, :) ! (iim + 1, jjm + 1) orography not smoothed
44      real, intent(out):: zmea(:, :) ! Mean orography      real, intent(out):: zmea(:, :) ! (iim + 1, jjm + 1) smoothed orography
45      real, intent(out):: zstd(:, :) ! Standard deviation      real, intent(out):: zstd(:, :) ! (iim + 1, jjm + 1) Standard deviation
46      REAL, intent(out):: zsig(:, :) ! Slope      REAL, intent(out):: zsig(:, :) ! (iim + 1, jjm + 1) Slope
47      real, intent(out):: zgam(:, :) ! Anisotropy      real, intent(out):: zgam(:, :) ! (iim + 1, jjm + 1) Anisotropy
48      real, intent(out):: zthe(:, :) ! Orientation of the small axis  
49      REAL, intent(out):: zpic(:, :) ! Maximum altitude      real, intent(out):: zthe(:, :) ! (iim + 1, jjm + 1)
50      real, intent(out):: zval(:, :) ! Minimum altitude      ! Orientation of the small axis
51    
52        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(:, :) ! 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 58  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
     REAL ztz(iim + 1, jjm + 1), zxtzx(iim + 1, jjm + 1)  
     REAL zytzy(iim + 1, jjm + 1), zxtzy(iim + 1, jjm + 1)  
     REAL weight(iim + 1, jjm + 1)  
68    
69      ! Correlations of US Navy orography gradients:      ! Correlations of US Navy 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 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.      real, parameter:: rad = 6371229.
79    
# Line 93  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"
96      zdeltay = 2. * pi / real(jusn) * rad      zdeltay = 2. * pi / real(jusn) * rad
97    
98      ! Extension of the US Navy database for computations at boundaries:      ! Extension of the US Navy database for computations at boundaries:
# Line 121  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 167  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 218  contains Line 217  contains
217         stop 1         stop 1
218      end if      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 243  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 253  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 264  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 met un seuil à 10  
     ! % de terre car en dessous les paramètres de surface n'ont pas de  
     ! sens.  
     mask_tmp = 0.  
     WHERE (mask >= 0.1) mask_tmp = 1.  
   
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 279  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 289  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 = 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)  
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 323  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 357  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.76  
changed lines
  Added in v.147
  ViewVC Help
Powered by ViewVC 1.1.21