source: LMDZ.3.3/branches/LF/libf/dyn3d/fxhyp.F

Last change on this file was 2, checked in by lmdz, 25 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 9.2 KB
Line 
1       SUBROUTINE fxhyp ( xzoomdeg,grossism,dzoom,tau ,
2     , rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025)
3
4c      Auteur :  P. Le Van
5
6       IMPLICIT NONE
7
8c    Calcule les longitudes et derivees dans la grille du GCM pour une
9c     fonction f(x) a tangente  hyperbolique  .
10c
11c     grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.)
12c     dzoom  etant  la distance totale de la zone du zoom
13c     tau  la transition ,   normalement  = 1  .
14c
15
16       INTEGER nmax, nmax2
17       PARAMETER (  nmax = 50000, nmax2 = 2*nmax )
18
19#include "dimensions.h"
20#include "paramet.h"
21       
22c     ......  arguments  d'entree   .......
23c
24       REAL xzoomdeg,dzoom,tau,grossism
25       REAL rlonm025(iip1),xprimm025(iip1),rlonv(iip1),xprimv(iip1),
26     ,  rlonu(iip1),xprimu(iip1),rlonp025(iip1),xprimp025(iip1)
27
28c    ......   arguments  de  sortie  ......
29c
30       REAL xlon(iip1),xprimm(iip1),xuv
31       REAL xtild(0:nmax2)
32       REAL fhyp(0:nmax),ffdx(0:nmax),beta,Xprimt(0:nmax2)
33       REAL Xf(0:nmax2),xxpr(0:nmax2)
34       REAL xvrai(iip1),xxprim(iip1)
35       REAL pi,depi,epsilon,xzoom
36       INTEGER i,it,ik,iter,ii,idif
37       REAL xi,xo1,xint,xmoy,xlon2,fxm,Xprimin
38       REAL champmin,champmax
39       INTEGER is2
40       SAVE is2
41       REAL dlon1(iip1),dlon2(iip1),dlon3(iip1)
42
43       pi       = 2. * ASIN(1.)
44       depi     = 2. * pi
45       epsilon  = 1.e-6
46       xzoom    = xzoomdeg * pi/180.
47
48
49
50       DO i = 0, nmax2
51         xtild(i) = FLOAT(i) /nmax2
52        IF( xtild(i).EQ. 0.5 )  xtild(i) = xtild(i) + 1.e-6
53       ENDDO
54
55       DO i = 1, nmax
56        fhyp(i) = TANH ( ( xtild(i) - 0.5*(1.- dzoom) )          /
57     ,                 ( tau * xtild(i) * ( 0.5 -xtild(i))) )
58       ENDDO
59
60        fhyp(   0  ) = - 1.
61        fhyp( nmax ) =   1.
62
63cc  ....  Calcul  de  beta  ....
64c
65       ffdx( 0 ) = 0.
66
67       DO i = 1, nmax
68        xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
69        fxm     = TANH ( ( xmoy - 0.5 * ( 1. - dzoom ) )    /
70     ,                 ( tau * xmoy * ( 0.5 -xmoy)) )
71        ffdx(i) = ffdx(i-1) + fxm * ( xtild(i) - xtild(i-1) )
72       ENDDO
73
74        beta  = ( grossism * ffdx(nmax) - 0.5 ) / ( ffdx(nmax) - 0.5 )
75c
76c   .....  calcul  de  Xprimt   .....
77c
78       
79       DO i = 0, nmax
80        Xprimt(i) = beta  + ( grossism - beta ) * fhyp(i)
81       ENDDO
82c   
83       DO i = 0, nmax
84        Xprimt( nmax2 - i ) = Xprimt( i )
85       ENDDO
86c
87
88c   .....  Calcul  de  Xf     ........
89
90        Xf(0) = 0.
91       DO i = 1, nmax
92        xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
93        fxm     = TANH ( ( xmoy - 0.5 * ( 1. - dzoom ) )    /
94     ,                 ( tau * xmoy * ( 0.5 -xmoy)) )
95        xxpr(i)    = beta + ( grossism - beta ) * fxm
96       ENDDO
97
98       DO i = 1,nmax
99         xxpr(nmax2-i+1) = xxpr(i)
100       ENDDO
101
102        DO i=1,nmax2
103         Xf(i)   = Xf(i-1) + xxpr(i) * ( xtild(i) - xtild(i-1) )
104        ENDDO
105        do i=1,nmax2
106        xf(i)=xf(i)/xf(nmax2)
107        enddo
108
109
110       PRINT *,' XF ',xf(0),xf(nmax),xf(nmax2)
111
112
113c    *****************************************************************
114c
115
116c     .....  xuv = 0.   si  calcul  aux pts   scalaires   ........
117c     .....  xuv = 0.5  si  calcul  aux pts      U        ........
118c
119c
120      DO 5000  ik = 1, 4
121
122       IF( ik.EQ.1 )        THEN
123         xuv = - 0.25
124       ELSE IF ( ik.EQ.2 )  THEN
125         xuv =   0.
126       ELSE IF ( ik.EQ.3 )  THEN
127         xuv =   0.5
128       ELSE IF ( ik.EQ.4 )  THEN
129         xuv =   0.25
130       ENDIF
131
132
133      DO 1500 i = 1, iim
134
135      xlon2 = (  FLOAT(i) + xuv - 0.75) / FLOAT(iim)
136
137      xo1   = 0.
138      xi    = xlon2
139c
140      DO 500 iter = 1,300
141
142      DO 250 it =  nmax2,0,-1
143      IF( xi.GE.xtild(it))  GO TO 350
144250   CONTINUE
145
146      it = 0
147      xi = xtild(it)
148
149350   CONTINUE
150       IF(it.EQ.nmax2)  THEN
151        it       = nmax2 -1
152        xf(it+1) = 1.
153       ENDIF
154c  .................................................................
155c  ....  Interpolation entre  xi(it) et xi(it+1)   pour avoir X(xi) 
156c      .....           et   X'(xi)                             .....
157c  .................................................................
158
159       xint   = ( Xf(it+1)-Xf(it) ) / ( xtild(it+1)-xtild(it) )      *
160     +                    ( xi-xtild(it) )  +  Xf(it)
161      Xprimin = ( Xprimt(it+1)-Xprimt(it) )/ ( xtild(it+1)-xtild(it) ) *
162     +                    ( xi-xtild(it) )  +  Xprimt(it)
163
164       xi = xi - (xint-xlon2)/Xprimin
165
166      IF( ABS(xi-xo1).LE.epsilon) GO TO 550
167      xo1 = xi
168c
169500   CONTINUE
170      PRINT *,' ***   PAS DE SOLUTION  ****  ',i,xlon2
171        STOP 4
172550   CONTINUE
173
174       xxprim(i)   = depi/( FLOAT(iim) * Xprimin)
175       xvrai(i)    = depi *  (xi - 0.5) + xzoom
176
177
1781500   CONTINUE
179
180       DO i = 1 , iim
181        xlon  (i)   = xvrai(i)
182        xprimm(i)   = xxprim(i)
183cc        xxlon(i) = xlon(i)*180./pi
184       ENDDO
185       
186cc      PRINT *,' XLON avant '
187cc      PRINT 68,(xxlon(i),i=1,iim)
188       
189       DO i = 1, iim -1
190        IF( xvrai(i+1). LT. xvrai(i) )  THEN
191         PRINT *,' PBS.  avec  rlonu(',i+1,' plus petit que rlonu(',i,
192     ,   ')'
193         STOP
194        ENDIF
195       ENDDO
196c
197c   ... Reorganisation  des  longitudes  pour les avoir  entre - pi et pi ..
198c   ........................................................................
199
200       champmin =  1.e12
201       champmax = -1.e12
202       DO i = 1, iim
203        champmin = MIN( champmin,xvrai(i) )
204        champmax = MAX( champmax,xvrai(i) )
205       ENDDO
206
207       PRINT *,' LONGITUDES  min max ',champmin,champmax
208
209      IF(champmin .GE. - pi .AND. champmax.LE. pi )  THEN
210                GO TO 1600
211      ELSE
212        PRINT 18
213        PRINT *,'Reorganisation des longitudes pour avoir entre - pi ',
214     ,  ' et  pi '
215c
216        IF( xzoom.LE.0.)  THEN
217          IF( ik.EQ. 1 )  THEN
218          DO i = 1, iim
219           IF( xvrai(i).GE. - pi )  GO TO 80
220          ENDDO
221            PRINT *, ' PBS. 1 '
222            STOP
223 80       CONTINUE
224          is2 = i
225          ENDIF
226
227          IF( is2.NE. 1 )  THEN
228            DO ii = is2 , iim
229             xlon  (ii-is2+1) = xvrai(ii)
230             xprimm(ii-is2+1) = xxprim(ii)
231            ENDDO
232            DO ii = 1 , is2 -1
233             xlon  (ii+iim-is2+1) = xvrai(ii) + depi
234             xprimm(ii+iim-is2+1) = xxprim(ii)
235            ENDDO
236          ENDIF
237        ELSE
238          IF( ik.EQ.1 )  THEN
239           DO i = iim,1,-1
240            IF( xvrai(i).LE. pi )  GO TO 90
241           ENDDO
242             PRINT *,' PBS.  2 '
243              STOP
244 90        CONTINUE
245            is2 = i
246          ENDIF
247cc           PRINT *,' IS2 ',is2
248           idif = iim -is2
249           DO ii = 1, is2
250            xlon  (ii+idif) = xvrai(ii)
251            xprimm(ii+idif) = xxprim(ii)
252           ENDDO
253           DO ii = 1, idif
254            xlon (ii)  = xvrai (ii+is2) - depi
255            xprimm(ii) = xxprim(ii+is2)
256           ENDDO
257         ENDIF
258      ENDIF
259c
260c     .........   Fin  de la reorganisation   ............................
261
262 1600    CONTINUE
263
264
265         xlon  ( iip1)  = xlon(1) + depi
266         xprimm( iip1 ) = xprimm (1 )
267       
268         DO i = 1, iim+1
269         xvrai(i) = xlon(i)*180./pi
270         ENDDO
271
272
273         IF( ik.EQ.1 )  THEN
274         PRINT *, ' XLON aux pts. V-0.25   apres ( en  deg. ) '
275         PRINT 18
276         PRINT 68,xvrai
277         PRINT *,' XPRIM '
278         PRINT 68, xprimm
279           DO i = 1,iim + 1
280             rlonm025(i) = xlon( i )
281            xprimm025(i) = xprimm(i)
282           ENDDO
283         ELSE IF( ik.EQ.2 )  THEN
284         PRINT 18
285         PRINT *, ' XLON aux pts. V   apres ( en  deg. ) '
286         PRINT 68,xvrai
287         PRINT *,' XPRIM '
288         PRINT 68, xprimm
289           DO i = 1,iim + 1
290             rlonv(i) = xlon( i )
291            xprimv(i) = xprimm(i)
292           ENDDO
293         ELSE IF( ik.EQ.3 )  THEN
294         PRINT 18
295         PRINT *, ' XLON aux pts. U   apres ( en  deg. ) '
296         PRINT 68,xvrai
297         PRINT *,' XPRIM '
298         PRINT 68, xprimm
299           DO i = 1,iim + 1
300             rlonu(i) = xlon( i )
301            xprimu(i) = xprimm(i)
302           ENDDO
303         ELSE IF( ik.EQ.4 )  THEN
304         PRINT 18
305         PRINT *, ' XLON aux pts. V+0.25   apres ( en  deg. ) '
306         PRINT 68,xvrai
307         PRINT *,' XPRIM '
308         PRINT 68, xprimm
309           DO i = 1,iim + 1
310             rlonp025(i) = xlon( i )
311            xprimp025(i) = xprimm(i)
312           ENDDO
313         ENDIF
314
3155000    CONTINUE
316c
317c       ...........  fin  de la boucle  do 5000      ............
318
319c
320        DO i = 1, iim + 1
321         dlon1(i) = rlonm025(i) - rlonv(i)
322         dlon2(i) = rlonm025(i) - rlonp025(i)
323         dlon3(i) = rlonm025(i) - rlonu(i)
324        ENDDO
325
326        DO i = 1, iim + 1
327         rlonm025(i) = rlonm025(i) + dlon1(i)
328        ENDDO
329
330        DO i = 1, iim + 1
331         rlonv(i)    = rlonm025(i) - dlon1(i)
332         rlonp025(i) = rlonm025(i) - dlon2(i)
333         rlonu(i)    = rlonm025(i) - dlon3(i)
334        ENDDO
335
336        DO i = 1, iim
337         xprimu   (i) = rlonu(i+1) - rlonu(i)
338         xprimv   (i) = rlonv(i+1) - rlonv(i)
339         xprimm025(i) = rlonm025(i+1) - rlonm025(i)
340         xprimp025(i) = rlonp025(i+1) - rlonp025(i)
341        ENDDO
342         xprimu   (iip1) = xprimu   (1)
343         xprimv   (iip1) = xprimv   (1)
344         xprimm025(iip1) = xprimm025(1)
345         xprimp025(iip1) = xprimp025(1)
346
347
34818       FORMAT(/)
34968       FORMAT(1x,7f9.2)
350
351
352         RETURN
353         END
Note: See TracBrowser for help on using the repository browser.