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 fxyhyper_m, only: fxyhyper
    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 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

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