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

Last change on this file since 406 was 404, checked in by lmdzadmin, 22 years ago

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