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	module comgeom
2
3	use dimens_m, only: iim, jjm
4
5	implicit none
6
7	private iim, jjm
8
9	real cu_2d(iim + 1, jjm + 1), cv_2d(iim + 1, jjm) ! in m
10	real cu((iim + 1) * (jjm + 1)), cv((iim + 1) * jjm) ! in m
11	equivalence (cu, cu_2d), (cv, cv_2d)
12
13	real unscu2_2d(iim + 1, jjm + 1) ! in m-2
14	real unscu2((iim + 1) * (jjm + 1)) ! in m-2
15	equivalence (unscu2, unscu2_2d)
16
17	real unscv2_2d(iim + 1, jjm) ! in m-2
18	real unscv2((iim + 1) * jjm) ! in m-2
19	equivalence (unscv2, unscv2_2d)
20
21	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	equivalence (aire, aire_2d), (airesurg, airesurg_2d)
24
25	real aireu_2d(iim + 1, jjm + 1) ! in m2
26	real aireu((iim + 1) * (jjm + 1)) ! in m2
27	equivalence (aireu, aireu_2d)
28
29	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	equivalence (airev, airev_2d), (unsaire, unsaire_2d)
32
33	real apoln, apols ! in m2
34
35	real unsairez_2d(iim + 1, jjm)
36	real unsairez((iim + 1) * jjm)
37	equivalence (unsairez, unsairez_2d)
38
39	real alpha1_2d(iim + 1, jjm + 1)
40	real alpha1((iim + 1) * (jjm + 1))
41	equivalence (alpha1, alpha1_2d)
42
43	real alpha2_2d(iim + 1, jjm + 1)
44	real alpha2((iim + 1) * (jjm + 1))
45	equivalence (alpha2, alpha2_2d)
46
47	real alpha3_2d(iim + 1, jjm + 1), alpha4_2d(iim + 1, jjm + 1)
48	real alpha3((iim + 1) * (jjm + 1)), alpha4((iim + 1) * (jjm + 1))
49	equivalence (alpha3, alpha3_2d), (alpha4, alpha4_2d)
50
51	real alpha1p2_2d(iim + 1, jjm + 1)
52	real alpha1p2((iim + 1) * (jjm + 1))
53	equivalence (alpha1p2, alpha1p2_2d)
54
55	real alpha1p4_2d(iim + 1, jjm + 1), alpha2p3_2d(iim + 1, jjm + 1)
56	real alpha1p4((iim + 1) * (jjm + 1)), alpha2p3((iim + 1) * (jjm + 1))
57	equivalence (alpha1p4, alpha1p4_2d), (alpha2p3, alpha2p3_2d)
58
59	real alpha3p4((iim + 1) * (jjm + 1))
60	real alpha3p4_2d(iim + 1, jjm + 1)
61	equivalence (alpha3p4, alpha3p4_2d)
62
63	real fext_2d(iim + 1, jjm), constang_2d(iim + 1, jjm + 1)
64	real fext((iim + 1) * jjm), constang((iim + 1) * (jjm + 1))
65	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	real cuvsurcv_2d(iim + 1, jjm), cvsurcuv_2d(iim + 1, jjm) ! no dimension
79	real cuvsurcv((iim + 1) * jjm), cvsurcuv((iim + 1) * jjm) ! no dimension
80	equivalence (cuvsurcv, cuvsurcv_2d), (cvsurcuv, cvsurcuv_2d)
81
82	real cvusurcu_2d(iim + 1, jjm + 1), cusurcvu_2d(iim + 1, jjm + 1)
83	! no dimension
84	real cvusurcu((iim + 1) * (jjm + 1)), cusurcvu((iim + 1) * (jjm + 1))
85	! no dimension
86	equivalence (cvusurcu, cvusurcu_2d), (cusurcvu, cusurcvu_2d)
87
88	real cuvscvgam1_2d(iim + 1, jjm)
89	real cuvscvgam1((iim + 1) * jjm)
90	equivalence (cuvscvgam1, cuvscvgam1_2d)
91
92	real cuvscvgam2_2d(iim + 1, jjm), cvuscugam1_2d(iim + 1, jjm + 1)
93	real cuvscvgam2((iim + 1) * jjm), cvuscugam1((iim + 1) * (jjm + 1))
94	equivalence (cuvscvgam2, cuvscvgam2_2d), (cvuscugam1, cvuscugam1_2d)
95
96	real cvuscugam2_2d(iim + 1, jjm + 1), cvscuvgam_2d(iim + 1, jjm)
97	real cvuscugam2((iim + 1) * (jjm + 1)), cvscuvgam((iim + 1) * jjm)
98	equivalence (cvuscugam2, cvuscugam2_2d), (cvscuvgam, cvscuvgam_2d)
99
100	real cuscvugam((iim + 1) * (jjm + 1))
101	real cuscvugam_2d(iim + 1, jjm + 1)
102	equivalence (cuscvugam, cuscvugam_2d)
103
104	real unsapolnga1, unsapolnga2, unsapolsga1, unsapolsga2
105
106	real unsair_gam1_2d(iim + 1, jjm + 1), unsair_gam2_2d(iim + 1, jjm + 1)
107	real unsair_gam1((iim + 1) * (jjm + 1)), unsair_gam2((iim + 1) * (jjm + 1))
108	equivalence (unsair_gam1, unsair_gam1_2d), (unsair_gam2, unsair_gam2_2d)
109
110	real unsairz_gam_2d(iim + 1, jjm)
111	real unsairz_gam((iim + 1) * jjm)
112	equivalence (unsairz_gam, unsairz_gam_2d)
113
114	real xprimu(iim + 1), xprimv(iim + 1)
115
116	save
117
118	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	! Local:
172	INTEGER i, j, unit
173	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	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
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	end module comgeom