source: trunk/LMDZ.TITAN/libf/phytitan/grid_noro.F @ 1243

Last change on this file since 1243 was 1056, checked in by slebonnois, 11 years ago

SL: Titan runs ! see DOC/chantiers/commit_importants.log
File size: 13.8 KB
RevLine 
[1056]1!
2! $Id: grid_noro.F 1442 2010-10-18 08:31:31Z jghattas $
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)
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 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
54      IMPLICIT none
55     
56#include "dimensions.h"
57#include "YOMCST.h"
58
59      INTEGER imdep, jmdep
60      REAL xdata(imdep),ydata(jmdep)
61      REAL zdata(imdep,jmdep)
62c
63      INTEGER imar, jmar
64c parametres lies au fichier d entree... A documenter...
65      integer iext
66      parameter(iext=216)
67      REAL xusn(imdep+2*iext),yusn(jmdep+2)
68      REAL zusn(imdep+2*iext,jmdep+2)
69 
70c local var
71      real zdeltax,zdeltay,zlenx,zleny,xincr
72      real zbordnor,zbordsud,zbordest,zbordoue,weighx,weighy
73      real zllmmea,zllmstd,zllmsig,zllmgam,zllmpic,zllmval,zllmthe
74      real zminthe,xk,xl,xm,xp,xq,xw
75      real zmeanor,zmeasud,zstdnor,zstdsud,zsignor,zsigsud
76      real zweinor,zweisud,zpicnor,zpicsud,zvalnor,zvalsud
77      integer i,j,ii,jj
78
79C INTERMEDIATE FIELDS  (CORRELATIONS OF OROGRAPHY GRADIENT)
80
81      REAL ztz(iim+1,jjm+1),zxtzx(iim+1,jjm+1)
82      REAL zytzy(iim+1,jjm+1),zxtzy(iim+1,jjm+1)
83      REAL weight(iim+1,jjm+1)
84
85C CORRELATIONS OF USN OROGRAPHY GRADIENTS
86
87      REAL zxtzxusn(imdep+2*iext,jmdep+2)
88      REAL zytzyusn(imdep+2*iext,jmdep+2)
89      REAL zxtzyusn(imdep+2*iext,jmdep+2)
90      REAL x(imar+1),y(jmar),zphi(imar+1,jmar)
91      REAL zmea(imar+1,jmar),zstd(imar+1,jmar)
92      REAL zsig(imar+1,jmar),zgam(imar+1,jmar),zthe(imar+1,jmar)
93      REAL zpic(imar+1,jmar),zval(imar+1,jmar)
94      real num_tot(2200,1100),num_lan(2200,1100)
95c
96      REAL a(2200),b(2200),c(1100),d(1100)
97c
98      print *,' parametres de l orographie a l echelle sous maille'
99
100      zdeltay=2.*RPI/REAL(jmdep)*RA
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(imar+1.ne.iim+1.or.jmar.ne.jjm+1)THEN
113        print *,' imar or jmar bad dimensions:',imar,jmar
114        call abort
115      ENDIF
116
117
118c      print *,'xdata:',xdata
119c      print *,'ydata:',ydata
120c      print *,'x:',x
121c      print *,'y:',y
122c
123C  EXTENSION OF THE USN DATABASE TO POCEED COMPUTATIONS AT
124C  BOUNDARIES:
125c
126      DO j=1,jmdep
127        yusn(j+1)=ydata(j)
128      DO i=1,imdep
129        zusn(i+iext,j+1)=zdata(i,j)
130        xusn(i+iext)=xdata(i)
131      ENDDO
132      DO i=1,iext
133        zusn(i,j+1)=zdata(imdep-iext+i,j)
134        xusn(i)=xdata(imdep-iext+i)-2.*RPI
135        zusn(imdep+iext+i,j+1)=zdata(i,j)
136        xusn(imdep+iext+i)=xdata(i)+2.*RPI
137      ENDDO
138      ENDDO
139
140        yusn(1)=ydata(1)+(ydata(1)-ydata(2))
141        yusn(jmdep+2)=ydata(jmdep)+(ydata(jmdep)-ydata(jmdep-1))
142       DO i=1,imdep/2+iext
143        zusn(i,1)=zusn(i+imdep/2,2)
144        zusn(i+imdep/2+iext,1)=zusn(i,2)
145        zusn(i,jmdep+2)=zusn(i+imdep/2,jmdep+1)
146        zusn(i+imdep/2+iext,jmdep+2)=zusn(i,jmdep+1)
147       ENDDO
148c 
149c COMPUTE LIMITS OF MODEL GRIDPOINT AREA
150C     ( REGULAR GRID)
151c
152      a(1) = x(1) - (x(2)-x(1))/2.0
153      b(1) = (x(1)+x(2))/2.0
154      DO i = 2, imar
155         a(i) = b(i-1)
156         b(i) = (x(i)+x(i+1))/2.0
157      ENDDO
158      a(imar+1) = b(imar)
159      b(imar+1) = x(imar+1) + (x(imar+1)-x(imar))/2.0
160
161      c(1) = y(1) - (y(2)-y(1))/2.0
162      d(1) = (y(1)+y(2))/2.0
163      DO j = 2, jmar-1
164         c(j) = d(j-1)
165         d(j) = (y(j)+y(j+1))/2.0
166      ENDDO
167      c(jmar) = d(jmar-1)
168      d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0
169c
170c  initialisations:
171c
172      DO i = 1, imar+1
173      DO j = 1, jmar
174         weight(i,j) = 0.0
175         zxtzx(i,j)  = 0.0
176         zytzy(i,j)  = 0.0
177         zxtzy(i,j)  = 0.0
178         ztz(i,j)    = 0.0
179         zmea(i,j)   = 0.0
180         zpic(i,j)  =-1.E+10
181         zval(i,j)  = 1.E+10
182      ENDDO
183      ENDDO
184c
185c  COMPUTE SLOPES CORRELATIONS ON USN GRID
186c
187         DO j = 1,jmdep+2
188         DO i = 1, imdep+2*iext
189            zytzyusn(i,j)=0.0
190            zxtzxusn(i,j)=0.0
191            zxtzyusn(i,j)=0.0
192         ENDDO
193         ENDDO
194
195
196         DO j = 2,jmdep+1
197            zdeltax=zdeltay*cos(yusn(j))
198         DO i = 2, imdep+2*iext-1
199            zytzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))**2/zdeltay**2
200            zxtzxusn(i,j)=(zusn(i+1,j)-zusn(i-1,j))**2/zdeltax**2
201            zxtzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))/zdeltay
202     *                   *(zusn(i+1,j)-zusn(i-1,j))/zdeltax
203         ENDDO
204         ENDDO
205c
206c  SUMMATION OVER GRIDPOINT AREA
207c
208      zleny=RPI/REAL(jmdep)*RA
209      xincr=RPI/2./REAL(jmdep)
210       DO ii = 1, imar+1
211       DO jj = 1, jmar
212       num_tot(ii,jj)=0.
213       num_lan(ii,jj)=0.
214c        PRINT *,' iteration ii jj:',ii,jj
215         DO j = 2,jmdep+1
216c         DO j = 3,jmdep
217            zlenx=zleny*cos(yusn(j))
218            zdeltax=zdeltay*cos(yusn(j))
219            zbordnor=(c(jj)-yusn(j)+xincr)*RA
220            zbordsud=(yusn(j)-d(jj)+xincr)*RA
221            weighy=AMAX1(0.,
222     *             amin1(zbordnor,zbordsud,zleny))
223         IF(weighy.ne.0)THEN
224         DO i = 2, imdep+2*iext-1
225            zbordest=(xusn(i)-a(ii)+xincr)*RA*cos(yusn(j))
226            zbordoue=(b(ii)+xincr-xusn(i))*RA*cos(yusn(j))
227            weighx=AMAX1(0.,
228     *             amin1(zbordest,zbordoue,zlenx))
229            IF(weighx.ne.0)THEN
230            num_tot(ii,jj)=num_tot(ii,jj)+1.0
231            if(zusn(i,j).ge.1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
232            weight(ii,jj)=weight(ii,jj)+weighx*weighy
233            zxtzx(ii,jj)=zxtzx(ii,jj)+zxtzxusn(i,j)*weighx*weighy
234            zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)*weighx*weighy
235            zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)*weighx*weighy
236            ztz(ii,jj)  =ztz(ii,jj)  +zusn(i,j)*zusn(i,j)*weighx*weighy
237c mean
238            zmea(ii,jj) =zmea(ii,jj)+zusn(i,j)*weighx*weighy
239c peacks
240            zpic(ii,jj)=amax1(zpic(ii,jj),zusn(i,j))
241c valleys
242            zval(ii,jj)=amin1(zval(ii,jj),zusn(i,j))
243            ENDIF
244         ENDDO
245         ENDIF
246         ENDDO
247       ENDDO
248       ENDDO
249c
250c  COMPUTE PARAMETERS NEEDED BY THE LOTT & MILLER (1997) AND
251C  LOTT (1999) SSO SCHEME.
252c
253      zllmmea=0.
254      zllmstd=0.
255      zllmsig=0.
256      zllmgam=0.
257      zllmpic=0.
258      zllmval=0.
259      zllmthe=0.
260      zminthe=0.
261c     print 100,' '
262c100  format(1X,A1,'II JJ',4X,'H',8X,'SD',8X,'SI',3X,'GA',3X,'TH')
263       DO ii = 1, imar+1
264       DO jj = 1, jmar
265         IF (weight(ii,jj) .NE. 0.0) THEN
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)
301       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 slope:
317           zsig(ii,jj)=sqrt(xq)
318c isotropy:
319           zgam(ii,jj)=xp/xq
320c angle theta:
321           zthe(ii,jj)=57.29577951*atan2(xm,xl)/2.
322           zphi(ii,jj)=zmea(ii,jj)
323           !
324           zmea(ii,jj)=zmea(ii,jj)
325           zpic(ii,jj)=zpic(ii,jj)
326           zval(ii,jj)=zval(ii,jj)
327           zstd(ii,jj)=zstd(ii,jj)
328c          print 101,ii,jj,
329c    *           zmea(ii,jj),zstd(ii,jj),zsig(ii,jj),zgam(ii,jj),
330c    *           zthe(ii,jj)
331c101  format(1x,2(1x,i2),2(1x,f7.1),1x,f7.4,2x,f4.2,1x,f5.1)     
332         ELSE
333c           PRINT*, 'probleme,ii,jj=', ii,jj
334         ENDIF
335      zllmmea=AMAX1(zmea(ii,jj),zllmmea)
336      zllmstd=AMAX1(zstd(ii,jj),zllmstd)
337      zllmsig=AMAX1(zsig(ii,jj),zllmsig)
338      zllmgam=AMAX1(zgam(ii,jj),zllmgam)
339      zllmthe=AMAX1(zthe(ii,jj),zllmthe)
340      zminthe=amin1(zthe(ii,jj),zminthe)
341      zllmpic=AMAX1(zpic(ii,jj),zllmpic)
342      zllmval=AMAX1(zval(ii,jj),zllmval)
343       ENDDO
344       ENDDO
345      print *,'  MEAN ORO:',zllmmea
346      print *,'  ST. DEV.:',zllmstd
347      print *,'  PENTE:',zllmsig
348      print *,' ANISOTROP:',zllmgam
349      print *,'  ANGLE:',zminthe,zllmthe
350      print *,'  pic:',zllmpic
351      print *,'  val:',zllmval
352     
353C
354c gamma and theta a 1. and 0. at poles
355c
356      DO jj=1,jmar
357      zmea(imar+1,jj)=zmea(1,jj)
358      zphi(imar+1,jj)=zphi(1,jj)
359      zpic(imar+1,jj)=zpic(1,jj)
360      zval(imar+1,jj)=zval(1,jj)
361      zstd(imar+1,jj)=zstd(1,jj)
362      zsig(imar+1,jj)=zsig(1,jj)
363      zgam(imar+1,jj)=zgam(1,jj)
364      zthe(imar+1,jj)=zthe(1,jj)
365      ENDDO
366
367
368      zmeanor=0.0
369      zmeasud=0.0
370      zstdnor=0.0
371      zstdsud=0.0
372      zsignor=0.0
373      zsigsud=0.0
374      zweinor=0.0
375      zweisud=0.0
376      zpicnor=0.0
377      zpicsud=0.0                                   
378      zvalnor=0.0
379      zvalsud=0.0
380
381      DO ii=1,imar
382      zweinor=zweinor+              weight(ii,   1)
383      zweisud=zweisud+              weight(ii,jmar)
384      zmeanor=zmeanor+zmea(ii,   1)*weight(ii,   1)
385      zmeasud=zmeasud+zmea(ii,jmar)*weight(ii,jmar)
386      zstdnor=zstdnor+zstd(ii,   1)*weight(ii,   1)
387      zstdsud=zstdsud+zstd(ii,jmar)*weight(ii,jmar)
388      zsignor=zsignor+zsig(ii,   1)*weight(ii,   1)
389      zsigsud=zsigsud+zsig(ii,jmar)*weight(ii,jmar)
390      zpicnor=zpicnor+zpic(ii,   1)*weight(ii,   1)
391      zpicsud=zpicsud+zpic(ii,jmar)*weight(ii,jmar)
392      zvalnor=zvalnor+zval(ii,   1)*weight(ii,   1)
393      zvalsud=zvalsud+zval(ii,jmar)*weight(ii,jmar)
394      ENDDO
395
396      DO ii=1,imar+1
397      zmea(ii,   1)=zmeanor/zweinor
398      zmea(ii,jmar)=zmeasud/zweisud
399      zphi(ii,   1)=zmeanor/zweinor
400      zphi(ii,jmar)=zmeasud/zweisud
401      zpic(ii,   1)=zpicnor/zweinor
402      zpic(ii,jmar)=zpicsud/zweisud
403      zval(ii,   1)=zvalnor/zweinor
404      zval(ii,jmar)=zvalsud/zweisud
405      zstd(ii,   1)=zstdnor/zweinor
406      zstd(ii,jmar)=zstdsud/zweisud
407      zsig(ii,   1)=zsignor/zweinor
408      zsig(ii,jmar)=zsigsud/zweisud
409      zgam(ii,   1)=1.
410      zgam(ii,jmar)=1.
411      zthe(ii,   1)=0.
412      zthe(ii,jmar)=0.
413      ENDDO
414
415      RETURN
416      END
417
418      SUBROUTINE MVA9(X,IMAR,JMAR)
419
420C MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
421
422      REAL X(IMAR,JMAR),XF(IMAR,JMAR)
423      real WEIGHTpb(-1:1,-1:1)
424
425
426      SUM=0.
427      DO IS=-1,1
428        DO JS=-1,1
429          WEIGHTpb(IS,JS)=1./REAL((1+IS**2)*(1+JS**2))
430          SUM=SUM+WEIGHTpb(IS,JS)
431        ENDDO
432      ENDDO
433     
434c     WRITE(*,*) 'MVA9 ', IMAR, JMAR
435c     WRITE(*,*) 'MVA9 ', WEIGHTpb
436c     WRITE(*,*) 'MVA9 SUM ', SUM
437      DO IS=-1,1
438        DO JS=-1,1
439          WEIGHTpb(IS,JS)=WEIGHTpb(IS,JS)/SUM
440        ENDDO
441      ENDDO
442
443      DO J=2,JMAR-1
444        DO I=2,IMAR-1
445          XF(I,J)=0.
446          DO IS=-1,1
447            DO JS=-1,1
448              XF(I,J)=XF(I,J)+X(I+IS,J+JS)*WEIGHTpb(IS,JS)
449            ENDDO
450          ENDDO
451        ENDDO
452      ENDDO
453
454      DO J=2,JMAR-1
455        XF(1,J)=0.
456        IS=IMAR-1
457        DO JS=-1,1
458          XF(1,J)=XF(1,J)+X(IS,J+JS)*WEIGHTpb(-1,JS)
459        ENDDO
460        DO IS=0,1
461          DO JS=-1,1
462            XF(1,J)=XF(1,J)+X(1+IS,J+JS)*WEIGHTpb(IS,JS)
463          ENDDO
464        ENDDO
465        XF(IMAR,J)=XF(1,J)
466      ENDDO
467
468      DO I=1,IMAR
469        XF(I,1)=XF(I,2)
470        XF(I,JMAR)=XF(I,JMAR-1)
471      ENDDO
472
473      DO I=1,IMAR
474        DO J=1,JMAR
475          X(I,J)=XF(I,J)
476        ENDDO
477      ENDDO
478
479      RETURN
480      END
481
482
483
Note: See TracBrowser for help on using the repository browser.