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

Last change on this file since 5069 was 2598, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn serre.h into module serre_mod.F90
EM

  • 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.6 KB
RevLine 
[2218]1module fyhyp_m
[524]2
[2218]3  IMPLICIT NONE
[524]4
[2218]5contains
[524]6
[2218]7  SUBROUTINE fyhyp(rlatu, yyprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)
[524]8
[2218]9    ! From LMDZ4/libf/dyn3d/fyhyp.F, version 1.2, 2005/06/03 09:11:32
[524]10
[2218]11    ! Author: P. Le Van, from analysis by R. Sadourny
[524]12
[2218]13    ! Calcule les latitudes et dérivées dans la grille du GCM pour une
14    ! fonction f(y) à dérivée tangente hyperbolique.
[524]15
[2218]16    ! Il vaut mieux avoir : grossismy * dzoom < pi / 2
[524]17
[2218]18    use coefpoly_m, only: coefpoly
[2228]19    use nrtype, only: k8
[2598]20    use serre_mod, only: clat, grossismy, dzoomy, tauy
[524]21
[2218]22    include "dimensions.h"
23    ! for jjm
[524]24
[2218]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)
[524]28
[2218]29    ! Local:
[524]30
[2228]31    REAL(K8) champmin, champmax
[2218]32    INTEGER, PARAMETER:: nmax=30000, nmax2=2*nmax
33    REAL dzoom ! distance totale de la zone du zoom (en radians)
[2228]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
[2218]47    SAVE y00, deply
[524]48
[2218]49    INTEGER i, j, it, ik, iter, jlat
50    INTEGER jpn, jjpn
51    SAVE jpn
[2228]52    REAL(K8) a0, a1, a2, a3, yi2, heavyy0, heavyy0m
53    REAL(K8) fa(0:nmax2), fb(0:nmax2)
[2218]54    REAL y0min, y0max
[524]55
[2228]56    REAL(K8) heavyside
[524]57
[2218]58    !-------------------------------------------------------------------
[524]59
[2218]60    print *, "Call sequence information: fyhyp"
[524]61
[2218]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
[524]70
[2218]71    DO i = 0, nmax2
72       yt(i) = -pis2 + real(i)*pi/nmax2
73    END DO
[524]74
[2218]75    heavyy0m = heavyside(-y0)
76    heavyy0 = heavyside(y0)
77    y0min = 2.*y0*heavyy0m - pis2
78    y0max = 2.*y0*heavyy0 + pis2
[524]79
[2218]80    fa = 999.999
81    fb = 999.999
[524]82
[2218]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
[524]91
[2218]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
[524]99
[2218]100       IF (yt(i)==y0) fhyp(i) = 1.
101       IF (yt(i)==y0min .OR. yt(i)==y0max) fhyp(i) = -1.
102    END DO
[524]103
[2218]104    ! Calcul de beta
[524]105
[2218]106    ffdy = 0.
[524]107
[2218]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
[524]117
[2218]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
[524]129
[2218]130    beta = (grossismy*ffdy-pi)/(ffdy-pi)
[524]131
[2218]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
[524]138
[2218]139    ! calcul de Ytprim
[524]140
[2218]141    DO i = 0, nmax2
142       ytprim(i) = beta + (grossismy-beta)*fhyp(i)
143    END DO
[524]144
[2218]145    ! Calcul de Yf
[524]146
[2218]147    yf(0) = -pis2
148    DO i = 1, nmax2
149       yypr(i) = beta + (grossismy-beta)*fxm(i)
150    END DO
[524]151
[2218]152    DO i = 1, nmax2
153       yf(i) = yf(i-1) + yypr(i)*(yt(i)-yt(i-1))
154    END DO
[524]155
[2218]156    ! yuv = 0. si calcul des latitudes aux pts. U
157    ! yuv = 0.5 si calcul des latitudes aux pts. V
[524]158
[2218]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
[524]173
[2218]174       yo1 = 0.
175       DO j = 1, jlat
176          yo1 = 0.
177          ylon2 = -pis2 + pisjm*(real(j) + yuv-1.)
178          yfi = ylon2
[524]179
[2218]180          it = nmax2
181          DO while (it >= 1 .and. yfi < yf(it))
182             it = it - 1
183          END DO
[524]184
[2218]185          yi = yt(it)
186          IF (it==nmax2) THEN
187             it = nmax2 - 1
188             yf(it + 1) = pis2
189          END IF
[524]190
[2218]191          ! Interpolation entre yi(it) et yi(it + 1) pour avoir Y(yi)
192          ! et Y'(yi)
[524]193
[2218]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
[524]255       DO j = 1, jlat
[2218]256          ylat(j) = ylatt(jlat + 1-j)
257          yprim(j) = yprimm(jlat + 1-j)
258       END DO
[524]259
[2218]260       DO j = 1, jlat
261          yvrai(j) = ylat(j)*180./pi
262       END DO
[524]263
[2218]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
[524]285
[2218]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
[524]297
[2218]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
[524]303
[2218]304       IF (rlatu2(j) <= rlatu(j+1)) THEN
305          print *, 'Attention ! rlatu2 < rlatup1 ', rlatu2(j), rlatu(j+1), j
306          STOP 14
307       ENDIF
[524]308
[2218]309       IF (rlatu(j) <= rlatu1(j)) THEN
310          print *, ' Attention ! rlatu < rlatu1 ', rlatu(j), rlatu1(j), j
311          STOP 15
312       ENDIF
[524]313
[2218]314       IF (rlatv(j) <= rlatu2(j)) THEN
315          print *, ' Attention ! rlatv < rlatu2 ', rlatv(j), rlatu2(j), j
316          STOP 16
317       ENDIF
[524]318
[2218]319       IF (rlatv(j) >= rlatu1(j)) THEN
320          print *, ' Attention ! rlatv > rlatu1 ', rlatv(j), rlatu1(j), j
321          STOP 17
322       ENDIF
[524]323
[2218]324       IF (rlatv(j) >= rlatu(j)) THEN
325          print *, ' Attention ! rlatv > rlatu ', rlatv(j), rlatu(j), j
326          STOP 18
327       ENDIF
328    ENDDO
[524]329
[2218]330    print *, 'Latitudes'
331    print 3, champmin, champmax
[524]332
[2218]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 ')
[524]337
[2218]338  END SUBROUTINE fyhyp
[524]339
[2218]340end module fyhyp_m
Note: See TracBrowser for help on using the repository browser.