GCC Code Coverage Report

Directory: ./
File: dyn3d_common/fyhyp_m.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 130 149 87.2%
Branches: 103 118 87.3%
Line Branch Exec Source
1 module fyhyp_m
2
3 IMPLICIT NONE
4
5 contains
6
7 1 SUBROUTINE fyhyp(rlatu, yyprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)
8
9 ! From LMDZ4/libf/dyn3d/fyhyp.F, version 1.2, 2005/06/03 09:11:32
10
11 ! Author: P. Le Van, from analysis by R. Sadourny
12
13 ! Calcule les latitudes et dérivées dans la grille du GCM pour une
14 ! fonction f(y) à dérivée tangente hyperbolique.
15
16 ! Il vaut mieux avoir : grossismy * dzoom < pi / 2
17
18 use coefpoly_m, only: coefpoly
19 use nrtype, only: k8
20 use serre_mod, only: clat, grossismy, dzoomy, tauy
21
22 include "dimensions.h"
23 ! for jjm
24
25 REAL, intent(out):: rlatu(jjm + 1), yyprimu(jjm + 1)
26 REAL, intent(out):: rlatv(jjm)
27 real, intent(out):: rlatu2(jjm), yprimu2(jjm), rlatu1(jjm), yprimu1(jjm)
28
29 ! Local:
30
31 REAL(K8) champmin, champmax
32 INTEGER, PARAMETER:: nmax=30000, nmax2=2*nmax
33 REAL dzoom ! distance totale de la zone du zoom (en radians)
34 REAL(K8) ylat(jjm + 1), yprim(jjm + 1)
35 REAL(K8) yuv
36 REAL(K8), save:: yt(0:nmax2)
37 REAL(K8) fhyp(0:nmax2), beta
38 REAL(K8), save:: ytprim(0:nmax2)
39 REAL(K8) fxm(0:nmax2)
40 REAL(K8), save:: yf(0:nmax2)
41 REAL(K8) yypr(0:nmax2)
42 REAL(K8) yvrai(jjm + 1), yprimm(jjm + 1), ylatt(jjm + 1)
43 REAL(K8) pi, pis2, epsilon, y0, pisjm
44 REAL(K8) yo1, yi, ylon2, ymoy, yprimin
45 REAL(K8) yfi, yf1, ffdy
46 REAL(K8) ypn, deply, y00
47 SAVE y00, deply
48
49 INTEGER i, j, it, ik, iter, jlat
50 INTEGER jpn, jjpn
51 SAVE jpn
52 REAL(K8) a0, a1, a2, a3, yi2, heavyy0, heavyy0m
53 REAL(K8) fa(0:nmax2), fb(0:nmax2)
54 REAL y0min, y0max
55
56 REAL(K8) heavyside
57
58 !-------------------------------------------------------------------
59
60 1 print *, "Call sequence information: fyhyp"
61
62 pi = 2.*asin(1.)
63 pis2 = pi/2.
64 pisjm = pi/real(jjm)
65 epsilon = 1e-3
66 1 y0 = clat*pi/180.
67 1 dzoom = dzoomy*pi
68 1 print *, 'yzoom(rad), grossismy, tauy, dzoom (rad):'
69 1 print *, y0, grossismy, tauy, dzoom
70
71
2/2
✓ Branch 0 taken 60001 times.
✓ Branch 1 taken 1 times.
 60002 DO i = 0, nmax2
72 60002 yt(i) = -pis2 + real(i)*pi/nmax2
73 END DO
74
75 1 heavyy0m = heavyside(-y0)
76 1 heavyy0 = heavyside(y0)
77 1 y0min = 2.*y0*heavyy0m - pis2
78 1 y0max = 2.*y0*heavyy0 + pis2
79
80
2/2
✓ Branch 0 taken 60001 times.
✓ Branch 1 taken 1 times.
 60002 fa = 999.999
81
2/2
✓ Branch 0 taken 60001 times.
✓ Branch 1 taken 1 times.
 60002 fb = 999.999
82
83
2/2
✓ Branch 0 taken 60001 times.
✓ Branch 1 taken 1 times.
 60002 DO i = 0, nmax2
84
2/2
✓ Branch 0 taken 30000 times.
✓ Branch 1 taken 30001 times.
 60001 IF (yt(i)<y0) THEN
85 30000 fa(i) = tauy*(yt(i)-y0 + dzoom/2.)
86 30000 fb(i) = (yt(i)-2.*y0*heavyy0m + pis2)*(y0-yt(i))
87
2/2
✓ Branch 0 taken 30000 times.
✓ Branch 1 taken 1 times.
 30001 ELSE IF (yt(i)>y0) THEN
88 30000 fa(i) = tauy*(y0-yt(i) + dzoom/2.)
89 30000 fb(i) = (2.*y0*heavyy0-yt(i) + pis2)*(yt(i)-y0)
90 END IF
91
92
2/2
✓ Branch 0 taken 488 times.
✓ Branch 1 taken 59513 times.
 60001 IF (200.*fb(i)<-fa(i)) THEN
93 488 fhyp(i) = -1.
94
2/2
✓ Branch 0 taken 84 times.
✓ Branch 1 taken 59429 times.
 59513 ELSE IF (200.*fb(i)<fa(i)) THEN
95 84 fhyp(i) = 1.
96 ELSE
97 59429 fhyp(i) = tanh(fa(i)/fb(i))
98 END IF
99
100
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 60000 times.
 60001 IF (yt(i)==y0) fhyp(i) = 1.
101
4/4
✓ Branch 0 taken 60000 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 59999 times.
 60002 IF (yt(i)==y0min .OR. yt(i)==y0max) fhyp(i) = -1.
102 END DO
103
104 ! Calcul de beta
105
106 ffdy = 0.
107
108
2/2
✓ Branch 0 taken 60000 times.
✓ Branch 1 taken 1 times.
 60001 DO i = 1, nmax2
109 60000 ymoy = 0.5*(yt(i-1) + yt(i))
110
2/2
✓ Branch 0 taken 30000 times.
✓ Branch 1 taken 30000 times.
 60000 IF (ymoy<y0) THEN
111 30000 fa(i) = tauy*(ymoy-y0 + dzoom/2.)
112 30000 fb(i) = (ymoy-2.*y0*heavyy0m + pis2)*(y0-ymoy)
113
1/2
✓ Branch 0 taken 30000 times.
✗ Branch 1 not taken.
 30000 ELSE IF (ymoy>y0) THEN
114 30000 fa(i) = tauy*(y0-ymoy + dzoom/2.)
115 30000 fb(i) = (2.*y0*heavyy0-ymoy + pis2)*(ymoy-y0)
116 END IF
117
118
2/2
✓ Branch 0 taken 486 times.
✓ Branch 1 taken 59514 times.
 60000 IF (200.*fb(i)<-fa(i)) THEN
119 486 fxm(i) = -1.
120
2/2
✓ Branch 0 taken 86 times.
✓ Branch 1 taken 59428 times.
 59514 ELSE IF (200.*fb(i)<fa(i)) THEN
121 86 fxm(i) = 1.
122 ELSE
123 59428 fxm(i) = tanh(fa(i)/fb(i))
124 END IF
125
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 60000 times.
 60000 IF (ymoy==y0) fxm(i) = 1.
126
3/4
✓ Branch 0 taken 60000 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 59999 times.
 60000 IF (ymoy==y0min .OR. yt(i)==y0max) fxm(i) = -1.
127 60001 ffdy = ffdy + fxm(i)*(yt(i)-yt(i-1))
128 END DO
129
130 1 beta = (grossismy*ffdy-pi)/(ffdy-pi)
131
132
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 IF (2. * beta - grossismy <= 0.) THEN
133 print *, 'Attention ! La valeur beta calculee dans la routine fyhyp ' &
134 // 'est mauvaise. Modifier les valeurs de grossismy, tauy ou ' &
135 // 'dzoomy et relancer.'
136 ✗ STOP 1
137 END IF
138
139 ! calcul de Ytprim
140
141
2/2
✓ Branch 0 taken 60001 times.
✓ Branch 1 taken 1 times.
 60002 DO i = 0, nmax2
142 60002 ytprim(i) = beta + (grossismy-beta)*fhyp(i)
143 END DO
144
145 ! Calcul de Yf
146
147 1 yf(0) = -pis2
148
2/2
✓ Branch 0 taken 60000 times.
✓ Branch 1 taken 1 times.
 60001 DO i = 1, nmax2
149 60001 yypr(i) = beta + (grossismy-beta)*fxm(i)
150 END DO
151
152
2/2
✓ Branch 0 taken 60000 times.
✓ Branch 1 taken 1 times.
 60001 DO i = 1, nmax2
153 60001 yf(i) = yf(i-1) + yypr(i)*(yt(i)-yt(i-1))
154 END DO
155
156 ! yuv = 0. si calcul des latitudes aux pts. U
157 ! yuv = 0.5 si calcul des latitudes aux pts. V
158
159
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
 5 loop_ik: DO ik = 1, 4
160
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
 4 IF (ik==1) THEN
161 yuv = 0.
162 jlat = jjm + 1
163
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
 3 ELSE IF (ik==2) THEN
164 yuv = 0.5
165 jlat = jjm
166
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 ELSE IF (ik==3) THEN
167 yuv = 0.25
168 jlat = jjm
169 ELSE IF (ik==4) THEN
170 yuv = 0.75
171 jlat = jjm
172 END IF
173
174 yo1 = 0.
175
2/2
✓ Branch 0 taken 129 times.
✓ Branch 1 taken 4 times.
 133 DO j = 1, jlat
176 yo1 = 0.
177 129 ylon2 = -pis2 + pisjm*(real(j) + yuv-1.)
178 yfi = ylon2
179
180 it = nmax2
181
4/4
✓ Branch 0 taken 3870178 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 128 times.
✓ Branch 3 taken 3870050 times.
 3870179 DO while (it >= 1 .and. yfi < yf(it))
182 3870178 it = it - 1
183 END DO
184
185 129 yi = yt(it)
186
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 128 times.
 129 IF (it==nmax2) THEN
187 it = nmax2 - 1
188 1 yf(it + 1) = pis2
189 END IF
190
191 ! Interpolation entre yi(it) et yi(it + 1) pour avoir Y(yi)
192 ! et Y'(yi)
193
194 CALL coefpoly(yf(it), yf(it + 1), ytprim(it), ytprim(it + 1), &
195 129 yt(it), yt(it + 1), a0, a1, a2, a3)
196
197 129 yf1 = yf(it)
198 129 yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi
199
200 iter = 1
201 128 DO
202 257 yi = yi - (yf1-yfi)/yprimin
203
2/2
✓ Branch 0 taken 128 times.
✓ Branch 1 taken 129 times.
 257 IF (abs(yi-yo1)<=epsilon .or. iter == 300) exit
204 yo1 = yi
205 128 yi2 = yi*yi
206 128 yf1 = a0 + a1*yi + a2*yi2 + a3*yi2*yi
207 128 yprimin = a1 + 2.*a2*yi + 3.*a3*yi2
208 END DO
209
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 129 times.
 129 if (abs(yi-yo1) > epsilon) then
210 ✗ print *, 'Pas de solution.', j, ylon2
211 ✗ STOP 1
212 end if
213
214 129 yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi
215 129 yprim(j) = pi/(jjm*yprimin)
216 262 yvrai(j) = yi
217 END DO
218
219
2/2
✓ Branch 0 taken 125 times.
✓ Branch 1 taken 4 times.
 129 DO j = 1, jlat - 1
220
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 125 times.
 129 IF (yvrai(j + 1)<yvrai(j)) THEN
221 ✗ print *, 'Problème avec rlat(', j + 1, ') plus petit que rlat(', &
222 ✗ j, ')'
223 ✗ STOP 1
224 END IF
225 END DO
226
227 4 print *, 'Reorganisation des latitudes pour avoir entre - pi/2 et pi/2'
228
229
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3 times.
 4 IF (ik==1) THEN
230 ypn = pis2
231
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 DO j = jjm + 1, 1, -1
232
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 IF (yvrai(j)<=ypn) exit
233 END DO
234
235 1 jpn = j
236 1 y00 = yvrai(jpn)
237 1 deply = pis2 - y00
238 END IF
239
240
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 DO j = 1, jjm + 1 - jpn
241 ✗ ylatt(j) = -pis2 - y00 + yvrai(jpn + j-1)
242 4 yprimm(j) = yprim(jpn + j-1)
243 END DO
244
245 jjpn = jpn
246
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
 4 IF (jlat==jjm) jjpn = jpn - 1
247
248
2/2
✓ Branch 0 taken 129 times.
✓ Branch 1 taken 4 times.
 133 DO j = 1, jjpn
249 129 ylatt(j + jjm + 1-jpn) = yvrai(j) + deply
250 133 yprimm(j + jjm + 1-jpn) = yprim(j)
251 END DO
252
253 ! Fin de la reorganisation
254
255
2/2
✓ Branch 0 taken 129 times.
✓ Branch 1 taken 4 times.
 133 DO j = 1, jlat
256 129 ylat(j) = ylatt(jlat + 1-j)
257 133 yprim(j) = yprimm(jlat + 1-j)
258 END DO
259
260
2/2
✓ Branch 0 taken 129 times.
✓ Branch 1 taken 4 times.
 133 DO j = 1, jlat
261 133 yvrai(j) = ylat(j)*180./pi
262 END DO
263
264
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3 times.
 5 IF (ik==1) THEN
265
2/2
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 1 times.
 34 DO j = 1, jjm + 1
266 33 rlatu(j) = ylat(j)
267 34 yyprimu(j) = yprim(j)
268 END DO
269
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
 3 ELSE IF (ik==2) THEN
270
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 32 times.
 33 DO j = 1, jjm
271 33 rlatv(j) = ylat(j)
272 END DO
273
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 ELSE IF (ik==3) THEN
274
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 32 times.
 33 DO j = 1, jjm
275 32 rlatu2(j) = ylat(j)
276 33 yprimu2(j) = yprim(j)
277 END DO
278 ELSE IF (ik==4) THEN
279
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 32 times.
 33 DO j = 1, jjm
280 32 rlatu1(j) = ylat(j)
281 33 yprimu1(j) = yprim(j)
282 END DO
283 END IF
284 END DO loop_ik
285
286
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 times.
 33 DO j = 1, jjm
287 33 ylat(j) = rlatu(j) - rlatu(j + 1)
288 END DO
289 1 champmin = 1e12
290 1 champmax = -1e12
291
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 times.
 33 DO j = 1, jjm
292 32 champmin = min(champmin, ylat(j))
293 33 champmax = max(champmax, ylat(j))
294 END DO
295 1 champmin = champmin*180./pi
296 1 champmax = champmax*180./pi
297
298
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 times.
 33 DO j = 1, jjm
299
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
 32 IF (rlatu1(j) <= rlatu2(j)) THEN
300 ✗ print *, 'Attention ! rlatu1 < rlatu2 ', rlatu1(j), rlatu2(j), j
301 ✗ STOP 13
302 ENDIF
303
304
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
 32 IF (rlatu2(j) <= rlatu(j+1)) THEN
305 ✗ print *, 'Attention ! rlatu2 < rlatup1 ', rlatu2(j), rlatu(j+1), j
306 ✗ STOP 14
307 ENDIF
308
309
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
 32 IF (rlatu(j) <= rlatu1(j)) THEN
310 ✗ print *, ' Attention ! rlatu < rlatu1 ', rlatu(j), rlatu1(j), j
311 ✗ STOP 15
312 ENDIF
313
314
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
 32 IF (rlatv(j) <= rlatu2(j)) THEN
315 ✗ print *, ' Attention ! rlatv < rlatu2 ', rlatv(j), rlatu2(j), j
316 ✗ STOP 16
317 ENDIF
318
319
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
 32 IF (rlatv(j) >= rlatu1(j)) THEN
320 ✗ print *, ' Attention ! rlatv > rlatu1 ', rlatv(j), rlatu1(j), j
321 ✗ STOP 17
322 ENDIF
323
324
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
 33 IF (rlatv(j) >= rlatu(j)) THEN
325 ✗ print *, ' Attention ! rlatv > rlatu ', rlatv(j), rlatu(j), j
326 ✗ STOP 18
327 ENDIF
328 ENDDO
329
330 1 print *, 'Latitudes'
331 1 print 3, champmin, champmax
332
333 3 Format(1x, ' Au centre du zoom, la longueur de la maille est', &
334 ' d environ ', f0.2, ' degres ', /, &
335 ' alors que la maille en dehors de la zone du zoom est ', &
336 "d'environ ", f0.2, ' degres ')
337
338 1 END SUBROUTINE fyhyp
339
340 end module fyhyp_m
341