source: trunk/LMDZ.GENERIC/libf/dyn3d/grid_noro.F @ 1243

Last change on this file since 1243 was 135, checked in by aslmd, 14 years ago

CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET

File size: 13.9 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
88
89c
90c  quelques tests de dimensions:
91c   
92c
93      if(iim.ne.imar) STOP 'Problem dim. x'
94      if(jjm.ne.jmar-1) STOP 'Problem dim. y'
95      IF (imar.GT.2200 .OR. jmar.GT.1100) THEN
96         PRINT*, 'imar or jmar too big', imar, jmar
97         CALL ABORT
98      ENDIF
99
100      IF(imdep.ne.iusn.or.jmdep.ne.jusn)then
101         print *,' imdep or jmdep bad dimensions:',imdep,jmdep
102         call abort
103      ENDIF
104
105      IF(imar+1.ne.iim+1.or.jmar.ne.jjm+1)THEN
106        print *,' imar or jmar bad dimensions:',imar,jmar
107        call abort
108      ENDIF
109
110
111c      print *,'xdata:',xdata
112c      print *,'ydata:',ydata
113c      print *,'x:',x
114c      print *,'y:',y
115c
116C  EXTENSION OF THE USN DATABASE TO POCEED COMPUTATIONS AT
117C  BOUNDARIES:
118c
119      DO j=1,jusn
120        yusn(j+1)=ydata(j)
121      DO i=1,iusn
122        zusn(i+iext,j+1)=zdata(i,j)
123        xusn(i+iext)=xdata(i)
124      ENDDO
125      DO i=1,iext
126        zusn(i,j+1)=zdata(iusn-iext+i,j)
127        xusn(i)=xdata(iusn-iext+i)-2.*xpi
128        zusn(iusn+iext+i,j+1)=zdata(i,j)
129        xusn(iusn+iext+i)=xdata(i)+2.*xpi
130      ENDDO
131      ENDDO
132
133        yusn(1)=ydata(1)+(ydata(1)-ydata(2))
134        yusn(jusn+2)=ydata(jusn)+(ydata(jusn)-ydata(jusn-1))
135       DO i=1,iusn/2+iext
136        zusn(i,1)=zusn(i+iusn/2,2)
137        zusn(i+iusn/2+iext,1)=zusn(i,2)
138        zusn(i,jusn+2)=zusn(i+iusn/2,jusn+1)
139        zusn(i+iusn/2+iext,jusn+2)=zusn(i,jusn+1)
140       ENDDO
141
142c COMPUTE LIMITS OF MODEL GRIDPOINT AREA
143C     ( REGULAR GRID)
144c
145      a(1) = x(1) - (x(2)-x(1))/2.0
146      b(1) = (x(1)+x(2))/2.0
147      DO i = 2, imar
148         a(i) = b(i-1)
149         b(i) = (x(i)+x(i+1))/2.0
150      ENDDO
151      a(imar+1) = b(imar)
152      b(imar+1) = x(imar+1) + (x(imar+1)-x(imar))/2.0
153
154      c(1) = y(1) - (y(2)-y(1))/2.0
155      d(1) = (y(1)+y(2))/2.0
156      DO j = 2, jmar-1
157         c(j) = d(j-1)
158         d(j) = (y(j)+y(j+1))/2.0
159      ENDDO
160      c(jmar) = d(jmar-1)
161      d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0
162
163c
164c  initialisations:
165c
166
167c     write (*,*) "zpic: ",zpic
168
169      DO i = 1, imar+1
170         DO j = 1, jmar
171            weight(i,j) = 0.0
172            zxtzx(i,j)  = 0.0
173            zytzy(i,j)  = 0.0
174            zxtzy(i,j)  = 0.0
175            ztz(i,j)    = 0.0
176            zmea(i,j)   = 0.0
177            zpic(i,j)  = -1.E+10
178            zval(i,j)  = 1.E+10
179         ENDDO
180      ENDDO
181c
182c  COMPUTE SLOPES CORRELATIONS ON USN GRID
183c
184         DO j = 1,jusn+2
185         DO i = 1, iusn+2*iext
186            zytzyusn(i,j)=0.0
187            zxtzxusn(i,j)=0.0
188            zxtzyusn(i,j)=0.0
189         ENDDO
190         ENDDO
191
192
193         DO j = 2,jusn+1
194            zdeltax=zdeltay*cos(yusn(j))
195         DO i = 2, iusn+2*iext-1
196            zytzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))**2/zdeltay**2
197            zxtzxusn(i,j)=(zusn(i+1,j)-zusn(i-1,j))**2/zdeltax**2
198            zxtzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))/zdeltay
199     *                   *(zusn(i+1,j)-zusn(i-1,j))/zdeltax
200         ENDDO
201         ENDDO
202c
203c  SUMMATION OVER GRIDPOINT AREA
204c
205      zleny=xpi/float(jusn)*rad
206      xincr=xpi/2./float(jusn)
207       DO ii = 1, imar+1
208       DO jj = 1, jmar
209       num_tot(ii,jj)=0.
210       num_lan(ii,jj)=0.
211c        PRINT *,' iteration ii jj:',ii,jj
212         DO j = 2,jusn+1
213c         DO j = 3,jusn
214            zlenx=zleny*cos(yusn(j))
215            zdeltax=zdeltay*cos(yusn(j))
216            zbordnor=(c(jj)-yusn(j)+xincr)*rad
217            zbordsud=(yusn(j)-d(jj)+xincr)*rad
218            weighy=AMAX1(0.,
219     *             amin1(zbordnor,zbordsud,zleny))
220         IF(weighy.ne.0)THEN
221         DO i = 2, iusn+2*iext-1
222            zbordest=(xusn(i)-a(ii)+xincr)*rad*cos(yusn(j))
223            zbordoue=(b(ii)+xincr-xusn(i))*rad*cos(yusn(j))
224            weighx=AMAX1(0.,
225     *             amin1(zbordest,zbordoue,zlenx))
226            IF(weighx.ne.0)THEN
227            num_tot(ii,jj)=num_tot(ii,jj)+1.0
228            if(zusn(i,j).ge.1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
229            weight(ii,jj)=weight(ii,jj)+weighx*weighy
230            zxtzx(ii,jj)=zxtzx(ii,jj)+zxtzxusn(i,j)*weighx*weighy
231            zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)*weighx*weighy
232            zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)*weighx*weighy
233            ztz(ii,jj)  =ztz(ii,jj)  +zusn(i,j)*zusn(i,j)*weighx*weighy
234c mean
235            zmea(ii,jj) =zmea(ii,jj)+zusn(i,j)*weighx*weighy
236c peacks
237            zpic(ii,jj)=amax1(zpic(ii,jj),zusn(i,j))
238c valleys
239            zval(ii,jj)=amin1(zval(ii,jj),zusn(i,j))
240            ENDIF
241         ENDDO
242         ENDIF
243         ENDDO
244       ENDDO
245       ENDDO
246c
247c  COMPUTE PARAMETERS NEEDED BY THE LOTT & MILLER (1997) AND
248C  LOTT (1999) SSO SCHEME.
249c
250      zllmmea=0.
251      zllmstd=0.
252      zllmsig=0.
253      zllmgam=0.
254      zllmpic=0.
255      zllmval=0.
256      zllmthe=0.
257      zminthe=0.
258c     print 100,' '
259c100  format(1X,A1,'II JJ',4X,'H',8X,'SD',8X,'SI',3X,'GA',3X,'TH')
260       DO ii = 1, imar+1
261       DO jj = 1, jmar
262         IF (weight(ii,jj) .NE. 0.0) THEN
263c  Mask
264           if(num_lan(ii,jj)/num_tot(ii,jj).ge.0.5)then
265             mask(ii,jj)=1
266           else
267             mask(ii,jj)=0
268           endif
269c  Mean Orography:
270           zmea (ii,jj)=zmea (ii,jj)/weight(ii,jj)
271           zxtzx(ii,jj)=zxtzx(ii,jj)/weight(ii,jj)
272           zytzy(ii,jj)=zytzy(ii,jj)/weight(ii,jj)
273           zxtzy(ii,jj)=zxtzy(ii,jj)/weight(ii,jj)
274           ztz(ii,jj)  =ztz(ii,jj)/weight(ii,jj)
275c  Standard deviation:
276           zstd(ii,jj)=sqrt(AMAX1(0.,ztz(ii,jj)-zmea(ii,jj)**2))
277         ELSE
278            PRINT*, 'probleme,ii,jj=', ii,jj
279         ENDIF
280       ENDDO
281       ENDDO
282
283C CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:
284
285       DO ii = 1, imar+1
286         zxtzx(ii,1)=zxtzx(ii,2)
287         zxtzx(ii,jmar)=zxtzx(ii,jmar-1)
288         zxtzy(ii,1)=zxtzy(ii,2)
289         zxtzy(ii,jmar)=zxtzy(ii,jmar-1)
290         zytzy(ii,1)=zytzy(ii,2)
291         zytzy(ii,jmar)=zytzy(ii,jmar-1)
292       ENDDO
293
294C  FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.
295
296C  FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.
297
298       CALL MVA9(zmea,iim+1,jjm+1)
299       CALL MVA9(zstd,iim+1,jjm+1)
300       CALL MVA9(zpic,iim+1,jjm+1)
301       CALL MVA9(zval,iim+1,jjm+1)
302       CALL MVA9(zxtzx,iim+1,jjm+1)
303       CALL MVA9(zxtzy,iim+1,jjm+1)
304       CALL MVA9(zytzy,iim+1,jjm+1)
305
306       DO ii = 1, imar
307       DO jj = 1, jmar
308         IF (weight(ii,jj) .NE. 0.0) THEN
309c  Coefficients K, L et M:
310           xk=(zxtzx(ii,jj)+zytzy(ii,jj))/2.
311           xl=(zxtzx(ii,jj)-zytzy(ii,jj))/2.
312           xm=zxtzy(ii,jj)
313           xp=xk-sqrt(xl**2+xm**2)
314           xq=xk+sqrt(xl**2+xm**2)
315           xw=1.e-8
316           if(xp.le.xw) xp=0.
317           if(xq.le.xw) xq=xw
318           if(abs(xm).le.xw) xm=xw*sign(1.,xm)
319c slope:
320           zsig(ii,jj)=sqrt(xq)*mask(ii,jj)
321c isotropy:
322           zgam(ii,jj)=xp/xq*mask(ii,jj)
323c angle theta:
324           zthe(ii,jj)=57.29577951*atan2(xm,xl)/2.*mask(ii,jj)
325           zphi(ii,jj)=zmea(ii,jj)*mask(ii,jj)
326           zmea(ii,jj)=zmea(ii,jj)*mask(ii,jj)
327           zpic(ii,jj)=zpic(ii,jj)*mask(ii,jj)
328           zval(ii,jj)=zval(ii,jj)*mask(ii,jj)
329           zstd(ii,jj)=zstd(ii,jj)*mask(ii,jj)
330
331c          print 101,ii,jj,
332c    *           zmea(ii,jj),zstd(ii,jj),zsig(ii,jj),zgam(ii,jj),
333c    *           zthe(ii,jj)
334c101  format(1x,2(1x,i2),2(1x,f7.1),1x,f7.4,2x,f4.2,1x,f5.1)     
335         ELSE
336c           PRINT*, 'probleme,ii,jj=', ii,jj
337         ENDIF
338      zllmmea=AMAX1(zmea(ii,jj),zllmmea)
339      zllmstd=AMAX1(zstd(ii,jj),zllmstd)
340      zllmsig=AMAX1(zsig(ii,jj),zllmsig)
341      zllmgam=AMAX1(zgam(ii,jj),zllmgam)
342      zllmthe=AMAX1(zthe(ii,jj),zllmthe)
343      zminthe=amin1(zthe(ii,jj),zminthe)
344      zllmpic=AMAX1(zpic(ii,jj),zllmpic)
345      zllmval=AMAX1(zval(ii,jj),zllmval)
346       ENDDO
347       ENDDO
348
349      print *,'  MEAN ORO:',zllmmea
350      print *,'  ST. DEV.:',zllmstd
351      print *,'  PENTE:',zllmsig
352      print *,' ANISOTROP:',zllmgam
353      print *,'  ANGLE:',zminthe,zllmthe       
354      print *,'  pic:',zllmpic
355      print *,'  val:',zllmval
356     
357C
358c gamma and theta a 1. and 0. at poles
359c
360      DO jj=1,jmar
361      zmea(imar+1,jj)=zmea(1,jj)
362      zphi(imar+1,jj)=zphi(1,jj)
363      zpic(imar+1,jj)=zpic(1,jj)
364      zval(imar+1,jj)=zval(1,jj)
365      zstd(imar+1,jj)=zstd(1,jj)
366      zsig(imar+1,jj)=zsig(1,jj)
367      zgam(imar+1,jj)=zgam(1,jj)
368      zthe(imar+1,jj)=zthe(1,jj)
369      ENDDO
370
371
372      zmeanor=0.0
373      zmeasud=0.0
374      zstdnor=0.0
375      zstdsud=0.0
376      zsignor=0.0
377      zsigsud=0.0
378      zweinor=0.0
379      zweisud=0.0
380      zpicnor=0.0
381      zpicsud=0.0                                   
382      zvalnor=0.0
383      zvalsud=0.0
384
385      DO ii=1,imar
386      zweinor=zweinor+              weight(ii,   1)
387      zweisud=zweisud+              weight(ii,jmar)
388      zmeanor=zmeanor+zmea(ii,   1)*weight(ii,   1)
389      zmeasud=zmeasud+zmea(ii,jmar)*weight(ii,jmar)
390      zstdnor=zstdnor+zstd(ii,   1)*weight(ii,   1)
391      zstdsud=zstdsud+zstd(ii,jmar)*weight(ii,jmar)
392      zsignor=zsignor+zsig(ii,   1)*weight(ii,   1)
393      zsigsud=zsigsud+zsig(ii,jmar)*weight(ii,jmar)
394      zpicnor=zpicnor+zpic(ii,   1)*weight(ii,   1)
395      zpicsud=zpicsud+zpic(ii,jmar)*weight(ii,jmar)
396      zvalnor=zvalnor+zval(ii,   1)*weight(ii,   1)
397      zvalsud=zvalsud+zval(ii,jmar)*weight(ii,jmar)
398      ENDDO
399
400      DO ii=1,imar+1
401      zmea(ii,   1)=zmeanor/zweinor
402      zmea(ii,jmar)=zmeasud/zweisud
403      zphi(ii,   1)=zmeanor/zweinor
404      zphi(ii,jmar)=zmeasud/zweisud
405      zpic(ii,   1)=zpicnor/zweinor
406      zpic(ii,jmar)=zpicsud/zweisud
407      zval(ii,   1)=zvalnor/zweinor
408      zval(ii,jmar)=zvalsud/zweisud
409      zstd(ii,   1)=zstdnor/zweinor
410      zstd(ii,jmar)=zstdsud/zweisud
411      zsig(ii,   1)=zsignor/zweinor
412      zsig(ii,jmar)=zsigsud/zweisud
413      zgam(ii,   1)=1.
414      zgam(ii,jmar)=1.
415      zthe(ii,   1)=0.
416      zthe(ii,jmar)=0.
417      ENDDO
418
419      RETURN
420      END
421
422      SUBROUTINE MVA9(X,IMAR,JMAR)
423
424C MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
425
426      PARAMETER (ISMo=300,JSMo=200)
427      REAL X(IMAR,JMAR),XF(ISMo,JSMo)
428      real weight(-1:1,-1:1)
429
430      if(imar.gt.ismo) stop'surdimensionner ismo dans mva9 (grid_noro)'
431      if(jmar.gt.jsmo) stop'surdimensionner jsmo dans mva9 (grid_noro)'
432     
433      SUM=0.
434      DO IS=-1,1
435      DO JS=-1,1
436      WEIGHT(IS,JS)=1./FLOAT((1+IS**2)*(1+JS**2))
437      SUM=SUM+WEIGHT(IS,JS)
438      ENDDO
439      ENDDO
440
441      DO IS=-1,1
442      DO JS=-1,1
443      WEIGHT(IS,JS)=WEIGHT(IS,JS)/SUM
444      ENDDO
445      ENDDO
446
447      DO J=2,JMAR-1
448      DO I=2,IMAR-1
449      XF(I,J)=0.
450      DO IS=-1,1
451      DO JS=-1,1
452         XF(I,J)=XF(I,J)+X(I+IS,J+JS)*WEIGHT(IS,JS)
453      ENDDO
454      ENDDO
455      ENDDO
456      ENDDO
457
458      DO J=2,JMAR-1
459         XF(1,J)=0.
460      IS=IMAR-1
461      DO JS=-1,1
462         XF(1,J)=XF(1,J)+X(IS,J+JS)*WEIGHT(-1,JS)
463      ENDDO
464      DO IS=0,1
465      DO JS=-1,1
466         XF(1,J)=XF(1,J)+X(1+IS,J+JS)*WEIGHT(IS,JS)
467      ENDDO
468      ENDDO
469         XF(IMAR,J)=XF(1,J)
470      ENDDO
471
472      DO I=1,IMAR
473      XF(I,1)=XF(I,2)
474      XF(I,JMAR)=XF(I,JMAR-1)
475      ENDDO
476
477      DO I=1,IMAR
478      DO J=1,JMAR
479         X(I,J)=XF(I,J)
480      ENDDO
481      ENDDO
482
483      RETURN
484      END
485
486
487
Note: See TracBrowser for help on using the repository browser.