source: LMDZ.3.3/branches/rel-LF/libf/dyn3d/fxhyp.F @ 242

Last change on this file since 242 was 242, checked in by lmdzadmin, 23 years ago

Pb d'initialisation surle Nec
LF

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 10.5 KB
Line 
1       SUBROUTINE fxhyp ( xzoomdeg,grossism,dzoom,tau ,
2     , rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025,
3     , champmin,champmax                                               )
4
5c      Auteur :  P. Le Van
6
7       IMPLICIT NONE
8
9c    Calcule les longitudes et derivees dans la grille du GCM pour une
10c     fonction f(x) a tangente  hyperbolique  .
11c
12c     grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.)
13c     dzoom  etant  la distance totale de la zone du zoom
14c     tau  la raideur de la transition de l'interieur a l'exterieur du zoom
15c
16c    On doit avoir grossism x dzoom <  pi ( radians )   , en longitude.
17c   ********************************************************************
18
19
20       INTEGER nmax, nmax2
21       PARAMETER (  nmax = 30000, nmax2 = 2*nmax )
22
23#include "dimensions.h"
24#include "paramet.h"
25       
26c     ......  arguments  d'entree   .......
27c
28       REAL xzoomdeg,dzoom,tau,grossism
29
30c    ......   arguments  de  sortie  ......
31
32       REAL rlonm025(iip1),xprimm025(iip1),rlonv(iip1),xprimv(iip1),
33     ,  rlonu(iip1),xprimu(iip1),rlonp025(iip1),xprimp025(iip1)
34
35c     .... variables locales  ....
36c
37       REAL*8 xlon(iip1),xprimm(iip1),xuv
38       REAL*8 xtild(0:nmax2)
39       REAL*8 fhyp(0:nmax2),ffdx,beta,Xprimt(0:nmax2)
40       REAL*8 Xf(0:nmax2),xxpr(0:nmax2)
41       REAL*8 xvrai(iip1),xxprim(iip1)
42       REAL*8 pi,depi,epsilon,xzoom,fa,fb
43       REAL*8 Xf1, Xfi , a0,a1,a2,a3,xi2
44       INTEGER i,it,ik,iter,ii,idif
45       REAL*8 xi,xo1,xmoy,xlon2,fxm,Xprimin
46       REAL*8 champmin,champmax
47       INTEGER is2
48       SAVE is2
49
50       REAL*8 heavyside
51       EXTERNAL coefpoly,heavyside
52
53       pi       = 2. * ASIN(1.)
54       depi     = 2. * pi
55       epsilon  = 1.e-3
56       xzoom    = xzoomdeg * pi/180.
57
58       IF( dzoom.LT.1.)  THEN
59         dzoom = dzoom * depi
60       ELSEIF( dzoom.LT. 25. ) THEN
61         WRITE(6,*) ' Le param. dzoomy pour fxhyp est trop petit ! L aug
62     ,menter et relancer ! '
63         STOP 1
64       ELSE
65         dzoom = dzoom * pi/180.
66       ENDIF
67
68       WRITE(6,*) ' xzoom( rad.),grossism,tau,dzoom (radians)'
69       WRITE(6,24) xzoom,grossism,tau,dzoom
70
71       DO i = 0, nmax2
72        xtild(i) = - pi + FLOAT(i) * depi /nmax2
73       ENDDO
74
75       DO i = 0, nmax2
76
77       fa  = tau*  ( dzoom/2.  - xtild(i) )
78       fb  = xtild(i) *  ( pi - xtild(i) )
79
80
81       if ( xtild(i) .ne. 0. .and. xtild(i) .ne. pi) then
82         IF( 200.* fb .LT. - fa )   THEN
83           fhyp ( i) = - 1.
84         ELSEIF( 200. * fb .LT. fa ) THEN
85           fhyp ( i) =   1.
86         ELSE
87           fhyp(i) =  TANH ( fa/fb )
88         ENDIF
89       endif
90
91       IF ( xtild(i).EQ. 0. )  fhyp(i) =  1.
92       IF ( xtild(i).EQ. pi )  fhyp(i) = -1.
93
94       ENDDO
95
96cc  ....  Calcul  de  beta  ....
97c   ............................
98
99       ffdx = 0.
100
101       DO i = nmax +1,nmax2
102
103       xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
104       fa  = tau*  ( dzoom/2.  - xmoy )
105       fb  = xmoy *  ( pi - xmoy )
106
107       IF( 200.* fb .LT. - fa )   THEN
108         fxm = - 1.
109       ELSEIF( 200. * fb .LT. fa ) THEN
110         fxm =   1.
111       ELSE
112         fxm =  TANH ( fa/fb )
113       ENDIF
114
115       IF ( xmoy.EQ. 0. )  fhyp(i) =  1.
116       IF ( xmoy.EQ. pi )  fhyp(i) = -1.
117       ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) )
118
119       ENDDO
120
121        beta  = ( grossism * ffdx - pi ) / ( ffdx - pi )
122
123       IF( 2.*beta - grossism.LE. 0.)  THEN
124        WRITE(6,*) ' **  Attention ! La valeur beta calculee dans la rou
125     ,tine fxhyp est mauvaise ! '
126        WRITE(6,*)'Modifier les valeurs de  grossismx ,tau ou dzoomx ',
127     , ' et relancer ! ***  '
128        CALL ABORT
129       ENDIF
130c
131c   .....  calcul  de  Xprimt   .....
132c
133       
134       DO i = nmax, nmax2
135        Xprimt(i) = beta  + ( grossism - beta ) * fhyp(i)
136       ENDDO
137c   
138       DO i =  nmax+1, nmax2
139        Xprimt( nmax2 - i ) = Xprimt( i )
140       ENDDO
141c
142
143c   .....  Calcul  de  Xf     ........
144
145       Xf(0) = - pi
146
147       DO i =  nmax +1, nmax2
148
149       xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
150       fa  = tau*  ( dzoom/2.  - xmoy )
151       fb  = xmoy *  ( pi - xmoy )
152
153       IF( 200.* fb .LT. - fa )   THEN
154         fxm = - 1.
155       ELSEIF( 200. * fb .LT. fa ) THEN
156         fxm =   1.
157       ELSE
158         fxm =  TANH ( fa/fb )
159       ENDIF
160
161       IF ( xmoy.EQ. 0. )  fxm =  1.
162       IF ( xmoy.EQ. pi )  fxm = -1.
163       xxpr(i)    = beta + ( grossism - beta ) * fxm
164
165       ENDDO
166
167       DO i = nmax+1, nmax2
168        xxpr(nmax2-i+1) = xxpr(i)
169       ENDDO
170
171        DO i=1,nmax2
172         Xf(i)   = Xf(i-1) + xxpr(i) * ( xtild(i) - xtild(i-1) )
173        ENDDO
174
175
176c    *****************************************************************
177c
178
179c     .....  xuv = 0.   si  calcul  aux pts   scalaires   ........
180c     .....  xuv = 0.5  si  calcul  aux pts      U        ........
181c
182      WRITE(6,18)
183c
184      DO 5000  ik = 1, 4
185
186       IF( ik.EQ.1 )        THEN
187         xuv = - 0.25
188       ELSE IF ( ik.EQ.2 )  THEN
189         xuv =   0.
190       ELSE IF ( ik.EQ.3 )  THEN
191         xuv =   0.5
192       ELSE IF ( ik.EQ.4 )  THEN
193         xuv =   0.25
194       ENDIF
195
196      xo1   = 0.
197
198      DO 1500 i = 1, iim
199
200      xlon2 = - pi + (FLOAT(i) + xuv - 0.75) * depi / FLOAT(iim)
201
202      Xfi    = xlon2
203c
204      DO 250 it =  nmax2,0,-1
205      IF( Xfi.GE.Xf(it))  GO TO 350
206250   CONTINUE
207
208      it = 0
209
210350   CONTINUE
211
212c    ......  Calcul de   Xf(xi)    ......
213c
214      xi  = xtild(it)
215
216      IF(it.EQ.nmax2)  THEN
217       it       = nmax2 -1
218       Xf(it+1) = pi
219      ENDIF
220c  .....................................................................
221c
222c   Appel de la routine qui calcule les coefficients a0,a1,a2,a3 d'un
223c   polynome de degre 3  qui passe  par les points (Xf(it),xtild(it) )
224c          et (Xf(it+1),xtild(it+1) )
225
226       CALL coefpoly ( Xf(it),Xf(it+1),Xprimt(it),Xprimt(it+1),
227     ,                xtild(it),xtild(it+1),  a0, a1, a2, a3  )
228
229       Xf1     = Xf(it)
230       Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi
231
232       DO 500 iter = 1,300
233        xi = xi - ( Xf1 - Xfi )/ Xprimin
234
235        IF( ABS(xi-xo1).LE.epsilon)  GO TO 550
236         xo1      = xi
237         xi2      = xi * xi
238         Xf1      = a0 +  a1 * xi +     a2 * xi2  +     a3 * xi2 * xi
239         Xprimin  =       a1      + 2.* a2 *  xi  + 3.* a3 * xi2
240500   CONTINUE
241        WRITE(6,*) ' Pas de solution ***** ',i,xlon2,iter
242          STOP 6
243550   CONTINUE
244
245       xxprim(i) = depi/ ( FLOAT(iim) * Xprimin )
246       xvrai(i)  =  xi + xzoom
247
2481500   CONTINUE
249
250       DO i = 1 , iim
251        xlon(i)     = xvrai(i)
252        xprimm(i)   = xxprim(i)
253       ENDDO
254       
255       DO i = 1, iim -1
256        IF( xvrai(i+1). LT. xvrai(i) )  THEN
257         WRITE(6,*) ' PBS. avec rlonu(',i+1,') plus petit que rlonu(',i,
258     ,  ')'
259        STOP 7
260        ENDIF
261       ENDDO
262c
263c   ... Reorganisation  des  longitudes  pour les avoir  entre - pi et pi ..
264c   ........................................................................
265
266       champmin =  1.e12
267       champmax = -1.e12
268       DO i = 1, iim
269        champmin = MIN( champmin,xvrai(i) )
270        champmax = MAX( champmax,xvrai(i) )
271       ENDDO
272
273      IF(champmin .GE. - pi .AND. champmax.LE. pi )  THEN
274                GO TO 1600
275      ELSE
276       WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi',
277     ,  ' et pi '
278c
279        IF( xzoom.LE.0.)  THEN
280          IF( ik.EQ. 1 )  THEN
281          DO i = 1, iim
282           IF( xvrai(i).GE. - pi )  GO TO 80
283          ENDDO
284            WRITE(6,*)  ' PBS. 1 !  Xvrai plus petit que  - pi ! '
285            STOP 8
286 80       CONTINUE
287          is2 = i
288          ENDIF
289
290          IF( is2.NE. 1 )  THEN
291            DO ii = is2 , iim
292             xlon  (ii-is2+1) = xvrai(ii)
293             xprimm(ii-is2+1) = xxprim(ii)
294            ENDDO
295            DO ii = 1 , is2 -1
296             xlon  (ii+iim-is2+1) = xvrai(ii) + depi
297             xprimm(ii+iim-is2+1) = xxprim(ii)
298            ENDDO
299          ENDIF
300        ELSE
301          IF( ik.EQ.1 )  THEN
302           DO i = iim,1,-1
303            IF( xvrai(i).LE. pi )  GO TO 90
304           ENDDO
305             WRITE(6,*) ' PBS.  2 ! Xvrai plus grand  que   pi ! '
306              STOP 9
307 90        CONTINUE
308            is2 = i
309          ENDIF
310           idif = iim -is2
311           DO ii = 1, is2
312            xlon  (ii+idif) = xvrai(ii)
313            xprimm(ii+idif) = xxprim(ii)
314           ENDDO
315           DO ii = 1, idif
316            xlon (ii)  = xvrai (ii+is2) - depi
317            xprimm(ii) = xxprim(ii+is2)
318           ENDDO
319         ENDIF
320      ENDIF
321c
322c     .........   Fin  de la reorganisation   ............................
323
324 1600    CONTINUE
325
326
327         xlon  ( iip1)  = xlon(1) + depi
328         xprimm( iip1 ) = xprimm (1 )
329       
330         DO i = 1, iim+1
331         xvrai(i) = xlon(i)*180./pi
332         ENDDO
333
334         IF( ik.EQ.1 )  THEN
335c          WRITE(6,*)  ' XLON aux pts. V-0.25   apres ( en  deg. ) '
336c          WRITE(6,18)
337c          WRITE(6,68) xvrai
338ccc          WRITE(6,*) ' XPRIM k ',ik
339ccc          WRITE(6,566)  xprimm
340
341           DO i = 1,iim + 1
342             rlonm025(i) = xlon( i )
343            xprimm025(i) = xprimm(i)
344           ENDDO
345
346         ELSE IF( ik.EQ.2 )  THEN
347c          WRITE(6,18)
348c          WRITE(6,*)  ' XLON aux pts. V   apres ( en  deg. ) '
349c          WRITE(6,68) xvrai
350ccc          WRITE(6,*) ' XPRIM k ',ik
351ccc          WRITE(6,566)  xprimm
352
353           DO i = 1,iim + 1
354             rlonv(i) = xlon( i )
355            xprimv(i) = xprimm(i)
356           ENDDO
357
358         ELSE IF( ik.EQ.3)   THEN
359c          WRITE(6,18)
360c          WRITE(6,*)  ' XLON aux pts. U   apres ( en  deg. ) '
361c          WRITE(6,68) xvrai
362ccc          WRITE(6,*) ' XPRIM ik ',ik
363ccc          WRITE(6,566)  xprimm
364
365           DO i = 1,iim + 1
366             rlonu(i) = xlon( i )
367            xprimu(i) = xprimm(i)
368           ENDDO
369
370         ELSE IF( ik.EQ.4 )  THEN
371c          WRITE(6,18)
372c          WRITE(6,*)  ' XLON aux pts. V+0.25   apres ( en  deg. ) '
373c          WRITE(6,68) xvrai
374ccc          WRITE(6,*) ' XPRIM ik ',ik
375ccc          WRITE(6,566)  xprimm
376
377           DO i = 1,iim + 1
378             rlonp025(i) = xlon( i )
379            xprimp025(i) = xprimm(i)
380           ENDDO
381
382         ENDIF
383
3845000    CONTINUE
385c
386       WRITE(6,18)
387c
388c    ...........  fin  de la boucle  do 5000      ............
389
390        DO i = 1, iim
391         xlon(i) = rlonv(i+1) - rlonv(i)
392        ENDDO
393        champmin =  1.e12
394        champmax = -1.e12
395        DO i = 1, iim
396         champmin = MIN( champmin, xlon(i) )
397         champmax = MAX( champmax, xlon(i) )
398        ENDDO
399         champmin = champmin * 180./pi
400         champmax = champmax * 180./pi
401
40218     FORMAT(/)
40324     FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3)
40468     FORMAT(1x,7f9.2)
405
406       RETURN
407       END
Note: See TracBrowser for help on using the repository browser.