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 | !----------------------------------------------------------------------- |
---|
31 | ! INCLUDE 'dimensions.h' |
---|
32 | ! |
---|
33 | ! dimensions.h contient les dimensions du modele |
---|
34 | ! ndm est tel que iim=2**ndm |
---|
35 | !----------------------------------------------------------------------- |
---|
36 | |
---|
37 | INTEGER iim,jjm,llm,ndm |
---|
38 | |
---|
39 | PARAMETER (iim= 128,jjm=96,llm=64,ndm=1) |
---|
40 | |
---|
41 | !----------------------------------------------------------------------- |
---|
42 | ! |
---|
43 | ! $Header$ |
---|
44 | ! |
---|
45 | ! |
---|
46 | ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre |
---|
47 | ! veillez n'utiliser que des ! pour les commentaires |
---|
48 | ! et bien positionner les & des lignes de continuation |
---|
49 | ! (les placer en colonne 6 et en colonne 73) |
---|
50 | ! |
---|
51 | ! |
---|
52 | !----------------------------------------------------------------------- |
---|
53 | ! INCLUDE 'paramet.h' |
---|
54 | |
---|
55 | INTEGER iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1 |
---|
56 | INTEGER kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm |
---|
57 | INTEGER ijmllm,mvar |
---|
58 | INTEGER jcfil,jcfllm |
---|
59 | |
---|
60 | PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3 & |
---|
61 | & ,jjp1=jjm+1-1/jjm) |
---|
62 | PARAMETER( llmp1 = llm+1, llmp2 = llm+2, llmm1 = llm-1 ) |
---|
63 | PARAMETER( kftd = iim/2 -ndm ) |
---|
64 | PARAMETER( ip1jm = iip1*jjm, ip1jmp1= iip1*jjp1 ) |
---|
65 | PARAMETER( ip1jmi1= ip1jm - iip1 ) |
---|
66 | PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm ) |
---|
67 | PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm ) |
---|
68 | PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm ) |
---|
69 | |
---|
70 | !----------------------------------------------------------------------- |
---|
71 | |
---|
72 | INTEGER nmax , nmax2 |
---|
73 | PARAMETER ( nmax = 30000, nmax2 = 2*nmax ) |
---|
74 | c |
---|
75 | c |
---|
76 | c ....... arguments d'entree ....... |
---|
77 | c |
---|
78 | REAL yzoomdeg, grossism,dzooma,tau |
---|
79 | c ( rentres par run.def ) |
---|
80 | |
---|
81 | c ....... arguments de sortie ....... |
---|
82 | c |
---|
83 | REAL rrlatu(jjp1), yyprimu(jjp1),rrlatv(jjm), yyprimv(jjm), |
---|
84 | , rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm) |
---|
85 | |
---|
86 | c |
---|
87 | c ..... champs locaux ..... |
---|
88 | c |
---|
89 | |
---|
90 | REAL dzoom |
---|
91 | REAL(KIND=8) ylat(jjp1), yprim(jjp1) |
---|
92 | REAL(KIND=8) yuv |
---|
93 | REAL(KIND=8) yt(0:nmax2) |
---|
94 | REAL(KIND=8) fhyp(0:nmax2),beta,Ytprim(0:nmax2),fxm(0:nmax2) |
---|
95 | SAVE Ytprim, yt,Yf |
---|
96 | REAL(KIND=8) Yf(0:nmax2),yypr(0:nmax2) |
---|
97 | REAL(KIND=8) yvrai(jjp1), yprimm(jjp1),ylatt(jjp1) |
---|
98 | REAL(KIND=8) pi,depi,pis2,epsilon,y0,pisjm |
---|
99 | REAL(KIND=8) yo1,yi,ylon2,ymoy,Yprimin,champmin,champmax |
---|
100 | REAL(KIND=8) yfi,Yf1,ffdy |
---|
101 | REAL(KIND=8) ypn,deply,y00 |
---|
102 | SAVE y00, deply |
---|
103 | |
---|
104 | INTEGER i,j,it,ik,iter,jlat |
---|
105 | INTEGER jpn,jjpn |
---|
106 | SAVE jpn |
---|
107 | REAL(KIND=8) a0,a1,a2,a3,yi2,heavyy0,heavyy0m |
---|
108 | REAL(KIND=8) fa(0:nmax2),fb(0:nmax2) |
---|
109 | REAL y0min,y0max |
---|
110 | |
---|
111 | REAL(KIND=8) heavyside |
---|
112 | |
---|
113 | pi = 2. * ASIN(1.) |
---|
114 | depi = 2. * pi |
---|
115 | pis2 = pi/2. |
---|
116 | pisjm = pi/ REAL(jjm) |
---|
117 | epsilon = 1.e-3 |
---|
118 | y0 = yzoomdeg * pi/180. |
---|
119 | |
---|
120 | IF( dzooma.LT.1.) THEN |
---|
121 | dzoom = dzooma * pi |
---|
122 | ELSEIF( dzooma.LT. 12. ) THEN |
---|
123 | WRITE(6,*) ' Le param. dzoomy pour fyhyp est trop petit ! L aug |
---|
124 | ,menter et relancer ! ' |
---|
125 | STOP 1 |
---|
126 | ELSE |
---|
127 | dzoom = dzooma * pi/180. |
---|
128 | ENDIF |
---|
129 | |
---|
130 | WRITE(6,18) |
---|
131 | WRITE(6,*) ' yzoom( rad.),grossism,tau,dzoom (radians)' |
---|
132 | WRITE(6,24) y0,grossism,tau,dzoom |
---|
133 | |
---|
134 | DO i = 0, nmax2 |
---|
135 | yt(i) = - pis2 + REAL(i)* pi /nmax2 |
---|
136 | ENDDO |
---|
137 | |
---|
138 | heavyy0m = heavyside( -y0 ) |
---|
139 | heavyy0 = heavyside( y0 ) |
---|
140 | y0min = 2.*y0*heavyy0m - pis2 |
---|
141 | y0max = 2.*y0*heavyy0 + pis2 |
---|
142 | |
---|
143 | fa = 999.999 |
---|
144 | fb = 999.999 |
---|
145 | |
---|
146 | DO i = 0, nmax2 |
---|
147 | IF( yt(i).LT.y0 ) THEN |
---|
148 | fa (i) = tau* (yt(i)-y0+dzoom/2. ) |
---|
149 | fb(i) = (yt(i)-2.*y0*heavyy0m +pis2) * ( y0 - yt(i) ) |
---|
150 | ELSEIF ( yt(i).GT.y0 ) THEN |
---|
151 | fa(i) = tau *(y0-yt(i)+dzoom/2. ) |
---|
152 | fb(i) = (2.*y0*heavyy0 -yt(i)+pis2) * ( yt(i) - y0 ) |
---|
153 | ENDIF |
---|
154 | |
---|
155 | IF( 200.* fb(i) .LT. - fa(i) ) THEN |
---|
156 | fhyp ( i) = - 1. |
---|
157 | ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN |
---|
158 | fhyp ( i) = 1. |
---|
159 | ELSE |
---|
160 | fhyp(i) = TANH ( fa(i)/fb(i) ) |
---|
161 | ENDIF |
---|
162 | |
---|
163 | IF( yt(i).EQ.y0 ) fhyp(i) = 1. |
---|
164 | IF(yt(i).EQ. y0min. OR.yt(i).EQ. y0max ) fhyp(i) = -1. |
---|
165 | |
---|
166 | ENDDO |
---|
167 | |
---|
168 | cc .... Calcul de beta .... |
---|
169 | c |
---|
170 | ffdy = 0. |
---|
171 | |
---|
172 | DO i = 1, nmax2 |
---|
173 | ymoy = 0.5 * ( yt(i-1) + yt( i ) ) |
---|
174 | IF( ymoy.LT.y0 ) THEN |
---|
175 | fa(i)= tau * ( ymoy-y0+dzoom/2.) |
---|
176 | fb(i) = (ymoy-2.*y0*heavyy0m +pis2) * ( y0 - ymoy ) |
---|
177 | ELSEIF ( ymoy.GT.y0 ) THEN |
---|
178 | fa(i)= tau * ( y0-ymoy+dzoom/2. ) |
---|
179 | fb(i) = (2.*y0*heavyy0 -ymoy+pis2) * ( ymoy - y0 ) |
---|
180 | ENDIF |
---|
181 | |
---|
182 | IF( 200.* fb(i) .LT. - fa(i) ) THEN |
---|
183 | fxm ( i) = - 1. |
---|
184 | ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN |
---|
185 | fxm ( i) = 1. |
---|
186 | ELSE |
---|
187 | fxm(i) = TANH ( fa(i)/fb(i) ) |
---|
188 | ENDIF |
---|
189 | IF( ymoy.EQ.y0 ) fxm(i) = 1. |
---|
190 | IF (ymoy.EQ. y0min. OR.yt(i).EQ. y0max ) fxm(i) = -1. |
---|
191 | ffdy = ffdy + fxm(i) * ( yt(i) - yt(i-1) ) |
---|
192 | |
---|
193 | ENDDO |
---|
194 | |
---|
195 | beta = ( grossism * ffdy - pi ) / ( ffdy - pi ) |
---|
196 | |
---|
197 | IF( 2.*beta - grossism.LE. 0.) THEN |
---|
198 | |
---|
199 | WRITE(6,*) ' ** Attention ! La valeur beta calculee dans la rou |
---|
200 | ,tine fyhyp est mauvaise ! ' |
---|
201 | WRITE(6,*)'Modifier les valeurs de grossismy ,tauy ou dzoomy', |
---|
202 | , ' et relancer ! *** ' |
---|
203 | CALL ABORT_GCM("FYHYP", "", 1) |
---|
204 | |
---|
205 | ENDIF |
---|
206 | c |
---|
207 | c ..... calcul de Ytprim ..... |
---|
208 | c |
---|
209 | |
---|
210 | DO i = 0, nmax2 |
---|
211 | Ytprim(i) = beta + ( grossism - beta ) * fhyp(i) |
---|
212 | ENDDO |
---|
213 | |
---|
214 | c ..... Calcul de Yf ........ |
---|
215 | |
---|
216 | Yf(0) = - pis2 |
---|
217 | DO i = 1, nmax2 |
---|
218 | yypr(i) = beta + ( grossism - beta ) * fxm(i) |
---|
219 | ENDDO |
---|
220 | |
---|
221 | DO i=1,nmax2 |
---|
222 | Yf(i) = Yf(i-1) + yypr(i) * ( yt(i) - yt(i-1) ) |
---|
223 | ENDDO |
---|
224 | |
---|
225 | c **************************************************************** |
---|
226 | c |
---|
227 | c ..... yuv = 0. si calcul des latitudes aux pts. U ..... |
---|
228 | c ..... yuv = 0.5 si calcul des latitudes aux pts. V ..... |
---|
229 | c |
---|
230 | WRITE(6,18) |
---|
231 | c |
---|
232 | DO 5000 ik = 1,4 |
---|
233 | |
---|
234 | IF( ik.EQ.1 ) THEN |
---|
235 | yuv = 0. |
---|
236 | jlat = jjm + 1 |
---|
237 | ELSE IF ( ik.EQ.2 ) THEN |
---|
238 | yuv = 0.5 |
---|
239 | jlat = jjm |
---|
240 | ELSE IF ( ik.EQ.3 ) THEN |
---|
241 | yuv = 0.25 |
---|
242 | jlat = jjm |
---|
243 | ELSE IF ( ik.EQ.4 ) THEN |
---|
244 | yuv = 0.75 |
---|
245 | jlat = jjm |
---|
246 | ENDIF |
---|
247 | c |
---|
248 | yo1 = 0. |
---|
249 | DO 1500 j = 1,jlat |
---|
250 | yo1 = 0. |
---|
251 | ylon2 = - pis2 + pisjm * ( REAL(j) + yuv -1.) |
---|
252 | yfi = ylon2 |
---|
253 | c |
---|
254 | DO 250 it = nmax2,0,-1 |
---|
255 | IF( yfi.GE.Yf(it)) GO TO 350 |
---|
256 | 250 CONTINUE |
---|
257 | it = 0 |
---|
258 | 350 CONTINUE |
---|
259 | |
---|
260 | yi = yt(it) |
---|
261 | IF(it.EQ.nmax2) THEN |
---|
262 | it = nmax2 -1 |
---|
263 | Yf(it+1) = pis2 |
---|
264 | ENDIF |
---|
265 | c ................................................................. |
---|
266 | c .... Interpolation entre yi(it) et yi(it+1) pour avoir Y(yi) |
---|
267 | c ..... et Y'(yi) ..... |
---|
268 | c ................................................................. |
---|
269 | |
---|
270 | CALL coefpoly ( Yf(it),Yf(it+1),Ytprim(it), Ytprim(it+1), |
---|
271 | , yt(it),yt(it+1) , a0,a1,a2,a3 ) |
---|
272 | |
---|
273 | Yf1 = Yf(it) |
---|
274 | Yprimin = a1 + 2.* a2 * yi + 3.*a3 * yi *yi |
---|
275 | |
---|
276 | DO 500 iter = 1,300 |
---|
277 | yi = yi - ( Yf1 - yfi )/ Yprimin |
---|
278 | |
---|
279 | IF( ABS(yi-yo1).LE.epsilon) GO TO 550 |
---|
280 | yo1 = yi |
---|
281 | yi2 = yi * yi |
---|
282 | Yf1 = a0 + a1 * yi + a2 * yi2 + a3 * yi2 * yi |
---|
283 | Yprimin = a1 + 2.* a2 * yi + 3.* a3 * yi2 |
---|
284 | 500 CONTINUE |
---|
285 | WRITE(6,*) ' Pas de solution ***** ',j,ylon2,iter |
---|
286 | STOP 2 |
---|
287 | 550 CONTINUE |
---|
288 | c |
---|
289 | Yprimin = a1 + 2.* a2 * yi + 3.* a3 * yi* yi |
---|
290 | yprim(j) = pi / ( jjm * Yprimin ) |
---|
291 | yvrai(j) = yi |
---|
292 | |
---|
293 | 1500 CONTINUE |
---|
294 | |
---|
295 | DO j = 1, jlat -1 |
---|
296 | IF( yvrai(j+1). LT. yvrai(j) ) THEN |
---|
297 | WRITE(6,*) ' PBS. avec rlat(',j+1,') plus petit que rlat(',j, |
---|
298 | , ')' |
---|
299 | STOP 3 |
---|
300 | ENDIF |
---|
301 | ENDDO |
---|
302 | |
---|
303 | WRITE(6,*) 'Reorganisation des latitudes pour avoir entre - pi/2' |
---|
304 | , ,' et pi/2 ' |
---|
305 | c |
---|
306 | IF( ik.EQ.1 ) THEN |
---|
307 | ypn = pis2 |
---|
308 | DO j = jlat,1,-1 |
---|
309 | IF( yvrai(j).LE. ypn ) GO TO 1502 |
---|
310 | ENDDO |
---|
311 | 1502 CONTINUE |
---|
312 | |
---|
313 | jpn = j |
---|
314 | y00 = yvrai(jpn) |
---|
315 | deply = pis2 - y00 |
---|
316 | ENDIF |
---|
317 | |
---|
318 | DO j = 1, jjm +1 - jpn |
---|
319 | ylatt (j) = -pis2 - y00 + yvrai(jpn+j-1) |
---|
320 | yprimm(j) = yprim(jpn+j-1) |
---|
321 | ENDDO |
---|
322 | |
---|
323 | jjpn = jpn |
---|
324 | IF( jlat.EQ. jjm ) jjpn = jpn -1 |
---|
325 | |
---|
326 | DO j = 1,jjpn |
---|
327 | ylatt (j + jjm+1 -jpn) = yvrai(j) + deply |
---|
328 | yprimm(j + jjm+1 -jpn) = yprim(j) |
---|
329 | ENDDO |
---|
330 | |
---|
331 | c *********** Fin de la reorganisation ************* |
---|
332 | c |
---|
333 | 1600 CONTINUE |
---|
334 | |
---|
335 | DO j = 1, jlat |
---|
336 | ylat(j) = ylatt( jlat +1 -j ) |
---|
337 | yprim(j) = yprimm( jlat +1 -j ) |
---|
338 | ENDDO |
---|
339 | |
---|
340 | DO j = 1, jlat |
---|
341 | yvrai(j) = ylat(j)*180./pi |
---|
342 | ENDDO |
---|
343 | |
---|
344 | IF( ik.EQ.1 ) THEN |
---|
345 | c WRITE(6,18) |
---|
346 | c WRITE(6,*) ' YLAT en U apres ( en deg. ) ' |
---|
347 | c WRITE(6,68) (yvrai(j),j=1,jlat) |
---|
348 | cc WRITE(6,*) ' YPRIM ' |
---|
349 | cc WRITE(6,445) ( yprim(j),j=1,jlat) |
---|
350 | |
---|
351 | DO j = 1, jlat |
---|
352 | rrlatu(j) = ylat( j ) |
---|
353 | yyprimu(j) = yprim( j ) |
---|
354 | ENDDO |
---|
355 | |
---|
356 | ELSE IF ( ik.EQ. 2 ) THEN |
---|
357 | c WRITE(6,18) |
---|
358 | c WRITE(6,*) ' YLAT en V apres ( en deg. ) ' |
---|
359 | c WRITE(6,68) (yvrai(j),j=1,jlat) |
---|
360 | cc WRITE(6,*)' YPRIM ' |
---|
361 | cc WRITE(6,445) ( yprim(j),j=1,jlat) |
---|
362 | |
---|
363 | DO j = 1, jlat |
---|
364 | rrlatv(j) = ylat( j ) |
---|
365 | yyprimv(j) = yprim( j ) |
---|
366 | ENDDO |
---|
367 | |
---|
368 | ELSE IF ( ik.EQ. 3 ) THEN |
---|
369 | c WRITE(6,18) |
---|
370 | c WRITE(6,*) ' YLAT en U + 0.75 apres ( en deg. ) ' |
---|
371 | c WRITE(6,68) (yvrai(j),j=1,jlat) |
---|
372 | cc WRITE(6,*) ' YPRIM ' |
---|
373 | cc WRITE(6,445) ( yprim(j),j=1,jlat) |
---|
374 | |
---|
375 | DO j = 1, jlat |
---|
376 | rlatu2(j) = ylat( j ) |
---|
377 | yprimu2(j) = yprim( j ) |
---|
378 | ENDDO |
---|
379 | |
---|
380 | ELSE IF ( ik.EQ. 4 ) THEN |
---|
381 | c WRITE(6,18) |
---|
382 | c WRITE(6,*) ' YLAT en U + 0.25 apres ( en deg. ) ' |
---|
383 | c WRITE(6,68)(yvrai(j),j=1,jlat) |
---|
384 | cc WRITE(6,*) ' YPRIM ' |
---|
385 | cc WRITE(6,68) ( yprim(j),j=1,jlat) |
---|
386 | |
---|
387 | DO j = 1, jlat |
---|
388 | rlatu1(j) = ylat( j ) |
---|
389 | yprimu1(j) = yprim( j ) |
---|
390 | ENDDO |
---|
391 | |
---|
392 | ENDIF |
---|
393 | |
---|
394 | 5000 CONTINUE |
---|
395 | c |
---|
396 | WRITE(6,18) |
---|
397 | c |
---|
398 | c ..... fin de la boucle do 5000 ..... |
---|
399 | |
---|
400 | DO j = 1, jjm |
---|
401 | ylat(j) = rrlatu(j) - rrlatu(j+1) |
---|
402 | ENDDO |
---|
403 | champmin = 1.e12 |
---|
404 | champmax = -1.e12 |
---|
405 | DO j = 1, jjm |
---|
406 | champmin = MIN( champmin, ylat(j) ) |
---|
407 | champmax = MAX( champmax, ylat(j) ) |
---|
408 | ENDDO |
---|
409 | champmin = champmin * 180./pi |
---|
410 | champmax = champmax * 180./pi |
---|
411 | |
---|
412 | 24 FORMAT(2x,'Parametres yzoom,gross,tau ,dzoom pour fyhyp ',4f8.3) |
---|
413 | 18 FORMAT(/) |
---|
414 | 68 FORMAT(1x,7f9.2) |
---|
415 | |
---|
416 | RETURN |
---|
417 | END |
---|