source: LMDZ.3.3/branches/rel-LF/libf/dyn3d/grid_noro.F @ 163

Last change on this file since 163 was 99, checked in by lmdzadmin, 24 years ago

Version de travail de l'interface avec les surfaces. LF

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 14.6 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)
79c$$$ PB     integer mask(imar+1,jmar)
80      real mask(imar+1,jmar)
81      real num_tot(2200,1100),num_lan(2200,1100)
82c
83      REAL a(2200),b(2200),c(1100),d(1100)
84c
85      print *,' parametres de l orographie a l echelle sous maille'
86      xpi=acos(-1.)
87      rad    = 6 371 229.
88      zdeltay=2.*xpi/float(jusn)*rad
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
162c
163c  initialisations:
164c
165      DO i = 1, imar+1
166      DO j = 1, jmar
167         weight(i,j) = 0.0
168         zxtzx(i,j)  = 0.0
169         zytzy(i,j)  = 0.0
170         zxtzy(i,j)  = 0.0
171         ztz(i,j)    = 0.0
172         zmea(i,j)   = 0.0
173         zpic(i,j)  =-1.E+10
174         zval(i,j)  = 1.E+10
175      ENDDO
176      ENDDO
177c
178c  COMPUTE SLOPES CORRELATIONS ON USN GRID
179c
180         DO j = 1,jusn+2
181         DO i = 1, iusn+2*iext
182            zytzyusn(i,j)=0.0
183            zxtzxusn(i,j)=0.0
184            zxtzyusn(i,j)=0.0
185         ENDDO
186         ENDDO
187
188
189         DO j = 2,jusn+1
190            zdeltax=zdeltay*cos(yusn(j))
191         DO i = 2, iusn+2*iext-1
192            zytzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))**2/zdeltay**2
193            zxtzxusn(i,j)=(zusn(i+1,j)-zusn(i-1,j))**2/zdeltax**2
194            zxtzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))/zdeltay
195     *                   *(zusn(i+1,j)-zusn(i-1,j))/zdeltax
196         ENDDO
197         ENDDO
198c
199c  SUMMATION OVER GRIDPOINT AREA
200c
201      zleny=xpi/float(jusn)*rad
202      xincr=xpi/2./float(jusn)
203       DO ii = 1, imar+1
204       DO jj = 1, jmar
205       num_tot(ii,jj)=0.
206       num_lan(ii,jj)=0.
207c        PRINT *,' iteration ii jj:',ii,jj
208         DO j = 2,jusn+1
209c         DO j = 3,jusn
210            zlenx=zleny*cos(yusn(j))
211            zdeltax=zdeltay*cos(yusn(j))
212            zbordnor=(c(jj)-yusn(j)+xincr)*rad
213            zbordsud=(yusn(j)-d(jj)+xincr)*rad
214            weighy=AMAX1(0.,
215     *             amin1(zbordnor,zbordsud,zleny))
216         IF(weighy.ne.0)THEN
217         DO i = 2, iusn+2*iext-1
218            zbordest=(xusn(i)-a(ii)+xincr)*rad*cos(yusn(j))
219            zbordoue=(b(ii)+xincr-xusn(i))*rad*cos(yusn(j))
220            weighx=AMAX1(0.,
221     *             amin1(zbordest,zbordoue,zlenx))
222            IF(weighx.ne.0)THEN
223            num_tot(ii,jj)=num_tot(ii,jj)+1.0
224            if(zusn(i,j).ge.1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
225            weight(ii,jj)=weight(ii,jj)+weighx*weighy
226            zxtzx(ii,jj)=zxtzx(ii,jj)+zxtzxusn(i,j)*weighx*weighy
227            zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)*weighx*weighy
228            zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)*weighx*weighy
229            ztz(ii,jj)  =ztz(ii,jj)  +zusn(i,j)*zusn(i,j)*weighx*weighy
230c mean
231            zmea(ii,jj) =zmea(ii,jj)+zusn(i,j)*weighx*weighy
232c peacks
233            zpic(ii,jj)=amax1(zpic(ii,jj),zusn(i,j))
234c valleys
235            zval(ii,jj)=amin1(zval(ii,jj),zusn(i,j))
236            ENDIF
237         ENDDO
238         ENDIF
239         ENDDO
240       ENDDO
241       ENDDO
242c
243c  COMPUTE PARAMETERS NEEDED BY THE LOTT & MILLER (1997) AND
244C  LOTT (1999) SSO SCHEME.
245c
246      zllmmea=0.
247      zllmstd=0.
248      zllmsig=0.
249      zllmgam=0.
250      zllmpic=0.
251      zllmval=0.
252      zllmthe=0.
253      zminthe=0.
254c     print 100,' '
255c100  format(1X,A1,'II JJ',4X,'H',8X,'SD',8X,'SI',3X,'GA',3X,'TH')
256       DO ii = 1, imar+1
257       DO jj = 1, jmar
258         IF (weight(ii,jj) .NE. 0.0) THEN
259c  Mask
260c$$$           if(num_lan(ii,jj)/num_tot(ii,jj).ge.0.5)then
261c$$$             mask(ii,jj)=1
262c$$$           else
263c$$$             mask(ii,jj)=0
264c$$$           ENDIF
265             mask(ii,jj) = num_lan(ii,jj)/num_tot(ii,jj)
266c  Mean Orography:
267           zmea (ii,jj)=zmea (ii,jj)/weight(ii,jj)
268           zxtzx(ii,jj)=zxtzx(ii,jj)/weight(ii,jj)
269           zytzy(ii,jj)=zytzy(ii,jj)/weight(ii,jj)
270           zxtzy(ii,jj)=zxtzy(ii,jj)/weight(ii,jj)
271           ztz(ii,jj)  =ztz(ii,jj)/weight(ii,jj)
272c  Standard deviation:
273           zstd(ii,jj)=sqrt(AMAX1(0.,ztz(ii,jj)-zmea(ii,jj)**2))
274         ELSE
275            PRINT*, 'probleme,ii,jj=', ii,jj
276         ENDIF
277       ENDDO
278       ENDDO
279
280C CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:
281
282       DO ii = 1, imar+1
283         zxtzx(ii,1)=zxtzx(ii,2)
284         zxtzx(ii,jmar)=zxtzx(ii,jmar-1)
285         zxtzy(ii,1)=zxtzy(ii,2)
286         zxtzy(ii,jmar)=zxtzy(ii,jmar-1)
287         zytzy(ii,1)=zytzy(ii,2)
288         zytzy(ii,jmar)=zytzy(ii,jmar-1)
289       ENDDO
290
291C  FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.
292
293C  FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.
294
295       CALL MVA9(zmea,iim+1,jjm+1)
296       CALL MVA9(zstd,iim+1,jjm+1)
297       CALL MVA9(zpic,iim+1,jjm+1)
298       CALL MVA9(zval,iim+1,jjm+1)
299       CALL MVA9(zxtzx,iim+1,jjm+1)
300       CALL MVA9(zxtzy,iim+1,jjm+1)
301c$$$       CALL MVA9(zytzy,iim+1,jjm+1)
302
303       DO ii = 1, imar
304       DO jj = 1, jmar
305         IF (weight(ii,jj) .NE. 0.0) THEN
306c  Coefficients K, L et M:
307           xk=(zxtzx(ii,jj)+zytzy(ii,jj))/2.
308           xl=(zxtzx(ii,jj)-zytzy(ii,jj))/2.
309           xm=zxtzy(ii,jj)
310           xp=xk-sqrt(xl**2+xm**2)
311           xq=xk+sqrt(xl**2+xm**2)
312           xw=1.e-8
313           if(xp.le.xw) xp=0.
314           if(xq.le.xw) xq=xw
315           if(abs(xm).le.xw) xm=xw*sign(1.,xm)
316c$$$c slope:
317c$$$           zsig(ii,jj)=sqrt(xq)*mask(ii,jj)
318c$$$c isotropy:
319c$$$           zgam(ii,jj)=xp/xq*mask(ii,jj)
320c$$$c angle theta:
321c$$$           zthe(ii,jj)=57.29577951*atan2(xm,xl)/2.*mask(ii,jj)
322c$$$           zphi(ii,jj)=zmea(ii,jj)*mask(ii,jj)
323c$$$           zmea(ii,jj)=zmea(ii,jj)*mask(ii,jj)
324c$$$           zpic(ii,jj)=zpic(ii,jj)*mask(ii,jj)
325c$$$           zval(ii,jj)=zval(ii,jj)*mask(ii,jj)
326c$$$           zstd(ii,jj)=zstd(ii,jj)*mask(ii,jj)
327C$$* PB modif pour maque de terre fractionnaire
328c slope:
329           zsig(ii,jj)=sqrt(xq)
330c isotropy:
331           zgam(ii,jj)=xp/xq
332c angle theta:
333           zthe(ii,jj)=57.29577951*atan2(xm,xl)/2.
334c$$$           zphi(ii,jj)=zmea(ii,jj)
335c$$$           zmea(ii,jj)=zmea(ii,jj)
336c$$$           zpic(ii,jj)=zpic(ii,jj)
337c$$$           zval(ii,jj)=zval(ii,jj)
338c$$$           zstd(ii,jj)=zstd(ii,jj)
339c          print 101,ii,jj,
340c    *           zmea(ii,jj),zstd(ii,jj),zsig(ii,jj),zgam(ii,jj),
341c    *           zthe(ii,jj)
342c101  format(1x,2(1x,i2),2(1x,f7.1),1x,f7.4,2x,f4.2,1x,f5.1)     
343         ELSE
344c           PRINT*, 'probleme,ii,jj=', ii,jj
345         ENDIF
346      zllmmea=AMAX1(zmea(ii,jj),zllmmea)
347      zllmstd=AMAX1(zstd(ii,jj),zllmstd)
348      zllmsig=AMAX1(zsig(ii,jj),zllmsig)
349      zllmgam=AMAX1(zgam(ii,jj),zllmgam)
350      zllmthe=AMAX1(zthe(ii,jj),zllmthe)
351      zminthe=amin1(zthe(ii,jj),zminthe)
352      zllmpic=AMAX1(zpic(ii,jj),zllmpic)
353      zllmval=AMAX1(zval(ii,jj),zllmval)
354       ENDDO
355       ENDDO
356
357      print *,'  MEAN ORO:',zllmmea
358      print *,'  ST. DEV.:',zllmstd
359      print *,'  PENTE:',zllmsig
360      print *,' ANISOTROP:',zllmgam
361      print *,'  ANGLE:',zminthe,zllmthe       
362      print *,'  pic:',zllmpic
363      print *,'  val:',zllmval
364     
365C
366c gamma and theta a 1. and 0. at poles
367c
368      DO jj=1,jmar
369      zmea(imar+1,jj)=zmea(1,jj)
370      zphi(imar+1,jj)=zphi(1,jj)
371      zpic(imar+1,jj)=zpic(1,jj)
372      zval(imar+1,jj)=zval(1,jj)
373      zstd(imar+1,jj)=zstd(1,jj)
374      zsig(imar+1,jj)=zsig(1,jj)
375      zgam(imar+1,jj)=zgam(1,jj)
376      zthe(imar+1,jj)=zthe(1,jj)
377      ENDDO
378
379
380      zmeanor=0.0
381      zmeasud=0.0
382      zstdnor=0.0
383      zstdsud=0.0
384      zsignor=0.0
385      zsigsud=0.0
386      zweinor=0.0
387      zweisud=0.0
388      zpicnor=0.0
389      zpicsud=0.0                                   
390      zvalnor=0.0
391      zvalsud=0.0
392
393      DO ii=1,imar
394      zweinor=zweinor+              weight(ii,   1)
395      zweisud=zweisud+              weight(ii,jmar)
396      zmeanor=zmeanor+zmea(ii,   1)*weight(ii,   1)
397      zmeasud=zmeasud+zmea(ii,jmar)*weight(ii,jmar)
398      zstdnor=zstdnor+zstd(ii,   1)*weight(ii,   1)
399      zstdsud=zstdsud+zstd(ii,jmar)*weight(ii,jmar)
400      zsignor=zsignor+zsig(ii,   1)*weight(ii,   1)
401      zsigsud=zsigsud+zsig(ii,jmar)*weight(ii,jmar)
402      zpicnor=zpicnor+zpic(ii,   1)*weight(ii,   1)
403      zpicsud=zpicsud+zpic(ii,jmar)*weight(ii,jmar)
404      zvalnor=zvalnor+zval(ii,   1)*weight(ii,   1)
405      zvalsud=zvalsud+zval(ii,jmar)*weight(ii,jmar)
406      ENDDO
407
408      DO ii=1,imar+1
409      zmea(ii,   1)=zmeanor/zweinor
410      zmea(ii,jmar)=zmeasud/zweisud
411      zphi(ii,   1)=zmeanor/zweinor
412      zphi(ii,jmar)=zmeasud/zweisud
413      zpic(ii,   1)=zpicnor/zweinor
414      zpic(ii,jmar)=zpicsud/zweisud
415      zval(ii,   1)=zvalnor/zweinor
416      zval(ii,jmar)=zvalsud/zweisud
417      zstd(ii,   1)=zstdnor/zweinor
418      zstd(ii,jmar)=zstdsud/zweisud
419      zsig(ii,   1)=zsignor/zweinor
420      zsig(ii,jmar)=zsigsud/zweisud
421      zgam(ii,   1)=1.
422      zgam(ii,jmar)=1.
423      zthe(ii,   1)=0.
424      zthe(ii,jmar)=0.
425      ENDDO
426
427      RETURN
428      END
429
430      SUBROUTINE MVA9(X,IMAR,JMAR)
431
432C MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
433
434      PARAMETER (ISMo=300,JSMo=200)
435      REAL X(IMAR,JMAR),XF(ISMo,JSMo)
436      real WEIGHTpb(-1:1,-1:1)
437
438      if(imar.gt.ismo) stop'surdimensionner ismo dans mva9 (grid_noro)'
439      if(jmar.gt.jsmo) stop'surdimensionner jsmo dans mva9 (grid_noro)'
440     
441      SUM=0.
442      DO IS=-1,1
443        DO JS=-1,1
444          WEIGHTpb(IS,JS)=1./FLOAT((1+IS**2)*(1+JS**2))
445          SUM=SUM+WEIGHTpb(IS,JS)
446        ENDDO
447      ENDDO
448     
449      WRITE(*,*) 'MVA9 ', IMAR, JMAR
450      WRITE(*,*) 'MVA9 ', WEIGHTpb
451      WRITE(*,*) 'MVA9 SUM ', SUM
452      DO IS=-1,1
453        DO JS=-1,1
454          WEIGHTpb(IS,JS)=WEIGHTpb(IS,JS)/SUM
455        ENDDO
456      ENDDO
457
458      DO J=2,JMAR-1
459        DO I=2,IMAR-1
460          XF(I,J)=0.
461          DO IS=-1,1
462            DO JS=-1,1
463              XF(I,J)=XF(I,J)+X(I+IS,J+JS)*WEIGHTpb(IS,JS)
464            ENDDO
465          ENDDO
466        ENDDO
467      ENDDO
468
469      DO J=2,JMAR-1
470        XF(1,J)=0.
471        IS=IMAR-1
472        DO JS=-1,1
473          XF(1,J)=XF(1,J)+X(IS,J+JS)*WEIGHTpb(-1,JS)
474        ENDDO
475        DO IS=0,1
476          DO JS=-1,1
477            XF(1,J)=XF(1,J)+X(1+IS,J+JS)*WEIGHTpb(IS,JS)
478          ENDDO
479        ENDDO
480        XF(IMAR,J)=XF(1,J)
481      ENDDO
482
483      DO I=1,IMAR
484        XF(I,1)=XF(I,2)
485        XF(I,JMAR)=XF(I,JMAR-1)
486      ENDDO
487
488      DO I=1,IMAR
489        DO J=1,JMAR
490          X(I,J)=XF(I,J)
491        ENDDO
492      ENDDO
493
494      RETURN
495      END
496
497
498
Note: See TracBrowser for help on using the repository browser.