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