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 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 / (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

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