source: LMDZ.3.3/branches/LF/libf/dyn3d/grid_noro.F @ 5456

Last change on this file since 5456 was 2, checked in by lmdz, 25 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 13.8 KB
Line 
1
2      SUBROUTINE grid_noro(imdep, jmdep, xdata, ydata, zdata,
3     .             imar, jmar, x, y,
4     .             zphi,zmea,zstd,zsig,zgam,zthe,
5     .             zpic,zval,mask)
6c=======================================================================
7c (F. Lott) (voir aussi z.x. Li, A. Harzallah et L. Fairhead)
8c
9c      Compute the Parameters of the SSO scheme as described in
10c      LOTT & MILLER (1997) and LOTT(1999).
11c      Target points are on a rectangular grid:
12c      iim+1 latitudes including North and South Poles;
13c      jjm+1 longitudes, with periodicity jjm+1=1.
14c      aux poles.  At the poles the fields value is repeated
15c      jjm+1 time.
16c      The parameters a,b,c,d represent the limite of the target
17c      gridpoint region. The means over this region are calculated
18c      from USN data, ponderated by a weight proportional to the
19c      surface occupated by the data inside the model gridpoint area.
20c      In most circumstances, this weight is the ratio between the
21c      surface of the USN gridpoint area and the surface of the
22c      model gridpoint area.
23c
24c           (c)
25c        ----d-----
26c        | . . . .|
27c        |        |
28c     (b)a . * . .b(a)
29c        |        |
30c        | . . . .|
31c        ----c-----
32c           (d)
33C=======================================================================
34c INPUT:
35c        imdep, jmdep: dimensions X and Y input field
36c        xdata, ydata: coordinates X and Y input field
37c        zdata: Input field
38c        In this version it is assumed that the entry data come from
39c        the USNavy dataset: imdep=iusn=2160, jmdep=jusn=1080.
40c OUTPUT:
41c        imar, jmar: dimensions X and Y Output field
42c        x, y: ccordinates  X and Y Output field.
43c             zmea:  Mean orographie   
44c             zstd:  Standard deviation
45c             zsig:  Slope
46c             zgam:  Anisotropy
47c             zthe:  Orientation of the small axis
48c             zpic:  Maximum altitude
49c             zval:  Minimum altitude
50C=======================================================================
51      IMPLICIT INTEGER (I,J)
52      IMPLICIT REAL(X,Z)
53     
54          parameter(iusn=2160,jusn=1080,iext=216)
55#include "dimensions.h"
56          REAL xusn(iusn+2*iext),yusn(jusn+2)   
57      REAL zusn(iusn+2*iext,jusn+2)
58
59      INTEGER imdep, jmdep
60      REAL xdata(imdep),ydata(jmdep)
61      REAL zdata(imdep,jmdep)
62c
63      INTEGER imar, jmar
64 
65C INTERMEDIATE FIELDS  (CORRELATIONS OF OROGRAPHY GRADIENT)
66
67      REAL ztz(iim+1,jjm+1),zxtzx(iim+1,jjm+1)
68      REAL zytzy(iim+1,jjm+1),zxtzy(iim+1,jjm+1)
69      REAL weight(iim+1,jjm+1)
70
71C CORRELATIONS OF USN OROGRAPHY GRADIENTS
72
73      REAL zxtzxusn(iusn+2*iext,jusn+2),zytzyusn(iusn+2*iext,jusn+2)
74      REAL zxtzyusn(iusn+2*iext,jusn+2)
75      REAL x(imar+1),y(jmar),zphi(imar+1,jmar)
76      REAL zmea(imar+1,jmar),zstd(imar+1,jmar)
77      REAL zsig(imar+1,jmar),zgam(imar+1,jmar),zthe(imar+1,jmar)
78      REAL zpic(imar+1,jmar),zval(imar+1,jmar)
79      integer mask(imar+1,jmar)
80      real num_tot(2200,1100),num_lan(2200,1100)
81c
82      REAL a(2200),b(2200),c(1100),d(1100)
83c
84      print *,' parametres de l orographie a l echelle sous maille'
85      xpi=acos(-1.)
86      rad    = 6 371 229.
87      zdeltay=2.*xpi/float(jusn)*rad
88c
89c  quelques tests de dimensions:
90c   
91c
92      if(iim.ne.imar) STOP 'Problem dim. x'
93      if(jjm.ne.jmar-1) STOP 'Problem dim. y'
94      IF (imar.GT.2200 .OR. jmar.GT.1100) THEN
95         PRINT*, 'imar or jmar too big', imar, jmar
96         CALL ABORT
97      ENDIF
98
99      IF(imdep.ne.iusn.or.jmdep.ne.jusn)then
100         print *,' imdep or jmdep bad dimensions:',imdep,jmdep
101         call abort
102      ENDIF
103
104      IF(imar+1.ne.iim+1.or.jmar.ne.jjm+1)THEN
105        print *,' imar or jmar bad dimensions:',imar,jmar
106        call abort
107      ENDIF
108
109
110c      print *,'xdata:',xdata
111c      print *,'ydata:',ydata
112c      print *,'x:',x
113c      print *,'y:',y
114c
115C  EXTENSION OF THE USN DATABASE TO POCEED COMPUTATIONS AT
116C  BOUNDARIES:
117c
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.*xpi
127        zusn(iusn+iext+i,j+1)=zdata(i,j)
128        xusn(iusn+iext+i)=xdata(i)+2.*xpi
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
141c COMPUTE LIMITS OF MODEL GRIDPOINT AREA
142C     ( REGULAR GRID)
143c
144      a(1) = x(1) - (x(2)-x(1))/2.0
145      b(1) = (x(1)+x(2))/2.0
146      DO i = 2, imar
147         a(i) = b(i-1)
148         b(i) = (x(i)+x(i+1))/2.0
149      ENDDO
150      a(imar+1) = b(imar)
151      b(imar+1) = x(imar+1) + (x(imar+1)-x(imar))/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, jmar-1
156         c(j) = d(j-1)
157         d(j) = (y(j)+y(j+1))/2.0
158      ENDDO
159      c(jmar) = d(jmar-1)
160      d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0
161c
162c  initialisations:
163c
164      DO i = 1, imar+1
165      DO j = 1, jmar
166         weight(i,j) = 0.0
167         zxtzx(i,j)  = 0.0
168         zytzy(i,j)  = 0.0
169         zxtzy(i,j)  = 0.0
170         ztz(i,j)    = 0.0
171         zmea(i,j)   = 0.0
172         zpic(i,j)  =-1.E+10
173         zval(i,j)  = 1.E+10
174      ENDDO
175      ENDDO
176c
177c  COMPUTE SLOPES CORRELATIONS ON USN GRID
178c
179         DO j = 1,jusn+2
180         DO i = 1, iusn+2*iext
181            zytzyusn(i,j)=0.0
182            zxtzxusn(i,j)=0.0
183            zxtzyusn(i,j)=0.0
184         ENDDO
185         ENDDO
186
187
188         DO j = 2,jusn+1
189            zdeltax=zdeltay*cos(yusn(j))
190         DO i = 2, iusn+2*iext-1
191            zytzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))**2/zdeltay**2
192            zxtzxusn(i,j)=(zusn(i+1,j)-zusn(i-1,j))**2/zdeltax**2
193            zxtzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))/zdeltay
194     *                   *(zusn(i+1,j)-zusn(i-1,j))/zdeltax
195         ENDDO
196         ENDDO
197c
198c  SUMMATION OVER GRIDPOINT AREA
199c
200      zleny=xpi/float(jusn)*rad
201      xincr=xpi/2./float(jusn)
202       DO ii = 1, imar+1
203       DO jj = 1, jmar
204       num_tot(ii,jj)=0.
205       num_lan(ii,jj)=0.
206c        PRINT *,' iteration ii jj:',ii,jj
207         DO j = 2,jusn+1
208c         DO j = 3,jusn
209            zlenx=zleny*cos(yusn(j))
210            zdeltax=zdeltay*cos(yusn(j))
211            zbordnor=(c(jj)-yusn(j)+xincr)*rad
212            zbordsud=(yusn(j)-d(jj)+xincr)*rad
213            weighy=AMAX1(0.,
214     *             amin1(zbordnor,zbordsud,zleny))
215         IF(weighy.ne.0)THEN
216         DO i = 2, iusn+2*iext-1
217            zbordest=(xusn(i)-a(ii)+xincr)*rad*cos(yusn(j))
218            zbordoue=(b(ii)+xincr-xusn(i))*rad*cos(yusn(j))
219            weighx=AMAX1(0.,
220     *             amin1(zbordest,zbordoue,zlenx))
221            IF(weighx.ne.0)THEN
222            num_tot(ii,jj)=num_tot(ii,jj)+1.0
223            if(zusn(i,j).ge.1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
224            weight(ii,jj)=weight(ii,jj)+weighx*weighy
225            zxtzx(ii,jj)=zxtzx(ii,jj)+zxtzxusn(i,j)*weighx*weighy
226            zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)*weighx*weighy
227            zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)*weighx*weighy
228            ztz(ii,jj)  =ztz(ii,jj)  +zusn(i,j)*zusn(i,j)*weighx*weighy
229c mean
230            zmea(ii,jj) =zmea(ii,jj)+zusn(i,j)*weighx*weighy
231c peacks
232            zpic(ii,jj)=amax1(zpic(ii,jj),zusn(i,j))
233c valleys
234            zval(ii,jj)=amin1(zval(ii,jj),zusn(i,j))
235            ENDIF
236         ENDDO
237         ENDIF
238         ENDDO
239       ENDDO
240       ENDDO
241c
242c  COMPUTE PARAMETERS NEEDED BY THE LOTT & MILLER (1997) AND
243C  LOTT (1999) SSO SCHEME.
244c
245      zllmmea=0.
246      zllmstd=0.
247      zllmsig=0.
248      zllmgam=0.
249      zllmpic=0.
250      zllmval=0.
251      zllmthe=0.
252      zminthe=0.
253c     print 100,' '
254c100  format(1X,A1,'II JJ',4X,'H',8X,'SD',8X,'SI',3X,'GA',3X,'TH')
255       DO ii = 1, imar+1
256       DO jj = 1, jmar
257         IF (weight(ii,jj) .NE. 0.0) THEN
258c  Mask
259           if(num_lan(ii,jj)/num_tot(ii,jj).ge.0.5)then
260             mask(ii,jj)=1
261           else
262             mask(ii,jj)=0
263           endif
264c  Mean Orography:
265           zmea (ii,jj)=zmea (ii,jj)/weight(ii,jj)
266           zxtzx(ii,jj)=zxtzx(ii,jj)/weight(ii,jj)
267           zytzy(ii,jj)=zytzy(ii,jj)/weight(ii,jj)
268           zxtzy(ii,jj)=zxtzy(ii,jj)/weight(ii,jj)
269           ztz(ii,jj)  =ztz(ii,jj)/weight(ii,jj)
270c  Standard deviation:
271           zstd(ii,jj)=sqrt(AMAX1(0.,ztz(ii,jj)-zmea(ii,jj)**2))
272         ELSE
273            PRINT*, 'probleme,ii,jj=', ii,jj
274         ENDIF
275       ENDDO
276       ENDDO
277
278C CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:
279
280       DO ii = 1, imar+1
281         zxtzx(ii,1)=zxtzx(ii,2)
282         zxtzx(ii,jmar)=zxtzx(ii,jmar-1)
283         zxtzy(ii,1)=zxtzy(ii,2)
284         zxtzy(ii,jmar)=zxtzy(ii,jmar-1)
285         zytzy(ii,1)=zytzy(ii,2)
286         zytzy(ii,jmar)=zytzy(ii,jmar-1)
287       ENDDO
288
289C  FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.
290
291C  FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.
292
293       CALL MVA9(zmea,iim+1,jjm+1)
294       CALL MVA9(zstd,iim+1,jjm+1)
295       CALL MVA9(zpic,iim+1,jjm+1)
296       CALL MVA9(zval,iim+1,jjm+1)
297       CALL MVA9(zxtzx,iim+1,jjm+1)
298       CALL MVA9(zxtzy,iim+1,jjm+1)
299       CALL MVA9(zytzy,iim+1,jjm+1)
300
301       DO ii = 1, imar
302       DO jj = 1, jmar
303         IF (weight(ii,jj) .NE. 0.0) THEN
304c  Coefficients K, L et M:
305           xk=(zxtzx(ii,jj)+zytzy(ii,jj))/2.
306           xl=(zxtzx(ii,jj)-zytzy(ii,jj))/2.
307           xm=zxtzy(ii,jj)
308           xp=xk-sqrt(xl**2+xm**2)
309           xq=xk+sqrt(xl**2+xm**2)
310           xw=1.e-8
311           if(xp.le.xw) xp=0.
312           if(xq.le.xw) xq=xw
313           if(abs(xm).le.xw) xm=xw*sign(1.,xm)
314c slope:
315           zsig(ii,jj)=sqrt(xq)*mask(ii,jj)
316c isotropy:
317           zgam(ii,jj)=xp/xq*mask(ii,jj)
318c angle theta:
319           zthe(ii,jj)=57.29577951*atan2(xm,xl)/2.*mask(ii,jj)
320           zphi(ii,jj)=zmea(ii,jj)*mask(ii,jj)
321           zmea(ii,jj)=zmea(ii,jj)*mask(ii,jj)
322           zpic(ii,jj)=zpic(ii,jj)*mask(ii,jj)
323           zval(ii,jj)=zval(ii,jj)*mask(ii,jj)
324           zstd(ii,jj)=zstd(ii,jj)*mask(ii,jj)
325
326c          print 101,ii,jj,
327c    *           zmea(ii,jj),zstd(ii,jj),zsig(ii,jj),zgam(ii,jj),
328c    *           zthe(ii,jj)
329c101  format(1x,2(1x,i2),2(1x,f7.1),1x,f7.4,2x,f4.2,1x,f5.1)     
330         ELSE
331c           PRINT*, 'probleme,ii,jj=', ii,jj
332         ENDIF
333      zllmmea=AMAX1(zmea(ii,jj),zllmmea)
334      zllmstd=AMAX1(zstd(ii,jj),zllmstd)
335      zllmsig=AMAX1(zsig(ii,jj),zllmsig)
336      zllmgam=AMAX1(zgam(ii,jj),zllmgam)
337      zllmthe=AMAX1(zthe(ii,jj),zllmthe)
338      zminthe=amin1(zthe(ii,jj),zminthe)
339      zllmpic=AMAX1(zpic(ii,jj),zllmpic)
340      zllmval=AMAX1(zval(ii,jj),zllmval)
341       ENDDO
342       ENDDO
343
344      print *,'  MEAN ORO:',zllmmea
345      print *,'  ST. DEV.:',zllmstd
346      print *,'  PENTE:',zllmsig
347      print *,' ANISOTROP:',zllmgam
348      print *,'  ANGLE:',zminthe,zllmthe       
349      print *,'  pic:',zllmpic
350      print *,'  val:',zllmval
351     
352C
353c gamma and theta a 1. and 0. at poles
354c
355      DO jj=1,jmar
356      zmea(imar+1,jj)=zmea(1,jj)
357      zpic(imar+1,jj)=zpic(1,jj)
358      zval(imar+1,jj)=zval(1,jj)
359      zstd(imar+1,jj)=zstd(1,jj)
360      zsig(imar+1,jj)=zsig(1,jj)
361      zgam(imar+1,jj)=zgam(1,jj)
362      zthe(imar+1,jj)=zthe(1,jj)
363      ENDDO
364
365
366      zmeanor=0.0
367      zmeasud=0.0
368      zstdnor=0.0
369      zstdsud=0.0
370      zsignor=0.0
371      zsigsud=0.0
372      zweinor=0.0
373      zweisud=0.0
374      zpicnor=0.0
375      zpicsud=0.0                                   
376      zvalnor=0.0
377      zvalsud=0.0
378
379      DO ii=1,imar
380      zweinor=zweinor+              weight(ii,   1)
381      zweisud=zweisud+              weight(ii,jmar)
382      zmeanor=zmeanor+zmea(ii,   1)*weight(ii,   1)
383      zmeasud=zmeasud+zmea(ii,jmar)*weight(ii,jmar)
384      zstdnor=zstdnor+zstd(ii,   1)*weight(ii,   1)
385      zstdsud=zstdsud+zstd(ii,jmar)*weight(ii,jmar)
386      zsignor=zsignor+zsig(ii,   1)*weight(ii,   1)
387      zsigsud=zsigsud+zsig(ii,jmar)*weight(ii,jmar)
388      zpicnor=zpicnor+zpic(ii,   1)*weight(ii,   1)
389      zpicsud=zpicsud+zpic(ii,jmar)*weight(ii,jmar)
390      zvalnor=zvalnor+zval(ii,   1)*weight(ii,   1)
391      zvalsud=zvalsud+zval(ii,jmar)*weight(ii,jmar)
392      ENDDO
393
394      DO ii=1,imar+1
395      zmea(ii,   1)=zmeanor/zweinor
396      zmea(ii,jmar)=zmeasud/zweisud
397      zphi(ii,   1)=zmeanor/zweinor
398      zphi(ii,jmar)=zmeasud/zweisud
399      zpic(ii,   1)=zpicnor/zweinor
400      zpic(ii,jmar)=zpicsud/zweisud
401      zval(ii,   1)=zvalnor/zweinor
402      zval(ii,jmar)=zvalsud/zweisud
403      zstd(ii,   1)=zstdnor/zweinor
404      zstd(ii,jmar)=zstdsud/zweisud
405      zsig(ii,   1)=zsignor/zweinor
406      zsig(ii,jmar)=zsigsud/zweisud
407      zgam(ii,   1)=1.
408      zgam(ii,jmar)=1.
409      zthe(ii,   1)=0.
410      zthe(ii,jmar)=0.
411      ENDDO
412
413      RETURN
414      END
415
416      SUBROUTINE MVA9(X,IMAR,JMAR)
417
418C MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
419
420      PARAMETER (ISMo=300,JSMo=200)
421      REAL X(IMAR,JMAR),XF(ISMo,JSMo)
422      real weight(-1:1,-1:1)
423
424      if(imar.gt.ismo) stop'surdimensionner ismo dans mva9 (grid_noro)'
425      if(jmar.gt.jsmo) stop'surdimensionner jsmo dans mva9 (grid_noro)'
426     
427      SUM=0.
428      DO IS=-1,1
429      DO JS=-1,1
430      WEIGHT(IS,JS)=1./FLOAT((1+IS**2)*(1+JS**2))
431      SUM=SUM+WEIGHT(IS,JS)
432      ENDDO
433      ENDDO
434
435      DO IS=-1,1
436      DO JS=-1,1
437      WEIGHT(IS,JS)=WEIGHT(IS,JS)/SUM
438      ENDDO
439      ENDDO
440
441      DO J=2,JMAR-1
442      DO I=2,IMAR-1
443      XF(I,J)=0.
444      DO IS=-1,1
445      DO JS=-1,1
446         XF(I,J)=XF(I,J)+X(I+IS,J+JS)*WEIGHT(IS,JS)
447      ENDDO
448      ENDDO
449      ENDDO
450      ENDDO
451
452      DO J=2,JMAR-1
453         XF(1,J)=0.
454      IS=IMAR-1
455      DO JS=-1,1
456         XF(1,J)=XF(1,J)+X(IS,J+JS)*WEIGHT(-1,JS)
457      ENDDO
458      DO IS=0,1
459      DO JS=-1,1
460         XF(1,J)=XF(1,J)+X(1+IS,J+JS)*WEIGHT(IS,JS)
461      ENDDO
462      ENDDO
463         XF(IMAR,J)=XF(1,J)
464      ENDDO
465
466      DO I=1,IMAR
467      XF(I,1)=XF(I,2)
468      XF(I,JMAR)=XF(I,JMAR-1)
469      ENDDO
470
471      DO I=1,IMAR
472      DO J=1,JMAR
473         X(I,J)=XF(I,J)
474      ENDDO
475      ENDDO
476
477      RETURN
478      END
479
480
481
Note: See TracBrowser for help on using the repository browser.