source: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/grid_noro1.F @ 3302

Last change on this file since 3302 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).

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