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

Diff of /trunk/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 70 by guez, Mon Jun 24 15:39:52 2013 UTC trunk/phylmd/Orography/grid_noro_m.f revision 265 by guez, Tue Mar 20 09:35:59 2018 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 28  contains Line 28  contains
28      ! the US Navy 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 dimensions, 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, in m
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:
     REAL, intent(out):: zphi(:, :)  
     real, intent(out):: zmea(:, :) ! Mean 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  
43    
44      real, intent(out):: mask(:, :) ! fraction of land      REAL, intent(out):: zphi(:, :) ! (iim + 1, jjm + 1)
45        ! geoptential height of orography, not smoothed, in m
46        
47        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    
52        real, intent(out):: zthe(:, :) ! (iim + 1, jjm + 1)
53        ! Orientation of the small axis
54    
55        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      ! Variables local to the procedure:      ! Local:
61    
62      ! In this version it is assumed that the input data come from      ! In this version it is assumed that the input data come from
63      ! the US Navy dataset:      ! the US Navy dataset:
64      integer, parameter:: iusn = 2160, jusn = 1080      integer, parameter:: iusn = 2160, jusn = 1080
65      integer, parameter:: iext = 216      integer, parameter:: iext = 216
66      REAL xusn(iusn + 2 * iext), yusn(jusn + 2)      REAL xusn(iusn + 2 * iext), yusn(jusn + 2)
67      REAL zusn(iusn + 2 * iext, jusn + 2)      REAL zusn(iusn + 2 * iext, jusn + 2) ! in m
68    
69      ! Intermediate fields (correlations of orography gradient)      ! Intermediate fields (correlations of orography gradient)
70        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)  
71    
72      ! Correlations of US Navy orography gradients:      ! Correlations of US Navy orography gradients:
73      REAL, dimension(iusn + 2 * iext, jusn + 2):: zxtzxusn, zytzyusn, zxtzyusn      REAL, dimension(iusn + 2 * iext, jusn + 2):: zxtzxusn, zytzyusn, zxtzyusn
74    
75      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)  
   
76      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)
77      real weighx, weighy, xincr, xk, xp, xm, xw, xq, xl      real weighx, weighy, xincr, xk, xp, xm, xw, xq, xl
78      real zbordnor, zdeltax, zbordsud, zdeltay, zbordoue, zlenx, zleny, zmeasud      real zbordnor, zdeltax, zbordsud, zdeltay, zbordoue, zlenx, zleny, zmeasud
79      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  
80      integer ii, i, jj, j      integer ii, i, jj, j
81      real, parameter:: rad = 6371229.      real, parameter:: rad = 6371229.
82    
# Line 93  contains Line 95  contains
95           size(zsig, 2), size(zgam, 2), size(zthe, 2), size(zpic, 2), &           size(zsig, 2), size(zgam, 2), size(zthe, 2), size(zpic, 2), &
96           size(zval, 2), size(mask, 2)/) == jjm + 1, "grid_noro jjm")           size(zval, 2), size(mask, 2)/) == jjm + 1, "grid_noro jjm")
97    
98      print *, "Paramètres de l'orographie à l'échelle sous-maille"      print *, "Parameters of subgrid-scale orography"
99      zdeltay = 2. * pi / real(jusn) * rad      zdeltay = 2. * pi / real(jusn) * rad
100    
101      ! Extension of the US Navy database for computations at boundaries:      ! Extension of the US Navy database for computations at boundaries:
# Line 121  contains Line 123  contains
123         zusn(i + iusn / 2 + iext, jusn + 2) = zusn(i, jusn + 1)         zusn(i + iusn / 2 + iext, jusn + 2) = zusn(i, jusn + 1)
124      ENDDO      ENDDO
125    
126      ! COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)      ! Compute limits of model gridpoint area (regular grid)
127    
128      a(1) = x(1) - (x(2) - x(1)) / 2.0      a(1) = x(1) - (x(2) - x(1)) / 2.0
129      b(1) = (x(1) + x(2)) / 2.0      b(1) = (x(1) + x(2)) / 2.0
# Line 167  contains Line 169  contains
169         ENDDO         ENDDO
170      ENDDO      ENDDO
171    
172      ! SUMMATION OVER GRIDPOINT AREA      ! Summation over gridpoint area
173    
174      zleny = pi / real(jusn) * rad      zleny = pi / real(jusn) * rad
175      xincr = pi / 2. / real(jusn)      xincr = pi / 2. / real(jusn)
# Line 218  contains Line 220  contains
220         stop 1         stop 1
221      end if      end if
222    
223      ! COMPUTE PARAMETERS NEEDED BY THE LOTT & MILLER (1997) AND      ! Compute parameters needed by the Lott & Miller (1997) and Lott
224      ! LOTT (1999) SSO SCHEME.      ! (1999) subgrid-scale orographic scheme.
225    
     zllmmea = 0.  
     zllmstd = 0.  
     zllmsig = 0.  
     zllmgam = 0.  
     zllmpic = 0.  
     zllmval = 0.  
     zllmthe = 0.  
     zminthe = 0.  
226      DO ii = 1, iim + 1      DO ii = 1, iim + 1
227         DO jj = 1, jjm + 1         DO jj = 1, jjm + 1
228            mask(ii, jj) = num_lan(ii, jj) / num_tot(ii, jj)            mask(ii, jj) = num_lan(ii, jj) / num_tot(ii, jj)
229            ! Mean Orography:            ! Mean orography:
230            zmea (ii, jj) = zmea (ii, jj) / weight(ii, jj)            zmea (ii, jj) = zmea (ii, jj) / weight(ii, jj)
231            zxtzx(ii, jj) = zxtzx(ii, jj) / weight(ii, jj)            zxtzx(ii, jj) = zxtzx(ii, jj) / weight(ii, jj)
232            zytzy(ii, jj) = zytzy(ii, jj) / weight(ii, jj)            zytzy(ii, jj) = zytzy(ii, jj) / weight(ii, jj)
# Line 243  contains Line 237  contains
237         ENDDO         ENDDO
238      ENDDO      ENDDO
239    
240      ! CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:      ! Correct values of horizontal slope near the poles:
241      DO ii = 1, iim + 1      DO ii = 1, iim + 1
242         zxtzx(ii, 1) = zxtzx(ii, 2)         zxtzx(ii, 1) = zxtzx(ii, 2)
243         zxtzx(ii, jjm + 1) = zxtzx(ii, jjm)         zxtzx(ii, jjm + 1) = zxtzx(ii, jjm)
# Line 253  contains Line 247  contains
247         zytzy(ii, jjm + 1) = zytzy(ii, jjm)         zytzy(ii, jjm + 1) = zytzy(ii, jjm)
248      ENDDO      ENDDO
249    
250      ! FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.      ! 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    
255        zphi(:iim, :) = zmea(:iim, :) * mask_tmp(:iim, :)
256        ! (zmea is not yet smoothed)
257    
258      ! FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.      ! Filters to smooth out fields for input into subgrid-scale
259        ! orographic scheme.
260    
261        ! First filter, moving average over 9 points.
262      CALL MVA9(zmea)      CALL MVA9(zmea)
263      CALL MVA9(zstd)      CALL MVA9(zstd)
264      CALL MVA9(zpic)      CALL MVA9(zpic)
# Line 264  contains Line 267  contains
267      CALL MVA9(zxtzy)      CALL MVA9(zxtzy)
268      CALL MVA9(zytzy)      CALL MVA9(zytzy)
269    
     ! 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.  
   
270      DO ii = 1, iim      DO ii = 1, iim
271         DO jj = 1, jjm + 1         DO jj = 1, jjm + 1
272            ! Coefficients K, L et M:            ! Coefficients K, L et M:
# Line 279  contains Line 276  contains
276            xp = xk - sqrt(xl**2 + xm**2)            xp = xk - sqrt(xl**2 + xm**2)
277            xq = xk + sqrt(xl**2 + xm**2)            xq = xk + sqrt(xl**2 + xm**2)
278            xw = 1e-8            xw = 1e-8
279            if(xp.le.xw) xp = 0.            if (xp <= xw) xp = 0.
280            if(xq.le.xw) xq = xw            if (xq <= xw) xq = xw
281            if(abs(xm).le.xw) xm = xw * sign(1., xm)            if (abs(xm) <= xw) xm = xw * sign(1., xm)
282            ! modification pour masque de terre fractionnaire            ! modification pour masque de terre fractionnaire
283            ! slope:            ! slope:
284            zsig(ii, jj) = sqrt(xq) * mask_tmp(ii, jj)            zsig(ii, jj) = sqrt(xq) * mask_tmp(ii, jj)
# Line 289  contains Line 286  contains
286            zgam(ii, jj) = xp / xq * mask_tmp(ii, jj)            zgam(ii, jj) = xp / xq * mask_tmp(ii, jj)
287            ! angle theta:            ! angle theta:
288            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)  
289         ENDDO         ENDDO
290      ENDDO      ENDDO
291    
292      print *, 'MEAN ORO: ', zllmmea      zmea(:iim, :) = zmea(:iim, :) * mask_tmp(:iim, :)
293      print *, 'ST. DEV.: ', zllmstd      zpic(:iim, :) = zpic(:iim, :) * mask_tmp(:iim, :)
294      print *, 'PENTE: ', zllmsig      zval(:iim, :) = zval(:iim, :) * mask_tmp(:iim, :)
295      print *, 'ANISOTROP: ', zllmgam      zstd(:iim, :) = zstd(:iim, :) * mask_tmp(:iim, :)
296      print *, 'ANGLE: ', zminthe, zllmthe  
297      print *, 'pic: ', zllmpic      print *, 'MEAN ORO: ', MAXVAL(zmea(:iim, :))
298      print *, 'val: ', zllmval      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      ! gamma and theta at 1. and 0. at poles      ! gamma and theta at 1. and 0. at poles
306      zmea(iim + 1, :) = zmea(1, :)      zmea(iim + 1, :) = zmea(1, :)
# Line 323  contains Line 312  contains
312      zgam(iim + 1, :) = zgam(1, :)      zgam(iim + 1, :) = zgam(1, :)
313      zthe(iim + 1, :) = zthe(1, :)      zthe(iim + 1, :) = zthe(1, :)
314    
315      zmeanor = 0.      zweinor = sum(weight(:iim, 1))
316      zmeasud = 0.      zweisud = sum(weight(:iim, jjm + 1))
317      zstdnor = 0.      zmeanor = sum(zmea(:iim, 1) * weight(:iim, 1))
318      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  
319    
320      zmea(:, 1) = zmeanor / zweinor      zmea(:, 1) = zmeanor / zweinor
321      zmea(:, jjm + 1) = zmeasud / zweisud      zmea(:, jjm + 1) = zmeasud / zweisud
# Line 357  contains Line 323  contains
323      zphi(:, 1) = zmeanor / zweinor      zphi(:, 1) = zmeanor / zweinor
324      zphi(:, jjm + 1) = zmeasud / zweisud      zphi(:, jjm + 1) = zmeasud / zweisud
325    
326      zpic(:, 1) = zpicnor / zweinor      zpic(:, 1) = sum(zpic(:iim, 1) * weight(:iim, 1)) / zweinor
327      zpic(:, jjm + 1) = zpicsud / zweisud      zpic(:, jjm + 1) = sum(zpic(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
328             / zweisud
329      zval(:, 1) = zvalnor / zweinor  
330      zval(:, jjm + 1) = zvalsud / zweisud      zval(:, 1) = sum(zval(:iim, 1) * weight(:iim, 1)) / zweinor
331        zval(:, jjm + 1) = sum(zval(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
332      zstd(:, 1) = zstdnor / zweinor           / zweisud
333      zstd(:, jjm + 1) = zstdsud / zweisud  
334        zstd(:, 1) = sum(zstd(:iim, 1) * weight(:iim, 1)) / zweinor
335      zsig(:, 1) = zsignor / zweinor      zstd(:, jjm + 1) = sum(zstd(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
336      zsig(:, jjm + 1) = zsigsud / zweisud           / zweisud
337    
338        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    
342      zgam(:, 1) = 1.      zgam(:, 1) = 1.
343      zgam(:, jjm + 1) = 1.      zgam(:, jjm + 1) = 1.
Legend:
Removed from v.70  
changed lines
  Added in v.265
  ViewVC Help
Powered by ViewVC 1.1.21