GCC Code Coverage Report

Directory: ./
File: dyn3d_common/inigeom.f
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 254 280 90.7%
Branches: 76 90 84.4%
Line Branch Exec Source
1 !
2 ! $Id: inigeom.F 2603 2016-07-25 09:31:56Z emillour $
3 !
4 c
5 c
6 5 SUBROUTINE inigeom
7 c
8 c Auteur : P. Le Van
9 c
10 c ............ Version du 01/04/2001 ........................
11 c
12 c Calcul des elongations cuij1,.cuij4 , cvij1,..cvij4 aux memes en-
13 c endroits que les aires aireij1,..aireij4 .
14
15 c Choix entre f(y) a derivee sinusoid. ou a derivee tangente hyperbol.
16 c
17 c
18 use fxhyp_m, only: fxhyp
19 use fyhyp_m, only: fyhyp
20 USE comconst_mod, ONLY: pi, g, omeg, rad
21 USE logic_mod, ONLY: fxyhypb, ysinus
22 USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy,
23 & alphax,alphay,taux,tauy,transx,transy,pxo,pyo
24 IMPLICIT NONE
25 c
26 include "dimensions.h"
27 include "paramet.h"
28 include "comgeom2.h"
29 include "comdissnew.h"
30
31 c-----------------------------------------------------------------------
32 c .... Variables locales ....
33 c
34 INTEGER i,j,itmax,itmay,iter
35 REAL cvu(iip1,jjp1),cuv(iip1,jjm)
36 REAL ai14,ai23,airez,rlatp,rlatm,xprm,xprp,un4rad2,yprp,yprm
37 REAL eps,x1,xo1,f,df,xdm,y1,yo1,ydm
38 REAL coslatm,coslatp,radclatm,radclatp
39 REAL cuij1(iip1,jjp1),cuij2(iip1,jjp1),cuij3(iip1,jjp1),
40 * cuij4(iip1,jjp1)
41 REAL cvij1(iip1,jjp1),cvij2(iip1,jjp1),cvij3(iip1,jjp1),
42 * cvij4(iip1,jjp1)
43 REAL rlonvv(iip1),rlatuu(jjp1)
44 REAL rlatu1(jjm),yprimu1(jjm),rlatu2(jjm),yprimu2(jjm) ,
45 * yprimv(jjm),yprimu(jjp1)
46 REAL gamdi_gdiv, gamdi_grot, gamdi_h
47
48 REAL rlonm025(iip1),xprimm025(iip1), rlonp025(iip1),
49 , xprimp025(iip1)
50 SAVE rlatu1,yprimu1,rlatu2,yprimu2,yprimv,yprimu
51 SAVE rlonm025,xprimm025,rlonp025,xprimp025
52
53 REAL SSUM
54 c
55 c
56 c ------------------------------------------------------------------
57 c - -
58 c - calcul des coeff. ( cu, cv , 1./cu**2, 1./cv**2 ) -
59 c - -
60 c ------------------------------------------------------------------
61 c
62 c les coef. ( cu, cv ) permettent de passer des vitesses naturelles
63 c aux vitesses covariantes et contravariantes , ou vice-versa ...
64 c
65 c
66 c on a : u (covariant) = cu * u (naturel) , u(contrav)= u(nat)/cu
67 c v (covariant) = cv * v (naturel) , v(contrav)= v(nat)/cv
68 c
69 c on en tire : u(covariant) = cu * cu * u(contravariant)
70 c v(covariant) = cv * cv * v(contravariant)
71 c
72 c
73 c on a l'application ( x(X) , y(Y) ) avec - im/2 +1 < X < im/2
74 c = =
75 c et - jm/2 < Y < jm/2
76 c = =
77 c
78 c ...................................................
79 c ...................................................
80 c . x est la longitude du point en radians .
81 c . y est la latitude du point en radians .
82 c . .
83 c . on a : cu(i,j) = rad * COS(y) * dx/dX .
84 c . cv( j ) = rad * dy/dY .
85 c . aire(i,j) = cu(i,j) * cv(j) .
86 c . .
87 c . y, dx/dX, dy/dY calcules aux points concernes .
88 c . .
89 c ...................................................
90 c ...................................................
91 c
92 c
93 c
94 c ,
95 c cv , bien que dependant de j uniquement,sera ici indice aussi en i
96 c pour un adressage plus facile en ij .
97 c
98 c
99 c
100 c ************** aux points u et v , *****************
101 c xprimu et xprimv sont respectivement les valeurs de dx/dX
102 c yprimu et yprimv . . . . . . . . . . . dy/dY
103 c rlatu et rlatv . . . . . . . . . . .la latitude
104 c cvu et cv . . . . . . . . . . . cv
105 c
106 c ************** aux points u, v, scalaires, et z ****************
107 c cu, cuv, cuscal, cuz sont respectiv. les valeurs de cu
108 c
109 c
110 c
111 c Exemple de distribution de variables sur la grille dans le
112 c domaine de travail ( X,Y ) .
113 c ................................................................
114 c DX=DY= 1
115 c
116 c
117 c + represente un point scalaire ( p.exp la pression )
118 c > represente la composante zonale du vent
119 c V represente la composante meridienne du vent
120 c o represente la vorticite
121 c
122 c ---- , car aux poles , les comp.zonales covariantes sont nulles
123 c
124 c
125 c
126 c i ->
127 c
128 c 1 2 3 4 5 6 7 8
129 c j
130 c v 1 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --
131 c
132 c V o V o V o V o V o V o V o V o
133 c
134 c 2 + > + > + > + > + > + > + > + >
135 c
136 c V o V o V o V o V o V o V o V o
137 c
138 c 3 + > + > + > + > + > + > + > + >
139 c
140 c V o V o V o V o V o V o V o V o
141 c
142 c 4 + > + > + > + > + > + > + > + >
143 c
144 c V o V o V o V o V o V o V o V o
145 c
146 c 5 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --
147 c
148 c
149 c Ci-dessus, on voit que le nombre de pts.en longitude est egal
150 c a IM = 8
151 c De meme , le nombre d'intervalles entre les 2 poles est egal
152 c a JM = 4
153 c
154 c Les points scalaires ( + ) correspondent donc a des valeurs
155 c entieres de i ( 1 a IM ) et de j ( 1 a JM +1 ) .
156 c
157 c Les vents U ( > ) correspondent a des valeurs semi-
158 c entieres de i ( 1+ 0.5 a IM+ 0.5) et entieres de j ( 1 a JM+1)
159 c
160 c Les vents V ( V ) correspondent a des valeurs entieres
161 c de i ( 1 a IM ) et semi-entieres de j ( 1 +0.5 a JM +0.5)
162 c
163 c
164 c
165 1 WRITE(6,3)
166 3 FORMAT( // 10x,' .... INIGEOM date du 01/06/98 ..... ',
167 * //5x,' Calcul des elongations cu et cv comme sommes des 4 ' /
168 * 5x,' elong. cuij1, .. 4 , cvij1,.. 4 qui les entourent , aux
169 * '/ 5x,' memes endroits que les aires aireij1,...j4 . ' / )
170 c
171 c
172
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 IF( nitergdiv.NE.2 ) THEN
173 1 gamdi_gdiv = coefdis/ ( REAL(nitergdiv) -2. )
174 ELSE
175 ✗ gamdi_gdiv = 0.
176 ENDIF
177
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 IF( nitergrot.NE.2 ) THEN
178 ✗ gamdi_grot = coefdis/ ( REAL(nitergrot) -2. )
179 ELSE
180 1 gamdi_grot = 0.
181 ENDIF
182
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 IF( niterh.NE.2 ) THEN
183 ✗ gamdi_h = coefdis/ ( REAL(niterh) -2. )
184 ELSE
185 1 gamdi_h = 0.
186 ENDIF
187
188 1 WRITE(6,*) ' gamdi_gd ',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis,
189 2 * nitergdiv,nitergrot,niterh
190 c
191 1 pi = 2.* ASIN(1.)
192 c
193 1 WRITE(6,990)
194
195 c ----------------------------------------------------------------
196 c
197
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 IF( .NOT.fxyhypb ) THEN
198 c
199 c
200 ✗ IF( ysinus ) THEN
201 c
202 ✗ WRITE(6,*) ' *** Inigeom , Y = Sinus ( Latitude ) *** '
203 c
204 c .... utilisation de f(x,y ) avec y = sinus de la latitude .....
205
206 CALL fxysinus (rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,
207 , rlatu2,yprimu2,
208 ✗ , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)
209
210 ELSE
211 c
212 ✗ WRITE(6,*) '*** Inigeom , Y = Latitude , der. sinusoid . ***'
213
214 c .... utilisation de f(x,y) a tangente sinusoidale , y etant la latit. ...
215 c
216
217 ✗ pxo = clon *pi /180.
218 ✗ pyo = 2.* clat* pi /180.
219 c
220 c .... determination de transx ( pour le zoom ) par Newton-Raphson ...
221 c
222 itmax = 10
223 eps = .1e-7
224 c
225 xo1 = 0.
226 ✗ DO 10 iter = 1, itmax
227 x1 = xo1
228 ✗ f = x1+ alphax *SIN(x1-pxo)
229 ✗ df = 1.+ alphax *COS(x1-pxo)
230 ✗ x1 = x1 - f/df
231 ✗ xdm = ABS( x1- xo1 )
232 ✗ IF( xdm.LE.eps )GO TO 11
233 xo1 = x1
234 ✗ 10 CONTINUE
235 11 CONTINUE
236 c
237 ✗ transx = xo1
238
239 itmay = 10
240 eps = .1e-7
241 C
242 yo1 = 0.
243 ✗ DO 15 iter = 1,itmay
244 y1 = yo1
245 ✗ f = y1 + alphay* SIN(y1-pyo)
246 ✗ df = 1. + alphay* COS(y1-pyo)
247 ✗ y1 = y1 -f/df
248 ✗ ydm = ABS(y1-yo1)
249 ✗ IF(ydm.LE.eps) GO TO 17
250 yo1 = y1
251 ✗ 15 CONTINUE
252 c
253 17 CONTINUE
254 ✗ transy = yo1
255
256 CALL fxy ( rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,
257 , rlatu2,yprimu2,
258 ✗ , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)
259
260 ENDIF
261 c
262 ELSE
263 c
264 c .... Utilisation de fxyhyper , f(x,y) a derivee tangente hyperbol.
265 c .....................................................................
266
267 1 WRITE(6,*)'*** Inigeom , Y = Latitude , der.tg. hyperbolique ***'
268
269 1 CALL fyhyp(rlatu, yprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)
270 1 CALL fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
271
272 ENDIF
273 c
274 c -------------------------------------------------------------------
275
276 c
277 1 rlatu(1) = ASIN(1.)
278 1 rlatu(jjp1) = - rlatu(1)
279 c
280 c
281 c .... calcul aux poles ....
282 c
283 1 yprimu(1) = 0.
284 1 yprimu(jjp1) = 0.
285 c
286 c
287 1 un4rad2 = 0.25 * rad * rad
288 c
289 c --------------------------------------------------------------------
290 c --------------------------------------------------------------------
291 c - -
292 c - calcul des aires ( aire,aireu,airev, 1./aire, 1./airez ) -
293 c - et de fext , force de coriolis extensive . -
294 c - -
295 c --------------------------------------------------------------------
296 c --------------------------------------------------------------------
297 c
298 c
299 c
300 c A 1 point scalaire P (i,j) de la grille, reguliere en (X,Y) , sont
301 c affectees 4 aires entourant P , calculees respectivement aux points
302 c ( i + 1/4, j - 1/4 ) : aireij1 (i,j)
303 c ( i + 1/4, j + 1/4 ) : aireij2 (i,j)
304 c ( i - 1/4, j + 1/4 ) : aireij3 (i,j)
305 c ( i - 1/4, j - 1/4 ) : aireij4 (i,j)
306 c
307 c ,
308 c Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X,Y).
309 c Chaque aire centree en 1 point scalaire P(i,j) est egale a la somme
310 c des 4 aires aireij1,aireij2,aireij3,aireij4 qui sont affectees au
311 c point (i,j) .
312 c On definit en outre les coefficients alpha comme etant egaux a
313 c (aireij / aire), c.a.d par exp. alpha1(i,j)=aireij1(i,j)/aire(i,j)
314 c
315 c De meme, toute aire centree en 1 point U est egale a la somme des
316 c 4 aires aireij1,aireij2,aireij3,aireij4 entourant le point U .
317 c Idem pour airev, airez .
318 c
319 c On a ,pour chaque maille : dX = dY = 1
320 c
321 c
322 c . V
323 c
324 c aireij4 . . aireij1
325 c
326 c U . . P . U
327 c
328 c aireij3 . . aireij2
329 c
330 c . V
331 c
332 c
333 c
334 c
335 c
336 c ....................................................................
337 c
338 c Calcul des 4 aires elementaires aireij1,aireij2,aireij3,aireij4
339 c qui entourent chaque aire(i,j) , ainsi que les 4 elongations elemen
340 c taires cuij et les 4 elongat. cvij qui sont calculees aux memes
341 c endroits que les aireij .
342 c
343 c ....................................................................
344 c
345 c ....... do 35 : boucle sur les jjm + 1 latitudes .....
346 c
347 c
348
2/2
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 1 times.
 34 DO 35 j = 1, jjp1
349 c
350
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 32 times.
 33 IF ( j. eq. 1 ) THEN
351 c
352 1 yprm = yprimu1(j)
353 1 rlatm = rlatu1(j)
354 c
355 1 coslatm = COS( rlatm )
356 1 radclatm = 0.5* rad * coslatm
357 c
358
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 times.
 33 DO 30 i = 1, iim
359 32 xprp = xprimp025( i )
360 32 xprm = xprimm025( i )
361 32 aireij2( i,1 ) = un4rad2 * coslatm * xprp * yprm
362 32 aireij3( i,1 ) = un4rad2 * coslatm * xprm * yprm
363 32 cuij2 ( i,1 ) = radclatm * xprp
364 32 cuij3 ( i,1 ) = radclatm * xprm
365 32 cvij2 ( i,1 ) = 0.5* rad * yprm
366 32 cvij3 ( i,1 ) = cvij2(i,1)
367 1 30 CONTINUE
368 c
369
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 times.
 33 DO i = 1, iim
370 32 aireij1( i,1 ) = 0.
371 32 aireij4( i,1 ) = 0.
372 32 cuij1 ( i,1 ) = 0.
373 32 cuij4 ( i,1 ) = 0.
374 32 cvij1 ( i,1 ) = 0.
375 33 cvij4 ( i,1 ) = 0.
376 ENDDO
377 c
378 END IF
379 c
380
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 32 times.
 33 IF ( j. eq. jjp1 ) THEN
381 1 yprp = yprimu2(j-1)
382 1 rlatp = rlatu2 (j-1)
383 ccc yprp = fyprim( REAL(j) - 0.25 )
384 ccc rlatp = fy ( REAL(j) - 0.25 )
385 c
386 1 coslatp = COS( rlatp )
387 1 radclatp = 0.5* rad * coslatp
388 c
389
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 times.
 33 DO 31 i = 1,iim
390 32 xprp = xprimp025( i )
391 32 xprm = xprimm025( i )
392 32 aireij1( i,jjp1 ) = un4rad2 * coslatp * xprp * yprp
393 32 aireij4( i,jjp1 ) = un4rad2 * coslatp * xprm * yprp
394 32 cuij1(i,jjp1) = radclatp * xprp
395 32 cuij4(i,jjp1) = radclatp * xprm
396 32 cvij1(i,jjp1) = 0.5 * rad* yprp
397 32 cvij4(i,jjp1) = cvij1(i,jjp1)
398 1 31 CONTINUE
399 c
400
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 times.
 33 DO i = 1, iim
401 32 aireij2( i,jjp1 ) = 0.
402 32 aireij3( i,jjp1 ) = 0.
403 32 cvij2 ( i,jjp1 ) = 0.
404 32 cvij3 ( i,jjp1 ) = 0.
405 32 cuij2 ( i,jjp1 ) = 0.
406 33 cuij3 ( i,jjp1 ) = 0.
407 ENDDO
408 c
409 END IF
410 c
411
412
2/2
✓ Branch 0 taken 31 times.
✓ Branch 1 taken 2 times.
 33 IF ( j .gt. 1 .AND. j .lt. jjp1 ) THEN
413 c
414 31 rlatp = rlatu2 ( j-1 )
415 31 yprp = yprimu2( j-1 )
416 31 rlatm = rlatu1 ( j )
417 31 yprm = yprimu1( j )
418 cc rlatp = fy ( REAL(j) - 0.25 )
419 cc yprp = fyprim( REAL(j) - 0.25 )
420 cc rlatm = fy ( REAL(j) + 0.25 )
421 cc yprm = fyprim( REAL(j) + 0.25 )
422
423 31 coslatm = COS( rlatm )
424 31 coslatp = COS( rlatp )
425 31 radclatp = 0.5* rad * coslatp
426 31 radclatm = 0.5* rad * coslatm
427 c
428 31 ai14 = un4rad2 * coslatp * yprp
429 31 ai23 = un4rad2 * coslatm * yprm
430
2/2
✓ Branch 0 taken 992 times.
✓ Branch 1 taken 31 times.
 1023 DO 32 i = 1,iim
431 992 xprp = xprimp025( i )
432 992 xprm = xprimm025( i )
433
434 992 aireij1 ( i,j ) = ai14 * xprp
435 992 aireij2 ( i,j ) = ai23 * xprp
436 992 aireij3 ( i,j ) = ai23 * xprm
437 992 aireij4 ( i,j ) = ai14 * xprm
438 992 cuij1 ( i,j ) = radclatp * xprp
439 992 cuij2 ( i,j ) = radclatm * xprp
440 992 cuij3 ( i,j ) = radclatm * xprm
441 992 cuij4 ( i,j ) = radclatp * xprm
442 992 cvij1 ( i,j ) = 0.5* rad * yprp
443 992 cvij2 ( i,j ) = 0.5* rad * yprm
444 992 cvij3 ( i,j ) = cvij2(i,j)
445 992 cvij4 ( i,j ) = cvij1(i,j)
446 31 32 CONTINUE
447 c
448 END IF
449 c
450 c ........ periodicite ............
451 c
452 33 cvij1 (iip1,j) = cvij1 (1,j)
453 33 cvij2 (iip1,j) = cvij2 (1,j)
454 33 cvij3 (iip1,j) = cvij3 (1,j)
455 33 cvij4 (iip1,j) = cvij4 (1,j)
456 33 cuij1 (iip1,j) = cuij1 (1,j)
457 33 cuij2 (iip1,j) = cuij2 (1,j)
458 33 cuij3 (iip1,j) = cuij3 (1,j)
459 33 cuij4 (iip1,j) = cuij4 (1,j)
460 33 aireij1 (iip1,j) = aireij1 (1,j )
461 33 aireij2 (iip1,j) = aireij2 (1,j )
462 33 aireij3 (iip1,j) = aireij3 (1,j )
463 33 aireij4 (iip1,j) = aireij4 (1,j )
464
465 1 35 CONTINUE
466 c
467 c ..............................................................
468 c
469
2/2
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 1 times.
 34 DO 37 j = 1, jjp1
470
2/2
✓ Branch 0 taken 1056 times.
✓ Branch 1 taken 33 times.
 1089 DO 36 i = 1, iim
471 aire ( i,j ) = aireij1(i,j) + aireij2(i,j) + aireij3(i,j) +
472 1056 * aireij4(i,j)
473 1056 alpha1 ( i,j ) = aireij1(i,j) / aire(i,j)
474 1056 alpha2 ( i,j ) = aireij2(i,j) / aire(i,j)
475 1056 alpha3 ( i,j ) = aireij3(i,j) / aire(i,j)
476 1056 alpha4 ( i,j ) = aireij4(i,j) / aire(i,j)
477 1056 alpha1p2( i,j ) = alpha1 (i,j) + alpha2 (i,j)
478 1056 alpha1p4( i,j ) = alpha1 (i,j) + alpha4 (i,j)
479 1056 alpha2p3( i,j ) = alpha2 (i,j) + alpha3 (i,j)
480 1056 alpha3p4( i,j ) = alpha3 (i,j) + alpha4 (i,j)
481 33 36 CONTINUE
482 c
483 c
484 33 aire (iip1,j) = aire (1,j)
485 33 alpha1 (iip1,j) = alpha1 (1,j)
486 33 alpha2 (iip1,j) = alpha2 (1,j)
487 33 alpha3 (iip1,j) = alpha3 (1,j)
488 33 alpha4 (iip1,j) = alpha4 (1,j)
489 33 alpha1p2(iip1,j) = alpha1p2(1,j)
490 33 alpha1p4(iip1,j) = alpha1p4(1,j)
491 33 alpha2p3(iip1,j) = alpha2p3(1,j)
492 33 alpha3p4(iip1,j) = alpha3p4(1,j)
493 1 37 CONTINUE
494 c
495
496
2/2
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 1 times.
 34 DO 42 j = 1,jjp1
497
2/2
✓ Branch 0 taken 1056 times.
✓ Branch 1 taken 33 times.
 1089 DO 41 i = 1,iim
498 aireu (i,j)= aireij1(i,j) + aireij2(i,j) + aireij4(i+1,j) +
499 1056 * aireij3(i+1,j)
500 1056 unsaire ( i,j)= 1./ aire(i,j)
501 1056 unsair_gam1( i,j)= unsaire(i,j)** ( - gamdi_gdiv )
502 1056 unsair_gam2( i,j)= unsaire(i,j)** ( - gamdi_h )
503 1056 airesurg ( i,j)= aire(i,j)/ g
504 33 41 CONTINUE
505 33 aireu (iip1,j) = aireu (1,j)
506 33 unsaire (iip1,j) = unsaire(1,j)
507 33 unsair_gam1(iip1,j) = unsair_gam1(1,j)
508 33 unsair_gam2(iip1,j) = unsair_gam2(1,j)
509 33 airesurg (iip1,j) = airesurg(1,j)
510 1 42 CONTINUE
511 c
512 c
513
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 times.
 33 DO 48 j = 1,jjm
514 c
515
2/2
✓ Branch 0 taken 1024 times.
✓ Branch 1 taken 32 times.
 1056 DO i=1,iim
516 airev (i,j) = aireij2(i,j)+ aireij3(i,j)+ aireij1(i,j+1) +
517 1056 * aireij4(i,j+1)
518 ENDDO
519
2/2
✓ Branch 0 taken 1024 times.
✓ Branch 1 taken 32 times.
 1056 DO i=1,iim
520 airez = aireij2(i,j)+aireij1(i,j+1)+aireij3(i+1,j) +
521 1024 * aireij4(i+1,j+1)
522 1024 unsairez(i,j) = 1./ airez
523 1024 unsairz_gam(i,j)= unsairez(i,j)** ( - gamdi_grot )
524 1056 fext (i,j) = airez * SIN(rlatv(j))* 2.* omeg
525 ENDDO
526 32 airev (iip1,j) = airev(1,j)
527 32 unsairez (iip1,j) = unsairez(1,j)
528 32 fext (iip1,j) = fext(1,j)
529 32 unsairz_gam(iip1,j) = unsairz_gam(1,j)
530 c
531 1 48 CONTINUE
532 c
533 c
534 c ..... Calcul des elongations cu,cv, cvu .........
535 c
536
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 times.
 33 DO j = 1, jjm
537
2/2
✓ Branch 0 taken 1024 times.
✓ Branch 1 taken 32 times.
 1056 DO i = 1, iim
538 1024 cv(i,j) = 0.5 *( cvij2(i,j)+cvij3(i,j)+cvij1(i,j+1)+cvij4(i,j+1))
539 1024 cvu(i,j)= 0.5 *( cvij1(i,j)+cvij4(i,j)+cvij2(i,j) +cvij3(i,j) )
540 1024 cuv(i,j)= 0.5 *( cuij2(i,j)+cuij3(i,j)+cuij1(i,j+1)+cuij4(i,j+1))
541 1056 unscv2(i,j) = 1./ ( cv(i,j)*cv(i,j) )
542 ENDDO
543
2/2
✓ Branch 0 taken 1024 times.
✓ Branch 1 taken 32 times.
 1056 DO i = 1, iim
544 1024 cuvsurcv (i,j) = airev(i,j) * unscv2(i,j)
545 1024 cvsurcuv (i,j) = 1./cuvsurcv(i,j)
546 1024 cuvscvgam1(i,j) = cuvsurcv (i,j) ** ( - gamdi_gdiv )
547 1024 cuvscvgam2(i,j) = cuvsurcv (i,j) ** ( - gamdi_h )
548 1056 cvscuvgam(i,j) = cvsurcuv (i,j) ** ( - gamdi_grot )
549 ENDDO
550 32 cv (iip1,j) = cv (1,j)
551 32 cvu (iip1,j) = cvu (1,j)
552 32 unscv2 (iip1,j) = unscv2 (1,j)
553 32 cuv (iip1,j) = cuv (1,j)
554 32 cuvsurcv (iip1,j) = cuvsurcv (1,j)
555 32 cvsurcuv (iip1,j) = cvsurcuv (1,j)
556 32 cuvscvgam1(iip1,j) = cuvscvgam1(1,j)
557 32 cuvscvgam2(iip1,j) = cuvscvgam2(1,j)
558 33 cvscuvgam(iip1,j) = cvscuvgam(1,j)
559 ENDDO
560
561
2/2
✓ Branch 0 taken 31 times.
✓ Branch 1 taken 1 times.
 32 DO j = 2, jjm
562
2/2
✓ Branch 0 taken 992 times.
✓ Branch 1 taken 31 times.
 1023 DO i = 1, iim
563 992 cu(i,j) = 0.5*(cuij1(i,j)+cuij4(i+1,j)+cuij2(i,j)+cuij3(i+1,j))
564 992 unscu2 (i,j) = 1./ ( cu(i,j) * cu(i,j) )
565 992 cvusurcu (i,j) = aireu(i,j) * unscu2(i,j)
566 992 cusurcvu (i,j) = 1./ cvusurcu(i,j)
567 992 cvuscugam1 (i,j) = cvusurcu(i,j) ** ( - gamdi_gdiv )
568 992 cvuscugam2 (i,j) = cvusurcu(i,j) ** ( - gamdi_h )
569 1023 cuscvugam (i,j) = cusurcvu(i,j) ** ( - gamdi_grot )
570 ENDDO
571 31 cu (iip1,j) = cu(1,j)
572 31 unscu2 (iip1,j) = unscu2(1,j)
573 31 cvusurcu (iip1,j) = cvusurcu(1,j)
574 31 cusurcvu (iip1,j) = cusurcvu(1,j)
575 31 cvuscugam1(iip1,j) = cvuscugam1(1,j)
576 31 cvuscugam2(iip1,j) = cvuscugam2(1,j)
577 32 cuscvugam (iip1,j) = cuscvugam(1,j)
578 ENDDO
579
580 c
581 c .... calcul aux poles ....
582 c
583
2/2
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 1 times.
 34 DO i = 1, iip1
584 33 cu ( i, 1 ) = 0.
585 33 unscu2( i, 1 ) = 0.
586 cvu ( i, 1 ) = 0.
587 c
588 33 cu (i, jjp1) = 0.
589 33 unscu2(i, jjp1) = 0.
590 1 cvu (i, jjp1) = 0.
591 ENDDO
592 c
593 c ..............................................................
594 c
595
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 times.
 33 DO j = 1, jjm
596
2/2
✓ Branch 0 taken 1024 times.
✓ Branch 1 taken 32 times.
 1056 DO i= 1, iim
597 1024 airvscu2 (i,j) = airev(i,j)/ ( cuv(i,j) * cuv(i,j) )
598 1056 aivscu2gam(i,j) = airvscu2(i,j)** ( - gamdi_grot )
599 ENDDO
600 32 airvscu2 (iip1,j) = airvscu2(1,j)
601 33 aivscu2gam(iip1,j) = aivscu2gam(1,j)
602 ENDDO
603
604
2/2
✓ Branch 0 taken 31 times.
✓ Branch 1 taken 1 times.
 32 DO j=2,jjm
605
2/2
✓ Branch 0 taken 992 times.
✓ Branch 1 taken 31 times.
 1023 DO i=1,iim
606 992 airuscv2 (i,j) = aireu(i,j)/ ( cvu(i,j) * cvu(i,j) )
607 1023 aiuscv2gam (i,j) = airuscv2(i,j)** ( - gamdi_grot )
608 ENDDO
609 31 airuscv2 (iip1,j) = airuscv2 (1,j)
610 32 aiuscv2gam(iip1,j) = aiuscv2gam(1,j)
611 ENDDO
612
613 c
614 c calcul des aires aux poles :
615 c -----------------------------
616 c
617 1 apoln = SSUM(iim,aire(1,1),1)
618 1 apols = SSUM(iim,aire(1,jjp1),1)
619 1 unsapolnga1 = 1./ ( apoln ** ( - gamdi_gdiv ) )
620 1 unsapolsga1 = 1./ ( apols ** ( - gamdi_gdiv ) )
621 1 unsapolnga2 = 1./ ( apoln ** ( - gamdi_h ) )
622 1 unsapolsga2 = 1./ ( apols ** ( - gamdi_h ) )
623 c
624 c-----------------------------------------------------------------------
625 c gtitre='Coriolis version ancienne'
626 c gfichier='fext1'
627 c CALL writestd(fext,iip1*jjm)
628 c
629 c changement F. Hourdin calcul conservatif pour fext
630 c constang contient le produit a * cos ( latitude ) * omega
631 c
632
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 times.
 33 DO i=1,iim
633 33 constang(i,1) = 0.
634 ENDDO
635
2/2
✓ Branch 0 taken 31 times.
✓ Branch 1 taken 1 times.
 32 DO j=1,jjm-1
636
2/2
✓ Branch 0 taken 992 times.
✓ Branch 1 taken 31 times.
 1024 DO i=1,iim
637 1023 constang(i,j+1) = rad*omeg*cu(i,j+1)*COS(rlatu(j+1))
638 ENDDO
639 ENDDO
640
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 times.
 33 DO i=1,iim
641 33 constang(i,jjp1) = 0.
642 ENDDO
643 c
644 c periodicite en longitude
645 c
646
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 times.
 33 DO j=1,jjm
647 33 fext(iip1,j) = fext(1,j)
648 ENDDO
649
2/2
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 1 times.
 34 DO j=1,jjp1
650 34 constang(iip1,j) = constang(1,j)
651 ENDDO
652
653 c fin du changement
654
655 c
656 c-----------------------------------------------------------------------
657 c
658 1 WRITE(6,*) ' *** Coordonnees de la grille *** '
659 1 WRITE(6,995)
660 c
661 1 WRITE(6,*) ' LONGITUDES aux pts. V ( degres ) '
662 1 WRITE(6,995)
663
2/2
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 1 times.
 34 DO i=1,iip1
664 34 rlonvv(i) = rlonv(i)*180./pi
665 ENDDO
666 1 WRITE(6,400) rlonvv
667 c
668 1 WRITE(6,995)
669 1 WRITE(6,*) ' LATITUDES aux pts. V ( degres ) '
670 1 WRITE(6,995)
671
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 times.
 33 DO i=1,jjm
672 33 rlatuu(i)=rlatv(i)*180./pi
673 ENDDO
674 1 WRITE(6,400) (rlatuu(i),i=1,jjm)
675 c
676
2/2
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 1 times.
 34 DO i=1,iip1
677 34 rlonvv(i)=rlonu(i)*180./pi
678 ENDDO
679 1 WRITE(6,995)
680 1 WRITE(6,*) ' LONGITUDES aux pts. U ( degres ) '
681 1 WRITE(6,995)
682 1 WRITE(6,400) rlonvv
683 1 WRITE(6,995)
684
685 1 WRITE(6,*) ' LATITUDES aux pts. U ( degres ) '
686 1 WRITE(6,995)
687
2/2
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 1 times.
 34 DO i=1,jjp1
688 34 rlatuu(i)=rlatu(i)*180./pi
689 ENDDO
690 1 WRITE(6,400) (rlatuu(i),i=1,jjp1)
691 1 WRITE(6,995)
692 c
693 444 format(f10.3,f6.0)
694 400 FORMAT(1x,8f8.2)
695 990 FORMAT(//)
696 995 FORMAT(/)
697 c
698 1 RETURN
699 END
700