/[lmdze]/trunk/Sources/dyn3d/fyhyp.f
ViewVC logotype

Annotation of /trunk/Sources/dyn3d/fyhyp.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 139 - (hide annotations)
Tue May 26 17:46:03 2015 UTC (9 years ago) by guez
File size: 8914 byte(s) 
dynetat0 read rlonu, rlatu, rlonv, rlatv, cu_2d, cv_2d, aire_2d from
"start.nc" and then these variables were overwritten by
inigeom. Corrected this. Now, inigeom does not compute rlonu, rlatu,
rlonv and rlatv. Moreover, cu_2d, cv_2d, aire_2d are not written to
"restart.nc". Since xprimu, xprimv, xprimm025, xprimp025, rlatu1,
rlatu2, yprimu1, yprimu2 are computed at the same time as rlonu,
rlatu, rlonv, rlatv, and since it would not be convenient to separate
those computations, we decide to write xprimu, xprimv, xprimm025,
xprimp025, rlatu1, rlatu2, yprimu1, yprimu2 into "restart.nc", read
them from "start.nc" and not compute them in inigeom. So, in summary,
"start.nc" contains all the coordinates and their derivatives, and
inigeom only computes the 2D-variables.

Technical details:

Moved variables rlatu, rlonv, rlonu, rlatv, xprimu, xprimv from module
comgeom to module dynetat0_m. Upgraded local variables rlatu1,
yprimu1, rlatu2, yprimu2, xprimm025, xprimp025 of procedure inigeom to
variables of module dynetat0_m.

Removed unused local variable yprimu of procedure inigeom and
corresponding argument yyprimu of fyhyp.

Moved variables clat, clon, grossismx, grossismy, dzoomx, dzoomy,
taux, tauy from module serre to module dynetat0_m (since they are read
from "start.nc"). The default values are now defined in read_serre
instead of in the declarations. Changed name of module serre to
read_serre_m, no more module variable here.

The calls to fxhyp and fyhyp are moved from inigeom to etat0.

Side effects in programs other than gcm: etat0 and read_serre write
variables of module dynetat0; the programs test_fxyp and
test_inter_barxy need more source files.

Removed unused arguments len and nd of cv3_tracer. Removed unused
argument PPSOL of LWU.

Bug fix in test_inter_barxy: forgotten call to read_serre.
1 guez 97 module fyhyp_m
2 guez 3
3 guez 112 IMPLICIT NONE
4 guez 3
5 guez 97 contains
6 guez 3
7 guez 139 SUBROUTINE fyhyp(rlatu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)
8 guez 3
9 guez 97 ! From LMDZ4/libf/dyn3d/fyhyp.F, version 1.2, 2005/06/03 09:11:32
10 guez 3
11 guez 97 ! Author: P. Le Van, from analysis by R. Sadourny
12 guez 3
13 guez 97 ! Calcule les latitudes et dérivées dans la grille du GCM pour une
14 guez 119 ! fonction f(y) à dérivée tangente hyperbolique.
15 guez 3
16 guez 121 ! Il vaut mieux avoir : grossismy * dzoom < pi / 2
17 guez 3
18 guez 120 use coefpoly_m, only: coefpoly
19 guez 97 USE dimens_m, only: jjm
20 guez 139 use dynetat0_m, only: clat, grossismy, dzoomy, tauy
21 guez 131 use heavyside_m, only: heavyside
22 guez 3
23 guez 139 REAL, intent(out):: rlatu(jjm + 1)
24 guez 119 REAL, intent(out):: rlatv(jjm)
25 guez 118 real, intent(out):: rlatu2(jjm), yprimu2(jjm), rlatu1(jjm), yprimu1(jjm)
26 guez 3
27 guez 97 ! Local:
28 guez 3
29 guez 119 DOUBLE PRECISION champmin, champmax
30 guez 97 INTEGER, PARAMETER:: nmax=30000, nmax2=2*nmax
31     REAL dzoom ! distance totale de la zone du zoom (en radians)
32 guez 118 DOUBLE PRECISION ylat(jjm + 1), yprim(jjm + 1)
33 guez 97 DOUBLE PRECISION yuv
34     DOUBLE PRECISION, save:: yt(0:nmax2)
35     DOUBLE PRECISION fhyp(0:nmax2), beta
36     DOUBLE PRECISION, save:: ytprim(0:nmax2)
37     DOUBLE PRECISION fxm(0:nmax2)
38     DOUBLE PRECISION, save:: yf(0:nmax2)
39     DOUBLE PRECISION yypr(0:nmax2)
40 guez 118 DOUBLE PRECISION yvrai(jjm + 1), yprimm(jjm + 1), ylatt(jjm + 1)
41 guez 131 DOUBLE PRECISION pi, pis2, epsilon, pisjm
42 guez 97 DOUBLE PRECISION yo1, yi, ylon2, ymoy, yprimin
43     DOUBLE PRECISION yfi, yf1, ffdy
44     DOUBLE PRECISION ypn, deply, y00
45     SAVE y00, deply
46 guez 3
47 guez 97 INTEGER i, j, it, ik, iter, jlat
48     INTEGER jpn, jjpn
49     SAVE jpn
50     DOUBLE PRECISION a0, a1, a2, a3, yi2, heavyy0, heavyy0m
51     DOUBLE PRECISION fa(0:nmax2), fb(0:nmax2)
52     REAL y0min, y0max
53 guez 3
54 guez 97 !-------------------------------------------------------------------
55 guez 3
56 guez 120 print *, "Call sequence information: fyhyp"
57    
58 guez 97 pi = 2.*asin(1.)
59     pis2 = pi/2.
60     pisjm = pi/real(jjm)
61     epsilon = 1e-3
62 guez 121 dzoom = dzoomy*pi
63 guez 119 print *, 'yzoom(rad), grossismy, tauy, dzoom (rad):'
64 guez 131 print *, clat, grossismy, tauy, dzoom
65 guez 81
66 guez 97 DO i = 0, nmax2
67     yt(i) = -pis2 + real(i)*pi/nmax2
68     END DO
69 guez 81
70 guez 131 heavyy0m = heavyside(-clat)
71     heavyy0 = heavyside(clat)
72     y0min = 2.*clat*heavyy0m - pis2
73     y0max = 2.*clat*heavyy0 + pis2
74 guez 81
75 guez 97 fa = 999.999
76     fb = 999.999
77 guez 81
78 guez 97 DO i = 0, nmax2
79 guez 131 IF (yt(i)<clat) THEN
80     fa(i) = tauy*(yt(i)-clat + dzoom/2.)
81     fb(i) = (yt(i)-2.*clat*heavyy0m + pis2)*(clat-yt(i))
82     ELSE IF (yt(i)>clat) THEN
83     fa(i) = tauy*(clat-yt(i) + dzoom/2.)
84     fb(i) = (2.*clat*heavyy0-yt(i) + pis2)*(yt(i)-clat)
85 guez 97 END IF
86 guez 81
87 guez 97 IF (200.*fb(i)<-fa(i)) THEN
88     fhyp(i) = -1.
89     ELSE IF (200.*fb(i)<fa(i)) THEN
90     fhyp(i) = 1.
91     ELSE
92     fhyp(i) = tanh(fa(i)/fb(i))
93     END IF
94 guez 81
95 guez 131 IF (yt(i)==clat) fhyp(i) = 1.
96 guez 97 IF (yt(i)==y0min .OR. yt(i)==y0max) fhyp(i) = -1.
97     END DO
98 guez 81
99 guez 97 ! Calcul de beta
100 guez 81
101 guez 97 ffdy = 0.
102 guez 81
103 guez 97 DO i = 1, nmax2
104     ymoy = 0.5*(yt(i-1) + yt(i))
105 guez 131 IF (ymoy<clat) THEN
106     fa(i) = tauy*(ymoy-clat + dzoom/2.)
107     fb(i) = (ymoy-2.*clat*heavyy0m + pis2)*(clat-ymoy)
108     ELSE IF (ymoy>clat) THEN
109     fa(i) = tauy*(clat-ymoy + dzoom/2.)
110     fb(i) = (2.*clat*heavyy0-ymoy + pis2)*(ymoy-clat)
111 guez 97 END IF
112 guez 81
113 guez 97 IF (200.*fb(i)<-fa(i)) THEN
114     fxm(i) = -1.
115     ELSE IF (200.*fb(i)<fa(i)) THEN
116     fxm(i) = 1.
117     ELSE
118     fxm(i) = tanh(fa(i)/fb(i))
119     END IF
120 guez 131 IF (ymoy==clat) fxm(i) = 1.
121 guez 97 IF (ymoy==y0min .OR. yt(i)==y0max) fxm(i) = -1.
122     ffdy = ffdy + fxm(i)*(yt(i)-yt(i-1))
123     END DO
124 guez 81
125 guez 119 beta = (grossismy*ffdy-pi)/(ffdy-pi)
126 guez 81
127 guez 119 IF (2. * beta - grossismy <= 0.) THEN
128 guez 97 print *, 'Attention ! La valeur beta calculee dans la routine fyhyp ' &
129     // 'est mauvaise. Modifier les valeurs de grossismy, tauy ou ' &
130     // 'dzoomy et relancer.'
131     STOP 1
132 guez 81 END IF
133    
134 guez 97 ! calcul de Ytprim
135 guez 81
136 guez 97 DO i = 0, nmax2
137 guez 119 ytprim(i) = beta + (grossismy-beta)*fhyp(i)
138 guez 97 END DO
139 guez 81
140 guez 97 ! Calcul de Yf
141 guez 81
142 guez 97 yf(0) = -pis2
143     DO i = 1, nmax2
144 guez 119 yypr(i) = beta + (grossismy-beta)*fxm(i)
145 guez 81 END DO
146 guez 3
147 guez 97 DO i = 1, nmax2
148     yf(i) = yf(i-1) + yypr(i)*(yt(i)-yt(i-1))
149 guez 81 END DO
150 guez 3
151 guez 97 ! yuv = 0. si calcul des latitudes aux pts. U
152     ! yuv = 0.5 si calcul des latitudes aux pts. V
153 guez 3
154 guez 97 loop_ik: DO ik = 1, 4
155     IF (ik==1) THEN
156     yuv = 0.
157     jlat = jjm + 1
158     ELSE IF (ik==2) THEN
159     yuv = 0.5
160     jlat = jjm
161     ELSE IF (ik==3) THEN
162     yuv = 0.25
163     jlat = jjm
164     ELSE IF (ik==4) THEN
165     yuv = 0.75
166     jlat = jjm
167     END IF
168 guez 3
169 guez 97 yo1 = 0.
170     DO j = 1, jlat
171     yo1 = 0.
172     ylon2 = -pis2 + pisjm*(real(j) + yuv-1.)
173     yfi = ylon2
174 guez 3
175 guez 97 it = nmax2
176     DO while (it >= 1 .and. yfi < yf(it))
177     it = it - 1
178     END DO
179 guez 3
180 guez 97 yi = yt(it)
181     IF (it==nmax2) THEN
182     it = nmax2 - 1
183     yf(it + 1) = pis2
184     END IF
185 guez 3
186 guez 97 ! Interpolation entre yi(it) et yi(it + 1) pour avoir Y(yi)
187     ! et Y'(yi)
188 guez 3
189 guez 97 CALL coefpoly(yf(it), yf(it + 1), ytprim(it), ytprim(it + 1), &
190     yt(it), yt(it + 1), a0, a1, a2, a3)
191 guez 3
192 guez 97 yf1 = yf(it)
193     yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi
194 guez 3
195 guez 97 iter = 1
196     DO
197     yi = yi - (yf1-yfi)/yprimin
198     IF (abs(yi-yo1)<=epsilon .or. iter == 300) exit
199     yo1 = yi
200     yi2 = yi*yi
201     yf1 = a0 + a1*yi + a2*yi2 + a3*yi2*yi
202     yprimin = a1 + 2.*a2*yi + 3.*a3*yi2
203     END DO
204     if (abs(yi-yo1) > epsilon) then
205     print *, 'Pas de solution.', j, ylon2
206     STOP 1
207     end if
208 guez 3
209 guez 97 yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi
210     yprim(j) = pi/(jjm*yprimin)
211     yvrai(j) = yi
212     END DO
213 guez 3
214 guez 97 DO j = 1, jlat - 1
215     IF (yvrai(j + 1)<yvrai(j)) THEN
216     print *, 'Problème avec rlat(', j + 1, ') plus petit que rlat(', &
217     j, ')'
218     STOP 1
219     END IF
220     END DO
221 guez 3
222 guez 97 print *, 'Reorganisation des latitudes pour avoir entre - pi/2 et pi/2'
223 guez 3
224 guez 97 IF (ik==1) THEN
225     ypn = pis2
226 guez 118 DO j = jjm + 1, 1, -1
227 guez 97 IF (yvrai(j)<=ypn) exit
228     END DO
229 guez 3
230 guez 97 jpn = j
231     y00 = yvrai(jpn)
232     deply = pis2 - y00
233     END IF
234 guez 3
235 guez 97 DO j = 1, jjm + 1 - jpn
236     ylatt(j) = -pis2 - y00 + yvrai(jpn + j-1)
237     yprimm(j) = yprim(jpn + j-1)
238     END DO
239 guez 3
240 guez 97 jjpn = jpn
241     IF (jlat==jjm) jjpn = jpn - 1
242 guez 3
243 guez 97 DO j = 1, jjpn
244     ylatt(j + jjm + 1-jpn) = yvrai(j) + deply
245     yprimm(j + jjm + 1-jpn) = yprim(j)
246     END DO
247 guez 3
248 guez 97 ! Fin de la reorganisation
249 guez 3
250 guez 97 DO j = 1, jlat
251     ylat(j) = ylatt(jlat + 1-j)
252     yprim(j) = yprimm(jlat + 1-j)
253     END DO
254 guez 81
255 guez 97 DO j = 1, jlat
256     yvrai(j) = ylat(j)*180./pi
257     END DO
258 guez 81
259 guez 97 IF (ik==1) THEN
260 guez 118 DO j = 1, jjm + 1
261 guez 119 rlatu(j) = ylat(j)
262 guez 97 END DO
263     ELSE IF (ik==2) THEN
264 guez 118 DO j = 1, jjm
265 guez 119 rlatv(j) = ylat(j)
266 guez 97 END DO
267     ELSE IF (ik==3) THEN
268 guez 118 DO j = 1, jjm
269 guez 97 rlatu2(j) = ylat(j)
270     yprimu2(j) = yprim(j)
271     END DO
272     ELSE IF (ik==4) THEN
273 guez 118 DO j = 1, jjm
274 guez 97 rlatu1(j) = ylat(j)
275     yprimu1(j) = yprim(j)
276     END DO
277     END IF
278     END DO loop_ik
279 guez 81
280 guez 97 DO j = 1, jjm
281 guez 119 ylat(j) = rlatu(j) - rlatu(j + 1)
282 guez 97 END DO
283     champmin = 1e12
284     champmax = -1e12
285     DO j = 1, jjm
286     champmin = min(champmin, ylat(j))
287     champmax = max(champmax, ylat(j))
288     END DO
289     champmin = champmin*180./pi
290     champmax = champmax*180./pi
291 guez 81
292 guez 119 DO j = 1, jjm
293     IF (rlatu1(j) <= rlatu2(j)) THEN
294     print *, 'Attention ! rlatu1 < rlatu2 ', rlatu1(j), rlatu2(j), j
295     STOP 13
296     ENDIF
297    
298     IF (rlatu2(j) <= rlatu(j+1)) THEN
299     print *, 'Attention ! rlatu2 < rlatup1 ', rlatu2(j), rlatu(j+1), j
300     STOP 14
301     ENDIF
302    
303     IF (rlatu(j) <= rlatu1(j)) THEN
304     print *, ' Attention ! rlatu < rlatu1 ', rlatu(j), rlatu1(j), j
305     STOP 15
306     ENDIF
307    
308     IF (rlatv(j) <= rlatu2(j)) THEN
309     print *, ' Attention ! rlatv < rlatu2 ', rlatv(j), rlatu2(j), j
310     STOP 16
311     ENDIF
312    
313     IF (rlatv(j) >= rlatu1(j)) THEN
314     print *, ' Attention ! rlatv > rlatu1 ', rlatv(j), rlatu1(j), j
315     STOP 17
316     ENDIF
317    
318     IF (rlatv(j) >= rlatu(j)) THEN
319     print *, ' Attention ! rlatv > rlatu ', rlatv(j), rlatu(j), j
320     STOP 18
321     ENDIF
322     ENDDO
323    
324     print *, 'Latitudes'
325     print 3, champmin, champmax
326    
327     3 Format(1x, ' Au centre du zoom, la longueur de la maille est', &
328     ' d environ ', f0.2, ' degres ', /, &
329     ' alors que la maille en dehors de la zone du zoom est ', &
330 guez 120 "d'environ ", f0.2, ' degres ')
331 guez 119
332 guez 97 END SUBROUTINE fyhyp
333 guez 81
334 guez 97 end module fyhyp_m
  ViewVC Help
Powered by ViewVC 1.1.21