/[lmdze]/trunk/libf/phylmd/Orography/grid_noro_m.f90
ViewVC logotype

Annotation of /trunk/libf/phylmd/Orography/grid_noro_m.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 40 - (hide annotations)
Tue Feb 22 13:49:36 2011 UTC (13 years, 2 months ago) by guez
File size: 14654 byte(s)
"alpha" useless, always 0, in "exner_hyb".

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

  ViewVC Help
Powered by ViewVC 1.1.21