source: LMDZ6/trunk/libf/dyn3d_common/fyhyp_m.f90 @ 5408

Last change on this file since 5408 was 5271, checked in by abarral, 7 weeks ago

Move dimensions.h into a module
Nb: doesn't compile yet

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