source: LMDZ4/trunk/libf/dyn3d/grid_noro.F @ 633

Last change on this file since 633 was 524, checked in by lmdzadmin, 21 years ago

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