source: LMDZ5/trunk/libf/dyn3d_common/fyhyp_m.F90 @ 2525

Last change on this file since 2525 was 2228, checked in by lguez, 10 years ago

Correcting a problem from revision 2218. The type double precision
with option "-fdefault-real-8" of gfortran is promoted to 16-byte
precision and there is no specific procedure in arth with this
precision. Could not add a specific procedure in arth with double
precision because, with ifort, the option "-real-size 64" does not
promote the double precision, so that would make two identical
specific procedures in arth.

In module nrtype, replaced double precision by a parameterized real
kind so that the effective precision does not depend on a compiler
option.

In coefpoly, fxhyp, fyhyp and invert_zoom_x, use the parameterized
real kind defined in nrtype, instead of double precision.

Also, in module nrtype, removed unused derived types sprs2_sp and
sprs2_dp.

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