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

Annotation of /trunk/phylmd/Orography/grid_noro.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 227 - (hide annotations)
Thu Nov 2 15:47:03 2017 UTC (6 years, 7 months ago) by guez
Original Path: trunk/Sources/phylmd/Orography/grid_noro_m.f
File size: 12450 byte(s) 
Rename phisinit to phis in restart.nc: clearer, same name as Fortran variable.

In aaam_bud, use rlat and rlon from phyetat0_m instead of having these
module variables associated to actual arguments in physiq.

In clmain, too many wind variables make the procedure hard to
understand. Use yu(:knon, 1) and yv(:knon, 1) instead of u1lay(:knon)
and v1lay(:knon). Note that when yu(:knon, 1) and yv(:knon, 1) are
used as actual arguments, they are probably copied to new arrays since
the elements are not contiguous. Rename yu10m to wind10m because this
is the norm of wind vector, not its zonal component. Rename yustar to
ustar. Rename uzon and vmer to u1 and v1 since these are wind
components at first layer and u1 and v1 are the names of corresponding
dummy arguments in stdlevvar.

In clmain, rename yzlev to zlev.

In clmain, screenc, stdlevvar and coefcdrag, remove the code
corresponding to zxli true (not used in LMDZ either).

Subroutine ustarhb becomes a function. Simplifications using the fact
that zx_alf2 = 0 and zx_alf1 = 1 (discarding the possibility to change
this).

In procedure vdif_kcay, remove unused dummy argument plev. Remove
useless computations of sss and sssq.

In clouds_gno, exp(100.) would overflow in single precision. Set
maximum to exp(80.) instead.

In physiq, use u(:, 1) and v(:, 1) as arguments to phytrac instead of
creating ad hoc variables yu1 and yv1.

In stdlevvar, rename dummy argument u_10m to wind10m, following the
corresponding modification in clmain. Simplifications using the fact
that ok_pred = 0 and ok_corr = 1 (discarding the possibility to change
this).
1 guez 3 module grid_noro_m
2    
3     implicit none
4    
5     contains
6    
7     SUBROUTINE grid_noro(xdata, ydata, zdata, x, y, zphi, zmea, zstd, zsig, &
8     zgam, zthe, zpic, zval, mask)
9    
10     ! From dyn3d/grid_noro.F, version 1.1.1.1 2004/05/19 12:53:06
11    
12 guez 147 ! Authors: Fran\c{}cois Lott, Laurent Li, A. Harzallah and Laurent
13 guez 70 ! Fairhead
14 guez 3
15 guez 70 ! Compute the parameters of the sub-grid scale orography scheme as
16     ! described in Lott and Miller (1997) and Lott (1999).
17    
18 guez 3 ! Target points are on a rectangular grid:
19 guez 68 ! jjm + 1 latitudes including North and South Poles;
20     ! iim + 1 longitudes, with periodicity: longitude(iim + 1) = longitude(1)
21     ! At the poles the field value is repeated iim + 1 times.
22 guez 3
23     ! The parameters a, b, c, d represent the limite of the target
24 guez 68 ! gridpoint region. The means over this region are calculated from
25 guez 70 ! US Navy data, ponderated by a weight proportional to the surface
26 guez 68 ! occupied by the data inside the model gridpoint area. In most
27     ! circumstances, this weight is the ratio between the surface of
28 guez 70 ! the US Navy gridpoint area and the surface of the model gridpoint
29 guez 68 ! area. See "grid_noto.txt".
30 guez 3
31     use dimens_m, only: iim, jjm
32 guez 95 use mva9_m, only: mva9
33 guez 39 use nr_util, only: assert, pi
34 guez 3
35 guez 147 ! Coordinates of input field:
36     REAL, intent(in):: xdata(:) ! (iusn)
37     REAL, intent(in):: ydata(:) ! (jusn)
38    
39 guez 227 REAL, intent(in):: zdata(:, :) ! (iusn, jusn) input field, in m
40 guez 78 REAL, intent(in):: x(:), y(:) ! coordinates of output field
41 guez 3
42 guez 70 ! Correlations of US Navy orography gradients:
43 guez 227
44     REAL, intent(out):: zphi(:, :) ! (iim + 1, jjm + 1)
45     ! geoptential height of orography, not smoothed, in m
46    
47 guez 147 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 guez 3
52 guez 147 real, intent(out):: zthe(:, :) ! (iim + 1, jjm + 1)
53     ! Orientation of the small axis
54 guez 3
55 guez 147 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 guez 227 ! Local:
61 guez 3
62     ! In this version it is assumed that the input data come from
63     ! the US Navy dataset:
64 guez 68 integer, parameter:: iusn = 2160, jusn = 1080
65     integer, parameter:: iext = 216
66     REAL xusn(iusn + 2 * iext), yusn(jusn + 2)
67 guez 227 REAL zusn(iusn + 2 * iext, jusn + 2) ! in m
68 guez 3
69 guez 68 ! Intermediate fields (correlations of orography gradient)
70 guez 78 REAL, dimension(iim + 1, jjm + 1):: ztz, zxtzx, zytzy, zxtzy, weight
71 guez 3
72 guez 70 ! Correlations of US Navy orography gradients:
73 guez 68 REAL, dimension(iusn + 2 * iext, jusn + 2):: zxtzxusn, zytzyusn, zxtzyusn
74 guez 3
75 guez 147 real, dimension(iim + 1, jjm + 1):: mask_tmp, num_tot, num_lan
76 guez 68 REAL a(iim + 1), b(iim + 1), c(jjm + 1), d(jjm + 1)
77 guez 70 real weighx, weighy, xincr, xk, xp, xm, xw, xq, xl
78 guez 3 real zbordnor, zdeltax, zbordsud, zdeltay, zbordoue, zlenx, zleny, zmeasud
79 guez 147 real zweinor, zweisud, zmeanor, zbordest
80 guez 3 integer ii, i, jj, j
81 guez 70 real, parameter:: rad = 6371229.
82 guez 3
83     !-------------------------------
84    
85     print *, "Call sequence information: grid_noro"
86    
87     call assert((/size(xdata), size(zdata, 1)/) == iusn, "grid_noro iusn")
88     call assert((/size(ydata), size(zdata, 2)/) == jusn, "grid_noro jusn")
89    
90     call assert((/size(x), size(zphi, 1), size(zmea, 1), size(zstd, 1), &
91     size(zsig, 1), size(zgam, 1), size(zthe, 1), size(zpic, 1), &
92     size(zval, 1), size(mask, 1)/) == iim + 1, "grid_noro iim")
93    
94     call assert((/size(y), size(zphi, 2), size(zmea, 2), size(zstd, 2), &
95     size(zsig, 2), size(zgam, 2), size(zthe, 2), size(zpic, 2), &
96     size(zval, 2), size(mask, 2)/) == jjm + 1, "grid_noro jjm")
97    
98 guez 147 print *, "Parameters of subgrid-scale orography"
99 guez 3 zdeltay = 2. * pi / real(jusn) * rad
100    
101 guez 70 ! Extension of the US Navy database for computations at boundaries:
102 guez 3
103 guez 68 DO j = 1, jusn
104     yusn(j + 1) = ydata(j)
105     DO i = 1, iusn
106     zusn(i + iext, j + 1) = zdata(i, j)
107     xusn(i + iext) = xdata(i)
108 guez 3 ENDDO
109 guez 68 DO i = 1, iext
110     zusn(i, j + 1) = zdata(iusn - iext + i, j)
111     xusn(i) = xdata(iusn - iext + i) - 2. * pi
112     zusn(iusn + iext + i, j + 1) = zdata(i, j)
113     xusn(iusn + iext + i) = xdata(i) + 2. * pi
114 guez 3 ENDDO
115     ENDDO
116    
117 guez 68 yusn(1) = ydata(1) + (ydata(1) - ydata(2))
118     yusn(jusn + 2) = ydata(jusn) + (ydata(jusn) - ydata(jusn - 1))
119     DO i = 1, iusn / 2 + iext
120     zusn(i, 1) = zusn(i + iusn / 2, 2)
121     zusn(i + iusn / 2 + iext, 1) = zusn(i, 2)
122     zusn(i, jusn + 2) = zusn(i + iusn / 2, jusn + 1)
123     zusn(i + iusn / 2 + iext, jusn + 2) = zusn(i, jusn + 1)
124 guez 3 ENDDO
125 guez 15
126 guez 78 ! Compute limits of model gridpoint area (regular grid)
127 guez 3
128 guez 68 a(1) = x(1) - (x(2) - x(1)) / 2.0
129     b(1) = (x(1) + x(2)) / 2.0
130 guez 3 DO i = 2, iim
131 guez 68 a(i) = b(i - 1)
132     b(i) = (x(i) + x(i + 1)) / 2.0
133 guez 3 ENDDO
134 guez 68 a(iim + 1) = b(iim)
135     b(iim + 1) = x(iim + 1) + (x(iim + 1) - x(iim)) / 2.0
136 guez 3
137 guez 68 c(1) = y(1) - (y(2) - y(1)) / 2.0
138     d(1) = (y(1) + y(2)) / 2.0
139 guez 3 DO j = 2, jjm
140 guez 68 c(j) = d(j - 1)
141     d(j) = (y(j) + y(j + 1)) / 2.0
142 guez 3 ENDDO
143     c(jjm + 1) = d(jjm)
144 guez 68 d(jjm + 1) = y(jjm + 1) + (y(jjm + 1) - y(jjm)) / 2.0
145 guez 3
146     ! Initialisations :
147 guez 68 weight = 0.
148     zxtzx = 0.
149     zytzy = 0.
150     zxtzy = 0.
151     ztz = 0.
152     zmea = 0.
153     zpic = - 1E10
154     zval = 1E10
155 guez 3
156 guez 70 ! Compute slopes correlations on US Navy grid
157 guez 3
158 guez 68 zytzyusn = 0.
159     zxtzxusn = 0.
160     zxtzyusn = 0.
161 guez 3
162 guez 68 DO j = 2, jusn + 1
163     zdeltax = zdeltay * cos(yusn(j))
164     DO i = 2, iusn + 2 * iext - 1
165     zytzyusn(i, j) = (zusn(i, j + 1) - zusn(i, j - 1))**2 / zdeltay**2
166     zxtzxusn(i, j) = (zusn(i + 1, j) - zusn(i - 1, j))**2 / zdeltax**2
167     zxtzyusn(i, j) = (zusn(i, j + 1) - zusn(i, j - 1)) / zdeltay &
168     * (zusn(i + 1, j) - zusn(i - 1, j)) / zdeltax
169 guez 3 ENDDO
170     ENDDO
171    
172 guez 78 ! Summation over gridpoint area
173 guez 15
174 guez 68 zleny = pi / real(jusn) * rad
175     xincr = pi / 2. / real(jusn)
176     DO ii = 1, iim + 1
177 guez 3 DO jj = 1, jjm + 1
178 guez 68 num_tot(ii, jj) = 0.
179     num_lan(ii, jj) = 0.
180     DO j = 2, jusn + 1
181     zlenx = zleny * cos(yusn(j))
182     zdeltax = zdeltay * cos(yusn(j))
183     zbordnor = (c(jj) - yusn(j) + xincr) * rad
184     zbordsud = (yusn(j) - d(jj) + xincr) * rad
185 guez 70 weighy = MAX(0., min(zbordnor, zbordsud, zleny))
186 guez 3 IF (weighy /= 0) THEN
187 guez 68 DO i = 2, iusn + 2 * iext - 1
188     zbordest = (xusn(i) - a(ii) + xincr) * rad * cos(yusn(j))
189     zbordoue = (b(ii) + xincr - xusn(i)) * rad * cos(yusn(j))
190 guez 70 weighx = MAX(0., min(zbordest, zbordoue, zlenx))
191 guez 3 IF (weighx /= 0) THEN
192     num_tot(ii, jj) = num_tot(ii, jj) + 1.
193     if (zusn(i, j) >= 1.) then
194     num_lan(ii, jj) = num_lan(ii, jj) + 1.
195     end if
196     weight(ii, jj) = weight(ii, jj) + weighx * weighy
197 guez 68 zxtzx(ii, jj) = zxtzx(ii, jj) &
198     + zxtzxusn(i, j) * weighx * weighy
199     zytzy(ii, jj) = zytzy(ii, jj) &
200     + zytzyusn(i, j) * weighx * weighy
201     zxtzy(ii, jj) = zxtzy(ii, jj) &
202     + zxtzyusn(i, j) * weighx * weighy
203 guez 3 ztz(ii, jj) = ztz(ii, jj) &
204     + zusn(i, j) * zusn(i, j) * weighx * weighy
205     ! mean
206 guez 68 zmea(ii, jj) = zmea(ii, jj) + zusn(i, j) * weighx * weighy
207 guez 3 ! peacks
208 guez 70 zpic(ii, jj) = max(zpic(ii, jj), zusn(i, j))
209 guez 3 ! valleys
210 guez 70 zval(ii, jj) = min(zval(ii, jj), zusn(i, j))
211 guez 3 ENDIF
212     ENDDO
213     ENDIF
214     ENDDO
215     ENDDO
216     ENDDO
217    
218 guez 70 if (any(weight == 0.)) then
219     print *, "zero weight in grid_noro"
220     stop 1
221     end if
222 guez 3
223 guez 78 ! Compute parameters needed by the Lott & Miller (1997) and Lott
224     ! (1999) subgrid-scale orographic scheme.
225 guez 3
226 guez 68 DO ii = 1, iim + 1
227 guez 3 DO jj = 1, jjm + 1
228 guez 68 mask(ii, jj) = num_lan(ii, jj) / num_tot(ii, jj)
229 guez 78 ! Mean orography:
230 guez 68 zmea (ii, jj) = zmea (ii, jj) / weight(ii, jj)
231     zxtzx(ii, jj) = zxtzx(ii, jj) / weight(ii, jj)
232     zytzy(ii, jj) = zytzy(ii, jj) / weight(ii, jj)
233     zxtzy(ii, jj) = zxtzy(ii, jj) / weight(ii, jj)
234     ztz(ii, jj) = ztz(ii, jj) / weight(ii, jj)
235     ! Standard deviation:
236     zstd(ii, jj) = sqrt(MAX(0., ztz(ii, jj) - zmea(ii, jj)**2))
237 guez 3 ENDDO
238     ENDDO
239    
240 guez 78 ! Correct values of horizontal slope near the poles:
241 guez 68 DO ii = 1, iim + 1
242     zxtzx(ii, 1) = zxtzx(ii, 2)
243     zxtzx(ii, jjm + 1) = zxtzx(ii, jjm)
244     zxtzy(ii, 1) = zxtzy(ii, 2)
245     zxtzy(ii, jjm + 1) = zxtzy(ii, jjm)
246     zytzy(ii, 1) = zytzy(ii, 2)
247     zytzy(ii, jjm + 1) = zytzy(ii, jjm)
248 guez 3 ENDDO
249    
250 guez 147 ! 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 guez 3
255 guez 147 zphi(:iim, :) = zmea(:iim, :) * mask_tmp(:iim, :)
256     ! (zmea is not yet smoothed)
257    
258 guez 78 ! Filters to smooth out fields for input into subgrid-scale
259     ! orographic scheme.
260    
261     ! First filter, moving average over 9 points.
262 guez 47 CALL MVA9(zmea)
263     CALL MVA9(zstd)
264     CALL MVA9(zpic)
265     CALL MVA9(zval)
266     CALL MVA9(zxtzx)
267     CALL MVA9(zxtzy)
268     CALL MVA9(zytzy)
269 guez 3
270     DO ii = 1, iim
271     DO jj = 1, jjm + 1
272 guez 68 ! Coefficients K, L et M:
273     xk = (zxtzx(ii, jj) + zytzy(ii, jj)) / 2.
274     xl = (zxtzx(ii, jj) - zytzy(ii, jj)) / 2.
275     xm = zxtzy(ii, jj)
276     xp = xk - sqrt(xl**2 + xm**2)
277     xq = xk + sqrt(xl**2 + xm**2)
278     xw = 1e-8
279 guez 147 if (xp <= xw) xp = 0.
280     if (xq <= xw) xq = xw
281     if (abs(xm) <= xw) xm = xw * sign(1., xm)
282 guez 68 ! modification pour masque de terre fractionnaire
283     ! slope:
284     zsig(ii, jj) = sqrt(xq) * mask_tmp(ii, jj)
285     ! isotropy:
286     zgam(ii, jj) = xp / xq * mask_tmp(ii, jj)
287     ! angle theta:
288     zthe(ii, jj) = 57.29577951 * atan2(xm, xl) / 2. * mask_tmp(ii, jj)
289 guez 3 ENDDO
290     ENDDO
291 guez 68
292 guez 147 zmea(:iim, :) = zmea(:iim, :) * mask_tmp(:iim, :)
293     zpic(:iim, :) = zpic(:iim, :) * mask_tmp(:iim, :)
294     zval(:iim, :) = zval(:iim, :) * mask_tmp(:iim, :)
295     zstd(:iim, :) = zstd(:iim, :) * mask_tmp(:iim, :)
296 guez 3
297 guez 147 print *, 'MEAN ORO: ', MAXVAL(zmea(:iim, :))
298     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 guez 68 ! gamma and theta at 1. and 0. at poles
306     zmea(iim + 1, :) = zmea(1, :)
307     zphi(iim + 1, :) = zphi(1, :)
308     zpic(iim + 1, :) = zpic(1, :)
309     zval(iim + 1, :) = zval(1, :)
310     zstd(iim + 1, :) = zstd(1, :)
311     zsig(iim + 1, :) = zsig(1, :)
312     zgam(iim + 1, :) = zgam(1, :)
313     zthe(iim + 1, :) = zthe(1, :)
314 guez 3
315 guez 147 zweinor = sum(weight(:iim, 1))
316     zweisud = sum(weight(:iim, jjm + 1))
317     zmeanor = sum(zmea(:iim, 1) * weight(:iim, 1))
318     zmeasud = sum(zmea(:iim, jjm + 1) * weight(:iim, jjm + 1))
319 guez 3
320 guez 68 zmea(:, 1) = zmeanor / zweinor
321     zmea(:, jjm + 1) = zmeasud / zweisud
322 guez 3
323 guez 68 zphi(:, 1) = zmeanor / zweinor
324     zphi(:, jjm + 1) = zmeasud / zweisud
325 guez 3
326 guez 147 zpic(:, 1) = sum(zpic(:iim, 1) * weight(:iim, 1)) / zweinor
327     zpic(:, jjm + 1) = sum(zpic(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
328     / zweisud
329 guez 3
330 guez 147 zval(:, 1) = sum(zval(:iim, 1) * weight(:iim, 1)) / zweinor
331     zval(:, jjm + 1) = sum(zval(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
332     / zweisud
333 guez 3
334 guez 147 zstd(:, 1) = sum(zstd(:iim, 1) * weight(:iim, 1)) / zweinor
335     zstd(:, jjm + 1) = sum(zstd(:iim, jjm + 1) * weight(:iim, jjm + 1)) &
336     / zweisud
337 guez 3
338 guez 147 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 guez 3
342 guez 68 zgam(:, 1) = 1.
343     zgam(:, jjm + 1) = 1.
344 guez 3
345 guez 68 zthe(:, 1) = 0.
346     zthe(:, jjm + 1) = 0.
347 guez 3
348     END SUBROUTINE grid_noro
349    
350     end module grid_noro_m
  ViewVC Help
Powered by ViewVC 1.1.21