source: trunk/LMDZ.COMMON/libf/dyn3d_common/fyhyp_m.F90 @ 1549

Last change on this file since 1549 was 1443, checked in by emillour, 10 years ago

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