My Project
 All Classes Files Functions Variables Macros
inigeom.F
Go to the documentation of this file.
1 !
2 ! $Id: inigeom.F 1403 2010-07-01 09:02:53Z fairhead $
3 !
4 c
5 c
6  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  IMPLICIT NONE
19 c
20 #include "dimensions.h"
21 #include "paramet.h"
22 #include "comconst.h"
23 #include "comgeom2.h"
24 #include "serre.h"
25 #include "logic.h"
26 #include "comdissnew.h"
27 
28 c-----------------------------------------------------------------------
29 c .... Variables locales ....
30 c
31  INTEGER i,j,itmax,itmay,iter
32  REAL cvu(iip1,jjp1),cuv(iip1,jjm)
33  REAL ai14,ai23,airez,rlatp,rlatm,xprm,xprp,un4rad2,yprp,yprm
34  REAL eps,x1,xo1,f,df,xdm,y1,yo1,ydm
35  REAL coslatm,coslatp,radclatm,radclatp
36  REAL cuij1(iip1,jjp1),cuij2(iip1,jjp1),cuij3(iip1,jjp1),
37  * cuij4(iip1,jjp1)
38  REAL cvij1(iip1,jjp1),cvij2(iip1,jjp1),cvij3(iip1,jjp1),
39  * cvij4(iip1,jjp1)
40  REAL rlonvv(iip1),rlatuu(jjp1)
41  REAL rlatu1(jjm),yprimu1(jjm),rlatu2(jjm),yprimu2(jjm) ,
42  * yprimv(jjm),yprimu(jjp1)
43  REAL gamdi_gdiv, gamdi_grot, gamdi_h
44 
45  REAL rlonm025(iip1),xprimm025(iip1), rlonp025(iip1),
46  , xprimp025(iip1)
47  SAVE rlatu1,yprimu1,rlatu2,yprimu2,yprimv,yprimu
48  SAVE rlonm025,xprimm025,rlonp025,xprimp025
49 
50  REAL ssum
51 c
52 c
162  WRITE(6,3)
163  3 FORMAT( // 10x,' .... INIGEOM date du 01/06/98 ..... ',
164  * //5x,' Calcul des elongations cu et cv comme sommes des 4 ' /
165  * 5x,
166 ' elong. cuij1, .. 4 , cvij1,.. 4 qui les entourent , aux * '/ 5x,' memes endroits que les aires aireij1,...j4 . ' / )
167 c
168 c
169  IF( nitergdiv.NE.2 ) THEN
170  gamdi_gdiv = coefdis/ ( REAL(nitergdiv) -2. )
171  ELSE
172  gamdi_gdiv = 0.
173  ENDIF
174  IF( nitergrot.NE.2 ) THEN
175  gamdi_grot = coefdis/ ( REAL(nitergrot) -2. )
176  ELSE
177  gamdi_grot = 0.
178  ENDIF
179  IF( niterh.NE.2 ) THEN
180  gamdi_h = coefdis/ ( REAL(niterh) -2. )
181  ELSE
182  gamdi_h = 0.
183  ENDIF
184 
185  WRITE(6,*) ' gamdi_gd ',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis,
186  * nitergdiv,nitergrot,niterh
187 c
188  pi = 2.* asin(1.)
189 c
190  WRITE(6,990)
191 
192 c ----------------------------------------------------------------
193 c
194  IF( .NOT.fxyhypb ) THEN
195 c
196 c
197  IF( ysinus ) THEN
198 c
199  WRITE(6,*) ' *** Inigeom , Y = Sinus ( Latitude ) *** '
200 c
201 c .... utilisation de f(x,y ) avec y = sinus de la latitude .....
202 
203  CALL fxysinus(rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,
204  , rlatu2,yprimu2,
205  , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)
206 
207  ELSE
208 c
209  WRITE(6,*) '*** Inigeom , Y = Latitude , der. sinusoid . ***'
210 
211 c .... utilisation de f(x,y) a tangente sinusoidale , y etant la latit. ...
212 c
213 
214  pxo = clon *pi /180.
215  pyo = 2.* clat* pi /180.
216 c
217 c .... determination de transx ( pour le zoom ) par Newton-Raphson ...
218 c
219  itmax = 10
220  eps = .1e-7
221 c
222  xo1 = 0.
223  DO 10 iter = 1, itmax
224  x1 = xo1
225  f = x1+ alphax *sin(x1-pxo)
226  df = 1.+ alphax *cos(x1-pxo)
227  x1 = x1 - f/df
228  xdm = abs( x1- xo1 )
229  IF( xdm.LE.eps )go to 11
230  xo1 = x1
231  10 CONTINUE
232  11 CONTINUE
233 c
234  transx = xo1
235 
236  itmay = 10
237  eps = .1e-7
238 C
239  yo1 = 0.
240  DO 15 iter = 1,itmay
241  y1 = yo1
242  f = y1 + alphay* sin(y1-pyo)
243  df = 1. + alphay* cos(y1-pyo)
244  y1 = y1 -f/df
245  ydm = abs(y1-yo1)
246  IF(ydm.LE.eps) go to 17
247  yo1 = y1
248  15 CONTINUE
249 c
250  17 CONTINUE
251  transy = yo1
252 
253  CALL fxy( rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,
254  , rlatu2,yprimu2,
255  , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)
256 
257  ENDIF
258 c
259  ELSE
260 c
261 c .... Utilisation de fxyhyper , f(x,y) a derivee tangente hyperbol.
262 c .....................................................................
263 
264  WRITE(6,*)'*** Inigeom , Y = Latitude , der.tg. hyperbolique ***'
265 
266  CALL fxyhyper( clat, grossismy, dzoomy, tauy ,
267  , clon, grossismx, dzoomx, taux ,
268  , rlatu,yprimu,rlatv, yprimv,rlatu1, yprimu1,rlatu2,yprimu2 ,
269  , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025 )
270 
271 
272  ENDIF
273 c
274 c -------------------------------------------------------------------
275 
276 c
277  rlatu(1) = asin(1.)
278  rlatu(jjp1) = - rlatu(1)
279 c
280 c
281 c .... calcul aux poles ....
282 c
283  yprimu(1) = 0.
284  yprimu(jjp1) = 0.
285 c
286 c
287  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  DO 35 j = 1, jjp1
349 c
350  IF ( j. eq. 1 ) THEN
351 c
352  yprm = yprimu1(j)
353  rlatm = rlatu1(j)
354 c
355  coslatm = cos( rlatm )
356  radclatm = 0.5* rad * coslatm
357 c
358  DO 30 i = 1, iim
359  xprp = xprimp025( i )
360  xprm = xprimm025( i )
361  aireij2( i,1 ) = un4rad2 * coslatm * xprp * yprm
362  aireij3( i,1 ) = un4rad2 * coslatm * xprm * yprm
363  cuij2( i,1 ) = radclatm * xprp
364  cuij3( i,1 ) = radclatm * xprm
365  cvij2( i,1 ) = 0.5* rad * yprm
366  cvij3( i,1 ) = cvij2(i,1)
367  30 CONTINUE
368 c
369  DO i = 1, iim
370  aireij1( i,1 ) = 0.
371  aireij4( i,1 ) = 0.
372  cuij1( i,1 ) = 0.
373  cuij4( i,1 ) = 0.
374  cvij1( i,1 ) = 0.
375  cvij4( i,1 ) = 0.
376  ENDDO
377 c
378  END IF
379 c
380  IF ( j. eq. jjp1 ) THEN
381  yprp = yprimu2(j-1)
382  rlatp = rlatu2(j-1)
383 ccc yprp = fyprim( REAL(j) - 0.25 )
384 ccc rlatp = fy ( REAL(j) - 0.25 )
385 c
386  coslatp = cos( rlatp )
387  radclatp = 0.5* rad * coslatp
388 c
389  DO 31 i = 1,iim
390  xprp = xprimp025( i )
391  xprm = xprimm025( i )
392  aireij1( i,jjp1 ) = un4rad2 * coslatp * xprp * yprp
393  aireij4( i,jjp1 ) = un4rad2 * coslatp * xprm * yprp
394  cuij1(i,jjp1) = radclatp * xprp
395  cuij4(i,jjp1) = radclatp * xprm
396  cvij1(i,jjp1) = 0.5 * rad* yprp
397  cvij4(i,jjp1) = cvij1(i,jjp1)
398  31 CONTINUE
399 c
400  DO i = 1, iim
401  aireij2( i,jjp1 ) = 0.
402  aireij3( i,jjp1 ) = 0.
403  cvij2( i,jjp1 ) = 0.
404  cvij3( i,jjp1 ) = 0.
405  cuij2( i,jjp1 ) = 0.
406  cuij3( i,jjp1 ) = 0.
407  ENDDO
408 c
409  END IF
410 c
411 
412  IF ( j .gt. 1 .AND. j .lt. jjp1 ) THEN
413 c
414  rlatp = rlatu2( j-1 )
415  yprp = yprimu2( j-1 )
416  rlatm = rlatu1( j )
417  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  coslatm = cos( rlatm )
424  coslatp = cos( rlatp )
425  radclatp = 0.5* rad * coslatp
426  radclatm = 0.5* rad * coslatm
427 c
428  DO 32 i = 1,iim
429  xprp = xprimp025( i )
430  xprm = xprimm025( i )
431 
432  ai14 = un4rad2 * coslatp * yprp
433  ai23 = un4rad2 * coslatm * yprm
434  aireij1( i,j ) = ai14 * xprp
435  aireij2( i,j ) = ai23 * xprp
436  aireij3( i,j ) = ai23 * xprm
437  aireij4( i,j ) = ai14 * xprm
438  cuij1( i,j ) = radclatp * xprp
439  cuij2( i,j ) = radclatm * xprp
440  cuij3( i,j ) = radclatm * xprm
441  cuij4( i,j ) = radclatp * xprm
442  cvij1( i,j ) = 0.5* rad * yprp
443  cvij2( i,j ) = 0.5* rad * yprm
444  cvij3( i,j ) = cvij2(i,j)
445  cvij4( i,j ) = cvij1(i,j)
446  32 CONTINUE
447 c
448  END IF
449 c
450 c ........ periodicite ............
451 c
452  cvij1(iip1,j) = cvij1(1,j)
453  cvij2(iip1,j) = cvij2(1,j)
454  cvij3(iip1,j) = cvij3(1,j)
455  cvij4(iip1,j) = cvij4(1,j)
456  cuij1(iip1,j) = cuij1(1,j)
457  cuij2(iip1,j) = cuij2(1,j)
458  cuij3(iip1,j) = cuij3(1,j)
459  cuij4(iip1,j) = cuij4(1,j)
460  aireij1(iip1,j) = aireij1(1,j )
461  aireij2(iip1,j) = aireij2(1,j )
462  aireij3(iip1,j) = aireij3(1,j )
463  aireij4(iip1,j) = aireij4(1,j )
464 
465  35 CONTINUE
466 c
467 c ..............................................................
468 c
469  DO 37 j = 1, jjp1
470  DO 36 i = 1, iim
471  aire( i,j ) = aireij1(i,j) + aireij2(i,j) + aireij3(i,j) +
472  * aireij4(i,j)
473  alpha1( i,j ) = aireij1(i,j) / aire(i,j)
474  alpha2( i,j ) = aireij2(i,j) / aire(i,j)
475  alpha3( i,j ) = aireij3(i,j) / aire(i,j)
476  alpha4( i,j ) = aireij4(i,j) / aire(i,j)
477  alpha1p2( i,j ) = alpha1(i,j) + alpha2(i,j)
478  alpha1p4( i,j ) = alpha1(i,j) + alpha4(i,j)
479  alpha2p3( i,j ) = alpha2(i,j) + alpha3(i,j)
480  alpha3p4( i,j ) = alpha3(i,j) + alpha4(i,j)
481  36 CONTINUE
482 c
483 c
484  aire(iip1,j) = aire(1,j)
485  alpha1(iip1,j) = alpha1(1,j)
486  alpha2(iip1,j) = alpha2(1,j)
487  alpha3(iip1,j) = alpha3(1,j)
488  alpha4(iip1,j) = alpha4(1,j)
489  alpha1p2(iip1,j) = alpha1p2(1,j)
490  alpha1p4(iip1,j) = alpha1p4(1,j)
491  alpha2p3(iip1,j) = alpha2p3(1,j)
492  alpha3p4(iip1,j) = alpha3p4(1,j)
493  37 CONTINUE
494 c
495 
496  DO 42 j = 1,jjp1
497  DO 41 i = 1,iim
498  aireu(i,j)= aireij1(i,j) + aireij2(i,j) + aireij4(i+1,j) +
499  * aireij3(i+1,j)
500  unsaire( i,j)= 1./ aire(i,j)
501  unsair_gam1( i,j)= unsaire(i,j)** ( - gamdi_gdiv )
502  unsair_gam2( i,j)= unsaire(i,j)** ( - gamdi_h )
503  airesurg( i,j)= aire(i,j)/ g
504  41 CONTINUE
505  aireu(iip1,j) = aireu(1,j)
506  unsaire(iip1,j) = unsaire(1,j)
507  unsair_gam1(iip1,j) = unsair_gam1(1,j)
508  unsair_gam2(iip1,j) = unsair_gam2(1,j)
509  airesurg(iip1,j) = airesurg(1,j)
510  42 CONTINUE
511 c
512 c
513  DO 48 j = 1,jjm
514 c
515  DO i=1,iim
516  airev(i,j) = aireij2(i,j)+ aireij3(i,j)+ aireij1(i,j+1) +
517  * aireij4(i,j+1)
518  ENDDO
519  DO i=1,iim
520  airez = aireij2(i,j)+aireij1(i,j+1)+aireij3(i+1,j) +
521  * aireij4(i+1,j+1)
522  unsairez(i,j) = 1./ airez
523  unsairz_gam(i,j)= unsairez(i,j)** ( - gamdi_grot )
524  fext(i,j) = airez * sin(rlatv(j))* 2.* omeg
525  ENDDO
526  airev(iip1,j) = airev(1,j)
527  unsairez(iip1,j) = unsairez(1,j)
528  fext(iip1,j) = fext(1,j)
529  unsairz_gam(iip1,j) = unsairz_gam(1,j)
530 c
531  48 CONTINUE
532 c
533 c
534 c ..... Calcul des elongations cu,cv, cvu .........
535 c
536  DO j = 1, jjm
537  DO i = 1, iim
538  cv(i,j) = 0.5 *( cvij2(i,j)+cvij3(i,j)+cvij1(i,j+1)+cvij4(i,j+1))
539  cvu(i,j)= 0.5 *( cvij1(i,j)+cvij4(i,j)+cvij2(i,j) +cvij3(i,j) )
540  cuv(i,j)= 0.5 *( cuij2(i,j)+cuij3(i,j)+cuij1(i,j+1)+cuij4(i,j+1))
541  unscv2(i,j) = 1./ ( cv(i,j)*cv(i,j) )
542  ENDDO
543  DO i = 1, iim
544  cuvsurcv(i,j) = airev(i,j) * unscv2(i,j)
545  cvsurcuv(i,j) = 1./cuvsurcv(i,j)
546  cuvscvgam1(i,j) = cuvsurcv(i,j) ** ( - gamdi_gdiv )
547  cuvscvgam2(i,j) = cuvsurcv(i,j) ** ( - gamdi_h )
548  cvscuvgam(i,j) = cvsurcuv(i,j) ** ( - gamdi_grot )
549  ENDDO
550  cv(iip1,j) = cv(1,j)
551  cvu(iip1,j) = cvu(1,j)
552  unscv2(iip1,j) = unscv2(1,j)
553  cuv(iip1,j) = cuv(1,j)
554  cuvsurcv(iip1,j) = cuvsurcv(1,j)
555  cvsurcuv(iip1,j) = cvsurcuv(1,j)
556  cuvscvgam1(iip1,j) = cuvscvgam1(1,j)
557  cuvscvgam2(iip1,j) = cuvscvgam2(1,j)
558  cvscuvgam(iip1,j) = cvscuvgam(1,j)
559  ENDDO
560 
561  DO j = 2, jjm
562  DO i = 1, iim
563  cu(i,j) = 0.5*(cuij1(i,j)+cuij4(i+1,j)+cuij2(i,j)+cuij3(i+1,j))
564  unscu2(i,j) = 1./ ( cu(i,j) * cu(i,j) )
565  cvusurcu(i,j) = aireu(i,j) * unscu2(i,j)
566  cusurcvu(i,j) = 1./ cvusurcu(i,j)
567  cvuscugam1(i,j) = cvusurcu(i,j) ** ( - gamdi_gdiv )
568  cvuscugam2(i,j) = cvusurcu(i,j) ** ( - gamdi_h )
569  cuscvugam(i,j) = cusurcvu(i,j) ** ( - gamdi_grot )
570  ENDDO
571  cu(iip1,j) = cu(1,j)
572  unscu2(iip1,j) = unscu2(1,j)
573  cvusurcu(iip1,j) = cvusurcu(1,j)
574  cusurcvu(iip1,j) = cusurcvu(1,j)
575  cvuscugam1(iip1,j) = cvuscugam1(1,j)
576  cvuscugam2(iip1,j) = cvuscugam2(1,j)
577  cuscvugam(iip1,j) = cuscvugam(1,j)
578  ENDDO
579 
580 c
581 c .... calcul aux poles ....
582 c
583  DO i = 1, iip1
584  cu( i, 1 ) = 0.
585  unscu2( i, 1 ) = 0.
586  cvu( i, 1 ) = 0.
587 c
588  cu(i, jjp1) = 0.
589  unscu2(i, jjp1) = 0.
590  cvu(i, jjp1) = 0.
591  ENDDO
592 c
593 c ..............................................................
594 c
595  DO j = 1, jjm
596  DO i= 1, iim
597  airvscu2(i,j) = airev(i,j)/ ( cuv(i,j) * cuv(i,j) )
598  aivscu2gam(i,j) = airvscu2(i,j)** ( - gamdi_grot )
599  ENDDO
600  airvscu2(iip1,j) = airvscu2(1,j)
601  aivscu2gam(iip1,j) = aivscu2gam(1,j)
602  ENDDO
603 
604  DO j=2,jjm
605  DO i=1,iim
606  airuscv2(i,j) = aireu(i,j)/ ( cvu(i,j) * cvu(i,j) )
607  aiuscv2gam(i,j) = airuscv2(i,j)** ( - gamdi_grot )
608  ENDDO
609  airuscv2(iip1,j) = airuscv2(1,j)
610  aiuscv2gam(iip1,j) = aiuscv2gam(1,j)
611  ENDDO
612 
613 c
614 c calcul des aires aux poles :
615 c -----------------------------
616 c
617  apoln = ssum(iim,aire(1,1),1)
618  apols = ssum(iim,aire(1,jjp1),1)
619  unsapolnga1 = 1./ ( apoln ** ( - gamdi_gdiv ) )
620  unsapolsga1 = 1./ ( apols ** ( - gamdi_gdiv ) )
621  unsapolnga2 = 1./ ( apoln ** ( - gamdi_h ) )
622  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  DO i=1,iim
633  constang(i,1) = 0.
634  ENDDO
635  DO j=1,jjm-1
636  DO i=1,iim
637  constang(i,j+1) = rad*omeg*cu(i,j+1)*cos(rlatu(j+1))
638  ENDDO
639  ENDDO
640  DO i=1,iim
641  constang(i,jjp1) = 0.
642  ENDDO
643 c
644 c periodicite en longitude
645 c
646  DO j=1,jjm
647  fext(iip1,j) = fext(1,j)
648  ENDDO
649  DO j=1,jjp1
650  constang(iip1,j) = constang(1,j)
651  ENDDO
652 
653 c fin du changement
654 
655 c
656 c-----------------------------------------------------------------------
657 c
658  WRITE(6,*) ' *** Coordonnees de la grille *** '
659  WRITE(6,995)
660 c
661  WRITE(6,*) ' LONGITUDES aux pts. V ( degres ) '
662  WRITE(6,995)
663  DO i=1,iip1
664  rlonvv(i) = rlonv(i)*180./pi
665  ENDDO
666  WRITE(6,400) rlonvv
667 c
668  WRITE(6,995)
669  WRITE(6,*) ' LATITUDES aux pts. V ( degres ) '
670  WRITE(6,995)
671  DO i=1,jjm
672  rlatuu(i)=rlatv(i)*180./pi
673  ENDDO
674  WRITE(6,400) (rlatuu(i),i=1,jjm)
675 c
676  DO i=1,iip1
677  rlonvv(i)=rlonu(i)*180./pi
678  ENDDO
679  WRITE(6,995)
680  WRITE(6,*) ' LONGITUDES aux pts. U ( degres ) '
681  WRITE(6,995)
682  WRITE(6,400) rlonvv
683  WRITE(6,995)
684 
685  WRITE(6,*) ' LATITUDES aux pts. U ( degres ) '
686  WRITE(6,995)
687  DO i=1,jjp1
688  rlatuu(i)=rlatu(i)*180./pi
689  ENDDO
690  WRITE(6,400) (rlatuu(i),i=1,jjp1)
691  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  RETURN
699  END