/[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 47 - (hide annotations)
Fri Jul 1 15:00:48 2011 UTC (12 years, 10 months ago) by guez
File size: 12978 byte(s)
Split "thermcell.f" and "cv3_routines.f".
Removed copies of files that are now in "L_util".
Moved "mva9" and "diagetpq" to their own files.
Unified variable names across procedures.

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

  ViewVC Help
Powered by ViewVC 1.1.21