source: LMDZ.3.3/trunk/libf/dyn3d/grid_noro.F @ 3671

Last change on this file since 3671 was 7, checked in by lmdzadmin, 25 years ago

Initialisation zphi pour assurer periodicite

  • 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      zphi(imar+1,jj)=zphi(1,jj)
358      zpic(imar+1,jj)=zpic(1,jj)
359      zval(imar+1,jj)=zval(1,jj)
360      zstd(imar+1,jj)=zstd(1,jj)
361      zsig(imar+1,jj)=zsig(1,jj)
362      zgam(imar+1,jj)=zgam(1,jj)
363      zthe(imar+1,jj)=zthe(1,jj)
364      ENDDO
365
366
367      zmeanor=0.0
368      zmeasud=0.0
369      zstdnor=0.0
370      zstdsud=0.0
371      zsignor=0.0
372      zsigsud=0.0
373      zweinor=0.0
374      zweisud=0.0
375      zpicnor=0.0
376      zpicsud=0.0                                   
377      zvalnor=0.0
378      zvalsud=0.0
379
380      DO ii=1,imar
381      zweinor=zweinor+              weight(ii,   1)
382      zweisud=zweisud+              weight(ii,jmar)
383      zmeanor=zmeanor+zmea(ii,   1)*weight(ii,   1)
384      zmeasud=zmeasud+zmea(ii,jmar)*weight(ii,jmar)
385      zstdnor=zstdnor+zstd(ii,   1)*weight(ii,   1)
386      zstdsud=zstdsud+zstd(ii,jmar)*weight(ii,jmar)
387      zsignor=zsignor+zsig(ii,   1)*weight(ii,   1)
388      zsigsud=zsigsud+zsig(ii,jmar)*weight(ii,jmar)
389      zpicnor=zpicnor+zpic(ii,   1)*weight(ii,   1)
390      zpicsud=zpicsud+zpic(ii,jmar)*weight(ii,jmar)
391      zvalnor=zvalnor+zval(ii,   1)*weight(ii,   1)
392      zvalsud=zvalsud+zval(ii,jmar)*weight(ii,jmar)
393      ENDDO
394
395      DO ii=1,imar+1
396      zmea(ii,   1)=zmeanor/zweinor
397      zmea(ii,jmar)=zmeasud/zweisud
398      zphi(ii,   1)=zmeanor/zweinor
399      zphi(ii,jmar)=zmeasud/zweisud
400      zpic(ii,   1)=zpicnor/zweinor
401      zpic(ii,jmar)=zpicsud/zweisud
402      zval(ii,   1)=zvalnor/zweinor
403      zval(ii,jmar)=zvalsud/zweisud
404      zstd(ii,   1)=zstdnor/zweinor
405      zstd(ii,jmar)=zstdsud/zweisud
406      zsig(ii,   1)=zsignor/zweinor
407      zsig(ii,jmar)=zsigsud/zweisud
408      zgam(ii,   1)=1.
409      zgam(ii,jmar)=1.
410      zthe(ii,   1)=0.
411      zthe(ii,jmar)=0.
412      ENDDO
413
414      RETURN
415      END
416
417      SUBROUTINE MVA9(X,IMAR,JMAR)
418
419C MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
420
421      PARAMETER (ISMo=300,JSMo=200)
422      REAL X(IMAR,JMAR),XF(ISMo,JSMo)
423      real weight(-1:1,-1:1)
424
425      if(imar.gt.ismo) stop'surdimensionner ismo dans mva9 (grid_noro)'
426      if(jmar.gt.jsmo) stop'surdimensionner jsmo dans mva9 (grid_noro)'
427     
428      SUM=0.
429      DO IS=-1,1
430      DO JS=-1,1
431      WEIGHT(IS,JS)=1./FLOAT((1+IS**2)*(1+JS**2))
432      SUM=SUM+WEIGHT(IS,JS)
433      ENDDO
434      ENDDO
435
436      DO IS=-1,1
437      DO JS=-1,1
438      WEIGHT(IS,JS)=WEIGHT(IS,JS)/SUM
439      ENDDO
440      ENDDO
441
442      DO J=2,JMAR-1
443      DO I=2,IMAR-1
444      XF(I,J)=0.
445      DO IS=-1,1
446      DO JS=-1,1
447         XF(I,J)=XF(I,J)+X(I+IS,J+JS)*WEIGHT(IS,JS)
448      ENDDO
449      ENDDO
450      ENDDO
451      ENDDO
452
453      DO J=2,JMAR-1
454         XF(1,J)=0.
455      IS=IMAR-1
456      DO JS=-1,1
457         XF(1,J)=XF(1,J)+X(IS,J+JS)*WEIGHT(-1,JS)
458      ENDDO
459      DO IS=0,1
460      DO JS=-1,1
461         XF(1,J)=XF(1,J)+X(1+IS,J+JS)*WEIGHT(IS,JS)
462      ENDDO
463      ENDDO
464         XF(IMAR,J)=XF(1,J)
465      ENDDO
466
467      DO I=1,IMAR
468      XF(I,1)=XF(I,2)
469      XF(I,JMAR)=XF(I,JMAR-1)
470      ENDDO
471
472      DO I=1,IMAR
473      DO J=1,JMAR
474         X(I,J)=XF(I,J)
475      ENDDO
476      ENDDO
477
478      RETURN
479      END
480
481
482
Note: See TracBrowser for help on using the repository browser.