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(:), rlonv(:), xprimv(:) ! (iim + 1)
    real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)

    ! Local:
    real rlonm025(iim + 1), rlonp025(iim + 1), d_rlonv(iim)
    REAL delta, step
    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
       step = twopi / iim

       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
       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(xf, xtild, G, rlonm025(:iim), xprimm025(:iim), &
            xuv = - 0.25d0)
       call invert_zoom_x(xf, xtild, G, rlonv(:iim), xprimv(:iim), xuv = 0d0)
       call invert_zoom_x(xf, xtild, G, rlonu(:iim), xprimu(:iim), xuv = 0.5d0)
       call invert_zoom_x(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	119	REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
28			real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
29	guez	3
30	guez	91	! Local:
31	guez	146	real rlonm025(iim + 1), rlonp025(iim + 1), d_rlonv(iim)
32	guez	148	REAL delta, step
33	guez	147	DOUBLE PRECISION, dimension(0:nmax):: xtild, fhyp, G, Xf, ffdx
34			DOUBLE PRECISION beta
35	guez	124	INTEGER i, is2
36	guez	147	DOUBLE PRECISION xmoy(nmax), fxm(nmax)
37	guez	3
38	guez	91	!----------------------------------------------------------------------
39
40	guez	120	print *, "Call sequence information: fxhyp"
41
42	guez	148	if (grossismx == 1.) then
43	guez	126	step = twopi / iim
44	guez	78
45	guez	126	xprimm025(:iim) = step
46			xprimp025(:iim) = step
47			xprimv(:iim) = step
48			xprimu(:iim) = step
49
50	guez	127	rlonv(:iim) = arth(- pi + clon, step, iim)
51	guez	126	rlonm025(:iim) = rlonv(:iim) - 0.25 * step
52			rlonp025(:iim) = rlonv(:iim) + 0.25 * step
53			rlonu(:iim) = rlonv(:iim) + 0.5 * step
54	guez	148	else
55			delta = dzoomx * twopi_d
56	guez	146	xtild = arth(0d0, pi_d / nmax, nmax + 1)
57			forall (i = 1:nmax) xmoy(i) = 0.5d0 * (xtild(i-1) + xtild(i))
58	guez	126
59			! Compute fhyp:
60	guez	148	fhyp(1:nmax - 1) = tanh_cautious(taux * (delta / 2d0 &
61	guez	146	- xtild(1:nmax - 1)), xtild(1:nmax - 1) &
62			* (pi_d - xtild(1:nmax - 1)))
63			fhyp(0) = 1d0
64			fhyp(nmax) = -1d0
65	guez	126
66	guez	148	fxm = tanh_cautious(taux * (delta / 2d0 - xmoy), xmoy * (pi_d - xmoy))
67	guez	3
68	guez	147	! Compute \int_0 ^{\tilde x} F:
69	guez	3
70	guez	147	ffdx(0) = 0d0
71	guez	3
72	guez	146	DO i = 1, nmax
73	guez	147	ffdx(i) = ffdx(i - 1) + fxm(i) * (xtild(i) - xtild(i-1))
74	guez	126	END DO
75	guez	3
76	guez	147	print *, "ffdx(nmax) = ", ffdx(nmax)
77			beta = (pi_d - grossismx * ffdx(nmax)) / (pi_d - ffdx(nmax))
78	guez	126	print *, "beta = ", beta
79	guez	3
80	guez	147	IF (2d0 * beta - grossismx <= 0d0) THEN
81	guez	126	print *, 'Bad choice of grossismx, taux, dzoomx.'
82			print *, 'Decrease dzoomx or grossismx.'
83			STOP 1
84			END IF
85	guez	78
86	guez	146	G = beta + (grossismx - beta) * fhyp
87	guez	78
88	guez	147	Xf(:nmax - 1) = beta * xtild(:nmax - 1) + (grossismx - beta) &
89			* ffdx(:nmax - 1)
90	guez	146	Xf(nmax) = pi_d
91	guez	3
92	guez	146	call invert_zoom_x(xf, xtild, G, rlonm025(:iim), xprimm025(:iim), &
93	guez	127	xuv = - 0.25d0)
94	guez	146	call invert_zoom_x(xf, xtild, G, rlonv(:iim), xprimv(:iim), xuv = 0d0)
95			call invert_zoom_x(xf, xtild, G, rlonu(:iim), xprimu(:iim), xuv = 0.5d0)
96			call invert_zoom_x(xf, xtild, G, rlonp025(:iim), xprimp025(:iim), &
97	guez	127	xuv = 0.25d0)
98	guez	148	end if
99	guez	123
100	guez	124	is2 = 0
101	guez	3
102	guez	124	IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
103			.or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
104			IF (clon <= 0.) THEN
105			is2 = 1
106	guez	3
107	guez	124	do while (rlonm025(is2) < - pi .and. is2 < iim)
108			is2 = is2 + 1
109			end do
110
111			if (rlonm025(is2) < - pi) then
112			print *, 'Rlonm025 plus petit que - pi !'
113			STOP 1
114			end if
115			ELSE
116			is2 = iim
117
118			do while (rlonm025(is2) > pi .and. is2 > 1)
119			is2 = is2 - 1
120			end do
121
122			if (rlonm025(is2) > pi) then
123			print *, 'Rlonm025 plus grand que pi !'
124			STOP 1
125			end if
126			END IF
127			END IF
128
129			call principal_cshift(is2, rlonm025, xprimm025)
130			call principal_cshift(is2, rlonv, xprimv)
131			call principal_cshift(is2, rlonu, xprimu)
132			call principal_cshift(is2, rlonp025, xprimp025)
133
134			forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
135			print , "Minimum longitude step:", MINval(d_rlonv) 180. / pi, "degrees"
136			print , "Maximum longitude step:", MAXval(d_rlonv) 180. / pi, "degrees"
137
138	guez	128	! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
139	guez	119	DO i = 1, iim + 1
140			IF (rlonp025(i) < rlonv(i)) THEN
141	guez	121	print *, 'rlonp025(', i, ') = ', rlonp025(i)
142			print *, "< rlonv(", i, ") = ", rlonv(i)
143	guez	119	STOP 1
144			END IF
145
146			IF (rlonv(i) < rlonm025(i)) THEN
147	guez	121	print *, 'rlonv(', i, ') = ', rlonv(i)
148			print *, "< rlonm025(", i, ") = ", rlonm025(i)
149	guez	119	STOP 1
150			END IF
151
152			IF (rlonp025(i) > rlonu(i)) THEN
153	guez	120	print *, 'rlonp025(', i, ') = ', rlonp025(i)
154			print *, "> rlonu(", i, ") = ", rlonu(i)
155	guez	119	STOP 1
156			END IF
157			END DO
158
159	guez	78	END SUBROUTINE fxhyp
160
161			end module fxhyp_m