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