Sources/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.

    ! Fonction "f(y)" à 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 fxhyp_m, only: fxhyp
    use fyhyp_m, only: fyhyp
    use jumble, only: new_unit
    use nr_util, only: pi
    USE paramet_m, ONLY : iip1, jjp1

    ! Local:
    INTEGER i, j, unit
    REAL cvu(iip1, jjp1), cuv(iip1, jjm)
    REAL ai14, ai23, airez, un4rad2
    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 yprimu(jjp1)
    REAL gamdi_gdiv, gamdi_grot, gamdi_h
    REAL xprimm025(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 / (nitergdiv - 2)
    ELSE
       gamdi_gdiv = 0.
    END IF

    IF (nitergrot /= 2) THEN
       gamdi_grot = coefdis / (nitergrot - 2)
    ELSE
       gamdi_grot = 0.
    END IF

    IF (niterh /= 2) THEN
       gamdi_h = coefdis / (niterh - 2)
    ELSE
       gamdi_h = 0.
    END IF

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

    CALL fyhyp(rlatu, yprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)
    CALL fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)

    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	guez	121	! Fonction "f(y)" à dérivée tangente hyperbolique. Calcul des
128			! coefficients cu_2d, cv_2d, 1. / cu_2d2, 1. / cv_2d2. Les
129			! coefficients cu_2d et cv_2d permettent de passer des vitesses
130			! naturelles aux vitesses covariantes et contravariantes, ou
131			! vice-versa.
132	guez	78
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	guez	119	use fxhyp_m, only: fxhyp
167			use fyhyp_m, only: fyhyp
168	guez	78	use jumble, only: new_unit
169			use nr_util, only: pi
170			USE paramet_m, ONLY : iip1, jjp1
171
172	guez	96	! Local:
173	guez	113	INTEGER i, j, unit
174	guez	78	REAL cvu(iip1, jjp1), cuv(iip1, jjm)
175			REAL ai14, ai23, airez, un4rad2
176			REAL coslatm, coslatp, radclatm, radclatp
177			REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
178			REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4 ! in m
179			REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
180	guez	119	real yprimu(jjp1)
181	guez	78	REAL gamdi_gdiv, gamdi_grot, gamdi_h
182	guez	119	REAL xprimm025(iip1), xprimp025(iip1)
183	guez	78	real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &
184			aireij4_2d ! in m2
185			real airuscv2_2d(iim + 1, jjm)
186			real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)
187			real aivscu2gam_2d(iim + 1, jjm)
188
189			!------------------------------------------------------------------
190
191			PRINT *, 'Call sequence information: inigeom'
192
193	guez	125	IF (nitergdiv /= 2) THEN
194			gamdi_gdiv = coefdis / (nitergdiv - 2)
195	guez	78	ELSE
196			gamdi_gdiv = 0.
197			END IF
198	guez	121
199	guez	125	IF (nitergrot /= 2) THEN
200			gamdi_grot = coefdis / (nitergrot - 2)
201	guez	78	ELSE
202			gamdi_grot = 0.
203			END IF
204	guez	121
205	guez	125	IF (niterh /= 2) THEN
206			gamdi_h = coefdis / (niterh - 2)
207	guez	78	ELSE
208			gamdi_h = 0.
209			END IF
210
211			print *, 'gamdi_gdiv = ', gamdi_gdiv
212			print *, "gamdi_grot = ", gamdi_grot
213			print *, "gamdi_h = ", gamdi_h
214
215	guez	119	CALL fyhyp(rlatu, yprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)
216			CALL fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
217	guez	78
218			rlatu(1) = pi / 2.
219			rlatu(jjp1) = -rlatu(1)
220
221			! Calcul aux pôles
222
223			yprimu(1) = 0.
224			yprimu(jjp1) = 0.
225
226			un4rad2 = 0.25 * rad * rad
227
228			! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
229			! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
230			! chaque aire_2d(i, j), ainsi que les quatre élongations
231			! élémentaires cuij et les quatre élongations cvij qui sont
232			! calculées aux mêmes endroits que les aireij.
233
234			coslatm = cos(rlatu1(1))
235			radclatm = 0.5 * rad * coslatm
236
237			aireij1_2d(:iim, 1) = 0.
238			aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
239			aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
240			aireij4_2d(:iim, 1) = 0.
241
242			cuij1(:iim, 1) = 0.
243			cuij2(:iim, 1) = radclatm * xprimp025(:iim)
244			cuij3(:iim, 1) = radclatm * xprimm025(:iim)
245			cuij4(:iim, 1) = 0.
246
247			cvij1(:iim, 1) = 0.
248			cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
249			cvij3(:iim, 1) = cvij2(:iim, 1)
250			cvij4(:iim, 1) = 0.
251
252			do j = 2, jjm
253			coslatm = cos(rlatu1(j))
254			coslatp = cos(rlatu2(j-1))
255			radclatp = 0.5 * rad * coslatp
256			radclatm = 0.5 * rad * coslatm
257			ai14 = un4rad2 * coslatp * yprimu2(j-1)
258			ai23 = un4rad2 * coslatm * yprimu1(j)
259
260			aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
261			aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
262			aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
263			aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
264			cuij1(:iim, j) = radclatp * xprimp025(:iim)
265			cuij2(:iim, j) = radclatm * xprimp025(:iim)
266			cuij3(:iim, j) = radclatm * xprimm025(:iim)
267			cuij4(:iim, j) = radclatp * xprimm025(:iim)
268			cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
269			cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
270			cvij3(:iim, j) = cvij2(:iim, j)
271			cvij4(:iim, j) = cvij1(:iim, j)
272			end do
273
274			coslatp = cos(rlatu2(jjm))
275			radclatp = 0.5 * rad * coslatp
276
277			aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
278			aireij2_2d(:iim, jjp1) = 0.
279			aireij3_2d(:iim, jjp1) = 0.
280			aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)
281
282			cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
283			cuij2(:iim, jjp1) = 0.
284			cuij3(:iim, jjp1) = 0.
285			cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)
286
287			cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
288			cvij2(:iim, jjp1) = 0.
289			cvij3(:iim, jjp1) = 0.
290			cvij4(:iim, jjp1) = cvij1(:iim, jjp1)
291
292			! Périodicité :
293
294			cvij1(iip1, :) = cvij1(1, :)
295			cvij2(iip1, :) = cvij2(1, :)
296			cvij3(iip1, :) = cvij3(1, :)
297			cvij4(iip1, :) = cvij4(1, :)
298
299			cuij1(iip1, :) = cuij1(1, :)
300			cuij2(iip1, :) = cuij2(1, :)
301			cuij3(iip1, :) = cuij3(1, :)
302			cuij4(iip1, :) = cuij4(1, :)
303
304			aireij1_2d(iip1, :) = aireij1_2d(1, :)
305			aireij2_2d(iip1, :) = aireij2_2d(1, :)
306			aireij3_2d(iip1, :) = aireij3_2d(1, :)
307			aireij4_2d(iip1, :) = aireij4_2d(1, :)
308
309			DO j = 1, jjp1
310			DO i = 1, iim
311			aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
312			+ aireij3_2d(i, j) + aireij4_2d(i, j)
313			alpha1_2d(i, j) = aireij1_2d(i, j) / aire_2d(i, j)
314			alpha2_2d(i, j) = aireij2_2d(i, j) / aire_2d(i, j)
315			alpha3_2d(i, j) = aireij3_2d(i, j) / aire_2d(i, j)
316			alpha4_2d(i, j) = aireij4_2d(i, j) / aire_2d(i, j)
317			alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
318			alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
319			alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)
320			alpha3p4_2d(i, j) = alpha3_2d(i, j) + alpha4_2d(i, j)
321			END DO
322
323			aire_2d(iip1, j) = aire_2d(1, j)
324			alpha1_2d(iip1, j) = alpha1_2d(1, j)
325			alpha2_2d(iip1, j) = alpha2_2d(1, j)
326			alpha3_2d(iip1, j) = alpha3_2d(1, j)
327			alpha4_2d(iip1, j) = alpha4_2d(1, j)
328			alpha1p2_2d(iip1, j) = alpha1p2_2d(1, j)
329			alpha1p4_2d(iip1, j) = alpha1p4_2d(1, j)
330			alpha2p3_2d(iip1, j) = alpha2p3_2d(1, j)
331			alpha3p4_2d(iip1, j) = alpha3p4_2d(1, j)
332			END DO
333
334			DO j = 1, jjp1
335			DO i = 1, iim
336			aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
337			aireij4_2d(i + 1, j) + aireij3_2d(i + 1, j)
338			unsaire_2d(i, j) = 1. / aire_2d(i, j)
339			unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
340			unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
341			airesurg_2d(i, j) = aire_2d(i, j) / g
342			END DO
343			aireu_2d(iip1, j) = aireu_2d(1, j)
344			unsaire_2d(iip1, j) = unsaire_2d(1, j)
345			unsair_gam1_2d(iip1, j) = unsair_gam1_2d(1, j)
346			unsair_gam2_2d(iip1, j) = unsair_gam2_2d(1, j)
347			airesurg_2d(iip1, j) = airesurg_2d(1, j)
348			END DO
349
350			DO j = 1, jjm
351			DO i = 1, iim
352			airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
353			aireij1_2d(i, j + 1) + aireij4_2d(i, j + 1)
354			END DO
355			DO i = 1, iim
356			airez = aireij2_2d(i, j) + aireij1_2d(i, j + 1) &
357			+ aireij3_2d(i + 1, j) + aireij4_2d(i + 1, j + 1)
358			unsairez_2d(i, j) = 1. / airez
359			unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
360			fext_2d(i, j) = airez * sin(rlatv(j)) * 2. * omeg
361			END DO
362			airev_2d(iip1, j) = airev_2d(1, j)
363			unsairez_2d(iip1, j) = unsairez_2d(1, j)
364			fext_2d(iip1, j) = fext_2d(1, j)
365			unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
366			END DO
367
368			! Calcul des élongations cu_2d, cv_2d, cvu
369
370			DO j = 1, jjm
371			DO i = 1, iim
372			cv_2d(i, j) = 0.5 * &
373			(cvij2(i, j) + cvij3(i, j) + cvij1(i, j + 1) + cvij4(i, j + 1))
374			cvu(i, j) = 0.5 * (cvij1(i, j) + cvij4(i, j) + cvij2(i, j) &
375			+ cvij3(i, j))
376			cuv(i, j) = 0.5 * (cuij2(i, j) + cuij3(i, j) + cuij1(i, j + 1) &
377			+ cuij4(i, j + 1))
378			unscv2_2d(i, j) = 1. / cv_2d(i, j)**2
379			END DO
380			DO i = 1, iim
381			cuvsurcv_2d(i, j) = airev_2d(i, j) * unscv2_2d(i, j)
382			cvsurcuv_2d(i, j) = 1. / cuvsurcv_2d(i, j)
383			cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
384			cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
385			cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
386			END DO
387			cv_2d(iip1, j) = cv_2d(1, j)
388			cvu(iip1, j) = cvu(1, j)
389			unscv2_2d(iip1, j) = unscv2_2d(1, j)
390			cuv(iip1, j) = cuv(1, j)
391			cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
392			cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
393			cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
394			cuvscvgam2_2d(iip1, j) = cuvscvgam2_2d(1, j)
395			cvscuvgam_2d(iip1, j) = cvscuvgam_2d(1, j)
396			END DO
397
398			DO j = 2, jjm
399			DO i = 1, iim
400			cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i + 1, j) + cuij2(i, j) &
401			+ cuij3(i + 1, j))
402			unscu2_2d(i, j) = 1. / cu_2d(i, j)**2
403			cvusurcu_2d(i, j) = aireu_2d(i, j) * unscu2_2d(i, j)
404			cusurcvu_2d(i, j) = 1. / cvusurcu_2d(i, j)
405			cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
406			cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
407			cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
408			END DO
409			cu_2d(iip1, j) = cu_2d(1, j)
410			unscu2_2d(iip1, j) = unscu2_2d(1, j)
411			cvusurcu_2d(iip1, j) = cvusurcu_2d(1, j)
412			cusurcvu_2d(iip1, j) = cusurcvu_2d(1, j)
413			cvuscugam1_2d(iip1, j) = cvuscugam1_2d(1, j)
414			cvuscugam2_2d(iip1, j) = cvuscugam2_2d(1, j)
415			cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
416			END DO
417
418			! Calcul aux pôles
419
420			cu_2d(:, 1) = 0.
421			unscu2_2d(:, 1) = 0.
422			cvu(:, 1) = 0.
423
424			cu_2d(:, jjp1) = 0.
425			unscu2_2d(:, jjp1) = 0.
426			cvu(:, jjp1) = 0.
427
428			DO j = 1, jjm
429			DO i = 1, iim
430			airvscu2_2d(i, j) = airev_2d(i, j) / (cuv(i, j) * cuv(i, j))
431			aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)
432			END DO
433			airvscu2_2d(iip1, j) = airvscu2_2d(1, j)
434			aivscu2gam_2d(iip1, j) = aivscu2gam_2d(1, j)
435			END DO
436
437			DO j = 2, jjm
438			DO i = 1, iim
439			airuscv2_2d(i, j) = aireu_2d(i, j) / (cvu(i, j) * cvu(i, j))
440			aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)
441			END DO
442			airuscv2_2d(iip1, j) = airuscv2_2d(1, j)
443			aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
444			END DO
445
446			! Calcul des aires aux pôles :
447
448			apoln = sum(aire_2d(:iim, 1))
449			apols = sum(aire_2d(:iim, jjp1))
450			unsapolnga1 = 1. / (apoln**(-gamdi_gdiv))
451			unsapolsga1 = 1. / (apols**(-gamdi_gdiv))
452			unsapolnga2 = 1. / (apoln**(-gamdi_h))
453			unsapolsga2 = 1. / (apols**(-gamdi_h))
454
455			! Changement F. Hourdin calcul conservatif pour fext_2d
456			! constang_2d contient le produit a * cos (latitude) * omega
457
458			DO i = 1, iim
459			constang_2d(i, 1) = 0.
460			END DO
461			DO j = 1, jjm - 1
462			DO i = 1, iim
463			constang_2d(i, j + 1) = rad * omeg * cu_2d(i, j + 1) &
464			* cos(rlatu(j + 1))
465			END DO
466			END DO
467			DO i = 1, iim
468			constang_2d(i, jjp1) = 0.
469			END DO
470
471			! Périodicité en longitude
472			DO j = 1, jjp1
473			constang_2d(iip1, j) = constang_2d(1, j)
474			END DO
475
476			call new_unit(unit)
477			open(unit, file="longitude_latitude.txt", status="replace", action="write")
478			write(unit, fmt=) '"longitudes at V points (degrees)"', rlonv 180. / pi
479			write(unit, fmt=) '"latitudes at V points (degrees)"', rlatv 180. / pi
480			write(unit, fmt=) '"longitudes at U points (degrees)"', rlonu 180. / pi
481			write(unit, fmt=) '"latitudes at U points (degrees)"', rlatu 180. / pi
482			close(unit)
483
484			END SUBROUTINE inigeom
485
486	guez	3	end module comgeom