trunk/dyn3d/comgeom.f

module comgeom

  use dimens_m, only: iim, jjm

  implicit none

  private iim, jjm

  real cu_2d(iim + 1, jjm + 1), cv_2d(iim + 1, jjm) ! in m
  real cu((iim + 1) * (jjm + 1)), cv((iim + 1) * jjm) ! in m
  equivalence (cu, cu_2d), (cv, cv_2d)

  real unscu2_2d(iim + 1, jjm + 1) ! in m-2
  real unscu2((iim + 1) * (jjm + 1)) ! in m-2
  equivalence (unscu2, unscu2_2d)

  real unscv2_2d(iim + 1, jjm) ! in m-2
  real unscv2((iim + 1) * jjm) ! in m-2
  equivalence (unscv2, unscv2_2d)

  real aire((iim + 1) * (jjm + 1)), aire_2d(iim + 1, jjm + 1) ! in m2
  real airesurg_2d(iim + 1, jjm + 1), airesurg((iim + 1) * (jjm + 1))
  equivalence (aire, aire_2d), (airesurg, airesurg_2d)

  real aireu_2d(iim + 1, jjm + 1) ! in m2
  real aireu((iim + 1) * (jjm + 1)) ! in m2
  equivalence (aireu, aireu_2d)

  real airev((iim + 1) * jjm), airev_2d(iim + 1, jjm) ! in m2
  real unsaire((iim + 1) * (jjm + 1)), unsaire_2d(iim + 1, jjm + 1) ! in m-2
  equivalence (airev, airev_2d), (unsaire, unsaire_2d)

  real apoln, apols ! in m2

  real unsairez_2d(iim + 1, jjm)
  real unsairez((iim + 1) * jjm)
  equivalence (unsairez, unsairez_2d)

  real alpha1_2d(iim + 1, jjm + 1)
  real alpha1((iim + 1) * (jjm + 1))
  equivalence (alpha1, alpha1_2d)

  real alpha2_2d(iim + 1, jjm + 1)         
  real alpha2((iim + 1) * (jjm + 1))
  equivalence (alpha2, alpha2_2d)

  real alpha3_2d(iim + 1, jjm + 1), alpha4_2d(iim + 1, jjm + 1)
  real alpha3((iim + 1) * (jjm + 1)), alpha4((iim + 1) * (jjm + 1))
  equivalence (alpha3, alpha3_2d), (alpha4, alpha4_2d)

  real alpha1p2_2d(iim + 1, jjm + 1)        
  real alpha1p2((iim + 1) * (jjm + 1))
  equivalence (alpha1p2, alpha1p2_2d)

  real alpha1p4_2d(iim + 1, jjm + 1), alpha2p3_2d(iim + 1, jjm + 1)
  real alpha1p4((iim + 1) * (jjm + 1)), alpha2p3((iim + 1) * (jjm + 1))
  equivalence (alpha1p4, alpha1p4_2d), (alpha2p3, alpha2p3_2d)

  real alpha3p4((iim + 1) * (jjm + 1))
  real alpha3p4_2d(iim + 1, jjm + 1)    
  equivalence (alpha3p4, alpha3p4_2d)

  real fext_2d(iim + 1, jjm), constang_2d(iim + 1, jjm + 1)
  real fext((iim + 1) * jjm), constang((iim + 1) * (jjm + 1))
  equivalence (fext, fext_2d), (constang, constang_2d)

  real rlatu(jjm + 1)
  ! (latitudes of points of the "scalar" and "u" grid, in rad)

  real rlatv(jjm) 
  ! (latitudes of points of the "v" grid, in rad, in decreasing order)

  real rlonu(iim + 1) ! longitudes of points of the "u" grid, in rad

  real rlonv(iim + 1)
  ! (longitudes of points of the "scalar" and "v" grid, in rad)

  real cuvsurcv_2d(iim + 1, jjm), cvsurcuv_2d(iim + 1, jjm) ! no dimension
  real cuvsurcv((iim + 1) * jjm), cvsurcuv((iim + 1) * jjm) ! no dimension
  equivalence (cuvsurcv, cuvsurcv_2d), (cvsurcuv, cvsurcuv_2d)

  real cvusurcu_2d(iim + 1, jjm + 1), cusurcvu_2d(iim + 1, jjm + 1)
  ! no dimension
  real cvusurcu((iim + 1) * (jjm + 1)), cusurcvu((iim + 1) * (jjm + 1))
  ! no dimension
  equivalence (cvusurcu, cvusurcu_2d), (cusurcvu, cusurcvu_2d)

  real cuvscvgam1_2d(iim + 1, jjm)
  real cuvscvgam1((iim + 1) * jjm)
  equivalence (cuvscvgam1, cuvscvgam1_2d)

  real cuvscvgam2_2d(iim + 1, jjm), cvuscugam1_2d(iim + 1, jjm + 1)
  real cuvscvgam2((iim + 1) * jjm), cvuscugam1((iim + 1) * (jjm + 1))
  equivalence (cuvscvgam2, cuvscvgam2_2d), (cvuscugam1, cvuscugam1_2d)

  real cvuscugam2_2d(iim + 1, jjm + 1), cvscuvgam_2d(iim + 1, jjm)
  real cvuscugam2((iim + 1) * (jjm + 1)), cvscuvgam((iim + 1) * jjm)
  equivalence (cvuscugam2, cvuscugam2_2d), (cvscuvgam, cvscuvgam_2d)

  real cuscvugam((iim + 1) * (jjm + 1))
  real cuscvugam_2d(iim + 1, jjm + 1) 
  equivalence (cuscvugam, cuscvugam_2d)

  real unsapolnga1, unsapolnga2, unsapolsga1, unsapolsga2                

  real unsair_gam1_2d(iim + 1, jjm + 1), unsair_gam2_2d(iim + 1, jjm + 1)
  real unsair_gam1((iim + 1) * (jjm + 1)), unsair_gam2((iim + 1) * (jjm + 1))
  equivalence (unsair_gam1, unsair_gam1_2d), (unsair_gam2, unsair_gam2_2d)

  real unsairz_gam_2d(iim + 1, jjm)
  real unsairz_gam((iim + 1) * jjm)
  equivalence (unsairz_gam, unsairz_gam_2d)

  real xprimu(iim + 1), xprimv(iim + 1)

  save

contains

  SUBROUTINE inigeom

    ! Auteur : P. Le Van

    ! 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 coefficients cu_2d et 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(contravariant) = u(naturel) / cu_2d
    ! v(covariant) = cv_2d * v(naturel), v(contravariant) = v(naturel) / 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 calculés aux points concernés. cv, bien que
    ! dépendant de j uniquement, sera ici indicé aussi en i pour un
    ! adressage plus facile en ij.

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

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

    USE comconst, ONLY : g, omeg, rad
    USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
    use conf_gcm_m, ONLY : fxyhypb, ysinus
    use fxy_m, only: fxy
    use fxyhyper_m, only: fxyhyper
    use jumble, only: new_unit
    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
    ! Modifiés pxo, pyo, transx, transy

    ! Local:
    INTEGER i, j, itmax, itmay, iter, unit
    REAL cvu(iip1, jjp1), cuv(iip1, jjm)
    REAL ai14, ai23, airez, un4rad2
    REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm
    REAL coslatm, coslatp, radclatm, radclatp
    REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
    REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4 ! in m
    REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
    real yprimv(jjm), yprimu(jjp1)
    REAL gamdi_gdiv, gamdi_grot, gamdi_h
    REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)
    real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &
         aireij4_2d ! in m2
    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

    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) = pi / 2.
    rlatu(jjp1) = -rlatu(1)

    ! Calcul aux pôles 

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

    un4rad2 = 0.25 * rad * rad

    ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
    ! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
    ! chaque aire_2d(i, j), ainsi que les quatre élongations
    ! élémentaires cuij et les quatre élongations cvij qui sont
    ! calculées aux mêmes endroits que les aireij.

    coslatm = cos(rlatu1(1))
    radclatm = 0.5 * rad * coslatm

    aireij1_2d(:iim, 1) = 0.
    aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
    aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
    aireij4_2d(:iim, 1) = 0.

    cuij1(:iim, 1) = 0.
    cuij2(:iim, 1) = radclatm * xprimp025(:iim)
    cuij3(:iim, 1) = radclatm * xprimm025(:iim)
    cuij4(:iim, 1) = 0.

    cvij1(:iim, 1) = 0.
    cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
    cvij3(:iim, 1) = cvij2(:iim, 1)
    cvij4(:iim, 1) = 0.

    do j = 2, jjm
       coslatm = cos(rlatu1(j))
       coslatp = cos(rlatu2(j-1))
       radclatp = 0.5 * rad * coslatp
       radclatm = 0.5 * rad * coslatm
       ai14 = un4rad2 * coslatp * yprimu2(j-1)
       ai23 = un4rad2 * coslatm * yprimu1(j)

       aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
       aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
       aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
       aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
       cuij1(:iim, j) = radclatp * xprimp025(:iim)
       cuij2(:iim, j) = radclatm * xprimp025(:iim)
       cuij3(:iim, j) = radclatm * xprimm025(:iim)
       cuij4(:iim, j) = radclatp * xprimm025(:iim)
       cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
       cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
       cvij3(:iim, j) = cvij2(:iim, j)
       cvij4(:iim, j) = cvij1(:iim, j)
    end do

    coslatp = cos(rlatu2(jjm))
    radclatp = 0.5 * rad * coslatp

    aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
    aireij2_2d(:iim, jjp1) = 0.
    aireij3_2d(:iim, jjp1) = 0.
    aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)

    cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
    cuij2(:iim, jjp1) = 0.
    cuij3(:iim, jjp1) = 0.
    cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)

    cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
    cvij2(:iim, jjp1) = 0.
    cvij3(:iim, jjp1) = 0.
    cvij4(:iim, jjp1) = cvij1(:iim, jjp1)

    ! Périodicité :

    cvij1(iip1, :) = cvij1(1, :)
    cvij2(iip1, :) = cvij2(1, :)
    cvij3(iip1, :) = cvij3(1, :)
    cvij4(iip1, :) = cvij4(1, :)

    cuij1(iip1, :) = cuij1(1, :)
    cuij2(iip1, :) = cuij2(1, :)
    cuij3(iip1, :) = cuij3(1, :)
    cuij4(iip1, :) = cuij4(1, :)

    aireij1_2d(iip1, :) = aireij1_2d(1, :)
    aireij2_2d(iip1, :) = aireij2_2d(1, :)
    aireij3_2d(iip1, :) = aireij3_2d(1, :)
    aireij4_2d(iip1, :) = aireij4_2d(1, :)

    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 élongations 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)**2
       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)**2
          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 pôles 

    cu_2d(:, 1) = 0.
    unscu2_2d(:, 1) = 0.
    cvu(:, 1) = 0.

    cu_2d(:, jjp1) = 0.
    unscu2_2d(:, jjp1) = 0.
    cvu(:, jjp1) = 0.

    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 pôles :

    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

    ! Périodicité en longitude
    DO j = 1, jjp1
       constang_2d(iip1, j) = constang_2d(1, j)
    END DO

    call new_unit(unit)
    open(unit, file="longitude_latitude.txt", status="replace", action="write")
    write(unit, fmt=*) '"longitudes at V points (degrees)"', rlonv * 180. / pi
    write(unit, fmt=*) '"latitudes at V points (degrees)"', rlatv * 180. / pi
    write(unit, fmt=*) '"longitudes at U points (degrees)"', rlonu * 180. / pi
    write(unit, fmt=*) '"latitudes at U points (degrees)"', rlatu * 180. / pi
    close(unit)

  END SUBROUTINE inigeom

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