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 f(x) à dérivée tangente hyperbolique.

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

    ! Le premier point scalaire pour une grille regulière (grossismx =
    ! 1., taux = 0., 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 dzoom, step
    DOUBLE PRECISION, dimension(0:nmax):: xtild, fhyp, G, Xf
    DOUBLE PRECISION ffdx, beta
    INTEGER i, is2
    DOUBLE PRECISION xxpr(nmax - 1), xmoy(nmax), fxm(nmax)

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

    print *, "Call sequence information: fxhyp"

    test_grossismx: 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 test_grossismx
       dzoom = 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 * (dzoom / 2. &
            - xtild(1:nmax - 1)), xtild(1:nmax - 1) &
            * (pi_d - xtild(1:nmax - 1)))
       fhyp(0) = 1d0
       fhyp(nmax) = -1d0

       fxm = tanh_cautious(taux * (dzoom / 2. - xmoy), xmoy * (pi_d - xmoy))

       ! Calcul de beta 

       ffdx = 0.

       DO i = 1, 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

       G = beta + (grossismx - beta) * fhyp

       ! Calcul de Xf

       xxpr = beta + (grossismx - beta) * fxm(:nmax - 1)
       Xf(0) = 0d0

       DO i = 1, nmax - 1
          Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
       END DO

       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 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

          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	119	! une fonction f(x) à dérivée tangente hyperbolique.
14	guez	3
15	guez	121	! Il vaut mieux avoir : grossismx \times dzoom < pi
16	guez	3
17	guez	120	! Le premier point scalaire pour une grille regulière (grossismx =
18	guez	145	! 1., taux = 0., 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	126	REAL dzoom, step
33	guez	146	DOUBLE PRECISION, dimension(0:nmax):: xtild, fhyp, G, Xf
34			DOUBLE PRECISION ffdx, beta
35	guez	124	INTEGER i, is2
36	guez	146	DOUBLE PRECISION xxpr(nmax - 1), xmoy(nmax), fxm(nmax)
37	guez	3
38	guez	91	!----------------------------------------------------------------------
39
40	guez	120	print *, "Call sequence information: fxhyp"
41
42	guez	126	test_grossismx: if (grossismx == 1.) then
43			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	132	else test_grossismx
55	guez	126	dzoom = 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	146	fhyp(1:nmax - 1) = tanh_cautious(taux * (dzoom / 2. &
61			- 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	144	fxm = tanh_cautious(taux * (dzoom / 2. - xmoy), xmoy * (pi_d - xmoy))
67	guez	3
68	guez	126	! Calcul de beta
69	guez	3
70	guez	126	ffdx = 0.
71	guez	3
72	guez	146	DO i = 1, nmax
73	guez	144	ffdx = ffdx + fxm(i) * (xtild(i) - xtild(i-1))
74	guez	126	END DO
75	guez	3
76	guez	126	print *, "ffdx = ", ffdx
77	guez	144	beta = (pi_d - grossismx * ffdx) / (pi_d - ffdx)
78	guez	126	print *, "beta = ", beta
79	guez	3
80	guez	126	IF (2. * beta - grossismx <= 0.) THEN
81			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	126	! Calcul de Xf
89	guez	78
90	guez	146	xxpr = beta + (grossismx - beta) * fxm(:nmax - 1)
91			Xf(0) = 0d0
92	guez	3
93	guez	146	DO i = 1, nmax - 1
94	guez	126	Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
95			END DO
96	guez	3
97	guez	146	Xf(nmax) = pi_d
98	guez	3
99	guez	146	call invert_zoom_x(xf, xtild, G, rlonm025(:iim), xprimm025(:iim), &
100	guez	127	xuv = - 0.25d0)
101	guez	146	call invert_zoom_x(xf, xtild, G, rlonv(:iim), xprimv(:iim), xuv = 0d0)
102			call invert_zoom_x(xf, xtild, G, rlonu(:iim), xprimu(:iim), xuv = 0.5d0)
103			call invert_zoom_x(xf, xtild, G, rlonp025(:iim), xprimp025(:iim), &
104	guez	127	xuv = 0.25d0)
105	guez	126	end if test_grossismx
106	guez	123
107	guez	124	is2 = 0
108	guez	3
109	guez	124	IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
110			.or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
111			IF (clon <= 0.) THEN
112			is2 = 1
113	guez	3
114	guez	124	do while (rlonm025(is2) < - pi .and. is2 < iim)
115			is2 = is2 + 1
116			end do
117
118			if (rlonm025(is2) < - pi) then
119			print *, 'Rlonm025 plus petit que - pi !'
120			STOP 1
121			end if
122			ELSE
123			is2 = iim
124
125			do while (rlonm025(is2) > pi .and. is2 > 1)
126			is2 = is2 - 1
127			end do
128
129			if (rlonm025(is2) > pi) then
130			print *, 'Rlonm025 plus grand que pi !'
131			STOP 1
132			end if
133			END IF
134			END IF
135
136			call principal_cshift(is2, rlonm025, xprimm025)
137			call principal_cshift(is2, rlonv, xprimv)
138			call principal_cshift(is2, rlonu, xprimu)
139			call principal_cshift(is2, rlonp025, xprimp025)
140
141			forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
142			print , "Minimum longitude step:", MINval(d_rlonv) 180. / pi, "degrees"
143			print , "Maximum longitude step:", MAXval(d_rlonv) 180. / pi, "degrees"
144
145	guez	128	! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
146	guez	119	DO i = 1, iim + 1
147			IF (rlonp025(i) < rlonv(i)) THEN
148	guez	121	print *, 'rlonp025(', i, ') = ', rlonp025(i)
149			print *, "< rlonv(", i, ") = ", rlonv(i)
150	guez	119	STOP 1
151			END IF
152
153			IF (rlonv(i) < rlonm025(i)) THEN
154	guez	121	print *, 'rlonv(', i, ') = ', rlonv(i)
155			print *, "< rlonm025(", i, ") = ", rlonm025(i)
156	guez	119	STOP 1
157			END IF
158
159			IF (rlonp025(i) > rlonu(i)) THEN
160	guez	120	print *, 'rlonp025(', i, ') = ', rlonp025(i)
161			print *, "> rlonu(", i, ") = ", rlonu(i)
162	guez	119	STOP 1
163			END IF
164			END DO
165
166	guez	78	END SUBROUTINE fxhyp
167
168			end module fxhyp_m