/[lmdze]/trunk/Sources/dyn3d/fxhyp.f

Diff of /trunk/Sources/dyn3d/fxhyp.f

Parent Directory | Revision Log | View Patch Patch

-trunk/dyn3d/fxhyp.f
revision 105 by guez,
Thu Sep  4 10:40:24 2014 UTC
+trunk/Sources/dyn3d/fxhyp.f
revision 144 by guez,
Wed Jun 10 16:46:46 2015 UTC
 Line 4 
 module fxhyp_m
  contains
-   SUBROUTINE fxhyp(xzoomdeg, grossism, dzooma, tau, rlonm025, xprimm025, &
+   SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
-        rlonv, xprimv, rlonu, xprimu, rlonp025, xprimp025, champmin, champmax)
      ! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32
-     ! Author: P. Le Van
+     ! Author: P. Le Van, from formulas by R. Sadourny
      ! Calcule les longitudes et dérivées dans la grille du GCM pour
-     ! une fonction f(x) à tangente hyperbolique.
+     ! une fonction f(x) à dérivée tangente hyperbolique.
-     ! On doit avoir grossism \times dzoom < pi (radians), en longitude.
+     ! Il vaut mieux avoir : grossismx \times dzoom < pi
-     USE dimens_m, ONLY: iim
+     ! Le premier point scalaire pour une grille regulière (grossismx =
-     USE paramet_m, ONLY: iip1
+     ! 1., taux=0., clon=0.) est à - 180 degrés.
-     REAL, intent(in):: xzoomdeg
-     REAL, intent(in):: grossism
-     ! grossissement (= 2 si 2 fois, = 3 si 3 fois, etc.)
-     REAL, intent(in):: dzooma ! distance totale de la zone du zoom
-     REAL, intent(in):: tau
-     ! raideur de la transition de l'intérieur à l'extérieur du zoom
-     ! arguments de sortie
+     USE dimens_m, ONLY: iim
+     use dynetat0_m, only: clon, grossismx, dzoomx, taux
-     REAL, dimension(iip1):: rlonm025, xprimm025, rlonv, xprimv
+     use invert_zoom_x_m, only: invert_zoom_x, nmax
-     real, dimension(iip1):: rlonu, xprimu, rlonp025, xprimp025
+     use nr_util, only: pi, pi_d, twopi, twopi_d, arth
+     use principal_cshift_m, only: principal_cshift
+     use tanh_cautious_m, only: tanh_cautious
-     DOUBLE PRECISION, intent(out):: champmin, champmax
+     REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
+     real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
      ! Local:
+     real rlonm025(iim + 1), rlonp025(iim + 1)
-     INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2*nmax
+     REAL dzoom, step
+     real d_rlonv(iim)
-     LOGICAL, PARAMETER:: scal180 = .TRUE.
+     DOUBLE PRECISION xtild(0:2 * nmax)
-     ! scal180 = .TRUE. si on veut avoir le premier point scalaire pour
+     DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
-     ! une grille reguliere (grossism = 1., tau=0., clon=0.) a
+     DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)
-     ! -180. degres. sinon scal180 = .FALSE.
+     INTEGER i, is2
+     DOUBLE PRECISION, dimension(nmax + 1:2 * nmax):: xmoy, fxm
-     REAL dzoom
-     DOUBLE PRECISION xlon(iip1), xprimm(iip1), xuv
-     DOUBLE PRECISION xtild(0:nmax2)
-     DOUBLE PRECISION fhyp(0:nmax2), ffdx, beta, Xprimt(0:nmax2)
-     DOUBLE PRECISION Xf(0:nmax2), xxpr(0:nmax2)
-     DOUBLE PRECISION xvrai(iip1), xxprim(iip1)
-     DOUBLE PRECISION pi, depi, epsilon, xzoom, fa, fb
-     DOUBLE PRECISION Xf1, Xfi, a0, a1, a2, a3, xi2
-     INTEGER i, it, ik, iter, ii, idif, ii1, ii2
-     DOUBLE PRECISION xi, xo1, xmoy, xlon2, fxm, Xprimin
-     DOUBLE PRECISION decalx
-     INTEGER is2
-     SAVE is2
      !----------------------------------------------------------------------
-     pi = 2. * ASIN(1.)
+     print *, "Call sequence information: fxhyp"
-     depi = 2. * pi
-     epsilon = 1.e-3
-     xzoom = xzoomdeg * pi/180.
-     decalx = .75
-     IF (grossism == 1. .AND. scal180) THEN
-        decalx = 1.
-     ENDIF
-     print *, 'FXHYP scal180, decalx', scal180, decalx
-     IF (dzooma.LT.1.) THEN
-        dzoom = dzooma * depi
-     ELSEIF (dzooma.LT. 25.) THEN
-        print *, "Le paramètre dzoomx pour fxhyp est trop petit. " &
-             // "L'augmenter et relancer."
-        STOP 1
-     ELSE
-        dzoom = dzooma * pi/180.
-     ENDIF
-     print *, ' xzoom(rad), grossism, tau, dzoom (rad):'
-     print *, xzoom, grossism, tau, dzoom
-     DO i = 0, nmax2
-        xtild(i) = - pi + FLOAT(i) * depi /nmax2
-     ENDDO
-     DO i = nmax, nmax2
-        fa = tau* (dzoom/2. - xtild(i))
-        fb = xtild(i) * (pi - xtild(i))
-        IF (200.* fb .LT. - fa) THEN
-           fhyp (i) = - 1.
-        ELSEIF (200. * fb .LT. fa) THEN
-           fhyp (i) = 1.
-        ELSE
-           IF (ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN
-              IF (200.*fb + fa.LT.1.e-10) THEN
-                 fhyp (i) = - 1.
-              ELSEIF (200.*fb - fa.LT.1.e-10) THEN
-                 fhyp (i) = 1.
-              ENDIF
-           ELSE
-              fhyp (i) = TANH (fa/fb)
-           ENDIF
-        ENDIF
-        IF (xtild(i) == 0.) fhyp(i) = 1.
-        IF (xtild(i) == pi) fhyp(i) = -1.
-     ENDDO
-     ! Calcul de beta
-     ffdx = 0.
-     DO i = nmax + 1, nmax2
-        xmoy = 0.5 * (xtild(i-1) + xtild(i))
-        fa = tau* (dzoom/2. - xmoy)
-        fb = xmoy * (pi - xmoy)
-        IF (200.* fb .LT. - fa) THEN
-           fxm = - 1.
-        ELSEIF (200. * fb .LT. fa) THEN
-           fxm = 1.
-        ELSE
-           IF (ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN
-              IF (200.*fb + fa.LT.1.e-10) THEN
-                 fxm = - 1.
-              ELSEIF (200.*fb - fa.LT.1.e-10) THEN
-                 fxm = 1.
-              ENDIF
-           ELSE
-              fxm = TANH (fa/fb)
-           ENDIF
-        ENDIF
-        IF (xmoy == 0.) fxm = 1.
-        IF (xmoy == pi) fxm = -1.
-        ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
-     ENDDO
-     beta = (grossism * ffdx - pi) / (ffdx - pi)
-     IF (2.*beta - grossism <= 0.) THEN
-        print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'
-        print *, 'Modifier les valeurs de grossismx, tau ou dzoomx et relancer.'
-        STOP 1
-     ENDIF
-     ! calcul de Xprimt
-     DO i = nmax, nmax2
-        Xprimt(i) = beta + (grossism - beta) * fhyp(i)
-     ENDDO
-     DO i = nmax + 1, nmax2
-        Xprimt(nmax2 - i) = Xprimt(i)
-     ENDDO
-     ! Calcul de Xf
-     Xf(0) = - pi
-     DO i = nmax + 1, nmax2
-        xmoy = 0.5 * (xtild(i-1) + xtild(i))
-        fa = tau* (dzoom/2. - xmoy)
-        fb = xmoy * (pi - xmoy)
-        IF (200.* fb .LT. - fa) THEN
-           fxm = - 1.
-        ELSEIF (200. * fb .LT. fa) THEN
-           fxm = 1.
-        ELSE
-           fxm = TANH (fa/fb)
-        ENDIF
-        IF (xmoy == 0.) fxm = 1.
+     test_grossismx: if (grossismx == 1.) then
-        IF (xmoy == pi) fxm = -1.
+        step = twopi / iim
-        xxpr(i) = beta + (grossism - beta) * fxm
-     ENDDO
-     DO i = nmax + 1, nmax2
-        xxpr(nmax2-i + 1) = xxpr(i)
-     ENDDO
-     DO i=1, nmax2
-        Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
-     ENDDO
-     ! xuv = 0. si calcul aux pts scalaires
-     ! xuv = 0.5 si calcul aux pts U
-     print *
-     DO ik = 1, 4
-        IF (ik == 1) THEN
-           xuv = -0.25
-        ELSE IF (ik == 2) THEN
-           xuv = 0.
-        ELSE IF (ik == 3) THEN
-           xuv = 0.50
-        ELSE IF (ik == 4) THEN
-           xuv = 0.25
-        ENDIF
-        xo1 = 0.
-        ii1=1
-        ii2=iim
-        IF (ik == 1.and.grossism == 1.) THEN
-           ii1 = 2
-           ii2 = iim + 1
-        ENDIF
-        DO i = ii1, ii2
-           xlon2 = - pi + (FLOAT(i) + xuv - decalx) * depi / FLOAT(iim)
-           Xfi = xlon2
-           it = nmax2
-           do while (xfi < xf(it) .and. it >= 1)
-              it = it - 1
-           end do
-           ! Calcul de Xf(xi)
+        xprimm025(:iim) = step
+        xprimp025(:iim) = step
+        xprimv(:iim) = step
+        xprimu(:iim) = step
+        rlonv(:iim) = arth(- pi + clon, step, iim)
+        rlonm025(:iim) = rlonv(:iim) - 0.25 * step
+        rlonp025(:iim) = rlonv(:iim) + 0.25 * step
+        rlonu(:iim) = rlonv(:iim) + 0.5 * step
+     else test_grossismx
+        dzoom = dzoomx * twopi_d
+        xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
+        forall (i = nmax + 1:2 * nmax) xmoy(i) = 0.5d0 * (xtild(i-1) + xtild(i))
+        ! Compute fhyp:
+        fhyp(nmax + 1:2 * nmax - 1) = tanh_cautious(taux * (dzoom / 2. &
+             - xtild(nmax + 1:2 * nmax - 1)), xtild(nmax + 1:2 * nmax - 1) &
+             * (pi_d - xtild(nmax + 1:2 * nmax - 1)))
+        fhyp(nmax) = 1d0
+        fhyp(2 * nmax) = -1d0
+        fxm = tanh_cautious(taux * (dzoom / 2. - xmoy), xmoy * (pi_d - xmoy))
+        ! Calcul de beta
+        ffdx = 0.
+        DO i = nmax + 1, 2 * nmax
+           ffdx = ffdx + fxm(i) * (xtild(i) - xtild(i-1))
+        END DO
+        print *, "ffdx = ", ffdx
+        beta = (pi_d - grossismx * ffdx) / (pi_d - ffdx)
+        print *, "beta = ", beta
+        IF (2. * beta - grossismx <= 0.) THEN
+           print *, 'Bad choice of grossismx, taux, dzoomx.'
+           print *, 'Decrease dzoomx or grossismx.'
+           STOP 1
+        END IF
+        ! calcul de Xprimt
+        Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
+        xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
+        ! Calcul de Xf
+        xxpr(nmax + 1:2 * nmax) = beta + (grossismx - beta) * fxm
+        xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
+        Xf(0) = - pi_d
+        DO i=1, 2 * nmax - 1
+           Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
+        END DO
+        Xf(2 * nmax) = pi_d
+        call invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), xprimm025(:iim), &
+             xuv = - 0.25d0)
+        call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
+             xuv = 0d0)
+        call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
+             xuv = 0.5d0)
+        call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), xprimp025(:iim), &
+             xuv = 0.25d0)
+     end if test_grossismx
+     is2 = 0
+     IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
+          .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
+        IF (clon <= 0.) THEN
+           is2 = 1
-           xi = xtild(it)
+           do while (rlonm025(is2) < - pi .and. is2 < iim)
+              is2 = is2 + 1
+           end do
-           IF (it == nmax2) THEN
+           if (rlonm025(is2) < - pi) then
-              it = nmax2 -1
+              print *, 'Rlonm025 plus petit que - pi !'
-              Xf(it + 1) = pi
-           ENDIF
-           ! Appel de la routine qui calcule les coefficients a0, a1,
-           ! a2, a3 d'un polynome de degre 3 qui passe par les points
-           ! (Xf(it), xtild(it)) et (Xf(it + 1), xtild(it + 1))
-           CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &
-                xtild(it), xtild(it + 1), a0, a1, a2, a3)
-           Xf1 = Xf(it)
-           Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi
-           iter = 1
-           do
-              xi = xi - (Xf1 - Xfi)/ Xprimin
-              IF (ABS(xi - xo1) <= epsilon .or. iter == 300) exit
-              xo1 = xi
-              xi2 = xi * xi
-              Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi
-              Xprimin = a1 + 2.* a2 * xi + 3.* a3 * xi2
-           end DO
-           if (ABS(xi - xo1) > epsilon) then
-              ! iter == 300
-              print *, 'Pas de solution.'
-              print *, i, xlon2
               STOP 1
            end if
+        ELSE
+           is2 = iim
+           do while (rlonm025(is2) > pi .and. is2 > 1)
+              is2 = is2 - 1
+           end do
-           xxprim(i) = depi/ (FLOAT(iim) * Xprimin)
+           if (rlonm025(is2) > pi) then
-           xvrai(i) = xi + xzoom
+              print *, 'Rlonm025 plus grand que pi !'
-        end DO
-        IF (ik == 1.and.grossism == 1.) THEN
-           xvrai(1) = xvrai(iip1)-depi
-           xxprim(1) = xxprim(iip1)
-        ENDIF
-        DO i = 1, iim
-           xlon(i) = xvrai(i)
-           xprimm(i) = xxprim(i)
-        ENDDO
-        DO i = 1, iim -1
-           IF (xvrai(i + 1).LT. xvrai(i)) THEN
-              print *, 'Problème avec rlonu(', i + 1, &
-                   ') plus petit que rlonu(', i, ')'
               STOP 1
-           ENDIF
+           end if
-        ENDDO
+        END IF
+     END IF
-        ! Reorganisation des longitudes pour les avoir entre - pi et pi
-        champmin = 1.e12
+     call principal_cshift(is2, rlonm025, xprimm025)
-        champmax = -1.e12
+     call principal_cshift(is2, rlonv, xprimv)
-        DO i = 1, iim
+     call principal_cshift(is2, rlonu, xprimu)
-           champmin = MIN(champmin, xvrai(i))
+     call principal_cshift(is2, rlonp025, xprimp025)
-           champmax = MAX(champmax, xvrai(i))
-        ENDDO
+     forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
+     print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
-        IF (.not. (champmin >= -pi-0.10.and.champmax <= pi + 0.10)) THEN
+     print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"
-           print *, 'Reorganisation des longitudes pour avoir entre - pi', &
-                ' et pi '
+     ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
+     DO i = 1, iim + 1
-           IF (xzoom <= 0.) THEN
+        IF (rlonp025(i) < rlonv(i)) THEN
-              IF (ik == 1) THEN
+           print *, 'rlonp025(', i, ') = ', rlonp025(i)
-                 i = 1
+           print *, "< rlonv(", i, ") = ", rlonv(i)
+           STOP 1
-                 do while (xvrai(i) < - pi .and. i < iim)
+        END IF
-                    i = i + 1
-                 end do
+        IF (rlonv(i) < rlonm025(i)) THEN
+           print *, 'rlonv(', i, ') = ', rlonv(i)
-                 if (xvrai(i) < - pi) then
+           print *, "< rlonm025(", i, ") = ", rlonm025(i)
-                    print *, ' PBS. 1 ! Xvrai plus petit que - pi ! '
+           STOP 1
-                    STOP 1
+        END IF
-                 end if
+        IF (rlonp025(i) > rlonu(i)) THEN
-                 is2 = i
+           print *, 'rlonp025(', i, ') = ', rlonp025(i)
-              ENDIF
+           print *, "> rlonu(", i, ") = ", rlonu(i)
+           STOP 1
-              IF (is2.NE. 1) THEN
+        END IF
-                 DO ii = is2, iim
+     END DO
-                    xlon (ii-is2 + 1) = xvrai(ii)
-                    xprimm(ii-is2 + 1) = xxprim(ii)
-                 ENDDO
-                 DO ii = 1, is2 -1
-                    xlon (ii + iim-is2 + 1) = xvrai(ii) + depi
-                    xprimm(ii + iim-is2 + 1) = xxprim(ii)
-                 ENDDO
-              ENDIF
-           ELSE
-              IF (ik == 1) THEN
-                 i = iim
-                 do while (xvrai(i) > pi .and. i > 1)
-                    i = i - 1
-                 end do
-                 if (xvrai(i) > pi) then
-                    print *, ' PBS. 2 ! Xvrai plus grand que pi ! '
-                    STOP 1
-                 end if
-                 is2 = i
-              ENDIF
-              idif = iim -is2
-              DO ii = 1, is2
-                 xlon (ii + idif) = xvrai(ii)
-                 xprimm(ii + idif) = xxprim(ii)
-              ENDDO
-              DO ii = 1, idif
-                 xlon (ii) = xvrai (ii + is2) - depi
-                 xprimm(ii) = xxprim(ii + is2)
-              ENDDO
-           ENDIF
-        ENDIF
-        ! Fin de la reorganisation
-        xlon (iip1) = xlon(1) + depi
-        xprimm(iip1) = xprimm (1)
-        DO i = 1, iim + 1
-           xvrai(i) = xlon(i)*180./pi
-        ENDDO
-        IF (ik == 1) THEN
-           DO i = 1, iim + 1
-              rlonm025(i) = xlon(i)
-              xprimm025(i) = xprimm(i)
-           ENDDO
-        ELSE IF (ik == 2) THEN
-           DO i = 1, iim + 1
-              rlonv(i) = xlon(i)
-              xprimv(i) = xprimm(i)
-           ENDDO
-        ELSE IF (ik == 3) THEN
-           DO i = 1, iim + 1
-              rlonu(i) = xlon(i)
-              xprimu(i) = xprimm(i)
-           ENDDO
-        ELSE IF (ik == 4) THEN
-           DO i = 1, iim + 1
-              rlonp025(i) = xlon(i)
-              xprimp025(i) = xprimm(i)
-           ENDDO
-        ENDIF
-     end DO
-     print *
-     DO i = 1, iim
-        xlon(i) = rlonv(i + 1) - rlonv(i)
-     ENDDO
-     champmin = 1.e12
-     champmax = -1.e12
-     DO i = 1, iim
-        champmin = MIN(champmin, xlon(i))
-        champmax = MAX(champmax, xlon(i))
-     ENDDO
-     champmin = champmin * 180./pi
-     champmax = champmax * 180./pi
    END SUBROUTINE fxhyp

 Legend:



Removed from v.105
 


changed lines


 
Added in v.144
 Legend:



Removed from v.105
 


changed lines


 
Added in v.144
-Removed from v.105
+Added in v.144

	ViewVC Help
Powered by ViewVC 1.1.21