source: LMDZ.3.3/trunk/libf/dyn3d/fyhyp.F @ 507

Last change on this file since 507 was 251, checked in by lmdz, 23 years ago

Initialisations pour eviter des divisions par 0
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.0 KB
Line 
1C
2c $Header$
3c
4       SUBROUTINE fyhyp ( yzoomdeg, grossism, dzoom,tau  , 
5     ,  rrlatu,yyprimu,rrlatv,yyprimv,rlatu2,yprimu2,rlatu1,yprimu1 ,
6     ,  champmin,champmax                                            )
7
8cc    ...  Version du 01/04/2001 ....
9
10       IMPLICIT NONE
11c
12c    ...   Auteur :  P. Le Van  ...
13c
14c    .......    d'apres  formulations  de R. Sadourny  .......
15c
16c     Calcule les latitudes et derivees dans la grille du GCM pour une
17c     fonction f(y) a tangente  hyperbolique  .
18c
19c     grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois , etc)
20c     dzoom  etant  la distance totale de la zone du zoom ( en radians )
21c     tau  la raideur de la transition de l'interieur a l'exterieur du zoom   
22c
23c
24c N.B : Il vaut mieux avoir : grossism * dzoom  <  pi/2  (radians) ,en lati.
25c      ********************************************************************
26c
27c
28#include "dimensions.h"
29#include "paramet.h"
30
31       INTEGER      nmax , nmax2
32       PARAMETER (  nmax = 30000, nmax2 = 2*nmax )
33c
34c
35c     .......  arguments  d'entree    .......
36c
37       REAL yzoomdeg, grossism,dzoom,tau
38c         ( rentres  par  run.def )
39
40c     .......  arguments  de sortie   .......
41c
42       REAL rrlatu(jjp1), yyprimu(jjp1),rrlatv(jjm), yyprimv(jjm),
43     , rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
44
45c
46c     .....     champs  locaux    .....
47c
48     
49       REAL*8 ylat(jjp1), yprim(jjp1)
50       REAL*8 yuv
51       REAL*8 yt(0:nmax2)
52       REAL*8 fhyp(0:nmax2),beta,Ytprim(0:nmax2),fxm(0:nmax2)
53       SAVE Ytprim, yt,Yf
54       REAL*8 Yf(0:nmax2),yypr(0:nmax2)
55       REAL*8 yvrai(jjp1), yprimm(jjp1),ylatt(jjp1)
56       REAL*8 pi,depi,pis2,epsilon,y0,pisjm
57       REAL*8 yo1,yi,ylon2,ymoy,Yprimin,champmin,champmax
58       REAL*8 yfi,Yf1,ffdy
59       REAL*8 ypn,deply,y00
60       SAVE y00, deply
61
62       INTEGER i,j,it,ik,iter,jlat
63       INTEGER jpn,jjpn
64       SAVE jpn
65       REAL*8 a0,a1,a2,a3,yi2,heavyy0,heavyy0m
66       REAL*8 fa(0:nmax2),fb(0:nmax2)
67       REAL y0min,y0max
68
69       REAL*8     heavyside
70       EXTERNAL   heavyside
71
72       pi       = 2. * ASIN(1.)
73       depi     = 2. * pi
74       pis2     = pi/2.
75       pisjm    = pi/ FLOAT(jjm)
76       epsilon  = 1.e-3
77       y0       =  yzoomdeg * pi/180.
78
79       IF( dzoom.LT.1.)  THEN
80         dzoom = dzoom * pi
81       ELSEIF( dzoom.LT. 12. ) THEN
82         WRITE(6,*) ' Le param. dzoomy pour fyhyp est trop petit ! L aug
83     ,menter et relancer ! '
84         STOP 1
85       ELSE
86         dzoom = dzoom * pi/180.
87       ENDIF
88
89       WRITE(6,18)
90       WRITE(6,*) ' yzoom( rad.),grossism,tau,dzoom (radians)'
91       WRITE(6,24) y0,grossism,tau,dzoom
92
93       DO i = 0, nmax2
94        yt(i) = - pis2  + FLOAT(i)* pi /nmax2
95       ENDDO
96
97       heavyy0m = heavyside( -y0 )
98       heavyy0  = heavyside(  y0 )
99       y0min    = 2.*y0*heavyy0m - pis2
100       y0max    = 2.*y0*heavyy0  + pis2
101
102       fa = 999.999
103       fb = 999.999
104
105       DO i = 0, nmax2
106        IF( yt(i).LT.y0 )  THEN
107         fa (i) = tau*  (yt(i)-y0+dzoom/2. )
108         fb(i) =   (yt(i)-2.*y0*heavyy0m +pis2) * ( y0 - yt(i) )
109        ELSEIF ( yt(i).GT.y0 )  THEN
110         fa(i) =   tau *(y0-yt(i)+dzoom/2. )
111         fb(i) = (2.*y0*heavyy0 -yt(i)+pis2) * ( yt(i) - y0 )
112       ENDIF
113       
114       IF( 200.* fb(i) .LT. - fa(i) )   THEN
115         fhyp ( i) = - 1.
116       ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN
117         fhyp ( i) =   1.
118       ELSE 
119         fhyp(i) =  TANH ( fa(i)/fb(i) )
120       ENDIF
121
122       IF( yt(i).EQ.y0 )  fhyp(i) = 1.
123       IF(yt(i).EQ. y0min. OR.yt(i).EQ. y0max ) fhyp(i) = -1.
124
125       ENDDO
126
127cc  ....  Calcul  de  beta  ....
128c
129       ffdy   = 0.
130
131       DO i = 1, nmax2
132        ymoy    = 0.5 * ( yt(i-1) + yt( i ) )
133        IF( ymoy.LT.y0 )  THEN
134         fa(i)= tau * ( ymoy-y0+dzoom/2.)
135         fb(i) = (ymoy-2.*y0*heavyy0m +pis2) * ( y0 - ymoy )
136        ELSEIF ( ymoy.GT.y0 )  THEN
137         fa(i)= tau * ( y0-ymoy+dzoom/2. )
138         fb(i) = (2.*y0*heavyy0 -ymoy+pis2) * ( ymoy - y0 )
139        ENDIF
140
141        IF( 200.* fb(i) .LT. - fa(i) )    THEN
142         fxm ( i) = - 1.
143        ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN
144         fxm ( i) =   1.
145        ELSE
146         fxm(i) =  TANH ( fa(i)/fb(i) )
147        ENDIF
148         IF( ymoy.EQ.y0 )  fxm(i) = 1.
149         IF (ymoy.EQ. y0min. OR.yt(i).EQ. y0max ) fxm(i) = -1.
150         ffdy = ffdy + fxm(i) * ( yt(i) - yt(i-1) )
151
152        ENDDO
153
154        beta  = ( grossism * ffdy - pi ) / ( ffdy - pi )
155
156       IF( 2.*beta - grossism.LE. 0.)  THEN
157
158        WRITE(6,*) ' **  Attention ! La valeur beta calculee dans la rou
159     ,tine fyhyp est mauvaise ! '
160        WRITE(6,*)'Modifier les valeurs de  grossismy ,tauy ou dzoomy',
161     , ' et relancer ! ***  '
162        CALL ABORT
163
164       ENDIF
165c
166c   .....  calcul  de  Ytprim   .....
167c
168       
169       DO i = 0, nmax2
170        Ytprim(i) = beta  + ( grossism - beta ) * fhyp(i)
171       ENDDO
172
173c   .....  Calcul  de  Yf     ........
174
175       Yf(0) = - pis2
176       DO i = 1, nmax2
177        yypr(i)    = beta + ( grossism - beta ) * fxm(i)
178       ENDDO
179
180       DO i=1,nmax2
181        Yf(i)   = Yf(i-1) + yypr(i) * ( yt(i) - yt(i-1) )
182       ENDDO
183
184c    ****************************************************************
185c
186c   .....   yuv  = 0.   si calcul des latitudes  aux pts.  U  .....
187c   .....   yuv  = 0.5  si calcul des latitudes  aux pts.  V  .....
188c
189      WRITE(6,18)
190c
191      DO 5000  ik = 1,4
192
193       IF( ik.EQ.1 )  THEN
194         yuv  = 0.
195         jlat = jjm + 1
196       ELSE IF ( ik.EQ.2 )  THEN
197         yuv  = 0.5
198         jlat = jjm
199       ELSE IF ( ik.EQ.3 )  THEN
200         yuv  = 0.25
201         jlat = jjm
202       ELSE IF ( ik.EQ.4 )  THEN
203         yuv  = 0.75
204         jlat = jjm
205       ENDIF
206c
207       yo1   = 0.
208       DO 1500 j =  1,jlat
209        yo1   = 0.
210        ylon2 =  - pis2 + pisjm * ( FLOAT(j)  + yuv  -1.) 
211        yfi    = ylon2
212c
213       DO 250 it =  nmax2,0,-1
214        IF( yfi.GE.Yf(it))  GO TO 350
215250    CONTINUE
216       it = 0
217350    CONTINUE
218
219       yi = yt(it)
220       IF(it.EQ.nmax2)  THEN
221        it       = nmax2 -1
222        Yf(it+1) = pis2
223       ENDIF
224c  .................................................................
225c  ....  Interpolation entre  yi(it) et yi(it+1)   pour avoir Y(yi) 
226c      .....           et   Y'(yi)                             .....
227c  .................................................................
228
229       CALL coefpoly ( Yf(it),Yf(it+1),Ytprim(it), Ytprim(it+1),   
230     ,                  yt(it),yt(it+1) ,   a0,a1,a2,a3   )     
231
232       Yf1     = Yf(it)
233       Yprimin = a1 + 2.* a2 * yi + 3.*a3 * yi *yi
234
235       DO 500 iter = 1,300
236         yi = yi - ( Yf1 - yfi )/ Yprimin
237
238        IF( ABS(yi-yo1).LE.epsilon)  GO TO 550
239         yo1      = yi
240         yi2      = yi * yi
241         Yf1      = a0 +  a1 * yi +     a2 * yi2  +     a3 * yi2 * yi
242         Yprimin  =       a1      + 2.* a2 *  yi  + 3.* a3 * yi2
243500   CONTINUE
244        WRITE(6,*) ' Pas de solution ***** ',j,ylon2,iter
245         STOP 2
246550   CONTINUE
247c
248       Yprimin   = a1  + 2.* a2 *  yi   + 3.* a3 * yi* yi
249       yprim(j)  = pi / ( jjm * Yprimin )
250       yvrai(j)  = yi
251
2521500    CONTINUE
253
254       DO j = 1, jlat -1
255        IF( yvrai(j+1). LT. yvrai(j) )  THEN
256         WRITE(6,*) ' PBS. avec  rlat(',j+1,') plus petit que rlat(',j,
257     ,  ')'
258         STOP 3
259        ENDIF
260       ENDDO
261
262       WRITE(6,*) 'Reorganisation des latitudes pour avoir entre - pi/2'
263     , ,' et  pi/2 '
264c
265        IF( ik.EQ.1 )   THEN
266           ypn = pis2
267          DO j = jlat,1,-1
268           IF( yvrai(j).LE. ypn ) GO TO 1502
269          ENDDO
2701502     CONTINUE
271
272         jpn   = j
273         y00   = yvrai(jpn)
274         deply = pis2 -  y00
275        ENDIF
276
277         DO  j = 1, jjm +1 - jpn
278           ylatt (j)  = -pis2 - y00  + yvrai(jpn+j-1)
279           yprimm(j)  = yprim(jpn+j-1)
280         ENDDO
281
282         jjpn  = jpn
283         IF( jlat.EQ. jjm ) jjpn = jpn -1
284
285         DO j = 1,jjpn
286          ylatt (j + jjm+1 -jpn) = yvrai(j) + deply
287          yprimm(j + jjm+1 -jpn) = yprim(j)
288         ENDDO
289
290c      ***********   Fin de la reorganisation     *************
291c
292 1600   CONTINUE
293
294       DO j = 1, jlat
295          ylat(j) =  ylatt( jlat +1 -j )
296         yprim(j) = yprimm( jlat +1 -j )
297       ENDDO
298 
299        DO j = 1, jlat
300         yvrai(j) = ylat(j)*180./pi
301        ENDDO
302
303        IF( ik.EQ.1 )  THEN
304c         WRITE(6,18)
305c         WRITE(6,*)  ' YLAT  en U   apres ( en  deg. ) '
306c         WRITE(6,68) (yvrai(j),j=1,jlat)
307cc         WRITE(6,*) ' YPRIM '
308cc         WRITE(6,445) ( yprim(j),j=1,jlat)
309
310          DO j = 1, jlat
311            rrlatu(j) =  ylat( j )
312           yyprimu(j) = yprim( j )
313          ENDDO
314
315        ELSE IF ( ik.EQ. 2 )  THEN
316c         WRITE(6,18)
317c         WRITE(6,*) ' YLAT   en V  apres ( en  deg. ) '
318c         WRITE(6,68) (yvrai(j),j=1,jlat)
319cc         WRITE(6,*)' YPRIM '
320cc         WRITE(6,445) ( yprim(j),j=1,jlat)
321
322          DO j = 1, jlat
323            rrlatv(j) =  ylat( j )
324           yyprimv(j) = yprim( j )
325          ENDDO
326
327        ELSE IF ( ik.EQ. 3 )  THEN
328c         WRITE(6,18)
329c         WRITE(6,*)  ' YLAT  en U + 0.75  apres ( en  deg. ) '
330c         WRITE(6,68) (yvrai(j),j=1,jlat)
331cc         WRITE(6,*) ' YPRIM '
332cc         WRITE(6,445) ( yprim(j),j=1,jlat)
333
334          DO j = 1, jlat
335            rlatu2(j) =  ylat( j )
336           yprimu2(j) = yprim( j )
337          ENDDO
338
339        ELSE IF ( ik.EQ. 4 )  THEN
340c         WRITE(6,18)
341c         WRITE(6,*)  ' YLAT en U + 0.25  apres ( en  deg. ) '
342c         WRITE(6,68)(yvrai(j),j=1,jlat)
343cc         WRITE(6,*) ' YPRIM '
344cc         WRITE(6,68) ( yprim(j),j=1,jlat)
345
346          DO j = 1, jlat
347            rlatu1(j) =  ylat( j )
348           yprimu1(j) = yprim( j )
349          ENDDO
350
351        ENDIF
352
3535000   CONTINUE
354c
355        WRITE(6,18)
356c
357c  .....     fin de la boucle  do 5000 .....
358
359        DO j = 1, jjm
360         ylat(j) = rrlatu(j) - rrlatu(j+1)
361        ENDDO
362        champmin =  1.e12
363        champmax = -1.e12
364        DO j = 1, jjm
365         champmin = MIN( champmin, ylat(j) )
366         champmax = MAX( champmax, ylat(j) )
367        ENDDO
368         champmin = champmin * 180./pi
369         champmax = champmax * 180./pi
370
37124     FORMAT(2x,'Parametres yzoom,gross,tau ,dzoom pour fyhyp ',4f8.3)
37218      FORMAT(/)
37368      FORMAT(1x,7f9.2)
374
375        RETURN
376        END
Note: See TracBrowser for help on using the repository browser.