source: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/grid_noro.F @ 1174

Last change on this file since 1174 was 773, checked in by Laurent Fairhead, 17 years ago

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