source: trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/grid_noro1.F @ 3937

Last change on this file since 3937 was 3893, checked in by gmilcareck, 3 months ago

Remove all "call abort" and "stop" statement in LMDZ.GENERIC and replacing them by call abort_physic().

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