/[lmdze]/trunk/dyn3d/fyhyp.f90
ViewVC logotype

Annotation of /trunk/dyn3d/fyhyp.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 81 - (hide annotations)
Wed Mar 5 14:38:41 2014 UTC (10 years, 2 months ago) by guez
File size: 9134 byte(s) 
 Converted to free source form files which were still in fixed source
form. The conversion was done using the polish mode of the NAG Fortran
Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

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