/[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/phylmd/Orography/grid_noro_m.f90 revision 78 by guez, Wed Feb 5 17:51:07 2014 UTC trunk/Sources/phylmd/Orography/grid_noro_m.f revision 227 by guez, Thu Nov 2 15:47:03 2017 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, in m
     REAL, intent(in):: zdata(:, :) ! input field  
40      REAL, intent(in):: x(:), y(:) ! coordinates of output field      REAL, intent(in):: x(:), y(:) ! coordinates of 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, dimension(iim + 1, jjm + 1):: ztz, zxtzx, zytzy, zxtzy, weight
# Line 63  contains Line 72  contains
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, dimension(iim + 1, jjm + 1):: mask_tmp, num_tot, num_lan, zmea0      real, dimension(iim + 1, jjm + 1):: mask_tmp, num_tot, num_lan
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 88  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 216  contains Line 223  contains
223      ! Compute parameters needed by the Lott & Miller (1997) and Lott      ! Compute parameters needed by the Lott & Miller (1997) and Lott
224      ! (1999) subgrid-scale orographic 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)
# Line 248  contains Line 247  contains
247         zytzy(ii, jjm + 1) = zytzy(ii, jjm)         zytzy(ii, jjm + 1) = zytzy(ii, jjm)
248      ENDDO      ENDDO
249    
250      zmea0 = zmea ! not smoothed      ! 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      ! Filters to smooth out fields for input into subgrid-scale      ! Filters to smooth out fields for input into subgrid-scale
259      ! orographic scheme.      ! orographic scheme.
# Line 262  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 = merge(1., 0., mask >= 0.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 276  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 286  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) = zmea0(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 320  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 354  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.78  
changed lines
  Added in v.227
  ViewVC Help
Powered by ViewVC 1.1.21