source: LMDZ5/branches/LF-private/libf/dyn3dpar/fxhyp.F @ 5440

Last change on this file since 5440 was 1674, checked in by Ehouarn Millour, 12 years ago

Modification pour activation du 2D latitude-pression
(pour pouvoir compiler en -d 1xjmxlm)
dyn3d/fxhyp.F : calcul des longitudes a la main pour iim=1
dyn3d/groupe.F : desactive si iim=1
dyn3d/paramet.h : iip1=iim+1 au lieu de iim+1-1/iim precedemment
phylmd/iophy.F90 : on enleve les -1/iim
phylmd/phyetat0.F90 : on enleve les -1/iim

Modification for activation of the 3D latitude-pressure version
(to be compiled with -d 1xjmxlm)
dyn3d/fxhyp.F : longitudes imposed for iim=1
dyn3d/groupe.F : desactived when iim=1
dyn3d/paramet.h : iip1=iim+1 instead of iim+1-1/iim previously
phylmd/iophy.F90 : -1/iim removed
phylmd/phyetat0.F90 : -1/iim removed

FH et EM

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