trunk/dyn3d/comgeom.f

module comgeom

  use dimens_m, only: iim, jjm

  implicit none

  private iim, jjm

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

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

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

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

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

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

  real apoln, apols ! in m2

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  real unsapolnga1, unsapolnga2, unsapolsga1, unsapolsga2                

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

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

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

  save

contains

  SUBROUTINE inigeom

    ! Auteur : P. Le Van

    ! Calcul des élongations cuij1, ..., cuij4, cvij1, ..., cvij4 aux mêmes
    ! endroits que les aires aireij1_2d, ..., aireij4_2d.

    ! Choix entre une fonction "f(y)" à dérivée sinusoïdale ou à
    ! dérivée tangente hyperbolique. Calcul des coefficients cu_2d,
    ! cv_2d, 1. / cu_2d**2, 1. / cv_2d**2. Les coefficients cu_2d et cv_2d
    ! permettent de passer des vitesses naturelles aux vitesses
    ! covariantes et contravariantes, ou vice-versa.

    ! On a :
    ! u(covariant) = cu_2d * u(naturel), u(contravariant) = u(naturel) / cu_2d
    ! v(covariant) = cv_2d * v(naturel), v(contravariant) = v(naturel) / cv_2d

    ! On en tire : 
    ! u(covariant) = cu_2d * cu_2d * u(contravariant)
    ! v(covariant) = cv_2d * cv_2d * v(contravariant)

    ! On a l'application (x(X), y(Y)) avec - im / 2 + 1 <= X <= im / 2
    ! et - jm / 2 <= Y <= jm / 2

    ! x est la longitude du point en radians.
    ! y est la latitude du point en radians.
    ! 
    ! On a : cu_2d(i, j) = rad * cos(y) * dx / dX
    ! cv(j) = rad * dy / dY
    ! aire_2d(i, j) = cu_2d(i, j) * cv(j)
    ! 
    ! y, dx / dX, dy / dY calculés aux points concernés. cv, bien que
    ! dépendant de j uniquement, sera ici indicé aussi en i pour un
    ! adressage plus facile en ij.

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

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

    USE comconst, ONLY : g, omeg, rad
    USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
    use conf_gcm_m, ONLY : fxyhypb, ysinus
    use fxy_m, only: fxy
    use fxyhyper_m, only: fxyhyper
    use fxysinus_m, only: fxysinus
    use jumble, only: new_unit
    use nr_util, only: pi
    USE paramet_m, ONLY : iip1, jjp1
    USE serre, ONLY : alphax, alphay, clat, clon, pxo, pyo, transx, transy
    ! Modifiés pxo, pyo, transx, transy

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

    !------------------------------------------------------------------

    PRINT *, 'Call sequence information: inigeom'

    IF (nitergdiv/=2) THEN
       gamdi_gdiv = coefdis / (real(nitergdiv)-2.)
    ELSE
       gamdi_gdiv = 0.
    END IF
    IF (nitergrot/=2) THEN
       gamdi_grot = coefdis / (real(nitergrot)-2.)
    ELSE
       gamdi_grot = 0.
    END IF
    IF (niterh/=2) THEN
       gamdi_h = coefdis / (real(niterh)-2.)
    ELSE
       gamdi_h = 0.
    END IF

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

    IF (fxyhypb) THEN
       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)
    ELSE
       IF (ysinus) THEN
          print *, 'inigeom: Y = sin(latitude)'
          ! Utilisation de f(x, y) avec y = sinus de la latitude 
          CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &
               rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &
               xprimm025, rlonp025, xprimp025)
       ELSE
          print *, 'Inigeom, Y = Latitude, der. sinusoid .'
          ! utilisation de f(x, y) a tangente sinusoidale, y etant la latit

          pxo = clon * pi / 180.
          pyo = 2. * clat * pi / 180.

          ! determination de transx (pour le zoom) par Newton-Raphson

          itmax = 10
          eps = .1E-7

          xo1 = 0.
          DO iter = 1, itmax
             x1 = xo1
             f = x1 + alphax * sin(x1-pxo)
             df = 1. + alphax * cos(x1-pxo)
             x1 = x1 - f / df
             xdm = abs(x1-xo1)
             IF (xdm<=eps) EXIT
             xo1 = x1
          END DO

          transx = xo1

          itmay = 10
          eps = .1E-7

          yo1 = 0.
          DO iter = 1, itmay
             y1 = yo1
             f = y1 + alphay * sin(y1-pyo)
             df = 1. + alphay * cos(y1-pyo)
             y1 = y1 - f / df
             ydm = abs(y1-yo1)
             IF (ydm<=eps) EXIT
             yo1 = y1
          END DO

          transy = yo1

          CALL fxy(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
               yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
               rlonp025, xprimp025)
       END IF
    END IF

    rlatu(1) = pi / 2.
    rlatu(jjp1) = -rlatu(1)

    ! Calcul aux pôles 

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

    un4rad2 = 0.25 * rad * rad

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

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

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

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

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

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

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

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

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

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

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

    ! Périodicité :

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

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

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

    DO j = 1, jjp1
       DO i = 1, iim
          aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
               + aireij3_2d(i, j) + aireij4_2d(i, j)
          alpha1_2d(i, j) = aireij1_2d(i, j) / aire_2d(i, j)
          alpha2_2d(i, j) = aireij2_2d(i, j) / aire_2d(i, j)
          alpha3_2d(i, j) = aireij3_2d(i, j) / aire_2d(i, j)
          alpha4_2d(i, j) = aireij4_2d(i, j) / aire_2d(i, j)
          alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
          alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
          alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)
          alpha3p4_2d(i, j) = alpha3_2d(i, j) + alpha4_2d(i, j)
       END DO

       aire_2d(iip1, j) = aire_2d(1, j)
       alpha1_2d(iip1, j) = alpha1_2d(1, j)
       alpha2_2d(iip1, j) = alpha2_2d(1, j)
       alpha3_2d(iip1, j) = alpha3_2d(1, j)
       alpha4_2d(iip1, j) = alpha4_2d(1, j)
       alpha1p2_2d(iip1, j) = alpha1p2_2d(1, j)
       alpha1p4_2d(iip1, j) = alpha1p4_2d(1, j)
       alpha2p3_2d(iip1, j) = alpha2p3_2d(1, j)
       alpha3p4_2d(iip1, j) = alpha3p4_2d(1, j)
    END DO

    DO j = 1, jjp1
       DO i = 1, iim
          aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
               aireij4_2d(i + 1, j) + aireij3_2d(i + 1, j)
          unsaire_2d(i, j) = 1. / aire_2d(i, j)
          unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
          unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
          airesurg_2d(i, j) = aire_2d(i, j) / g
       END DO
       aireu_2d(iip1, j) = aireu_2d(1, j)
       unsaire_2d(iip1, j) = unsaire_2d(1, j)
       unsair_gam1_2d(iip1, j) = unsair_gam1_2d(1, j)
       unsair_gam2_2d(iip1, j) = unsair_gam2_2d(1, j)
       airesurg_2d(iip1, j) = airesurg_2d(1, j)
    END DO

    DO j = 1, jjm
       DO i = 1, iim
          airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
               aireij1_2d(i, j + 1) + aireij4_2d(i, j + 1)
       END DO
       DO i = 1, iim
          airez = aireij2_2d(i, j) + aireij1_2d(i, j + 1) &
               + aireij3_2d(i + 1, j) + aireij4_2d(i + 1, j + 1)
          unsairez_2d(i, j) = 1. / airez
          unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
          fext_2d(i, j) = airez * sin(rlatv(j)) * 2. * omeg
       END DO
       airev_2d(iip1, j) = airev_2d(1, j)
       unsairez_2d(iip1, j) = unsairez_2d(1, j)
       fext_2d(iip1, j) = fext_2d(1, j)
       unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
    END DO

    ! Calcul des élongations cu_2d, cv_2d, cvu 

    DO j = 1, jjm
       DO i = 1, iim
          cv_2d(i, j) = 0.5 * &
               (cvij2(i, j) + cvij3(i, j) + cvij1(i, j + 1) + cvij4(i, j + 1))
          cvu(i, j) = 0.5 * (cvij1(i, j) + cvij4(i, j) + cvij2(i, j) &
               + cvij3(i, j))
          cuv(i, j) = 0.5 * (cuij2(i, j) + cuij3(i, j) + cuij1(i, j + 1) &
               + cuij4(i, j + 1))
          unscv2_2d(i, j) = 1. / cv_2d(i, j)**2
       END DO
       DO i = 1, iim
          cuvsurcv_2d(i, j) = airev_2d(i, j) * unscv2_2d(i, j)
          cvsurcuv_2d(i, j) = 1. / cuvsurcv_2d(i, j)
          cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
          cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
          cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
       END DO
       cv_2d(iip1, j) = cv_2d(1, j)
       cvu(iip1, j) = cvu(1, j)
       unscv2_2d(iip1, j) = unscv2_2d(1, j)
       cuv(iip1, j) = cuv(1, j)
       cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
       cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
       cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
       cuvscvgam2_2d(iip1, j) = cuvscvgam2_2d(1, j)
       cvscuvgam_2d(iip1, j) = cvscuvgam_2d(1, j)
    END DO

    DO j = 2, jjm
       DO i = 1, iim
          cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i + 1, j) + cuij2(i, j) &
               + cuij3(i + 1, j))
          unscu2_2d(i, j) = 1. / cu_2d(i, j)**2
          cvusurcu_2d(i, j) = aireu_2d(i, j) * unscu2_2d(i, j)
          cusurcvu_2d(i, j) = 1. / cvusurcu_2d(i, j)
          cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
          cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
          cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
       END DO
       cu_2d(iip1, j) = cu_2d(1, j)
       unscu2_2d(iip1, j) = unscu2_2d(1, j)
       cvusurcu_2d(iip1, j) = cvusurcu_2d(1, j)
       cusurcvu_2d(iip1, j) = cusurcvu_2d(1, j)
       cvuscugam1_2d(iip1, j) = cvuscugam1_2d(1, j)
       cvuscugam2_2d(iip1, j) = cvuscugam2_2d(1, j)
       cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
    END DO

    ! Calcul aux pôles 

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

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

    DO j = 1, jjm
       DO i = 1, iim
          airvscu2_2d(i, j) = airev_2d(i, j) / (cuv(i, j) * cuv(i, j))
          aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)
       END DO
       airvscu2_2d(iip1, j) = airvscu2_2d(1, j)
       aivscu2gam_2d(iip1, j) = aivscu2gam_2d(1, j)
    END DO

    DO j = 2, jjm
       DO i = 1, iim
          airuscv2_2d(i, j) = aireu_2d(i, j) / (cvu(i, j) * cvu(i, j))
          aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)
       END DO
       airuscv2_2d(iip1, j) = airuscv2_2d(1, j)
       aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
    END DO

    ! Calcul des aires aux pôles :

    apoln = sum(aire_2d(:iim, 1))
    apols = sum(aire_2d(:iim, jjp1))
    unsapolnga1 = 1. / (apoln**(-gamdi_gdiv))
    unsapolsga1 = 1. / (apols**(-gamdi_gdiv))
    unsapolnga2 = 1. / (apoln**(-gamdi_h))
    unsapolsga2 = 1. / (apols**(-gamdi_h))

    ! Changement F. Hourdin calcul conservatif pour fext_2d
    ! constang_2d contient le produit a * cos (latitude) * omega

    DO i = 1, iim
       constang_2d(i, 1) = 0.
    END DO
    DO j = 1, jjm - 1
       DO i = 1, iim
          constang_2d(i, j + 1) = rad * omeg * cu_2d(i, j + 1) &
               * cos(rlatu(j + 1))
       END DO
    END DO
    DO i = 1, iim
       constang_2d(i, jjp1) = 0.
    END DO

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

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

  END SUBROUTINE inigeom

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