[1441] | 1 | module fxhyp_m |
---|
[1] | 2 | |
---|
[1441] | 3 | IMPLICIT NONE |
---|
[1] | 4 | |
---|
[1441] | 5 | contains |
---|
[1] | 6 | |
---|
[1441] | 7 | SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025) |
---|
[1] | 8 | |
---|
[1441] | 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 |
---|
[1] | 11 | |
---|
[1441] | 12 | ! Calcule les longitudes et dérivées dans la grille du GCM pour |
---|
| 13 | ! une fonction f(x) à dérivée tangente hyperbolique. |
---|
[1] | 14 | |
---|
[1441] | 15 | ! Il vaut mieux avoir : grossismx \times dzoom < pi |
---|
[1] | 16 | |
---|
[1441] | 17 | ! Le premier point scalaire pour une grille regulière (grossismx = |
---|
| 18 | ! 1., taux=0., clon=0.) est à - 180 degrés. |
---|
[1] | 19 | |
---|
[1441] | 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 |
---|
[1443] | 24 | use serre_mod, only: clon,grossismx,dzoomx,taux |
---|
[1] | 25 | |
---|
[1441] | 26 | include "dimensions.h" |
---|
| 27 | ! for iim |
---|
[1] | 28 | |
---|
[1441] | 29 | REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1) |
---|
| 30 | real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1) |
---|
[1] | 31 | |
---|
[1441] | 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 |
---|
[841] | 42 | |
---|
[1441] | 43 | !---------------------------------------------------------------------- |
---|
[841] | 44 | |
---|
[1441] | 45 | print *, "Call sequence information: fxhyp" |
---|
[841] | 46 | |
---|
[1441] | 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 |
---|
[841] | 52 | |
---|
[1441] | 53 | xprimm025(:)=1. |
---|
| 54 | xprimv(:)=1. |
---|
| 55 | xprimu(:)=1. |
---|
| 56 | xprimp025(:)=1. |
---|
| 57 | else test_iim |
---|
| 58 | test_grossismx: if (grossismx == 1.) then |
---|
| 59 | step = twopi / iim |
---|
[1] | 60 | |
---|
[1441] | 61 | xprimm025(:iim) = step |
---|
| 62 | xprimp025(:iim) = step |
---|
| 63 | xprimv(:iim) = step |
---|
| 64 | xprimu(:iim) = step |
---|
[1] | 65 | |
---|
[1441] | 66 | rlonv(:iim) = arth(- pi + clon / 180. * pi, step, iim) |
---|
| 67 | rlonm025(:iim) = rlonv(:iim) - 0.25 * step |
---|
| 68 | rlonp025(:iim) = rlonv(:iim) + 0.25 * step |
---|
| 69 | 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) |
---|
[1] | 73 | |
---|
[1441] | 74 | ! Compute fhyp: |
---|
| 75 | DO i = nmax, 2 * nmax |
---|
| 76 | fa = taux * (dzoom / 2. - xtild(i)) |
---|
| 77 | fb = xtild(i) * (pi_d - xtild(i)) |
---|
[1] | 78 | |
---|
[1441] | 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 |
---|
[1] | 94 | |
---|
[1441] | 95 | IF (xtild(i) == 0.) fhyp(i) = 1. |
---|
| 96 | IF (xtild(i) == pi_d) fhyp(i) = -1. |
---|
| 97 | END DO |
---|
[1] | 98 | |
---|
[1441] | 99 | ! Calcul de beta |
---|
[1] | 100 | |
---|
[1441] | 101 | ffdx = 0. |
---|
[1] | 102 | |
---|
[1441] | 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) |
---|
[1] | 107 | |
---|
[1441] | 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 |
---|
[1] | 123 | |
---|
[1441] | 124 | IF (xmoy == 0.) fxm = 1. |
---|
| 125 | IF (xmoy == pi_d) fxm = -1. |
---|
[1] | 126 | |
---|
[1441] | 127 | ffdx = ffdx + fxm * (xtild(i) - xtild(i-1)) |
---|
| 128 | END DO |
---|
[1] | 129 | |
---|
[1441] | 130 | print *, "ffdx = ", ffdx |
---|
| 131 | beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d) |
---|
| 132 | print *, "beta = ", beta |
---|
[1] | 133 | |
---|
[1441] | 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 |
---|
[1] | 139 | |
---|
[1441] | 140 | ! calcul de Xprimt |
---|
| 141 | Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp |
---|
| 142 | xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1) |
---|
[1] | 143 | |
---|
[1441] | 144 | ! Calcul de Xf |
---|
[1] | 145 | |
---|
[1441] | 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) |
---|
[1] | 150 | |
---|
[1441] | 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 |
---|
[1] | 158 | |
---|
[1441] | 159 | IF (xmoy == 0.) fxm = 1. |
---|
| 160 | IF (xmoy == pi_d) fxm = -1. |
---|
| 161 | xxpr(i) = beta + (grossismx - beta) * fxm |
---|
| 162 | END DO |
---|
[1] | 163 | |
---|
[1441] | 164 | xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1) |
---|
[1] | 165 | |
---|
[1441] | 166 | Xf(0) = - pi_d |
---|
[1] | 167 | |
---|
[1441] | 168 | DO i=1, 2 * nmax - 1 |
---|
| 169 | Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1)) |
---|
| 170 | END DO |
---|
[1] | 171 | |
---|
[1441] | 172 | Xf(2 * nmax) = pi_d |
---|
[1] | 173 | |
---|
[1441] | 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 |
---|
[1] | 183 | |
---|
[1441] | 184 | is2 = 0 |
---|
[1] | 185 | |
---|
[1441] | 186 | IF (MINval(rlonm025(:iim)) < - pi - 0.1 & |
---|
| 187 | .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN |
---|
| 188 | IF (clon <= 0.) THEN |
---|
| 189 | is2 = 1 |
---|
[1] | 190 | |
---|
[1441] | 191 | do while (rlonm025(is2) < - pi .and. is2 < iim) |
---|
| 192 | is2 = is2 + 1 |
---|
| 193 | end do |
---|
[1] | 194 | |
---|
[1441] | 195 | if (rlonm025(is2) < - pi) then |
---|
| 196 | print *, 'Rlonm025 plus petit que - pi !' |
---|
| 197 | STOP 1 |
---|
| 198 | end if |
---|
| 199 | ELSE |
---|
| 200 | is2 = iim |
---|
[1] | 201 | |
---|
[1441] | 202 | do while (rlonm025(is2) > pi .and. is2 > 1) |
---|
| 203 | is2 = is2 - 1 |
---|
| 204 | end do |
---|
[1] | 205 | |
---|
[1441] | 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 |
---|
[1] | 212 | |
---|
[1441] | 213 | call principal_cshift(is2, rlonm025, xprimm025) |
---|
| 214 | call principal_cshift(is2, rlonv, xprimv) |
---|
| 215 | call principal_cshift(is2, rlonu, xprimu) |
---|
| 216 | call principal_cshift(is2, rlonp025, xprimp025) |
---|
[1] | 217 | |
---|
[1441] | 218 | forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i) |
---|
| 219 | print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, & |
---|
| 220 | "degrees" |
---|
| 221 | print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, & |
---|
| 222 | "degrees" |
---|
[1] | 223 | |
---|
[1441] | 224 | ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu: |
---|
| 225 | DO i = 1, iim + 1 |
---|
| 226 | IF (rlonp025(i) < rlonv(i)) THEN |
---|
| 227 | print *, 'rlonp025(', i, ') = ', rlonp025(i) |
---|
| 228 | print *, "< rlonv(", i, ") = ", rlonv(i) |
---|
| 229 | STOP 1 |
---|
| 230 | END IF |
---|
[1] | 231 | |
---|
[1441] | 232 | IF (rlonv(i) < rlonm025(i)) THEN |
---|
| 233 | print *, 'rlonv(', i, ') = ', rlonv(i) |
---|
| 234 | print *, "< rlonm025(", i, ") = ", rlonm025(i) |
---|
| 235 | STOP 1 |
---|
| 236 | END IF |
---|
[1] | 237 | |
---|
[1441] | 238 | 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 |
---|
[1] | 245 | |
---|
[1441] | 246 | END SUBROUTINE fxhyp |
---|
[1] | 247 | |
---|
[1441] | 248 | end module fxhyp_m |
---|