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

Annotation of /trunk/Sources/phylmd/Orography/grid_noro_m.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 52 - (hide annotations)
Fri Sep 23 12:28:01 2011 UTC (12 years, 9 months ago) by guez
Original Path: trunk/libf/phylmd/Orography/grid_noro_m.f90
File size: 12921 byte(s)
Split "conflx.f" into single-procedure files in directory "Conflx".

Split "cv_routines.f" into single-procedure files in directory
"CV_routines". Made module "cvparam" from included file
"cvparam.h". No included file other than "netcdf.inc" left in LMDZE.

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     ! Authors: F. Lott, Z. X. Li, A. Harzallah and L. Fairhead
13    
14     ! Compute the parameters of the SSO scheme as described in
15     ! Lott and Miller (1997) and Lott (1999).
16     ! Target points are on a rectangular grid:
17     ! jjm+1 latitudes including North and South Poles;
18     ! iim+1 longitudes, with periodicity: longitude(iim+1)=longitude(1)
19     ! At the poles the field value is repeated iim+1 times.
20    
21     ! The parameters a, b, c, d represent the limite of the target
22     ! gridpoint region. The means over this region are calculated
23     ! from USN data, ponderated by a weight proportional to the
24     ! surface occupied by the data inside the model gridpoint area.
25     ! In most circumstances, this weight is the ratio between the
26     ! surface of the USN gridpoint area and the surface of the
27     ! model gridpoint area.
28    
29     ! (c)
30     ! ----d-----
31     ! | . . . .|
32     ! | |
33     ! (b)a . * . .b(a)
34     ! | |
35     ! | . . . .|
36     ! ----c-----
37     ! (d)
38    
39     use dimens_m, only: iim, jjm
40 guez 39 use nr_util, only: assert, pi
41 guez 47 use mva9_m, only: mva9
42 guez 3
43     REAL, intent(in):: xdata(:), ydata(:) ! coordinates of input field
44     REAL, intent(in):: zdata(:, :) ! input field
45     REAL, intent(in):: x(:), y(:) ! ccordinates output field
46    
47     ! Correlations of USN orography gradients:
48    
49     REAL zphi(:, :)
50     real, intent(out):: zmea(:, :) ! Mean orography
51     real, intent(out):: zstd(:, :) ! Standard deviation
52     REAL zsig(:, :) ! Slope
53     real zgam(:, :) ! Anisotropy
54     real zthe(:, :) ! Orientation of the small axis
55     REAL, intent(out):: zpic(:, :) ! Maximum altitude
56     real, intent(out):: zval(:, :) ! Minimum altitude
57    
58 guez 15 real, intent(out):: mask(:, :) ! fraction of land
59 guez 3
60     ! Variables local to the procedure:
61    
62     ! In this version it is assumed that the input data come from
63     ! the US Navy dataset:
64     integer, parameter:: iusn=2160, jusn=1080
65    
66     integer, parameter:: iext=216
67     REAL xusn(iusn+2*iext), yusn(jusn+2)
68     REAL zusn(iusn+2*iext, jusn+2)
69    
70     ! Intermediate fields (correlations of orography gradient)
71    
72     REAL ztz(iim+1, jjm+1), zxtzx(iim+1, jjm+1)
73     REAL zytzy(iim+1, jjm+1), zxtzy(iim+1, jjm+1)
74     REAL weight(iim+1, jjm+1)
75    
76     ! Correlations of USN orography gradients:
77    
78     REAL zxtzxusn(iusn+2*iext, jusn+2), zytzyusn(iusn+2*iext, jusn+2)
79     REAL zxtzyusn(iusn+2*iext, jusn+2)
80    
81     real mask_tmp(size(x), size(y))
82     real num_tot(2200, 1100), num_lan(2200, 1100)
83    
84     REAL a(2200), b(2200), c(1100), d(1100)
85     real rad, weighx, weighy, xincr, xk, xp, xm, xw, xq, xl
86     real zbordnor, zdeltax, zbordsud, zdeltay, zbordoue, zlenx, zleny, zmeasud
87     real zllmpic, zllmmea, zllmgam, zllmthe, zllmstd, zllmsig, zllmval
88     real zpicnor, zminthe, zsigsud, zstdnor, zstdsud, zvalsud, zvalnor
89     real zweinor, zweisud, zsignor, zpicsud, zmeanor, zbordest
90     integer ii, i, jj, j
91    
92     !-------------------------------
93    
94     print *, "Call sequence information: grid_noro"
95    
96     call assert((/size(xdata), size(zdata, 1)/) == iusn, "grid_noro iusn")
97     call assert((/size(ydata), size(zdata, 2)/) == jusn, "grid_noro jusn")
98    
99     call assert((/size(x), size(zphi, 1), size(zmea, 1), size(zstd, 1), &
100     size(zsig, 1), size(zgam, 1), size(zthe, 1), size(zpic, 1), &
101     size(zval, 1), size(mask, 1)/) == iim + 1, "grid_noro iim")
102    
103     call assert((/size(y), size(zphi, 2), size(zmea, 2), size(zstd, 2), &
104     size(zsig, 2), size(zgam, 2), size(zthe, 2), size(zpic, 2), &
105     size(zval, 2), size(mask, 2)/) == jjm + 1, "grid_noro jjm")
106    
107     IF (iim > 2200 .OR. jjm > 1099) THEN
108     print *, "iim = ", iim, ", jjm = ", jjm
109     stop '"iim" or "jjm" is too big'
110     ENDIF
111    
112     print *, "Paramètres de l'orographie à l'échelle sous-maille"
113     rad = 6371229.
114     zdeltay = 2. * pi / real(jusn) * rad
115    
116     ! Extension of the USN database to POCEED computations at boundaries:
117    
118     DO j=1, jusn
119     yusn(j+1)=ydata(j)
120     DO i=1, iusn
121     zusn(i+iext, j+1)=zdata(i, j)
122     xusn(i+iext)=xdata(i)
123     ENDDO
124     DO i=1, iext
125     zusn(i, j+1)=zdata(iusn-iext+i, j)
126     xusn(i)=xdata(iusn-iext+i)-2.*pi
127     zusn(iusn+iext+i, j+1)=zdata(i, j)
128     xusn(iusn+iext+i)=xdata(i)+2.*pi
129     ENDDO
130     ENDDO
131    
132     yusn(1)=ydata(1)+(ydata(1)-ydata(2))
133     yusn(jusn+2)=ydata(jusn)+(ydata(jusn)-ydata(jusn-1))
134     DO i=1, iusn/2+iext
135     zusn(i, 1)=zusn(i+iusn/2, 2)
136     zusn(i+iusn/2+iext, 1)=zusn(i, 2)
137     zusn(i, jusn+2)=zusn(i+iusn/2, jusn+1)
138     zusn(i+iusn/2+iext, jusn+2)=zusn(i, jusn+1)
139     ENDDO
140 guez 15
141 guez 3 ! COMPUTE LIMITS OF MODEL GRIDPOINT AREA
142     ! ( REGULAR GRID)
143    
144     a(1) = x(1) - (x(2)-x(1))/2.0
145     b(1) = (x(1)+x(2))/2.0
146     DO i = 2, iim
147     a(i) = b(i-1)
148     b(i) = (x(i)+x(i+1))/2.0
149     ENDDO
150     a(iim+1) = b(iim)
151     b(iim+1) = x(iim+1) + (x(iim+1)-x(iim))/2.0
152    
153     c(1) = y(1) - (y(2)-y(1))/2.0
154     d(1) = (y(1)+y(2))/2.0
155     DO j = 2, jjm
156     c(j) = d(j-1)
157     d(j) = (y(j)+y(j+1))/2.0
158     ENDDO
159     c(jjm + 1) = d(jjm)
160     d(jjm + 1) = y(jjm + 1) + (y(jjm + 1)-y(jjm))/2.0
161    
162     ! Initialisations :
163     weight(:, :) = 0.
164     zxtzx(:, :) = 0.
165     zytzy(:, :) = 0.
166     zxtzy(:, :) = 0.
167     ztz(:, :) = 0.
168     zmea(:, :) = 0.
169     zpic(:, :) =-1.E+10
170     zval(:, :) = 1.E+10
171    
172     ! COMPUTE SLOPES CORRELATIONS ON USN GRID
173    
174     zytzyusn(:, :)=0.
175     zxtzxusn(:, :)=0.
176     zxtzyusn(:, :)=0.
177    
178     DO j = 2, jusn+1
179     zdeltax=zdeltay*cos(yusn(j))
180     DO i = 2, iusn+2*iext-1
181     zytzyusn(i, j)=(zusn(i, j+1)-zusn(i, j-1))**2/zdeltay**2
182     zxtzxusn(i, j)=(zusn(i+1, j)-zusn(i-1, j))**2/zdeltax**2
183     zxtzyusn(i, j)=(zusn(i, j+1)-zusn(i, j-1))/zdeltay &
184     *(zusn(i+1, j)-zusn(i-1, j))/zdeltax
185     ENDDO
186     ENDDO
187    
188     ! SUMMATION OVER GRIDPOINT AREA
189 guez 15
190 guez 3 zleny=pi/real(jusn)*rad
191     xincr=pi/2./real(jusn)
192     DO ii = 1, iim+1
193     DO jj = 1, jjm + 1
194     num_tot(ii, jj)=0.
195     num_lan(ii, jj)=0.
196     DO j = 2, jusn+1
197     zlenx=zleny*cos(yusn(j))
198     zdeltax=zdeltay*cos(yusn(j))
199     zbordnor=(c(jj)-yusn(j)+xincr)*rad
200     zbordsud=(yusn(j)-d(jj)+xincr)*rad
201     weighy=AMAX1(0., amin1(zbordnor, zbordsud, zleny))
202     IF (weighy /= 0) THEN
203     DO i = 2, iusn+2*iext-1
204     zbordest=(xusn(i)-a(ii)+xincr)*rad*cos(yusn(j))
205     zbordoue=(b(ii)+xincr-xusn(i))*rad*cos(yusn(j))
206     weighx=AMAX1(0., amin1(zbordest, zbordoue, zlenx))
207     IF (weighx /= 0) THEN
208     num_tot(ii, jj) = num_tot(ii, jj) + 1.
209     if (zusn(i, j) >= 1.) then
210     num_lan(ii, jj) = num_lan(ii, jj) + 1.
211     end if
212     weight(ii, jj) = weight(ii, jj) + weighx * weighy
213     zxtzx(ii, jj)=zxtzx(ii, jj)+zxtzxusn(i, j)*weighx*weighy
214     zytzy(ii, jj)=zytzy(ii, jj)+zytzyusn(i, j)*weighx*weighy
215     zxtzy(ii, jj)=zxtzy(ii, jj)+zxtzyusn(i, j)*weighx*weighy
216     ztz(ii, jj) = ztz(ii, jj) &
217     + zusn(i, j) * zusn(i, j) * weighx * weighy
218     ! mean
219     zmea(ii, jj) =zmea(ii, jj)+zusn(i, j)*weighx*weighy
220     ! peacks
221     zpic(ii, jj)=amax1(zpic(ii, jj), zusn(i, j))
222     ! valleys
223     zval(ii, jj)=amin1(zval(ii, jj), zusn(i, j))
224     ENDIF
225     ENDDO
226     ENDIF
227     ENDDO
228     ENDDO
229     ENDDO
230    
231     if (any(weight == 0.)) stop "zero weight in grid_noro"
232    
233     ! COMPUTE PARAMETERS NEEDED BY THE LOTT & MILLER (1997) AND
234     ! LOTT (1999) SSO SCHEME.
235    
236     zllmmea=0.
237     zllmstd=0.
238     zllmsig=0.
239     zllmgam=0.
240     zllmpic=0.
241     zllmval=0.
242     zllmthe=0.
243     zminthe=0.
244     DO ii = 1, iim+1
245     DO jj = 1, jjm + 1
246     mask(ii, jj) = num_lan(ii, jj)/num_tot(ii, jj)
247     ! Mean Orography:
248     zmea (ii, jj)=zmea (ii, jj)/weight(ii, jj)
249     zxtzx(ii, jj)=zxtzx(ii, jj)/weight(ii, jj)
250     zytzy(ii, jj)=zytzy(ii, jj)/weight(ii, jj)
251     zxtzy(ii, jj)=zxtzy(ii, jj)/weight(ii, jj)
252     ztz(ii, jj) =ztz(ii, jj)/weight(ii, jj)
253     ! Standard deviation:
254 guez 40 zstd(ii, jj)=sqrt(MAX(0., ztz(ii, jj) - zmea(ii, jj)**2))
255 guez 3 ENDDO
256     ENDDO
257    
258     ! CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:
259    
260     DO ii = 1, iim+1
261     zxtzx(ii, 1)=zxtzx(ii, 2)
262     zxtzx(ii, jjm + 1)=zxtzx(ii, jjm)
263     zxtzy(ii, 1)=zxtzy(ii, 2)
264     zxtzy(ii, jjm + 1)=zxtzy(ii, jjm)
265     zytzy(ii, 1)=zytzy(ii, 2)
266     zytzy(ii, jjm + 1)=zytzy(ii, jjm)
267     ENDDO
268    
269     ! FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.
270    
271     ! FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.
272    
273 guez 47 CALL MVA9(zmea)
274     CALL MVA9(zstd)
275     CALL MVA9(zpic)
276     CALL MVA9(zval)
277     CALL MVA9(zxtzx)
278     CALL MVA9(zxtzy)
279     CALL MVA9(zytzy)
280 guez 3
281     ! Masque prenant en compte maximum de terre
282     ! On seuil a 10% de terre de terre car en dessous les parametres
283     ! de surface n'ont pas de sens (PB)
284     mask_tmp= 0.
285     WHERE (mask >= 0.1) mask_tmp = 1.
286    
287     DO ii = 1, iim
288     DO jj = 1, jjm + 1
289     IF (weight(ii, jj) /= 0.) THEN
290     ! Coefficients K, L et M:
291     xk=(zxtzx(ii, jj)+zytzy(ii, jj))/2.
292     xl=(zxtzx(ii, jj)-zytzy(ii, jj))/2.
293     xm=zxtzy(ii, jj)
294     xp=xk-sqrt(xl**2+xm**2)
295     xq=xk+sqrt(xl**2+xm**2)
296     xw=1.e-8
297     if(xp.le.xw) xp=0.
298     if(xq.le.xw) xq=xw
299     if(abs(xm).le.xw) xm=xw*sign(1., xm)
300     !$$* PB modif pour maque de terre fractionnaire
301     ! slope:
302     zsig(ii, jj)=sqrt(xq)*mask_tmp(ii, jj)
303     ! isotropy:
304     zgam(ii, jj)=xp/xq*mask_tmp(ii, jj)
305     ! angle theta:
306     zthe(ii, jj)=57.29577951*atan2(xm, xl)/2.*mask_tmp(ii, jj)
307     zphi(ii, jj)=zmea(ii, jj)*mask_tmp(ii, jj)
308     zmea(ii, jj)=zmea(ii, jj)*mask_tmp(ii, jj)
309     zpic(ii, jj)=zpic(ii, jj)*mask_tmp(ii, jj)
310     zval(ii, jj)=zval(ii, jj)*mask_tmp(ii, jj)
311     zstd(ii, jj)=zstd(ii, jj)*mask_tmp(ii, jj)
312     ENDIF
313     zllmmea=AMAX1(zmea(ii, jj), zllmmea)
314     zllmstd=AMAX1(zstd(ii, jj), zllmstd)
315     zllmsig=AMAX1(zsig(ii, jj), zllmsig)
316     zllmgam=AMAX1(zgam(ii, jj), zllmgam)
317     zllmthe=AMAX1(zthe(ii, jj), zllmthe)
318     zminthe=amin1(zthe(ii, jj), zminthe)
319     zllmpic=AMAX1(zpic(ii, jj), zllmpic)
320     zllmval=AMAX1(zval(ii, jj), zllmval)
321     ENDDO
322     ENDDO
323 guez 15 print *, 'MEAN ORO: ', zllmmea
324     print *, 'ST. DEV.: ', zllmstd
325     print *, 'PENTE: ', zllmsig
326     print *, 'ANISOTROP: ', zllmgam
327     print *, 'ANGLE: ', zminthe, zllmthe
328     print *, 'pic: ', zllmpic
329     print *, 'val: ', zllmval
330 guez 3
331     ! gamma and theta a 1. and 0. at poles
332     zmea(iim+1, :)=zmea(1, :)
333     zphi(iim+1, :)=zphi(1, :)
334     zpic(iim+1, :)=zpic(1, :)
335     zval(iim+1, :)=zval(1, :)
336     zstd(iim+1, :)=zstd(1, :)
337     zsig(iim+1, :)=zsig(1, :)
338     zgam(iim+1, :)=zgam(1, :)
339     zthe(iim+1, :)=zthe(1, :)
340    
341     zmeanor=0.
342     zmeasud=0.
343     zstdnor=0.
344     zstdsud=0.
345     zsignor=0.
346     zsigsud=0.
347     zweinor=0.
348     zweisud=0.
349     zpicnor=0.
350     zpicsud=0.
351     zvalnor=0.
352     zvalsud=0.
353    
354     DO ii=1, iim
355     zweinor=zweinor+ weight(ii, 1)
356     zweisud=zweisud+ weight(ii, jjm + 1)
357     zmeanor=zmeanor+zmea(ii, 1)*weight(ii, 1)
358     zmeasud=zmeasud+zmea(ii, jjm + 1)*weight(ii, jjm + 1)
359     zstdnor=zstdnor+zstd(ii, 1)*weight(ii, 1)
360     zstdsud=zstdsud+zstd(ii, jjm + 1)*weight(ii, jjm + 1)
361     zsignor=zsignor+zsig(ii, 1)*weight(ii, 1)
362     zsigsud=zsigsud+zsig(ii, jjm + 1)*weight(ii, jjm + 1)
363     zpicnor=zpicnor+zpic(ii, 1)*weight(ii, 1)
364     zpicsud=zpicsud+zpic(ii, jjm + 1)*weight(ii, jjm + 1)
365     zvalnor=zvalnor+zval(ii, 1)*weight(ii, 1)
366     zvalsud=zvalsud+zval(ii, jjm + 1)*weight(ii, jjm + 1)
367     ENDDO
368    
369     zmea(:, 1)=zmeanor/zweinor
370     zmea(:, jjm + 1)=zmeasud/zweisud
371    
372     zphi(:, 1)=zmeanor/zweinor
373     zphi(:, jjm + 1)=zmeasud/zweisud
374    
375     zpic(:, 1)=zpicnor/zweinor
376     zpic(:, jjm + 1)=zpicsud/zweisud
377    
378     zval(:, 1)=zvalnor/zweinor
379     zval(:, jjm + 1)=zvalsud/zweisud
380    
381     zstd(:, 1)=zstdnor/zweinor
382     zstd(:, jjm + 1)=zstdsud/zweisud
383    
384     zsig(:, 1)=zsignor/zweinor
385     zsig(:, jjm + 1)=zsigsud/zweisud
386    
387     zgam(:, 1)=1.
388     zgam(:, jjm + 1)=1.
389    
390     zthe(:, 1)=0.
391     zthe(:, jjm + 1)=0.
392    
393     END SUBROUTINE grid_noro
394    
395     end module grid_noro_m

  ViewVC Help
Powered by ViewVC 1.1.21