Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	dyn3d_common/fxhyp_m.F90	Lines:	28	109	25.7 %
Date:	2023-06-30 12:56:34	Branches:	53	144	36.8 %


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 arth_m, only: arth
    use invert_zoom_x_m, only: invert_zoom_x, nmax
    use nrtype, only: pi, pi_d, twopi, twopi_d, k8
    use principal_cshift_m, only: principal_cshift
    use serre_mod, only: clon, grossismx, dzoomx, taux

    include "dimensions.h"
    ! for iim

    REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
    real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)

    ! Local:
    real rlonm025(iim + 1), rlonp025(iim + 1)
    REAL dzoom, step
    real d_rlonv(iim)
    REAL(K8) xtild(0:2 * nmax)
    REAL(K8) fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
    REAL(K8) Xf(0:2 * nmax), xxpr(2 * nmax)
    REAL(K8) fa, fb
    INTEGER i, is2
    REAL(K8) xmoy, fxm

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

    print *, "Call sequence information: fxhyp"

    test_iim: if (iim==1) then
       rlonv(1)=0.
       rlonu(1)=pi
       rlonv(2)=rlonv(1)+twopi
       rlonu(2)=rlonu(1)+twopi

       xprimm025(:)=1.
       xprimv(:)=1.
       xprimu(:)=1.
       xprimp025(:)=1.
    else test_iim
       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 / 180. * pi, 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)

          ! Compute fhyp:
          DO i = nmax, 2 * nmax
             fa = taux * (dzoom / 2. - xtild(i))
             fb = xtild(i) * (pi_d - xtild(i))

             IF (200. * fb < - fa) THEN
                fhyp(i) = - 1.
             ELSE IF (200. * fb < fa) THEN
                fhyp(i) = 1.
             ELSE
                IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
                   IF (200. * fb + fa < 1e-10) THEN
                      fhyp(i) = - 1.
                   ELSE IF (200. * fb - fa < 1e-10) THEN
                      fhyp(i) = 1.
                   END IF
                ELSE
                   fhyp(i) = TANH(fa / fb)
                END IF
             END IF

             IF (xtild(i) == 0.) fhyp(i) = 1.
             IF (xtild(i) == pi_d) fhyp(i) = -1.
          END DO

          ! Calcul de beta

          ffdx = 0.

          DO i = nmax + 1, 2 * nmax
             xmoy = 0.5 * (xtild(i-1) + xtild(i))
             fa = taux * (dzoom / 2. - xmoy)
             fb = xmoy * (pi_d - xmoy)

             IF (200. * fb < - fa) THEN
                fxm = - 1.
             ELSE IF (200. * fb < fa) THEN
                fxm = 1.
             ELSE
                IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
                   IF (200. * fb + fa < 1e-10) THEN
                      fxm = - 1.
                   ELSE IF (200. * fb - fa < 1e-10) THEN
                      fxm = 1.
                   END IF
                ELSE
                   fxm = TANH(fa / fb)
                END IF
             END IF

             IF (xmoy == 0.) fxm = 1.
             IF (xmoy == pi_d) fxm = -1.

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

          print *, "ffdx = ", ffdx
          beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
          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

          DO i = nmax + 1, 2 * nmax
             xmoy = 0.5 * (xtild(i-1) + xtild(i))
             fa = taux * (dzoom / 2. - xmoy)
             fb = xmoy * (pi_d - xmoy)

             IF (200. * fb < - fa) THEN
                fxm = - 1.
             ELSE IF (200. * fb < fa) THEN
                fxm = 1.
             ELSE
                fxm = TANH(fa / fb)
             END IF

             IF (xmoy == 0.) fxm = 1.
             IF (xmoy == pi_d) fxm = -1.
             xxpr(i) = beta + (grossismx - beta) * fxm
          END DO

          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.25_k8)
          call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
               xuv = 0._k8)
          call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
               xuv = 0.5_k8)
          call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), &
               xprimp025(:iim), xuv = 0.25_k8)
       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 if test_iim

  END SUBROUTINE fxhyp

end module fxhyp_m


Generated by: GCOVR (Version 4.2)

Line	Branch	Exec	Source
1			module fxhyp_m
2
3			IMPLICIT NONE
4
5			contains
6
7	✗✓✗✓ ✗✓✗✓ ✗✓✗✓	1	SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
8
9			! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32
10			! Author: P. Le Van, from formulas by R. Sadourny
11
12			! Calcule les longitudes et dérivées dans la grille du GCM pour
13			! une fonction f(x) à dérivée tangente hyperbolique.
14
15			! Il vaut mieux avoir : grossismx \times dzoom < pi
16
17			! Le premier point scalaire pour une grille regulière (grossismx =
18			! 1., taux=0., clon=0.) est à - 180 degrés.
19
20			use arth_m, only: arth
21			use invert_zoom_x_m, only: invert_zoom_x, nmax
22			use nrtype, only: pi, pi_d, twopi, twopi_d, k8
23			use principal_cshift_m, only: principal_cshift
24			use serre_mod, only: clon, grossismx, dzoomx, taux
25
26			include "dimensions.h"
27			! for iim
28
29			REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
30			real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
31
32			! Local:
33			real rlonm025(iim + 1), rlonp025(iim + 1)
34			REAL dzoom, step
35			real d_rlonv(iim)
36			REAL(K8) xtild(0:2 * nmax)
37			REAL(K8) fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
38			REAL(K8) Xf(0:2 * nmax), xxpr(2 * nmax)
39			REAL(K8) fa, fb
40			INTEGER i, is2
41			REAL(K8) xmoy, fxm
42
43			!----------------------------------------------------------------------
44
45		1	print *, "Call sequence information: fxhyp"
46
47			test_iim: if (iim==1) then
48			rlonv(1)=0.
49			rlonu(1)=pi
50			rlonv(2)=rlonv(1)+twopi
51			rlonu(2)=rlonu(1)+twopi
52
53			xprimm025(:)=1.
54			xprimv(:)=1.
55			xprimu(:)=1.
56			xprimp025(:)=1.
57			else test_iim
58	✓✗	1	test_grossismx: if (grossismx == 1.) then
59		1	step = twopi / iim
60
61	✓✓	33	xprimm025(:iim) = step
62	✓✓	33	xprimp025(:iim) = step
63	✓✓	33	xprimv(:iim) = step
64	✓✓	33	xprimu(:iim) = step
65
66		1	rlonv(:iim) = arth(- pi + clon / 180. * pi, step, iim)
67	✓✓	33	rlonm025(:iim) = rlonv(:iim) - 0.25 * step
68	✓✓	33	rlonp025(:iim) = rlonv(:iim) + 0.25 * step
69	✓✓	33	rlonu(:iim) = rlonv(:iim) + 0.5 * step
70			else test_grossismx
71			dzoom = dzoomx * twopi_d
72			xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
73
74			! Compute fhyp:
75			DO i = nmax, 2 * nmax
76			fa = taux * (dzoom / 2. - xtild(i))
77			fb = xtild(i) * (pi_d - xtild(i))
78
79			IF (200. * fb < - fa) THEN
80			fhyp(i) = - 1.
81			ELSE IF (200. * fb < fa) THEN
82			fhyp(i) = 1.
83			ELSE
84			IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
85			IF (200. * fb + fa < 1e-10) THEN
86			fhyp(i) = - 1.
87			ELSE IF (200. * fb - fa < 1e-10) THEN
88			fhyp(i) = 1.
89			END IF
90			ELSE
91			fhyp(i) = TANH(fa / fb)
92			END IF
93			END IF
94
95			IF (xtild(i) == 0.) fhyp(i) = 1.
96			IF (xtild(i) == pi_d) fhyp(i) = -1.
97			END DO
98
99			! Calcul de beta
100
101			ffdx = 0.
102
103			DO i = nmax + 1, 2 * nmax
104			xmoy = 0.5 * (xtild(i-1) + xtild(i))
105			fa = taux * (dzoom / 2. - xmoy)
106			fb = xmoy * (pi_d - xmoy)
107
108			IF (200. * fb < - fa) THEN
109			fxm = - 1.
110			ELSE IF (200. * fb < fa) THEN
111			fxm = 1.
112			ELSE
113			IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
114			IF (200. * fb + fa < 1e-10) THEN
115			fxm = - 1.
116			ELSE IF (200. * fb - fa < 1e-10) THEN
117			fxm = 1.
118			END IF
119			ELSE
120			fxm = TANH(fa / fb)
121			END IF
122			END IF
123
124			IF (xmoy == 0.) fxm = 1.
125			IF (xmoy == pi_d) fxm = -1.
126
127			ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
128			END DO
129
130			print *, "ffdx = ", ffdx
131			beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
132			print *, "beta = ", beta
133
134			IF (2. * beta - grossismx <= 0.) THEN
135			print *, 'Bad choice of grossismx, taux, dzoomx.'
136			print *, 'Decrease dzoomx or grossismx.'
137			STOP 1
138			END IF
139
140			! calcul de Xprimt
141			Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
142			xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
143
144			! Calcul de Xf
145
146			DO i = nmax + 1, 2 * nmax
147			xmoy = 0.5 * (xtild(i-1) + xtild(i))
148			fa = taux * (dzoom / 2. - xmoy)
149			fb = xmoy * (pi_d - xmoy)
150
151			IF (200. * fb < - fa) THEN
152			fxm = - 1.
153			ELSE IF (200. * fb < fa) THEN
154			fxm = 1.
155			ELSE
156			fxm = TANH(fa / fb)
157			END IF
158
159			IF (xmoy == 0.) fxm = 1.
160			IF (xmoy == pi_d) fxm = -1.
161			xxpr(i) = beta + (grossismx - beta) * fxm
162			END DO
163
164			xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
165
166			Xf(0) = - pi_d
167
168			DO i=1, 2 * nmax - 1
169			Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
170			END DO
171
172			Xf(2 * nmax) = pi_d
173
174			call invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), &
175			xprimm025(:iim), xuv = - 0.25_k8)
176			call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
177			xuv = 0._k8)
178			call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
179			xuv = 0.5_k8)
180			call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), &
181			xprimp025(:iim), xuv = 0.25_k8)
182			end if test_grossismx
183
184		1	is2 = 0
185
186			IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
187	✗✓✓✗ ✓✓✓✓ ✗✓✓✗ ✓✓✓✗ ✓✗✗✓	68	.or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
188			IF (clon <= 0.) THEN
189			is2 = 1
190
191			do while (rlonm025(is2) < - pi .and. is2 < iim)
192			is2 = is2 + 1
193			end do
194
195			if (rlonm025(is2) < - pi) then
196			print *, 'Rlonm025 plus petit que - pi !'
197			STOP 1
198			end if
199			ELSE
200			is2 = iim
201
202			do while (rlonm025(is2) > pi .and. is2 > 1)
203			is2 = is2 - 1
204			end do
205
206			if (rlonm025(is2) > pi) then
207			print *, 'Rlonm025 plus grand que pi !'
208			STOP 1
209			end if
210			END IF
211			END IF
212
213		1	call principal_cshift(is2, rlonm025, xprimm025)
214		1	call principal_cshift(is2, rlonv, xprimv)
215		1	call principal_cshift(is2, rlonu, xprimu)
216		1	call principal_cshift(is2, rlonp025, xprimp025)
217
218	✓✓	33	forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
219	✗✓✓✗ ✓✓✓✓	34	print , "Minimum longitude step:", MINval(d_rlonv) 180. / pi, &
220		2	"degrees"
221	✗✓✓✗ ✓✓✓✓	34	print , "Maximum longitude step:", MAXval(d_rlonv) 180. / pi, &
222		2	"degrees"
223
224			! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
225	✓✓	34	DO i = 1, iim + 1
226	✗✓	33	IF (rlonp025(i) < rlonv(i)) THEN
227			print *, 'rlonp025(', i, ') = ', rlonp025(i)
228			print *, "< rlonv(", i, ") = ", rlonv(i)
229			STOP 1
230			END IF
231
232	✗✓	33	IF (rlonv(i) < rlonm025(i)) THEN
233			print *, 'rlonv(', i, ') = ', rlonv(i)
234			print *, "< rlonm025(", i, ") = ", rlonm025(i)
235			STOP 1
236			END IF
237
238	✗✓	34	IF (rlonp025(i) > rlonu(i)) THEN
239			print *, 'rlonp025(', i, ') = ', rlonp025(i)
240			print *, "> rlonu(", i, ") = ", rlonu(i)
241			STOP 1
242			END IF
243			END DO
244			end if test_iim
245
246		1	END SUBROUTINE fxhyp
247
248			end module fxhyp_m