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

Diff of /trunk/phylmd/Orography/grid_noro.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/phylmd/Orography/grid_noro_m.f90 revision 78 by guez, Wed Feb 5 17:51:07 2014 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çois 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
# Line 32  contains Line 34  contains
34    
35      REAL, intent(in):: xdata(:), ydata(:) ! coordinates of input field      REAL, intent(in):: xdata(:), ydata(:) ! coordinates of input field
36      REAL, intent(in):: zdata(:, :) ! input field      REAL, intent(in):: zdata(:, :) ! input field
37      REAL, intent(in):: x(:), y(:) ! ccordinates output field      REAL, intent(in):: x(:), y(:) ! coordinates of output field
   
     ! Correlations of USN orography gradients:  
38    
39        ! Correlations of US Navy orography gradients:
40      REAL, intent(out):: zphi(:, :)      REAL, intent(out):: zphi(:, :)
41      real, intent(out):: zmea(:, :) ! Mean orography      real, intent(out):: zmea(:, :) ! Mean orography
42      real, intent(out):: zstd(:, :) ! Standard deviation      real, intent(out):: zstd(:, :) ! Standard deviation
43      REAL zsig(:, :) ! Slope      REAL, intent(out):: zsig(:, :) ! Slope
44      real zgam(:, :) ! Anisotropy      real, intent(out):: zgam(:, :) ! Anisotropy
45      real zthe(:, :) ! Orientation of the small axis      real, intent(out):: zthe(:, :) ! Orientation of the small axis
46      REAL, intent(out):: zpic(:, :) ! Maximum altitude      REAL, intent(out):: zpic(:, :) ! Maximum altitude
47      real, intent(out):: zval(:, :) ! Minimum altitude      real, intent(out):: zval(:, :) ! Minimum altitude
48    
# Line 57  contains Line 58  contains
58      REAL zusn(iusn + 2 * iext, jusn + 2)      REAL zusn(iusn + 2 * iext, jusn + 2)
59    
60      ! Intermediate fields (correlations of orography gradient)      ! Intermediate fields (correlations of orography gradient)
61        REAL, dimension(iim + 1, jjm + 1):: ztz, zxtzx, zytzy, zxtzy, weight
62    
63      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:  
64      REAL, dimension(iusn + 2 * iext, jusn + 2):: zxtzxusn, zytzyusn, zxtzyusn      REAL, dimension(iusn + 2 * iext, jusn + 2):: zxtzxusn, zytzyusn, zxtzyusn
65    
66      real mask_tmp(size(x), size(y))      real, dimension(iim + 1, jjm + 1):: mask_tmp, num_tot, num_lan, zmea0
     real num_tot(iim + 1, jjm + 1), num_lan(iim + 1, jjm + 1)  
   
67      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)
68      real rad, weighx, weighy, xincr, xk, xp, xm, xw, xq, xl      real weighx, weighy, xincr, xk, xp, xm, xw, xq, xl
69      real zbordnor, zdeltax, zbordsud, zdeltay, zbordoue, zlenx, zleny, zmeasud      real zbordnor, zdeltax, zbordsud, zdeltay, zbordoue, zlenx, zleny, zmeasud
70      real zllmpic, zllmmea, zllmgam, zllmthe, zllmstd, zllmsig, zllmval      real zllmpic, zllmmea, zllmgam, zllmthe, zllmstd, zllmsig, zllmval
71      real zpicnor, zminthe, zsigsud, zstdnor, zstdsud, zvalsud, zvalnor      real zpicnor, zminthe, zsigsud, zstdnor, zstdsud, zvalsud, zvalnor
72      real zweinor, zweisud, zsignor, zpicsud, zmeanor, zbordest      real zweinor, zweisud, zsignor, zpicsud, zmeanor, zbordest
73      integer ii, i, jj, j      integer ii, i, jj, j
74        real, parameter:: rad = 6371229.
75    
76      !-------------------------------      !-------------------------------
77    
# Line 92  contains Line 89  contains
89           size(zval, 2), size(mask, 2)/) == jjm + 1, "grid_noro jjm")           size(zval, 2), size(mask, 2)/) == jjm + 1, "grid_noro jjm")
90    
91      print *, "Paramètres de l'orographie à l'échelle sous-maille"      print *, "Paramètres de l'orographie à l'échelle sous-maille"
     rad = 6371229.  
92      zdeltay = 2. * pi / real(jusn) * rad      zdeltay = 2. * pi / real(jusn) * rad
93    
94      ! Extension of the USN database to POCEED computations at boundaries:      ! Extension of the US Navy database for computations at boundaries:
95    
96      DO j = 1, jusn      DO j = 1, jusn
97         yusn(j + 1) = ydata(j)         yusn(j + 1) = ydata(j)
# Line 120  contains Line 116  contains
116         zusn(i + iusn / 2 + iext, jusn + 2) = zusn(i, jusn + 1)         zusn(i + iusn / 2 + iext, jusn + 2) = zusn(i, jusn + 1)
117      ENDDO      ENDDO
118    
119      ! COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)      ! Compute limits of model gridpoint area (regular grid)
120    
121      a(1) = x(1) - (x(2) - x(1)) / 2.0      a(1) = x(1) - (x(2) - x(1)) / 2.0
122      b(1) = (x(1) + x(2)) / 2.0      b(1) = (x(1) + x(2)) / 2.0
# Line 150  contains Line 146  contains
146      zpic = - 1E10      zpic = - 1E10
147      zval = 1E10      zval = 1E10
148    
149      ! COMPUTE SLOPES CORRELATIONS ON USN GRID      ! Compute slopes correlations on US Navy grid
150    
151      zytzyusn = 0.      zytzyusn = 0.
152      zxtzxusn = 0.      zxtzxusn = 0.
# Line 166  contains Line 162  contains
162         ENDDO         ENDDO
163      ENDDO      ENDDO
164    
165      ! SUMMATION OVER GRIDPOINT AREA      ! Summation over gridpoint area
166    
167      zleny = pi / real(jusn) * rad      zleny = pi / real(jusn) * rad
168      xincr = pi / 2. / real(jusn)      xincr = pi / 2. / real(jusn)
# Line 179  contains Line 175  contains
175               zdeltax = zdeltay * cos(yusn(j))               zdeltax = zdeltay * cos(yusn(j))
176               zbordnor = (c(jj) - yusn(j) + xincr) * rad               zbordnor = (c(jj) - yusn(j) + xincr) * rad
177               zbordsud = (yusn(j) - d(jj) + xincr) * rad               zbordsud = (yusn(j) - d(jj) + xincr) * rad
178               weighy = AMAX1(0., amin1(zbordnor, zbordsud, zleny))               weighy = MAX(0., min(zbordnor, zbordsud, zleny))
179               IF (weighy /= 0) THEN               IF (weighy /= 0) THEN
180                  DO i = 2, iusn + 2 * iext - 1                  DO i = 2, iusn + 2 * iext - 1
181                     zbordest = (xusn(i) - a(ii) + xincr) * rad * cos(yusn(j))                     zbordest = (xusn(i) - a(ii) + xincr) * rad * cos(yusn(j))
182                     zbordoue = (b(ii) + xincr - xusn(i)) * rad * cos(yusn(j))                     zbordoue = (b(ii) + xincr - xusn(i)) * rad * cos(yusn(j))
183                     weighx = AMAX1(0., amin1(zbordest, zbordoue, zlenx))                     weighx = MAX(0., min(zbordest, zbordoue, zlenx))
184                     IF (weighx /= 0) THEN                     IF (weighx /= 0) THEN
185                        num_tot(ii, jj) = num_tot(ii, jj) + 1.                        num_tot(ii, jj) = num_tot(ii, jj) + 1.
186                        if (zusn(i, j) >= 1.) then                        if (zusn(i, j) >= 1.) then
# Line 202  contains Line 198  contains
198                        ! mean                        ! mean
199                        zmea(ii, jj) = zmea(ii, jj) + zusn(i, j) * weighx * weighy                        zmea(ii, jj) = zmea(ii, jj) + zusn(i, j) * weighx * weighy
200                        ! peacks                        ! peacks
201                        zpic(ii, jj) = amax1(zpic(ii, jj), zusn(i, j))                        zpic(ii, jj) = max(zpic(ii, jj), zusn(i, j))
202                        ! valleys                        ! valleys
203                        zval(ii, jj) = amin1(zval(ii, jj), zusn(i, j))                        zval(ii, jj) = min(zval(ii, jj), zusn(i, j))
204                     ENDIF                     ENDIF
205                  ENDDO                  ENDDO
206               ENDIF               ENDIF
# Line 212  contains Line 208  contains
208         ENDDO         ENDDO
209      ENDDO      ENDDO
210    
211      if (any(weight == 0.)) stop "zero weight in grid_noro"      if (any(weight == 0.)) then
212           print *, "zero weight in grid_noro"
213           stop 1
214        end if
215    
216      ! COMPUTE PARAMETERS NEEDED BY THE LOTT & MILLER (1997) AND      ! Compute parameters needed by the Lott & Miller (1997) and Lott
217      ! LOTT (1999) SSO SCHEME.      ! (1999) subgrid-scale orographic scheme.
218    
219      zllmmea = 0.      zllmmea = 0.
220      zllmstd = 0.      zllmstd = 0.
# Line 228  contains Line 227  contains
227      DO ii = 1, iim + 1      DO ii = 1, iim + 1
228         DO jj = 1, jjm + 1         DO jj = 1, jjm + 1
229            mask(ii, jj) = num_lan(ii, jj) / num_tot(ii, jj)            mask(ii, jj) = num_lan(ii, jj) / num_tot(ii, jj)
230            ! Mean Orography:            ! Mean orography:
231            zmea (ii, jj) = zmea (ii, jj) / weight(ii, jj)            zmea (ii, jj) = zmea (ii, jj) / weight(ii, jj)
232            zxtzx(ii, jj) = zxtzx(ii, jj) / weight(ii, jj)            zxtzx(ii, jj) = zxtzx(ii, jj) / weight(ii, jj)
233            zytzy(ii, jj) = zytzy(ii, jj) / weight(ii, jj)            zytzy(ii, jj) = zytzy(ii, jj) / weight(ii, jj)
# Line 239  contains Line 238  contains
238         ENDDO         ENDDO
239      ENDDO      ENDDO
240    
241      ! CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:      ! Correct values of horizontal slope near the poles:
242      DO ii = 1, iim + 1      DO ii = 1, iim + 1
243         zxtzx(ii, 1) = zxtzx(ii, 2)         zxtzx(ii, 1) = zxtzx(ii, 2)
244         zxtzx(ii, jjm + 1) = zxtzx(ii, jjm)         zxtzx(ii, jjm + 1) = zxtzx(ii, jjm)
# Line 249  contains Line 248  contains
248         zytzy(ii, jjm + 1) = zytzy(ii, jjm)         zytzy(ii, jjm + 1) = zytzy(ii, jjm)
249      ENDDO      ENDDO
250    
251      ! FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.      zmea0 = zmea ! not smoothed
252    
253        ! Filters to smooth out fields for input into subgrid-scale
254        ! orographic scheme.
255    
256      ! FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.      ! First filter, moving average over 9 points.
257      CALL MVA9(zmea)      CALL MVA9(zmea)
258      CALL MVA9(zstd)      CALL MVA9(zstd)
259      CALL MVA9(zpic)      CALL MVA9(zpic)
# Line 260  contains Line 262  contains
262      CALL MVA9(zxtzy)      CALL MVA9(zxtzy)
263      CALL MVA9(zytzy)      CALL MVA9(zytzy)
264    
265      ! Masque prenant en compte maximum de terre. On seuille à 10 % de      ! Masque prenant en compte maximum de terre. On met un seuil à 10
266      ! terre car en dessous les paramètres de surface n'ont pas de      ! % de terre car en dessous les paramètres de surface n'ont pas de
267      ! sens.      ! sens.
268      mask_tmp = 0.      mask_tmp = merge(1., 0., mask >= 0.1)
     WHERE (mask >= 0.1) mask_tmp = 1.  
269    
270      DO ii = 1, iim      DO ii = 1, iim
271         DO jj = 1, jjm + 1         DO jj = 1, jjm + 1
# Line 285  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)
289            zphi(ii, jj) = zmea(ii, jj) * mask_tmp(ii, jj)            zphi(ii, jj) = zmea0(ii, jj) * mask_tmp(ii, jj)
290            zmea(ii, jj) = zmea(ii, jj) * mask_tmp(ii, jj)            zmea(ii, jj) = zmea(ii, jj) * mask_tmp(ii, jj)
291            zpic(ii, jj) = zpic(ii, jj) * mask_tmp(ii, jj)            zpic(ii, jj) = zpic(ii, jj) * mask_tmp(ii, jj)
292            zval(ii, jj) = zval(ii, jj) * mask_tmp(ii, jj)            zval(ii, jj) = zval(ii, jj) * mask_tmp(ii, jj)
293            zstd(ii, jj) = zstd(ii, jj) * mask_tmp(ii, jj)            zstd(ii, jj) = zstd(ii, jj) * mask_tmp(ii, jj)
294            zllmmea = AMAX1(zmea(ii, jj), zllmmea)            zllmmea = MAX(zmea(ii, jj), zllmmea)
295            zllmstd = AMAX1(zstd(ii, jj), zllmstd)            zllmstd = MAX(zstd(ii, jj), zllmstd)
296            zllmsig = AMAX1(zsig(ii, jj), zllmsig)            zllmsig = MAX(zsig(ii, jj), zllmsig)
297            zllmgam = AMAX1(zgam(ii, jj), zllmgam)            zllmgam = MAX(zgam(ii, jj), zllmgam)
298            zllmthe = AMAX1(zthe(ii, jj), zllmthe)            zllmthe = MAX(zthe(ii, jj), zllmthe)
299            zminthe = amin1(zthe(ii, jj), zminthe)            zminthe = min(zthe(ii, jj), zminthe)
300            zllmpic = AMAX1(zpic(ii, jj), zllmpic)            zllmpic = MAX(zpic(ii, jj), zllmpic)
301            zllmval = AMAX1(zval(ii, jj), zllmval)            zllmval = MAX(zval(ii, jj), zllmval)
302         ENDDO         ENDDO
303      ENDDO      ENDDO
304    
Legend:
Removed from v.68  
changed lines
  Added in v.78
  ViewVC Help
Powered by ViewVC 1.1.21