libf/dyn3d/inigeom.f90

module inigeom_m

  IMPLICIT NONE

contains

  SUBROUTINE inigeom

    ! Auteur : P. Le Van
    ! Version du 01/04/2001

    ! Calcul des élongations cuij1, ..., cuij4, cvij1, ..., cvij4 aux mêmes
    ! endroits que les aires aireij1_2d, ..., aireij4_2d.

    ! Choix entre une fonction "f(y)" à dérivée sinusoïdale ou à dérivée
    ! tangente hyperbolique
    ! calcul des coefficients (cu_2d, cv_2d, 1./cu_2d**2, 1./cv_2d**2)

    ! les coef. (cu_2d, cv_2d) permettent de passer des vitesses naturelles
    ! aux vitesses covariantes et contravariantes, ou vice-versa

    ! on a :
    ! u (covariant) = cu_2d * u (naturel), u(contrav)= u(nat)/cu_2d
    ! v (covariant) = cv_2d * v (naturel), v(contrav)= v(nat)/cv_2d

    ! on en tire : 
    ! u(covariant) = cu_2d * cu_2d * u(contravariant)
    ! v(covariant) = cv_2d * cv_2d * v(contravariant)

    ! on a l'application (x(X), y(Y)) avec - im/2 +1 <= X <= im/2
    ! et - jm/2 <= Y <= jm/2

    ! x est la longitude du point en radians.
    ! y est la latitude du point en radians.
    ! 
    ! on a : cu_2d(i, j) = rad * cos(y) * dx/dX
    ! cv(j) = rad * dy/dY
    ! aire_2d(i, j) = cu_2d(i, j) * cv(j)
    ! 
    ! y, dx/dX, dy/dY calcules aux points concernes
    ! cv, bien que dependant de j uniquement, sera ici indice aussi en i
    ! pour un adressage plus facile en ij.

    ! aux points u et v, 
    ! xprimu et xprimv sont respectivement les valeurs de dx/dX
    ! yprimu et yprimv sont respectivement les valeurs de dy/dY
    ! rlatu et rlatv sont respectivement les valeurs de la latitude
    ! cvu et cv_2d sont respectivement les valeurs de cv_2d

    ! aux points u, v, scalaires, et z 
    ! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d
    ! Cf. "inigeom.txt".

    USE comconst, ONLY : g, omeg, rad
    USE comgeom, ONLY : airesurg_2d, aireu_2d, airev_2d, aire_2d, &
         alpha1p2_2d, alpha1p4_2d, alpha1_2d, &
         alpha2p3_2d, alpha2_2d, alpha3p4_2d, alpha3_2d, alpha4_2d, apoln, &
         apols, constang_2d, cuscvugam_2d, cusurcvu_2d, cuvscvgam1_2d, &
         cuvscvgam2_2d, cuvsurcv_2d, cu_2d, cvscuvgam_2d, cvsurcuv_2d, &
         cvuscugam1_2d, cvuscugam2_2d, cvusurcu_2d, cv_2d, fext_2d, rlatu, &
         rlatv, rlonu, rlonv, unsairez_2d, unsaire_2d, unsairz_gam_2d, &
         unsair_gam1_2d, unsair_gam2_2d, unsapolnga1, unsapolnga2, &
         unsapolsga1, unsapolsga2, unscu2_2d, unscv2_2d, xprimu, xprimv
    USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
    USE dimens_m, ONLY : iim, jjm
    USE logic, ONLY : fxyhypb, ysinus
    use nr_util, only: pi
    USE paramet_m, ONLY : iip1, jjp1
    USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &
         grossismy, pxo, pyo, taux, tauy, transx, transy

    ! Variables locales

    INTEGER i, j, itmax, itmay, iter
    REAL cvu(iip1, jjp1), cuv(iip1, jjm)
    REAL ai14, ai23, airez, rlatp, rlatm, xprm, xprp, un4rad2, yprp, yprm
    REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm
    REAL coslatm, coslatp, radclatm, radclatp
    REAL cuij1(iip1, jjp1), cuij2(iip1, jjp1), cuij3(iip1, jjp1), &
         cuij4(iip1, jjp1)
    REAL cvij1(iip1, jjp1), cvij2(iip1, jjp1), cvij3(iip1, jjp1), &
         cvij4(iip1, jjp1)
    REAL rlonvv(iip1), rlatuu(jjp1)
    REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm), yprimv(jjm), &
         yprimu(jjp1)
    REAL gamdi_gdiv, gamdi_grot, gamdi_h

    REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)
    SAVE rlatu1, yprimu1, rlatu2, yprimu2, yprimv, yprimu
    SAVE rlonm025, xprimm025, rlonp025, xprimp025

    real aireij1_2d(iim + 1, jjm + 1)
    real aireij2_2d(iim + 1, jjm + 1)
    real aireij3_2d(iim + 1, jjm + 1), aireij4_2d(iim + 1, jjm + 1) 
    real airuscv2_2d(iim + 1, jjm) 
    real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm) 
    real aivscu2gam_2d(iim + 1, jjm)

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

    PRINT *, 'Call sequence information: inigeom'

    IF (nitergdiv/=2) THEN
       gamdi_gdiv = coefdis/(real(nitergdiv)-2.)
    ELSE
       gamdi_gdiv = 0.
    END IF
    IF (nitergrot/=2) THEN
       gamdi_grot = coefdis/(real(nitergrot)-2.)
    ELSE
       gamdi_grot = 0.
    END IF
    IF (niterh/=2) THEN
       gamdi_h = coefdis/(real(niterh)-2.)
    ELSE
       gamdi_h = 0.
    END IF

    print *, 'gamdi_gdiv = ', gamdi_gdiv
    print *, "gamdi_grot = ", gamdi_grot
    print *, "gamdi_h = ", gamdi_h

    WRITE (6, 990)

    IF (.NOT. fxyhypb) THEN
       IF (ysinus) THEN
          print *, ' Inigeom, Y = Sinus (Latitude) '
          ! utilisation de f(x, y) avec y = sinus de la latitude 
          CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &
               rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &
               xprimm025, rlonp025, xprimp025)
       ELSE
          print *, 'Inigeom, Y = Latitude, der. sinusoid .'
          ! utilisation de f(x, y) a tangente sinusoidale, y etant la latit

          pxo = clon*pi/180.
          pyo = 2.*clat*pi/180.

          ! determination de transx (pour le zoom) par Newton-Raphson

          itmax = 10
          eps = .1E-7

          xo1 = 0.
          DO iter = 1, itmax
             x1 = xo1
             f = x1 + alphax*sin(x1-pxo)
             df = 1. + alphax*cos(x1-pxo)
             x1 = x1 - f/df
             xdm = abs(x1-xo1)
             IF (xdm<=eps) EXIT
             xo1 = x1
          END DO

          transx = xo1

          itmay = 10
          eps = .1E-7

          yo1 = 0.
          DO iter = 1, itmay
             y1 = yo1
             f = y1 + alphay*sin(y1-pyo)
             df = 1. + alphay*cos(y1-pyo)
             y1 = y1 - f/df
             ydm = abs(y1-yo1)
             IF (ydm<=eps) EXIT
             yo1 = y1
          END DO

          transy = yo1

          CALL fxy(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
               yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
               rlonp025, xprimp025)
       END IF
    ELSE
       ! Utilisation de fxyhyper, f(x, y) à dérivée tangente hyperbolique
       print *, 'Inigeom, Y = Latitude, dérivée tangente hyperbolique'
       CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &
            taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
            yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
            rlonp025, xprimp025)
    END IF

    rlatu(1) = asin(1.)
    rlatu(jjp1) = -rlatu(1)

    ! calcul aux poles 

    yprimu(1) = 0.
    yprimu(jjp1) = 0.

    un4rad2 = 0.25*rad*rad

    ! calcul des aires (aire_2d, aireu_2d, airev_2d, 1./aire_2d, 1./airez)
    ! - et de fext_2d, force de coriolis extensive

    ! A 1 point scalaire P (i, j) de la grille, reguliere en (X, Y), sont
    ! affectees 4 aires entourant P, calculees respectivement aux points
    ! (i + 1/4, j - 1/4) : aireij1_2d (i, j)
    ! (i + 1/4, j + 1/4) : aireij2_2d (i, j)
    ! (i - 1/4, j + 1/4) : aireij3_2d (i, j)
    ! (i - 1/4, j - 1/4) : aireij4_2d (i, j)

    !,
    ! Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X, Y).
    ! Chaque aire centree en 1 point scalaire P(i, j) est egale a la somme
    ! des 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui sont
    ! affectees au
    ! point (i, j).
    ! On definit en outre les coefficients alpha comme etant egaux a
    ! (aireij / aire_2d), c.a.d par exp.
    ! alpha1_2d(i, j)=aireij1_2d(i, j)/aire_2d(i, j)

    ! De meme, toute aire centree en 1 point U est egale a la somme des
    ! 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d entourant
    ! le point U.
    ! Idem pour airev_2d, airez.

    ! On a, pour chaque maille : dX = dY = 1

    ! V

    ! aireij4_2d . . aireij1_2d

    ! U . . P . U

    ! aireij3_2d . . aireij2_2d

    ! V

    ! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d,
    ! aireij3_2d, aireij4_2d
    ! qui entourent chaque aire_2d(i, j), ainsi que les 4 elongations
    ! elementaires
    ! cuij et les 4 elongat. cvij qui sont calculees aux memes
    ! endroits que les aireij.

    ! do 35 : boucle sur les jjm + 1 latitudes 

    DO j = 1, jjp1

       IF (j==1) THEN

          yprm = yprimu1(j)
          rlatm = rlatu1(j)

          coslatm = cos(rlatm)
          radclatm = 0.5*rad*coslatm

          DO i = 1, iim
             xprp = xprimp025(i)
             xprm = xprimm025(i)
             aireij2_2d(i, 1) = un4rad2*coslatm*xprp*yprm
             aireij3_2d(i, 1) = un4rad2*coslatm*xprm*yprm
             cuij2(i, 1) = radclatm*xprp
             cuij3(i, 1) = radclatm*xprm
             cvij2(i, 1) = 0.5*rad*yprm
             cvij3(i, 1) = cvij2(i, 1)
          END DO

          DO i = 1, iim
             aireij1_2d(i, 1) = 0.
             aireij4_2d(i, 1) = 0.
             cuij1(i, 1) = 0.
             cuij4(i, 1) = 0.
             cvij1(i, 1) = 0.
             cvij4(i, 1) = 0.
          END DO

       END IF

       IF (j==jjp1) THEN
          yprp = yprimu2(j-1)
          rlatp = rlatu2(j-1)

          coslatp = cos(rlatp)
          radclatp = 0.5*rad*coslatp

          DO i = 1, iim
             xprp = xprimp025(i)
             xprm = xprimm025(i)
             aireij1_2d(i, jjp1) = un4rad2*coslatp*xprp*yprp
             aireij4_2d(i, jjp1) = un4rad2*coslatp*xprm*yprp
             cuij1(i, jjp1) = radclatp*xprp
             cuij4(i, jjp1) = radclatp*xprm
             cvij1(i, jjp1) = 0.5*rad*yprp
             cvij4(i, jjp1) = cvij1(i, jjp1)
          END DO

          DO i = 1, iim
             aireij2_2d(i, jjp1) = 0.
             aireij3_2d(i, jjp1) = 0.
             cvij2(i, jjp1) = 0.
             cvij3(i, jjp1) = 0.
             cuij2(i, jjp1) = 0.
             cuij3(i, jjp1) = 0.
          END DO

       END IF

       IF (j>1 .AND. j<jjp1) THEN

          rlatp = rlatu2(j-1)
          yprp = yprimu2(j-1)
          rlatm = rlatu1(j)
          yprm = yprimu1(j)

          coslatm = cos(rlatm)
          coslatp = cos(rlatp)
          radclatp = 0.5*rad*coslatp
          radclatm = 0.5*rad*coslatm

          DO i = 1, iim
             xprp = xprimp025(i)
             xprm = xprimm025(i)

             ai14 = un4rad2*coslatp*yprp
             ai23 = un4rad2*coslatm*yprm
             aireij1_2d(i, j) = ai14*xprp
             aireij2_2d(i, j) = ai23*xprp
             aireij3_2d(i, j) = ai23*xprm
             aireij4_2d(i, j) = ai14*xprm
             cuij1(i, j) = radclatp*xprp
             cuij2(i, j) = radclatm*xprp
             cuij3(i, j) = radclatm*xprm
             cuij4(i, j) = radclatp*xprm
             cvij1(i, j) = 0.5*rad*yprp
             cvij2(i, j) = 0.5*rad*yprm
             cvij3(i, j) = cvij2(i, j)
             cvij4(i, j) = cvij1(i, j)
          END DO

       END IF

       ! periodicite 

       cvij1(iip1, j) = cvij1(1, j)
       cvij2(iip1, j) = cvij2(1, j)
       cvij3(iip1, j) = cvij3(1, j)
       cvij4(iip1, j) = cvij4(1, j)
       cuij1(iip1, j) = cuij1(1, j)
       cuij2(iip1, j) = cuij2(1, j)
       cuij3(iip1, j) = cuij3(1, j)
       cuij4(iip1, j) = cuij4(1, j)
       aireij1_2d(iip1, j) = aireij1_2d(1, j)
       aireij2_2d(iip1, j) = aireij2_2d(1, j)
       aireij3_2d(iip1, j) = aireij3_2d(1, j)
       aireij4_2d(iip1, j) = aireij4_2d(1, j)

    END DO

    DO j = 1, jjp1
       DO i = 1, iim
          aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
               + aireij3_2d(i, j) + aireij4_2d(i, j)
          alpha1_2d(i, j) = aireij1_2d(i, j)/aire_2d(i, j)
          alpha2_2d(i, j) = aireij2_2d(i, j)/aire_2d(i, j)
          alpha3_2d(i, j) = aireij3_2d(i, j)/aire_2d(i, j)
          alpha4_2d(i, j) = aireij4_2d(i, j)/aire_2d(i, j)
          alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
          alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
          alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)
          alpha3p4_2d(i, j) = alpha3_2d(i, j) + alpha4_2d(i, j)
       END DO

       aire_2d(iip1, j) = aire_2d(1, j)
       alpha1_2d(iip1, j) = alpha1_2d(1, j)
       alpha2_2d(iip1, j) = alpha2_2d(1, j)
       alpha3_2d(iip1, j) = alpha3_2d(1, j)
       alpha4_2d(iip1, j) = alpha4_2d(1, j)
       alpha1p2_2d(iip1, j) = alpha1p2_2d(1, j)
       alpha1p4_2d(iip1, j) = alpha1p4_2d(1, j)
       alpha2p3_2d(iip1, j) = alpha2p3_2d(1, j)
       alpha3p4_2d(iip1, j) = alpha3p4_2d(1, j)
    END DO

    DO j = 1, jjp1
       DO i = 1, iim
          aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
               aireij4_2d(i+1, j) + aireij3_2d(i+1, j)
          unsaire_2d(i, j) = 1./aire_2d(i, j)
          unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
          unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
          airesurg_2d(i, j) = aire_2d(i, j)/g
       END DO
       aireu_2d(iip1, j) = aireu_2d(1, j)
       unsaire_2d(iip1, j) = unsaire_2d(1, j)
       unsair_gam1_2d(iip1, j) = unsair_gam1_2d(1, j)
       unsair_gam2_2d(iip1, j) = unsair_gam2_2d(1, j)
       airesurg_2d(iip1, j) = airesurg_2d(1, j)
    END DO

    DO j = 1, jjm

       DO i = 1, iim
          airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
               aireij1_2d(i, j+1) + aireij4_2d(i, j+1)
       END DO
       DO i = 1, iim
          airez = aireij2_2d(i, j) + aireij1_2d(i, j+1) + aireij3_2d(i+1, j) &
               + aireij4_2d(i+1, j+1)
          unsairez_2d(i, j) = 1./airez
          unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
          fext_2d(i, j) = airez*sin(rlatv(j))*2.*omeg
       END DO
       airev_2d(iip1, j) = airev_2d(1, j)
       unsairez_2d(iip1, j) = unsairez_2d(1, j)
       fext_2d(iip1, j) = fext_2d(1, j)
       unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)

    END DO

    ! Calcul des elongations cu_2d, cv_2d, cvu 

    DO j = 1, jjm
       DO i = 1, iim
          cv_2d(i, j) = 0.5 * &
               (cvij2(i, j) + cvij3(i, j) + cvij1(i, j+1) + cvij4(i, j+1))
          cvu(i, j) = 0.5*(cvij1(i, j)+cvij4(i, j)+cvij2(i, j)+cvij3(i, j))
          cuv(i, j) = 0.5*(cuij2(i, j)+cuij3(i, j)+cuij1(i, j+1)+cuij4(i, j+1))
          unscv2_2d(i, j) = 1./(cv_2d(i, j)*cv_2d(i, j))
       END DO
       DO i = 1, iim
          cuvsurcv_2d(i, j) = airev_2d(i, j)*unscv2_2d(i, j)
          cvsurcuv_2d(i, j) = 1./cuvsurcv_2d(i, j)
          cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
          cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
          cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
       END DO
       cv_2d(iip1, j) = cv_2d(1, j)
       cvu(iip1, j) = cvu(1, j)
       unscv2_2d(iip1, j) = unscv2_2d(1, j)
       cuv(iip1, j) = cuv(1, j)
       cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
       cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
       cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
       cuvscvgam2_2d(iip1, j) = cuvscvgam2_2d(1, j)
       cvscuvgam_2d(iip1, j) = cvscuvgam_2d(1, j)
    END DO

    DO j = 2, jjm
       DO i = 1, iim
          cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i+1, j) + cuij2(i, j) &
               + cuij3(i+1, j))
          unscu2_2d(i, j) = 1./(cu_2d(i, j)*cu_2d(i, j))
          cvusurcu_2d(i, j) = aireu_2d(i, j)*unscu2_2d(i, j)
          cusurcvu_2d(i, j) = 1./cvusurcu_2d(i, j)
          cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
          cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
          cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
       END DO
       cu_2d(iip1, j) = cu_2d(1, j)
       unscu2_2d(iip1, j) = unscu2_2d(1, j)
       cvusurcu_2d(iip1, j) = cvusurcu_2d(1, j)
       cusurcvu_2d(iip1, j) = cusurcvu_2d(1, j)
       cvuscugam1_2d(iip1, j) = cvuscugam1_2d(1, j)
       cvuscugam2_2d(iip1, j) = cvuscugam2_2d(1, j)
       cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
    END DO

    ! calcul aux poles 

    DO i = 1, iip1
       cu_2d(i, 1) = 0.
       unscu2_2d(i, 1) = 0.
       cvu(i, 1) = 0.

       cu_2d(i, jjp1) = 0.
       unscu2_2d(i, jjp1) = 0.
       cvu(i, jjp1) = 0.
    END DO

    DO j = 1, jjm
       DO i = 1, iim
          airvscu2_2d(i, j) = airev_2d(i, j)/(cuv(i, j)*cuv(i, j))
          aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)
       END DO
       airvscu2_2d(iip1, j) = airvscu2_2d(1, j)
       aivscu2gam_2d(iip1, j) = aivscu2gam_2d(1, j)
    END DO

    DO j = 2, jjm
       DO i = 1, iim
          airuscv2_2d(i, j) = aireu_2d(i, j)/(cvu(i, j)*cvu(i, j))
          aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)
       END DO
       airuscv2_2d(iip1, j) = airuscv2_2d(1, j)
       aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
    END DO

    ! calcul des aires aux poles :

    apoln = sum(aire_2d(:iim, 1))
    apols = sum(aire_2d(:iim, jjp1))
    unsapolnga1 = 1./(apoln**(-gamdi_gdiv))
    unsapolsga1 = 1./(apols**(-gamdi_gdiv))
    unsapolnga2 = 1./(apoln**(-gamdi_h))
    unsapolsga2 = 1./(apols**(-gamdi_h))

    ! changement F. Hourdin calcul conservatif pour fext_2d
    ! constang_2d contient le produit a * cos (latitude) * omega

    DO i = 1, iim
       constang_2d(i, 1) = 0.
    END DO
    DO j = 1, jjm - 1
       DO i = 1, iim
          constang_2d(i, j+1) = rad*omeg*cu_2d(i, j+1)*cos(rlatu(j+1))
       END DO
    END DO
    DO i = 1, iim
       constang_2d(i, jjp1) = 0.
    END DO

    ! periodicite en longitude

    DO j = 1, jjm
       fext_2d(iip1, j) = fext_2d(1, j)
    END DO
    DO j = 1, jjp1
       constang_2d(iip1, j) = constang_2d(1, j)
    END DO

    ! fin du changement

    print *, ' Coordonnees de la grille '
    WRITE (6, 995)

    print *, ' LONGITUDES aux pts. V (degres) '
    WRITE (6, 995)
    DO i = 1, iip1
       rlonvv(i) = rlonv(i)*180./pi
    END DO
    WRITE (6, 400) rlonvv

    WRITE (6, 995)
    print *, ' LATITUDES aux pts. V (degres) '
    WRITE (6, 995)
    DO i = 1, jjm
       rlatuu(i) = rlatv(i)*180./pi
    END DO
    WRITE (6, 400) (rlatuu(i), i=1, jjm)

    DO i = 1, iip1
       rlonvv(i) = rlonu(i)*180./pi
    END DO
    WRITE (6, 995)
    print *, ' LONGITUDES aux pts. U (degres) '
    WRITE (6, 995)
    WRITE (6, 400) rlonvv
    WRITE (6, 995)

    print *, ' LATITUDES aux pts. U (degres) '
    WRITE (6, 995)
    DO i = 1, jjp1
       rlatuu(i) = rlatu(i)*180./pi
    END DO
    WRITE (6, 400) (rlatuu(i), i=1, jjp1)
    WRITE (6, 995)

400 FORMAT (1X, 8F8.2)
990 FORMAT (//)
995 FORMAT (/)

  END SUBROUTINE inigeom

end module inigeom_m
1	guez	25	module inigeom_m
2	guez	3
3	guez	7	IMPLICIT NONE
4	guez	3
5	guez	25	contains
6	guez	7
7	guez	25	SUBROUTINE inigeom
8	guez	7
9	guez	25	! Auteur : P. Le Van
10			! Version du 01/04/2001
11	guez	7
12	guez	25	! Calcul des élongations cuij1, ..., cuij4, cvij1, ..., cvij4 aux mêmes
13			! endroits que les aires aireij1_2d, ..., aireij4_2d.
14	guez	7
15	guez	25	! Choix entre une fonction "f(y)" à dérivée sinusoïdale ou à dérivée
16			! tangente hyperbolique
17			! calcul des coefficients (cu_2d, cv_2d, 1./cu_2d2, 1./cv_2d2)
18	guez	7
19	guez	38	! les coef. (cu_2d, cv_2d) permettent de passer des vitesses naturelles
20			! aux vitesses covariantes et contravariantes, ou vice-versa
21	guez	7
22	guez	38	! on a :
23			! u (covariant) = cu_2d * u (naturel), u(contrav)= u(nat)/cu_2d
24			! v (covariant) = cv_2d * v (naturel), v(contrav)= v(nat)/cv_2d
25	guez	7
26	guez	38	! on en tire :
27			! u(covariant) = cu_2d * cu_2d * u(contravariant)
28			! v(covariant) = cv_2d * cv_2d * v(contravariant)
29	guez	7
30	guez	38	! on a l'application (x(X), y(Y)) avec - im/2 +1 <= X <= im/2
31			! et - jm/2 <= Y <= jm/2
32	guez	7
33	guez	38	! x est la longitude du point en radians.
34			! y est la latitude du point en radians.
35			!
36			! on a : cu_2d(i, j) = rad * cos(y) * dx/dX
37			! cv(j) = rad * dy/dY
38			! aire_2d(i, j) = cu_2d(i, j) * cv(j)
39			!
40			! y, dx/dX, dy/dY calcules aux points concernes
41			! cv, bien que dependant de j uniquement, sera ici indice aussi en i
42			! pour un adressage plus facile en ij.
43	guez	7
44	guez	38	! aux points u et v,
45			! xprimu et xprimv sont respectivement les valeurs de dx/dX
46			! yprimu et yprimv sont respectivement les valeurs de dy/dY
47			! rlatu et rlatv sont respectivement les valeurs de la latitude
48			! cvu et cv_2d sont respectivement les valeurs de cv_2d
49	guez	7
50	guez	38	! aux points u, v, scalaires, et z
51			! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d
52			! Cf. "inigeom.txt".
53	guez	7
54	guez	39	USE comconst, ONLY : g, omeg, rad
55	guez	25	USE comgeom, ONLY : airesurg_2d, aireu_2d, airev_2d, aire_2d, &
56			alpha1p2_2d, alpha1p4_2d, alpha1_2d, &
57			alpha2p3_2d, alpha2_2d, alpha3p4_2d, alpha3_2d, alpha4_2d, apoln, &
58			apols, constang_2d, cuscvugam_2d, cusurcvu_2d, cuvscvgam1_2d, &
59			cuvscvgam2_2d, cuvsurcv_2d, cu_2d, cvscuvgam_2d, cvsurcuv_2d, &
60			cvuscugam1_2d, cvuscugam2_2d, cvusurcu_2d, cv_2d, fext_2d, rlatu, &
61			rlatv, rlonu, rlonv, unsairez_2d, unsaire_2d, unsairz_gam_2d, &
62			unsair_gam1_2d, unsair_gam2_2d, unsapolnga1, unsapolnga2, &
63			unsapolsga1, unsapolsga2, unscu2_2d, unscv2_2d, xprimu, xprimv
64	guez	39	USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
65			USE dimens_m, ONLY : iim, jjm
66			USE logic, ONLY : fxyhypb, ysinus
67			use nr_util, only: pi
68			USE paramet_m, ONLY : iip1, jjp1
69	guez	25	USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &
70			grossismy, pxo, pyo, taux, tauy, transx, transy
71	guez	7
72	guez	25	! Variables locales
73	guez	7
74	guez	25	INTEGER i, j, itmax, itmay, iter
75			REAL cvu(iip1, jjp1), cuv(iip1, jjm)
76			REAL ai14, ai23, airez, rlatp, rlatm, xprm, xprp, un4rad2, yprp, yprm
77			REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm
78			REAL coslatm, coslatp, radclatm, radclatp
79			REAL cuij1(iip1, jjp1), cuij2(iip1, jjp1), cuij3(iip1, jjp1), &
80			cuij4(iip1, jjp1)
81			REAL cvij1(iip1, jjp1), cvij2(iip1, jjp1), cvij3(iip1, jjp1), &
82			cvij4(iip1, jjp1)
83			REAL rlonvv(iip1), rlatuu(jjp1)
84			REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm), yprimv(jjm), &
85			yprimu(jjp1)
86			REAL gamdi_gdiv, gamdi_grot, gamdi_h
87	guez	7
88	guez	25	REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)
89			SAVE rlatu1, yprimu1, rlatu2, yprimu2, yprimv, yprimu
90			SAVE rlonm025, xprimm025, rlonp025, xprimp025
91	guez	7
92	guez	25	real aireij1_2d(iim + 1, jjm + 1)
93			real aireij2_2d(iim + 1, jjm + 1)
94			real aireij3_2d(iim + 1, jjm + 1), aireij4_2d(iim + 1, jjm + 1)
95			real airuscv2_2d(iim + 1, jjm)
96	guez	38	real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)
97	guez	25	real aivscu2gam_2d(iim + 1, jjm)
98	guez	7
99	guez	25	!------------------------------------------------------------------
100	guez	7
101	guez	25	PRINT *, 'Call sequence information: inigeom'
102	guez	7
103	guez	25	IF (nitergdiv/=2) THEN
104			gamdi_gdiv = coefdis/(real(nitergdiv)-2.)
105			ELSE
106			gamdi_gdiv = 0.
107			END IF
108			IF (nitergrot/=2) THEN
109			gamdi_grot = coefdis/(real(nitergrot)-2.)
110			ELSE
111			gamdi_grot = 0.
112			END IF
113			IF (niterh/=2) THEN
114			gamdi_h = coefdis/(real(niterh)-2.)
115			ELSE
116			gamdi_h = 0.
117			END IF
118	guez	7
119	guez	25	print *, 'gamdi_gdiv = ', gamdi_gdiv
120			print *, "gamdi_grot = ", gamdi_grot
121			print *, "gamdi_h = ", gamdi_h
122	guez	7
123	guez	25	WRITE (6, 990)
124	guez	7
125	guez	38	IF (.NOT. fxyhypb) THEN
126	guez	25	IF (ysinus) THEN
127	guez	38	print *, ' Inigeom, Y = Sinus (Latitude) '
128			! utilisation de f(x, y) avec y = sinus de la latitude
129	guez	25	CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &
130			rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &
131			xprimm025, rlonp025, xprimp025)
132			ELSE
133	guez	38	print *, 'Inigeom, Y = Latitude, der. sinusoid .'
134			! utilisation de f(x, y) a tangente sinusoidale, y etant la latit
135	guez	7
136	guez	25	pxo = clon*pi/180.
137			pyo = 2.clatpi/180.
138	guez	7
139	guez	38	! determination de transx (pour le zoom) par Newton-Raphson
140	guez	7
141	guez	25	itmax = 10
142			eps = .1E-7
143	guez	7
144	guez	25	xo1 = 0.
145			DO iter = 1, itmax
146			x1 = xo1
147			f = x1 + alphax*sin(x1-pxo)
148			df = 1. + alphax*cos(x1-pxo)
149			x1 = x1 - f/df
150			xdm = abs(x1-xo1)
151			IF (xdm<=eps) EXIT
152			xo1 = x1
153			END DO
154	guez	7
155	guez	25	transx = xo1
156	guez	7
157	guez	25	itmay = 10
158			eps = .1E-7
159	guez	7
160	guez	25	yo1 = 0.
161			DO iter = 1, itmay
162			y1 = yo1
163			f = y1 + alphay*sin(y1-pyo)
164			df = 1. + alphay*cos(y1-pyo)
165			y1 = y1 - f/df
166			ydm = abs(y1-yo1)
167			IF (ydm<=eps) EXIT
168			yo1 = y1
169			END DO
170	guez	7
171	guez	25	transy = yo1
172	guez	7
173	guez	25	CALL fxy(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
174			yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
175			rlonp025, xprimp025)
176			END IF
177			ELSE
178	guez	38	! Utilisation de fxyhyper, f(x, y) à dérivée tangente hyperbolique
179			print *, 'Inigeom, Y = Latitude, dérivée tangente hyperbolique'
180	guez	25	CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &
181			taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
182			yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
183			rlonp025, xprimp025)
184			END IF
185	guez	7
186	guez	25	rlatu(1) = asin(1.)
187			rlatu(jjp1) = -rlatu(1)
188	guez	7
189	guez	38	! calcul aux poles
190	guez	7
191	guez	25	yprimu(1) = 0.
192			yprimu(jjp1) = 0.
193	guez	7
194	guez	25	un4rad2 = 0.25radrad
195	guez	7
196	guez	38	! calcul des aires (aire_2d, aireu_2d, airev_2d, 1./aire_2d, 1./airez)
197			! - et de fext_2d, force de coriolis extensive
198	guez	7
199	guez	38	! A 1 point scalaire P (i, j) de la grille, reguliere en (X, Y), sont
200			! affectees 4 aires entourant P, calculees respectivement aux points
201			! (i + 1/4, j - 1/4) : aireij1_2d (i, j)
202			! (i + 1/4, j + 1/4) : aireij2_2d (i, j)
203			! (i - 1/4, j + 1/4) : aireij3_2d (i, j)
204			! (i - 1/4, j - 1/4) : aireij4_2d (i, j)
205	guez	7
206	guez	38	!,
207	guez	25	! Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X, Y).
208	guez	38	! Chaque aire centree en 1 point scalaire P(i, j) est egale a la somme
209	guez	25	! des 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui sont
210			! affectees au
211	guez	38	! point (i, j).
212			! On definit en outre les coefficients alpha comme etant egaux a
213	guez	25	! (aireij / aire_2d), c.a.d par exp.
214			! alpha1_2d(i, j)=aireij1_2d(i, j)/aire_2d(i, j)
215	guez	7
216	guez	38	! De meme, toute aire centree en 1 point U est egale a la somme des
217	guez	25	! 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d entourant
218			! le point U.
219	guez	38	! Idem pour airev_2d, airez.
220	guez	7
221	guez	38	! On a, pour chaque maille : dX = dY = 1
222	guez	7
223	guez	38	! V
224	guez	7
225	guez	38	! aireij4_2d . . aireij1_2d
226	guez	7
227	guez	38	! U . . P . U
228	guez	7
229	guez	38	! aireij3_2d . . aireij2_2d
230	guez	7
231	guez	38	! V
232	guez	3
233	guez	25	! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d,
234			! aireij3_2d, aireij4_2d
235	guez	38	! qui entourent chaque aire_2d(i, j), ainsi que les 4 elongations
236	guez	25	! elementaires
237			! cuij et les 4 elongat. cvij qui sont calculees aux memes
238	guez	38	! endroits que les aireij.
239	guez	7
240	guez	38	! do 35 : boucle sur les jjm + 1 latitudes
241	guez	7
242	guez	25	DO j = 1, jjp1
243	guez	3
244	guez	25	IF (j==1) THEN
245	guez	3
246	guez	25	yprm = yprimu1(j)
247			rlatm = rlatu1(j)
248	guez	3
249	guez	25	coslatm = cos(rlatm)
250			radclatm = 0.5radcoslatm
251	guez	3
252	guez	25	DO i = 1, iim
253			xprp = xprimp025(i)
254			xprm = xprimm025(i)
255			aireij2_2d(i, 1) = un4rad2coslatmxprp*yprm
256			aireij3_2d(i, 1) = un4rad2coslatmxprm*yprm
257			cuij2(i, 1) = radclatm*xprp
258			cuij3(i, 1) = radclatm*xprm
259			cvij2(i, 1) = 0.5radyprm
260			cvij3(i, 1) = cvij2(i, 1)
261			END DO
262	guez	3
263	guez	25	DO i = 1, iim
264			aireij1_2d(i, 1) = 0.
265			aireij4_2d(i, 1) = 0.
266			cuij1(i, 1) = 0.
267			cuij4(i, 1) = 0.
268			cvij1(i, 1) = 0.
269			cvij4(i, 1) = 0.
270			END DO
271	guez	3
272	guez	25	END IF
273	guez	3
274	guez	25	IF (j==jjp1) THEN
275			yprp = yprimu2(j-1)
276			rlatp = rlatu2(j-1)
277	guez	3
278	guez	25	coslatp = cos(rlatp)
279			radclatp = 0.5radcoslatp
280	guez	3
281	guez	25	DO i = 1, iim
282			xprp = xprimp025(i)
283			xprm = xprimm025(i)
284			aireij1_2d(i, jjp1) = un4rad2coslatpxprp*yprp
285			aireij4_2d(i, jjp1) = un4rad2coslatpxprm*yprp
286			cuij1(i, jjp1) = radclatp*xprp
287			cuij4(i, jjp1) = radclatp*xprm
288			cvij1(i, jjp1) = 0.5radyprp
289			cvij4(i, jjp1) = cvij1(i, jjp1)
290			END DO
291	guez	3
292	guez	25	DO i = 1, iim
293			aireij2_2d(i, jjp1) = 0.
294			aireij3_2d(i, jjp1) = 0.
295			cvij2(i, jjp1) = 0.
296			cvij3(i, jjp1) = 0.
297			cuij2(i, jjp1) = 0.
298			cuij3(i, jjp1) = 0.
299			END DO
300	guez	3
301	guez	25	END IF
302	guez	3
303	guez	25	IF (j>1 .AND. j<jjp1) THEN
304	guez	3
305	guez	25	rlatp = rlatu2(j-1)
306			yprp = yprimu2(j-1)
307			rlatm = rlatu1(j)
308			yprm = yprimu1(j)
309	guez	3
310	guez	25	coslatm = cos(rlatm)
311			coslatp = cos(rlatp)
312			radclatp = 0.5radcoslatp
313			radclatm = 0.5radcoslatm
314	guez	3
315	guez	25	DO i = 1, iim
316			xprp = xprimp025(i)
317			xprm = xprimm025(i)
318	guez	7
319	guez	25	ai14 = un4rad2coslatpyprp
320			ai23 = un4rad2coslatmyprm
321			aireij1_2d(i, j) = ai14*xprp
322			aireij2_2d(i, j) = ai23*xprp
323			aireij3_2d(i, j) = ai23*xprm
324			aireij4_2d(i, j) = ai14*xprm
325			cuij1(i, j) = radclatp*xprp
326			cuij2(i, j) = radclatm*xprp
327			cuij3(i, j) = radclatm*xprm
328			cuij4(i, j) = radclatp*xprm
329			cvij1(i, j) = 0.5radyprp
330			cvij2(i, j) = 0.5radyprm
331			cvij3(i, j) = cvij2(i, j)
332			cvij4(i, j) = cvij1(i, j)
333			END DO
334	guez	7
335	guez	25	END IF
336	guez	7
337	guez	38	! periodicite
338	guez	7
339	guez	25	cvij1(iip1, j) = cvij1(1, j)
340			cvij2(iip1, j) = cvij2(1, j)
341			cvij3(iip1, j) = cvij3(1, j)
342			cvij4(iip1, j) = cvij4(1, j)
343			cuij1(iip1, j) = cuij1(1, j)
344			cuij2(iip1, j) = cuij2(1, j)
345			cuij3(iip1, j) = cuij3(1, j)
346			cuij4(iip1, j) = cuij4(1, j)
347			aireij1_2d(iip1, j) = aireij1_2d(1, j)
348			aireij2_2d(iip1, j) = aireij2_2d(1, j)
349			aireij3_2d(iip1, j) = aireij3_2d(1, j)
350			aireij4_2d(iip1, j) = aireij4_2d(1, j)
351	guez	7
352	guez	25	END DO
353	guez	7
354	guez	25	DO j = 1, jjp1
355			DO i = 1, iim
356			aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
357			+ aireij3_2d(i, j) + aireij4_2d(i, j)
358			alpha1_2d(i, j) = aireij1_2d(i, j)/aire_2d(i, j)
359			alpha2_2d(i, j) = aireij2_2d(i, j)/aire_2d(i, j)
360			alpha3_2d(i, j) = aireij3_2d(i, j)/aire_2d(i, j)
361			alpha4_2d(i, j) = aireij4_2d(i, j)/aire_2d(i, j)
362			alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
363			alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
364			alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)
365			alpha3p4_2d(i, j) = alpha3_2d(i, j) + alpha4_2d(i, j)
366			END DO
367	guez	7
368	guez	25	aire_2d(iip1, j) = aire_2d(1, j)
369			alpha1_2d(iip1, j) = alpha1_2d(1, j)
370			alpha2_2d(iip1, j) = alpha2_2d(1, j)
371			alpha3_2d(iip1, j) = alpha3_2d(1, j)
372			alpha4_2d(iip1, j) = alpha4_2d(1, j)
373			alpha1p2_2d(iip1, j) = alpha1p2_2d(1, j)
374			alpha1p4_2d(iip1, j) = alpha1p4_2d(1, j)
375			alpha2p3_2d(iip1, j) = alpha2p3_2d(1, j)
376			alpha3p4_2d(iip1, j) = alpha3p4_2d(1, j)
377			END DO
378	guez	7
379	guez	25	DO j = 1, jjp1
380			DO i = 1, iim
381			aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
382			aireij4_2d(i+1, j) + aireij3_2d(i+1, j)
383			unsaire_2d(i, j) = 1./aire_2d(i, j)
384			unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
385			unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
386			airesurg_2d(i, j) = aire_2d(i, j)/g
387			END DO
388			aireu_2d(iip1, j) = aireu_2d(1, j)
389			unsaire_2d(iip1, j) = unsaire_2d(1, j)
390			unsair_gam1_2d(iip1, j) = unsair_gam1_2d(1, j)
391			unsair_gam2_2d(iip1, j) = unsair_gam2_2d(1, j)
392			airesurg_2d(iip1, j) = airesurg_2d(1, j)
393			END DO
394	guez	7
395	guez	25	DO j = 1, jjm
396	guez	7
397	guez	25	DO i = 1, iim
398			airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
399			aireij1_2d(i, j+1) + aireij4_2d(i, j+1)
400			END DO
401			DO i = 1, iim
402			airez = aireij2_2d(i, j) + aireij1_2d(i, j+1) + aireij3_2d(i+1, j) &
403			+ aireij4_2d(i+1, j+1)
404			unsairez_2d(i, j) = 1./airez
405			unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
406			fext_2d(i, j) = airezsin(rlatv(j))2.*omeg
407			END DO
408			airev_2d(iip1, j) = airev_2d(1, j)
409			unsairez_2d(iip1, j) = unsairez_2d(1, j)
410			fext_2d(iip1, j) = fext_2d(1, j)
411			unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
412	guez	7
413	guez	25	END DO
414	guez	7
415	guez	38	! Calcul des elongations cu_2d, cv_2d, cvu
416	guez	7
417	guez	25	DO j = 1, jjm
418			DO i = 1, iim
419			cv_2d(i, j) = 0.5 * &
420			(cvij2(i, j) + cvij3(i, j) + cvij1(i, j+1) + cvij4(i, j+1))
421			cvu(i, j) = 0.5*(cvij1(i, j)+cvij4(i, j)+cvij2(i, j)+cvij3(i, j))
422			cuv(i, j) = 0.5*(cuij2(i, j)+cuij3(i, j)+cuij1(i, j+1)+cuij4(i, j+1))
423			unscv2_2d(i, j) = 1./(cv_2d(i, j)*cv_2d(i, j))
424			END DO
425			DO i = 1, iim
426			cuvsurcv_2d(i, j) = airev_2d(i, j)*unscv2_2d(i, j)
427			cvsurcuv_2d(i, j) = 1./cuvsurcv_2d(i, j)
428			cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
429			cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
430			cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
431			END DO
432			cv_2d(iip1, j) = cv_2d(1, j)
433			cvu(iip1, j) = cvu(1, j)
434			unscv2_2d(iip1, j) = unscv2_2d(1, j)
435			cuv(iip1, j) = cuv(1, j)
436			cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
437			cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
438			cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
439			cuvscvgam2_2d(iip1, j) = cuvscvgam2_2d(1, j)
440			cvscuvgam_2d(iip1, j) = cvscuvgam_2d(1, j)
441			END DO
442	guez	7
443	guez	25	DO j = 2, jjm
444			DO i = 1, iim
445			cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i+1, j) + cuij2(i, j) &
446			+ cuij3(i+1, j))
447			unscu2_2d(i, j) = 1./(cu_2d(i, j)*cu_2d(i, j))
448			cvusurcu_2d(i, j) = aireu_2d(i, j)*unscu2_2d(i, j)
449			cusurcvu_2d(i, j) = 1./cvusurcu_2d(i, j)
450			cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
451			cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
452			cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
453			END DO
454			cu_2d(iip1, j) = cu_2d(1, j)
455			unscu2_2d(iip1, j) = unscu2_2d(1, j)
456			cvusurcu_2d(iip1, j) = cvusurcu_2d(1, j)
457			cusurcvu_2d(iip1, j) = cusurcvu_2d(1, j)
458			cvuscugam1_2d(iip1, j) = cvuscugam1_2d(1, j)
459			cvuscugam2_2d(iip1, j) = cvuscugam2_2d(1, j)
460			cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
461			END DO
462	guez	7
463	guez	38	! calcul aux poles
464	guez	7
465	guez	25	DO i = 1, iip1
466			cu_2d(i, 1) = 0.
467			unscu2_2d(i, 1) = 0.
468			cvu(i, 1) = 0.
469	guez	7
470	guez	25	cu_2d(i, jjp1) = 0.
471			unscu2_2d(i, jjp1) = 0.
472			cvu(i, jjp1) = 0.
473			END DO
474	guez	7
475	guez	25	DO j = 1, jjm
476			DO i = 1, iim
477			airvscu2_2d(i, j) = airev_2d(i, j)/(cuv(i, j)*cuv(i, j))
478			aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)
479			END DO
480			airvscu2_2d(iip1, j) = airvscu2_2d(1, j)
481			aivscu2gam_2d(iip1, j) = aivscu2gam_2d(1, j)
482			END DO
483	guez	7
484	guez	25	DO j = 2, jjm
485			DO i = 1, iim
486			airuscv2_2d(i, j) = aireu_2d(i, j)/(cvu(i, j)*cvu(i, j))
487			aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)
488			END DO
489			airuscv2_2d(iip1, j) = airuscv2_2d(1, j)
490			aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
491			END DO
492	guez	7
493	guez	38	! calcul des aires aux poles :
494	guez	7
495	guez	25	apoln = sum(aire_2d(:iim, 1))
496			apols = sum(aire_2d(:iim, jjp1))
497			unsapolnga1 = 1./(apoln**(-gamdi_gdiv))
498			unsapolsga1 = 1./(apols**(-gamdi_gdiv))
499			unsapolnga2 = 1./(apoln**(-gamdi_h))
500			unsapolsga2 = 1./(apols**(-gamdi_h))
501	guez	7
502	guez	38	! changement F. Hourdin calcul conservatif pour fext_2d
503			! constang_2d contient le produit a * cos (latitude) * omega
504	guez	7
505	guez	25	DO i = 1, iim
506			constang_2d(i, 1) = 0.
507			END DO
508			DO j = 1, jjm - 1
509			DO i = 1, iim
510			constang_2d(i, j+1) = radomegcu_2d(i, j+1)*cos(rlatu(j+1))
511			END DO
512			END DO
513			DO i = 1, iim
514			constang_2d(i, jjp1) = 0.
515			END DO
516	guez	7
517	guez	38	! periodicite en longitude
518	guez	7
519	guez	25	DO j = 1, jjm
520			fext_2d(iip1, j) = fext_2d(1, j)
521			END DO
522			DO j = 1, jjp1
523			constang_2d(iip1, j) = constang_2d(1, j)
524			END DO
525	guez	7
526	guez	25	! fin du changement
527	guez	7
528	guez	38	print *, ' Coordonnees de la grille '
529	guez	25	WRITE (6, 995)
530	guez	7
531	guez	38	print *, ' LONGITUDES aux pts. V (degres) '
532	guez	25	WRITE (6, 995)
533			DO i = 1, iip1
534			rlonvv(i) = rlonv(i)*180./pi
535			END DO
536			WRITE (6, 400) rlonvv
537	guez	7
538	guez	25	WRITE (6, 995)
539	guez	38	print *, ' LATITUDES aux pts. V (degres) '
540	guez	25	WRITE (6, 995)
541			DO i = 1, jjm
542			rlatuu(i) = rlatv(i)*180./pi
543			END DO
544			WRITE (6, 400) (rlatuu(i), i=1, jjm)
545	guez	7
546	guez	25	DO i = 1, iip1
547			rlonvv(i) = rlonu(i)*180./pi
548			END DO
549			WRITE (6, 995)
550	guez	38	print *, ' LONGITUDES aux pts. U (degres) '
551	guez	25	WRITE (6, 995)
552			WRITE (6, 400) rlonvv
553			WRITE (6, 995)
554	guez	7
555	guez	38	print *, ' LATITUDES aux pts. U (degres) '
556	guez	25	WRITE (6, 995)
557			DO i = 1, jjp1
558			rlatuu(i) = rlatu(i)*180./pi
559			END DO
560			WRITE (6, 400) (rlatuu(i), i=1, jjp1)
561			WRITE (6, 995)
562	guez	7
563	guez	25	400 FORMAT (1X, 8F8.2)
564	guez	7	990 FORMAT (//)
565			995 FORMAT (/)
566
567	guez	25	END SUBROUTINE inigeom
568
569			end module inigeom_m