trunk/dyn3d/fxhyp.f

module fxhyp_m

  IMPLICIT NONE

contains

  SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)

    ! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32
    ! Author: P. Le Van, from formulas by R. Sadourny

    ! Calcule les longitudes et dérivées dans la grille du GCM pour
    ! une fonction x_f(\tilde x) à dérivée tangente hyperbolique.

    ! Il vaut mieux avoir : grossismx \times delta < pi

    ! Le premier point scalaire pour une grille regulière (grossismx =
    ! 1) avec clon = 0 est à - 180 degrés.

    USE dimens_m, ONLY: iim
    use dynetat0_m, only: clon, grossismx, dzoomx, taux
    use invert_zoom_x_m, only: invert_zoom_x, nmax
    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

    REAL, intent(out):: xprimm025(:) ! (iim + 1)

    REAL, intent(out):: rlonv(:) ! (iim + 1)
    ! longitudes of points of the "scalar" and "v" grid, in rad

    REAL, intent(out):: xprimv(:) ! (iim + 1)
    ! 2 pi / iim * (derivative of the longitudinal zoom function)(rlonv)

    real, intent(out):: rlonu(:) ! (iim + 1)
    ! longitudes of points of the "u" grid, in rad

    real, intent(out):: xprimu(:) ! (iim + 1)
    ! 2 pi / iim * (derivative of the longitudinal zoom function)(rlonu)

    real, intent(out):: xprimp025(:) ! (iim + 1)

    ! Local:
    real rlonm025(iim + 1), rlonp025(iim + 1), d_rlonv(iim)
    REAL delta, h
    DOUBLE PRECISION, dimension(0:nmax):: xtild, fhyp, G, Xf, ffdx
    DOUBLE PRECISION beta
    INTEGER i, is2
    DOUBLE PRECISION xmoy(nmax), fxm(nmax)

    !----------------------------------------------------------------------

    print *, "Call sequence information: fxhyp"

    if (grossismx == 1.) then
       h = twopi / iim

       xprimm025(:iim) = h
       xprimp025(:iim) = h
       xprimv(:iim) = h
       xprimu(:iim) = h

       rlonv(:iim) = arth(- pi + clon, h, iim)
       rlonm025(:iim) = rlonv(:iim) - 0.25 * h
       rlonp025(:iim) = rlonv(:iim) + 0.25 * h
       rlonu(:iim) = rlonv(:iim) + 0.5 * h
    else
       delta = dzoomx * twopi_d
       xtild = arth(0d0, pi_d / nmax, nmax + 1)
       forall (i = 1:nmax) xmoy(i) = 0.5d0 * (xtild(i-1) + xtild(i))

       ! Compute fhyp:
       fhyp(1:nmax - 1) = tanh_cautious(taux * (delta / 2d0 &
            - xtild(1:nmax - 1)), xtild(1:nmax - 1) &
            * (pi_d - xtild(1:nmax - 1)))
       fhyp(0) = 1d0
       fhyp(nmax) = -1d0

       fxm = tanh_cautious(taux * (delta / 2d0 - xmoy), xmoy * (pi_d - xmoy))

       ! Compute \int_0 ^{\tilde x} F:

       ffdx(0) = 0d0

       DO i = 1, nmax
          ffdx(i) = ffdx(i - 1) + fxm(i) * (xtild(i) - xtild(i-1))
       END DO

       print *, "ffdx(nmax) = ", ffdx(nmax)
       beta = (pi_d - grossismx * ffdx(nmax)) / (pi_d - ffdx(nmax))
       print *, "beta = ", beta

       IF (2d0 * beta - grossismx <= 0d0) THEN
          print *, 'Bad choice of grossismx, taux, dzoomx.'
          print *, 'Decrease dzoomx or grossismx.'
          STOP 1
       END IF

       G = beta + (grossismx - beta) * fhyp

       Xf(:nmax - 1) = beta * xtild(:nmax - 1) + (grossismx - beta) &
            * ffdx(:nmax - 1)
       Xf(nmax) = pi_d

       call invert_zoom_x(beta, xf, xtild, G, rlonm025(:iim), xprimm025(:iim), &
            xuv = - 0.25d0)
       call invert_zoom_x(beta, xf, xtild, G, rlonv(:iim), xprimv(:iim), &
            xuv = 0d0)
       call invert_zoom_x(beta, xf, xtild, G, rlonu(:iim), xprimu(:iim), &
            xuv = 0.5d0)
       call invert_zoom_x(beta, xf, xtild, G, rlonp025(:iim), xprimp025(:iim), &
            xuv = 0.25d0)
    end if

    is2 = 0

    IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
         .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
       IF (clon <= 0.) THEN
          is2 = 1

          do while (rlonm025(is2) < - pi .and. is2 < iim)
             is2 = is2 + 1
          end do

          if (rlonm025(is2) < - pi) then
             print *, 'Rlonm025 plus petit que - pi !'
             STOP 1
          end if
       ELSE
          is2 = iim

          do while (rlonm025(is2) > pi .and. is2 > 1)
             is2 = is2 - 1
          end do

          if (rlonm025(is2) > pi) then
             print *, 'Rlonm025 plus grand que pi !'
             STOP 1
          end if
       END IF
    END IF

    call principal_cshift(is2, rlonm025, xprimm025)
    call principal_cshift(is2, rlonv, xprimv)
    call principal_cshift(is2, rlonu, xprimu)
    call principal_cshift(is2, rlonp025, xprimp025)

    forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
    print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
    print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"

    ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
    DO i = 1, iim + 1
       IF (rlonp025(i) < rlonv(i)) THEN
          print *, 'rlonp025(', i, ') = ', rlonp025(i)
          print *, "< rlonv(", i, ") = ", rlonv(i)
          STOP 1
       END IF

       IF (rlonv(i) < rlonm025(i)) THEN 
          print *, 'rlonv(', i, ') = ', rlonv(i)
          print *, "< rlonm025(", i, ") = ", rlonm025(i)
          STOP 1
       END IF

       IF (rlonp025(i) > rlonu(i)) THEN
          print *, 'rlonp025(', i, ') = ', rlonp025(i)
          print *, "> rlonu(", i, ") = ", rlonu(i)
          STOP 1
       END IF
    END DO

  END SUBROUTINE fxhyp

end module fxhyp_m
1	guez	78	module fxhyp_m
2	guez	3
3	guez	78	IMPLICIT NONE
4	guez	3
5	guez	78	contains
6	guez	3
7	guez	119	SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
8	guez	3
9	guez	91	! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32
10	guez	119	! Author: P. Le Van, from formulas by R. Sadourny
11	guez	3
12	guez	78	! Calcule les longitudes et dérivées dans la grille du GCM pour
13	guez	151	! une fonction x_f(\tilde x) à dérivée tangente hyperbolique.
14	guez	3
15	guez	148	! Il vaut mieux avoir : grossismx \times delta < pi
16	guez	3
17	guez	120	! Le premier point scalaire pour une grille regulière (grossismx =
18	guez	151	! 1) avec clon = 0 est à - 180 degrés.
19	guez	120
20	guez	78	USE dimens_m, ONLY: iim
21	guez	139	use dynetat0_m, only: clon, grossismx, dzoomx, taux
22	guez	131	use invert_zoom_x_m, only: invert_zoom_x, nmax
23	guez	126	use nr_util, only: pi, pi_d, twopi, twopi_d, arth
24	guez	124	use principal_cshift_m, only: principal_cshift
25	guez	144	use tanh_cautious_m, only: tanh_cautious
26	guez	3
27	guez	156	REAL, intent(out):: xprimm025(:) ! (iim + 1)
28	guez	3
29	guez	156	REAL, intent(out):: rlonv(:) ! (iim + 1)
30			! longitudes of points of the "scalar" and "v" grid, in rad
31
32			REAL, intent(out):: xprimv(:) ! (iim + 1)
33			! 2 pi / iim * (derivative of the longitudinal zoom function)(rlonv)
34
35			real, intent(out):: rlonu(:) ! (iim + 1)
36			! longitudes of points of the "u" grid, in rad
37
38			real, intent(out):: xprimu(:) ! (iim + 1)
39			! 2 pi / iim * (derivative of the longitudinal zoom function)(rlonu)
40
41			real, intent(out):: xprimp025(:) ! (iim + 1)
42
43	guez	91	! Local:
44	guez	146	real rlonm025(iim + 1), rlonp025(iim + 1), d_rlonv(iim)
45	guez	156	REAL delta, h
46	guez	147	DOUBLE PRECISION, dimension(0:nmax):: xtild, fhyp, G, Xf, ffdx
47			DOUBLE PRECISION beta
48	guez	124	INTEGER i, is2
49	guez	147	DOUBLE PRECISION xmoy(nmax), fxm(nmax)
50	guez	3
51	guez	91	!----------------------------------------------------------------------
52
53	guez	120	print *, "Call sequence information: fxhyp"
54
55	guez	148	if (grossismx == 1.) then
56	guez	156	h = twopi / iim
57	guez	78
58	guez	156	xprimm025(:iim) = h
59			xprimp025(:iim) = h
60			xprimv(:iim) = h
61			xprimu(:iim) = h
62	guez	126
63	guez	156	rlonv(:iim) = arth(- pi + clon, h, iim)
64			rlonm025(:iim) = rlonv(:iim) - 0.25 * h
65			rlonp025(:iim) = rlonv(:iim) + 0.25 * h
66			rlonu(:iim) = rlonv(:iim) + 0.5 * h
67	guez	148	else
68			delta = dzoomx * twopi_d
69	guez	146	xtild = arth(0d0, pi_d / nmax, nmax + 1)
70			forall (i = 1:nmax) xmoy(i) = 0.5d0 * (xtild(i-1) + xtild(i))
71	guez	126
72			! Compute fhyp:
73	guez	148	fhyp(1:nmax - 1) = tanh_cautious(taux * (delta / 2d0 &
74	guez	146	- xtild(1:nmax - 1)), xtild(1:nmax - 1) &
75			* (pi_d - xtild(1:nmax - 1)))
76			fhyp(0) = 1d0
77			fhyp(nmax) = -1d0
78	guez	126
79	guez	148	fxm = tanh_cautious(taux * (delta / 2d0 - xmoy), xmoy * (pi_d - xmoy))
80	guez	3
81	guez	147	! Compute \int_0 ^{\tilde x} F:
82	guez	3
83	guez	147	ffdx(0) = 0d0
84	guez	3
85	guez	146	DO i = 1, nmax
86	guez	147	ffdx(i) = ffdx(i - 1) + fxm(i) * (xtild(i) - xtild(i-1))
87	guez	126	END DO
88	guez	3
89	guez	147	print *, "ffdx(nmax) = ", ffdx(nmax)
90			beta = (pi_d - grossismx * ffdx(nmax)) / (pi_d - ffdx(nmax))
91	guez	126	print *, "beta = ", beta
92	guez	3
93	guez	147	IF (2d0 * beta - grossismx <= 0d0) THEN
94	guez	126	print *, 'Bad choice of grossismx, taux, dzoomx.'
95			print *, 'Decrease dzoomx or grossismx.'
96			STOP 1
97			END IF
98	guez	78
99	guez	146	G = beta + (grossismx - beta) * fhyp
100	guez	78
101	guez	147	Xf(:nmax - 1) = beta * xtild(:nmax - 1) + (grossismx - beta) &
102			* ffdx(:nmax - 1)
103	guez	146	Xf(nmax) = pi_d
104	guez	3
105	guez	167	call invert_zoom_x(beta, xf, xtild, G, rlonm025(:iim), xprimm025(:iim), &
106	guez	127	xuv = - 0.25d0)
107	guez	167	call invert_zoom_x(beta, xf, xtild, G, rlonv(:iim), xprimv(:iim), &
108			xuv = 0d0)
109			call invert_zoom_x(beta, xf, xtild, G, rlonu(:iim), xprimu(:iim), &
110			xuv = 0.5d0)
111			call invert_zoom_x(beta, xf, xtild, G, rlonp025(:iim), xprimp025(:iim), &
112	guez	127	xuv = 0.25d0)
113	guez	148	end if
114	guez	123
115	guez	124	is2 = 0
116	guez	3
117	guez	124	IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
118			.or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
119			IF (clon <= 0.) THEN
120			is2 = 1
121	guez	3
122	guez	124	do while (rlonm025(is2) < - pi .and. is2 < iim)
123			is2 = is2 + 1
124			end do
125
126			if (rlonm025(is2) < - pi) then
127			print *, 'Rlonm025 plus petit que - pi !'
128			STOP 1
129			end if
130			ELSE
131			is2 = iim
132
133			do while (rlonm025(is2) > pi .and. is2 > 1)
134			is2 = is2 - 1
135			end do
136
137			if (rlonm025(is2) > pi) then
138			print *, 'Rlonm025 plus grand que pi !'
139			STOP 1
140			end if
141			END IF
142			END IF
143
144			call principal_cshift(is2, rlonm025, xprimm025)
145			call principal_cshift(is2, rlonv, xprimv)
146			call principal_cshift(is2, rlonu, xprimu)
147			call principal_cshift(is2, rlonp025, xprimp025)
148
149			forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
150			print , "Minimum longitude step:", MINval(d_rlonv) 180. / pi, "degrees"
151			print , "Maximum longitude step:", MAXval(d_rlonv) 180. / pi, "degrees"
152
153	guez	128	! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
154	guez	119	DO i = 1, iim + 1
155			IF (rlonp025(i) < rlonv(i)) THEN
156	guez	121	print *, 'rlonp025(', i, ') = ', rlonp025(i)
157			print *, "< rlonv(", i, ") = ", rlonv(i)
158	guez	119	STOP 1
159			END IF
160
161			IF (rlonv(i) < rlonm025(i)) THEN
162	guez	121	print *, 'rlonv(', i, ') = ', rlonv(i)
163			print *, "< rlonm025(", i, ") = ", rlonm025(i)
164	guez	119	STOP 1
165			END IF
166
167			IF (rlonp025(i) > rlonu(i)) THEN
168	guez	120	print *, 'rlonp025(', i, ') = ', rlonp025(i)
169			print *, "> rlonu(", i, ") = ", rlonu(i)
170	guez	119	STOP 1
171			END IF
172			END DO
173
174	guez	78	END SUBROUTINE fxhyp
175
176			end module fxhyp_m