source: trunk/LMDZ.MARS/libf/dyn3d/grid_noro1.F @ 937

Last change on this file since 937 was 38, checked in by emillour, 14 years ago

Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs.
EM

File size: 13.3 KB
Line 
1      SUBROUTINE grid_noro1(imdep, jmdep, xdata, ydata, entree,
2     .                 imar, jmar, x, y, zmea,zstd,zsig,zgam,zthe)
3c=======================================================================
4c (F. Lott) (voir aussi z.x. Li, A. Harzallah et L. Fairhead)
5c
6c      Calcul des parametres de l'orographie sous-maille necessaires
7c      au nouveau shema de representation des montagnes meso-echelles
8c      dans le modele.  Les points sont mis sur une grille rectangulaire
9c      pseudo-physique.  Typiquement, il y a iim+1 latitudes incluant
10c      le pole nord et le pole sud.  Il y a jjm+1 longitudes, y compris
11c      aux poles.  Aux poles les champs peuvent ont une valeurs repetee
12c      jjm+1 fois.....  La valeur du champs en jjm+1 (jmar) est celle
13c      en j=1. 
14c      Les parametres a,b,c,d representent les limites de la region
15c      de point de grille correspondant a un point decrit precedemment.
16c      Les moyennes sur ces regions des valeurs calculees a partir de
17c      l'USN, sont ponderees par un poids, fonction de la surface
18c      occuppe par ces donnees a l'interieure de la grille du modele.
19c      Dans la plupart des cas ce poid est le rapport entre la surface
20c      de la region de point de grille USN et la surface de la region
21c      de point de grille du modele.
22c       
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 et Y pour depart
36c        xdata, ydata: coordonnees X et Y pour depart
37c        entree: champ d'entree a transformer
38c        dans ce programme, on assume que les donnees sont les altitudes
39c        de l'USNavy: imdep=iusn=2160, jmdep=jusn=1080.
40c OUTPUT:
41c        imar, jmar: dimensions X et Y d'arrivee
42c        x, y: coordonnees X et Y d'arrivee
43c        les champs de sorties sont sur une grille physique:
44c             zmea:  orographie moyenne
45c             zstd:  deviation standard de l'orographie sous-maille
46c             zsig:  pente de l'orographie sous-maille
47c             zgam:  anisotropy de l'orographie sous maille
48c             zthe:  orientation de l'axe oriente dans la direction
49c                    de plus grande pente de l'orographie sous maille
50C=======================================================================
51c     IMPLICIT INTEGER (I,J)
52c     IMPLICIT REAL(X,Z)
53       implicit none
54       integer iusn,jusn,iext
55       parameter(iusn=360,jusn=180,iext=40)
56c!-*-      include 'param1'
57c!-*-      include 'comcstfi.h'
58#include "dimensions.h"
59#include "comconst.h"
60c!-*-
61c!-*-      parameter(iim=cols,jjm=rows)
62      REAL xusn(iusn+2*iext),yusn(jusn+2)       
63      REAL zusn(iusn+2*iext,jusn+2),zusnfi(iusn+2*iext,jusn+2)
64
65c   modif declarations pour implicit none
66      real zmeanor,zmeasud,zstdnor,zstdsud,zsignor
67      real zsigsud,zweinor,zweisud
68      real xk,xl,xm,xw,xp,xq
69      real zmaxmea,zmaxstd,zmaxsig,zmaxgam,zmaxthe,zminthe
70      real zbordnor,zbordsud,zbordest,zbordoue,xpi
71      real zdeltax,zdeltay,zlenx,zleny,weighx,weighy,xincr
72      integer i,j,ii,jj,ideltax,ihalph
73
74      INTEGER imdep, jmdep
75      REAL xdata(imdep),ydata(jmdep)
76      REAL entree(imdep,jmdep)
77c
78      INTEGER imar, jmar
79 
80      REAL ztz(iim+1,jjm+1),zxtzx(iim+1,jjm+1)
81      REAL zytzy(iim+1,jjm+1),zxtzy(iim+1,jjm+1)
82      REAL zxtzxusn(iusn+2*iext,jusn+2),zytzyusn(iusn+2*iext,jusn+2)
83      REAL zxtzyusn(iusn+2*iext,jusn+2)
84      REAL weight(iim+1,jjm+1)
85      REAL x(imar+1),y(jmar)
86      REAL zmea(imar+1,jmar),zstd(imar+1,jmar)
87      REAL zsig(imar+1,jmar),zgam(imar+1,jmar),zthe(imar+1,jmar)
88c
89      REAL a(2200),b(2200),c(1100),d(1100)
90c
91c  quelques constantes:
92c
93      print *,' parametres de l orographie a l echelle sous maille'
94      print*,'rad =',rad
95      print*,'Long et lat entree'
96      print*,(x(i),i=1,imar+1)
97      print*,(y(j),j=1,jmar)
98       print*,'Long et lat donnees'
99      print*,(xdata(i),i=1,imdep)
100      print*,(ydata(j),j=1,jmdep)
101
102      xpi=acos(-1.)
103      zdeltay=2.*xpi/real(jusn)*rad
104c
105c  quelques tests de dimensions:
106c   
107      IF (imar.GT.2200 .OR. jmar.GT.1100) THEN
108         PRINT*, 'imar ou jmar trop grand', imar, jmar
109         CALL ABORT
110      ENDIF
111
112      IF(imdep.ne.iusn.or.jmdep.ne.jusn)then
113         print *,' imdep ou jmdep mal dimensionnes:',imdep,jmdep
114         call abort
115      ENDIF
116
117      IF(imar+1.gt.iim+1.or.jmar.gt.jjm+1)THEN
118        print *,' imar ou jmar mal dimensionnes:',imar,jmar
119        call abort
120      ENDIF
121c
122C  Extension de la base de donnee de l'USN pour faciliter
123C  les calculs ulterieurs:
124c
125      DO j=1,jusn
126        yusn(j+1)=ydata(j)
127      DO i=1,iusn
128        zusn(i+iext,j+1)=entree(i,j)
129        xusn(i+iext)=xdata(i)
130      ENDDO
131      DO i=1,iext
132        zusn(i,j+1)=entree(iusn-iext+i,j)
133        xusn(i)=xdata(iusn-iext+i)-2.*xpi
134        zusn(iusn+iext+i,j+1)=entree(i,j)
135        xusn(iusn+iext+i)=xdata(i)+2.*xpi
136      ENDDO
137      ENDDO
138
139        yusn(1)=ydata(1)+(ydata(1)-ydata(2))
140        yusn(jusn+2)=ydata(jusn)+(ydata(jusn)-ydata(jusn-1))
141       DO i=1,iusn/2+iext
142        zusn(i,1)=zusn(i+iusn/2,2)
143        zusn(i+iusn/2+iext,1)=zusn(i,2)
144        zusn(i,jusn+2)=zusn(i+iusn/2,jusn+1)
145        zusn(i+iusn/2+iext,jusn+2)=zusn(i,jusn+1)
146       ENDDO
147c
148c  Calcul d'une orographie filtree aux hautes latitudes
149c  pour permettre des calculs plus isotropiques sur la pente
150c  des montagnes
151c
152       DO i=1,IUSN+2*iext
153       DO J=1,JUSN+2
154          zusnfi(i,j)=0.0
155       ENDDO
156       ENDDO
157
158      DO j=1,jusn
159            ideltax=1./cos(yusn(j+1))
160            ideltax=min(iusn/2-1,ideltax)
161            IF(MOD(IDELTAX,2).EQ.0)THEN
162              IDELTAX=IDELTAX+1
163            ENDIF
164            IHALPH=(IDELTAX-1)/2
165c           print *,' ideltax=',ideltax
166         IF(ideltax.eq.1)THEN
167            DO i=1,iusn
168               zusnfi(i+iext,j+1)=entree(i,j)
169            ENDDO   
170         ELSE
171            DO i=1,ihalph
172               DO ii=1,i+ihalph
173               zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)+entree(ii,j)
174               ENDDO
175               DO ii=ihalph-i,0,-1
176               zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)+entree(iusn-ii,j)
177               ENDDO 
178               zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)/real(ideltax)
179            ENDDO   
180            DO i=iusn-ihalph+1,iusn
181               DO ii = i-ihalph,iusn
182               zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)+entree(ii,j)
183               ENDDO
184               DO ii = 1,ihalph+i-iusn
185               zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)+entree(ii,j)
186               ENDDO
187               zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)/real(ideltax)
188            ENDDO
189            DO i=ihalph+1,iusn-ihalph
190               DO ii=-ihalph,ihalph
191               zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)+entree(i+ii,j)
192               ENDDO
193               zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)/real(ideltax)
194            ENDDO
195         ENDIF
196            DO i=1,iext
197               zusnfi(i,j+1)=zusnfi(iusn-iext+i,j+1)
198               zusnfi(i+iusn+iext,j+1)=zusnfi(i,j+1)
199            ENDDO
200      ENDDO
201
202c Calculer les limites des zones des nouveaux points
203c
204      a(1) = x(1) - (x(2)-x(1))/2.0
205      b(1) = (x(1)+x(2))/2.0
206      DO i = 2, imar-1
207         a(i) = b(i-1)
208         b(i) = (x(i)+x(i+1))/2.0
209      ENDDO
210      a(imar) = b(imar-1)
211      b(imar) = x(imar) + (x(imar)-x(imar-1))/2.0
212
213      c(1) = y(1) - (y(2)-y(1))/2.0
214      d(1) = (y(1)+y(2))/2.0
215      DO j = 2, jmar-1
216         c(j) = d(j-1)
217         d(j) = (y(j)+y(j+1))/2.0
218      ENDDO
219      c(jmar) = d(jmar-1)
220      d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0
221c
222c      quelques initialisations:
223      print*,'OKM1'
224c
225      DO i = 1, imar
226      DO j = 1, jmar
227         weight(i,j) = 0.0
228         zxtzx(i,j) = 0.0
229         zytzy(i,j) = 0.0
230         zxtzy(i,j) = 0.0
231         ztz(i,j) = 0.0
232         zmea(i,j) = 0.0
233         zstd(i,j)=0.0
234      ENDDO
235      ENDDO
236c
237c  calculs des correlations de pentes sur la grille de l'USN.
238c
239         DO j = 2,jusn+1
240         DO i = 1, iusn+2*iext
241            zytzyusn(i,j)=0.0
242            zxtzxusn(i,j)=0.0
243            zxtzyusn(i,j)=0.0
244         ENDDO
245         ENDDO
246
247
248         DO j = 2,jusn+1
249            zdeltax=zdeltay*cos(yusn(j))
250         DO i = 2, iusn+2*iext-1
251            zytzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))**2/zdeltay**2
252            zxtzxusn(i,j)=(zusnfi(i+1,j)-zusnfi(i-1,j))**2/zdeltax**2
253            zxtzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))/zdeltay
254     *                   *(zusnfi(i+1,j)-zusnfi(i-1,j))/zdeltax
255         ENDDO
256
257         ENDDO
258
259 
260
261      print*,'OK0'
262c
263c  sommations des differentes quantites definies precedemment
264c  sur une grille du modele.
265c
266      zleny=xpi/real(jusn)*rad
267      xincr=xpi/2./real(jusn)
268       DO ii = 1, imar
269       DO jj = 1, jmar
270c        PRINT *,' iteration ii jj:',ii,jj
271         DO j = 2,jusn+1
272c         DO j = 3,jusn
273            zlenx=zleny*cos(yusn(j))
274            zdeltax=zdeltay*cos(yusn(j))
275            zbordnor=(c(jj)-yusn(j)+xincr)*rad
276            zbordsud=(yusn(j)-d(jj)+xincr)*rad
277            weighy=amax1(0.,
278     *             amin1(zbordnor,zbordsud,zleny))
279         IF(weighy.ne.0)THEN
280         DO i = 2, iusn+2*iext-1
281            zbordest=(xusn(i)-a(ii)+xincr)*rad*cos(yusn(j))
282            zbordoue=(b(ii)+xincr-xusn(i))*rad*cos(yusn(j))
283            weighx=amax1(0.,
284     *             amin1(zbordest,zbordoue,zlenx))
285            IF(weighx.ne.0)THEN
286            weight(ii,jj)=weight(ii,jj)+weighx*weighy
287            zxtzx(ii,jj)=zxtzx(ii,jj)+zxtzxusn(i,j)*weighx*weighy
288            zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)*weighx*weighy
289            zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)*weighx*weighy
290            ztz(ii,jj)  =ztz(ii,jj)  +zusn(i,j)*zusn(i,j)*weighx*weighy
291            zmea(ii,jj) =zmea(ii,jj)+zusn(i,j)*weighx*weighy
292            ENDIF
293         ENDDO
294         ENDIF
295         ENDDO
296       ENDDO
297       ENDDO
298c
299c  calculs des differents parametres necessaires au programme
300c  de parametrisation de l'orographie a l'echelle moyenne:
301c
302      zmaxmea=0.
303      zmaxstd=0.
304      zmaxsig=0.
305      zmaxgam=0.
306      zmaxthe=0.
307      zminthe=0.
308c     print 100,' '
309c100  format(1X,A1,'II JJ',4X,'H',8X,'SD',8X,'SI',3X,'GA',3X,'TH')
310       print*,'OK1'
311       DO ii = 1, imar
312       DO jj = 1, jmar
313c       print*,'ok0'
314         IF (weight(ii,jj) .NE. 0.0) THEN
315c  Orography moyenne:
316c         print*,'ok1'
317           zmea (ii,jj)=zmea (ii,jj)/weight(ii,jj)
318           zxtzx(ii,jj)=zxtzx(ii,jj)/weight(ii,jj)
319           zytzy(ii,jj)=zytzy(ii,jj)/weight(ii,jj)
320           zxtzy(ii,jj)=zxtzy(ii,jj)/weight(ii,jj)
321           ztz(ii,jj)  =ztz(ii,jj)/weight(ii,jj)
322c         print*,'ok2'
323c  Deviation standard:
324           zstd(ii,jj)=sqrt(amax1(0.,ztz(ii,jj)-zmea(ii,jj)**2))
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          print*,'ok3'
336c pente:
337           zsig(ii,jj)=sqrt(xq)
338c           zsig(ii,jj)=sqrt(2.*xk)
339c isotropy:
340           zgam(ii,jj)=xp/xq
341c angle theta:
342           zthe(ii,jj)=57.29577951*atan2(xm,xl)/2.
343
344c          print 101,ii,jj,
345c    *           zmea(ii,jj),zstd(ii,jj),zsig(ii,jj),zgam(ii,jj),
346c    *           zthe(ii,jj)
347c101  format(1x,2(1x,i2),2(1x,f7.1),1x,f7.4,2x,f4.2,1x,f5.1)     
348c          print*,'ok4'
349         ELSE
350c           PRINT*, 'probleme,ii,jj=', ii,jj
351c          print*,'ok1b'
352         ENDIF
353      zmaxmea=amax1(zmea(ii,jj),zmaxmea)
354c         print*,'oka'
355      zmaxstd=amax1(zstd(ii,jj),zmaxstd)
356c         print*,'okb'
357      zmaxsig=amax1(zsig(ii,jj),zmaxsig)
358c         print*,'okc'
359      zmaxgam=amax1(zgam(ii,jj),zmaxgam)
360c         print*,'okd'
361      zmaxthe=amax1(zthe(ii,jj),zmaxthe)
362c         print*,'oke'
363      zminthe=amin1(zthe(ii,jj),zminthe)
364c      print*,'ok5'
365       ENDDO
366       ENDDO
367
368      print *,'  MEAN ORO:',zmaxmea
369          print *,'  ST. DEV.:',zmaxstd
370      print *,'  PENTE:',zmaxsig
371      print *,' ANISOTROP:',zmaxgam
372      print *,'  ANGLE:',zminthe,zmaxthe       
373     
374C
375c  On passe ce donnees sur la grille dite physique....(?)
376c  On met gamma et theta a 1. et 0. aux poles ou ces quantites
377c  n'ont pas vraiment de sens
378c
379      DO jj=1,jmar
380      zmea(imar+1,jj)=zmea(1,jj)
381      zstd(imar+1,jj)=zstd(1,jj)
382      zsig(imar+1,jj)=zsig(1,jj)
383      zgam(imar+1,jj)=zgam(1,jj)
384      zthe(imar+1,jj)=zthe(1,jj)
385      ENDDO
386
387
388      zmeanor=0.0
389      zmeasud=0.0
390      zstdnor=0.0
391      zstdsud=0.0
392      zsignor=0.0
393      zsigsud=0.0
394      zweinor=0.0
395      zweisud=0.0
396
397      DO ii=1,imar
398      zweinor=zweinor+              weight(ii,   1)
399      zweisud=zweisud+              weight(ii,jmar)
400      zmeanor=zmeanor+zmea(ii,   1)*weight(ii,   1)
401      zmeasud=zmeasud+zmea(ii,jmar)*weight(ii,jmar)
402      zstdnor=zstdnor+zstd(ii,   1)*weight(ii,   1)
403      zstdsud=zstdsud+zstd(ii,jmar)*weight(ii,jmar)
404      zsignor=zsignor+zsig(ii,   1)*weight(ii,   1)
405      zsigsud=zsigsud+zsig(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      zstd(ii,   1)=zstdnor/zweinor
412      zstd(ii,jmar)=zstdsud/zweisud
413      zsig(ii,   1)=zsignor/zweinor
414      zsig(ii,jmar)=zsigsud/zweisud
415      zgam(ii,   1)=1.
416      zgam(ii,jmar)=1.
417      zthe(ii,   1)=0.
418      zthe(ii,jmar)=0.
419      ENDDO
420
421
422      RETURN
423      END
Note: See TracBrowser for help on using the repository browser.