trunk/dyn3d/comgeom.f

module comgeom

  use dimens_m, only: iim, jjm
  use paramet_m, only: ip1jmp1, ip1jm

  implicit none

  private iim, jjm, ip1jmp1, ip1jm

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

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

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

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

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

  real airev(ip1jm), airev_2d(iim + 1, jjm) ! in m2
  real unsaire(ip1jmp1), 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(ip1jm)
  equivalence (unsairez, unsairez_2d)

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

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

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

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

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

  real alpha3p4(ip1jmp1)
  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(ip1jm), constang(ip1jmp1)
  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(ip1jm), cvsurcuv(ip1jm) ! 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(ip1jmp1), cusurcvu(ip1jmp1) ! no dimension
  equivalence (cvusurcu, cvusurcu_2d), (cusurcvu, cusurcvu_2d)

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

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

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

  real cuscvugam(ip1jmp1)
  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(ip1jmp1), unsair_gam2(ip1jmp1)
  equivalence (unsair_gam1, unsair_gam1_2d), (unsair_gam2, unsair_gam2_2d)

  real unsairz_gam_2d(iim + 1, jjm)
  real unsairz_gam(ip1jm)
  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 dimens_m, ONLY : iim, jjm
    use fxy_m, only: fxy
    use fxyhyper_m, only: fxyhyper
    use jumble, only: new_unit
    use nr_util, only: pi
    USE paramet_m, ONLY : iip1, jjp1
    USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &
         grossismy, pxo, pyo, taux, tauy, transx, transy
    ! Modifies pxo, pyo, transx, transy

    ! Variables locales

    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 (.NOT. fxyhypb) THEN
       IF (ysinus) THEN
          print *, ' Inigeom, Y = Sinus (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
    ELSE
       ! Utilisation de fxyhyper, f(x, y) à dérivée tangente hyperbolique
       print *, 'Inigeom, Y = Latitude, dérivée tangente hyperbolique'
       CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &
            taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
            yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
            rlonp025, xprimp025)
    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, jjm
       fext_2d(iip1, j) = fext_2d(1, j)
    END DO
    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	use paramet_m, only: ip1jmp1, ip1jm
5
6	implicit none
7
8	private iim, jjm, ip1jmp1, ip1jm
9
10	real cu_2d(iim + 1, jjm + 1), cv_2d(iim + 1, jjm) ! in m
11	real cu(ip1jmp1), cv(ip1jm) ! in m
12	equivalence (cu, cu_2d), (cv, cv_2d)
13
14	real unscu2_2d(iim + 1, jjm + 1) ! in m-2
15	real unscu2(ip1jmp1) ! in m-2
16	equivalence (unscu2, unscu2_2d)
17
18	real unscv2_2d(iim + 1, jjm) ! in m-2
19	real unscv2(ip1jm) ! in m-2
20	equivalence (unscv2, unscv2_2d)
21
22	real aire(ip1jmp1), aire_2d(iim + 1, jjm + 1) ! in m2
23	real airesurg_2d(iim + 1, jjm + 1), airesurg(ip1jmp1)
24	equivalence (aire, aire_2d), (airesurg, airesurg_2d)
25
26	real aireu_2d(iim + 1, jjm + 1) ! in m2
27	real aireu(ip1jmp1) ! in m2
28	equivalence (aireu, aireu_2d)
29
30	real airev(ip1jm), airev_2d(iim + 1, jjm) ! in m2
31	real unsaire(ip1jmp1), unsaire_2d(iim + 1, jjm + 1) ! in m-2
32	equivalence (airev, airev_2d), (unsaire, unsaire_2d)
33
34	real apoln, apols ! in m2
35
36	real unsairez_2d(iim + 1, jjm)
37	real unsairez(ip1jm)
38	equivalence (unsairez, unsairez_2d)
39
40	real alpha1_2d(iim + 1, jjm + 1)
41	real alpha1(ip1jmp1)
42	equivalence (alpha1, alpha1_2d)
43
44	real alpha2_2d(iim + 1, jjm + 1)
45	real alpha2(ip1jmp1)
46	equivalence (alpha2, alpha2_2d)
47
48	real alpha3_2d(iim + 1, jjm + 1), alpha4_2d(iim + 1, jjm + 1)
49	real alpha3(ip1jmp1), alpha4(ip1jmp1)
50	equivalence (alpha3, alpha3_2d), (alpha4, alpha4_2d)
51
52	real alpha1p2_2d(iim + 1, jjm + 1)
53	real alpha1p2(ip1jmp1)
54	equivalence (alpha1p2, alpha1p2_2d)
55
56	real alpha1p4_2d(iim + 1, jjm + 1), alpha2p3_2d(iim + 1, jjm + 1)
57	real alpha1p4(ip1jmp1), alpha2p3(ip1jmp1)
58	equivalence (alpha1p4, alpha1p4_2d), (alpha2p3, alpha2p3_2d)
59
60	real alpha3p4(ip1jmp1)
61	real alpha3p4_2d(iim + 1, jjm + 1)
62	equivalence (alpha3p4, alpha3p4_2d)
63
64	real fext_2d(iim + 1, jjm), constang_2d(iim + 1, jjm + 1)
65	real fext(ip1jm), constang(ip1jmp1)
66	equivalence (fext, fext_2d), (constang, constang_2d)
67
68	real rlatu(jjm + 1)
69	! (latitudes of points of the "scalar" and "u" grid, in rad)
70
71	real rlatv(jjm)
72	! (latitudes of points of the "v" grid, in rad, in decreasing order)
73
74	real rlonu(iim + 1) ! longitudes of points of the "u" grid, in rad
75
76	real rlonv(iim + 1)
77	! (longitudes of points of the "scalar" and "v" grid, in rad)
78
79	real cuvsurcv_2d(iim + 1, jjm), cvsurcuv_2d(iim + 1, jjm) ! no dimension
80	real cuvsurcv(ip1jm), cvsurcuv(ip1jm) ! no dimension
81	equivalence (cuvsurcv, cuvsurcv_2d), (cvsurcuv, cvsurcuv_2d)
82
83	real cvusurcu_2d(iim + 1, jjm + 1), cusurcvu_2d(iim + 1, jjm + 1)
84	! no dimension
85	real cvusurcu(ip1jmp1), cusurcvu(ip1jmp1) ! no dimension
86	equivalence (cvusurcu, cvusurcu_2d), (cusurcvu, cusurcvu_2d)
87
88	real cuvscvgam1_2d(iim + 1, jjm)
89	real cuvscvgam1(ip1jm)
90	equivalence (cuvscvgam1, cuvscvgam1_2d)
91
92	real cuvscvgam2_2d(iim + 1, jjm), cvuscugam1_2d(iim + 1, jjm + 1)
93	real cuvscvgam2(ip1jm), cvuscugam1(ip1jmp1)
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(ip1jmp1), cvscuvgam(ip1jm)
98	equivalence (cvuscugam2, cvuscugam2_2d), (cvscuvgam, cvscuvgam_2d)
99
100	real cuscvugam(ip1jmp1)
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(ip1jmp1), unsair_gam2(ip1jmp1)
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(ip1jm)
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 dimens_m, ONLY : iim, jjm
168	use fxy_m, only: fxy
169	use fxyhyper_m, only: fxyhyper
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, dzoomx, dzoomy, grossismx, &
174	grossismy, pxo, pyo, taux, tauy, transx, transy
175	! Modifies pxo, pyo, transx, transy
176
177	! Variables locales
178
179	INTEGER i, j, itmax, itmay, iter, unit
180	REAL cvu(iip1, jjp1), cuv(iip1, jjm)
181	REAL ai14, ai23, airez, un4rad2
182	REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm
183	REAL coslatm, coslatp, radclatm, radclatp
184	REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
185	REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4 ! in m
186	REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
187	real yprimv(jjm), yprimu(jjp1)
188	REAL gamdi_gdiv, gamdi_grot, gamdi_h
189	REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)
190	real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &
191	aireij4_2d ! in m2
192	real airuscv2_2d(iim + 1, jjm)
193	real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)
194	real aivscu2gam_2d(iim + 1, jjm)
195
196	!------------------------------------------------------------------
197
198	PRINT *, 'Call sequence information: inigeom'
199
200	IF (nitergdiv/=2) THEN
201	gamdi_gdiv = coefdis / (real(nitergdiv)-2.)
202	ELSE
203	gamdi_gdiv = 0.
204	END IF
205	IF (nitergrot/=2) THEN
206	gamdi_grot = coefdis / (real(nitergrot)-2.)
207	ELSE
208	gamdi_grot = 0.
209	END IF
210	IF (niterh/=2) THEN
211	gamdi_h = coefdis / (real(niterh)-2.)
212	ELSE
213	gamdi_h = 0.
214	END IF
215
216	print *, 'gamdi_gdiv = ', gamdi_gdiv
217	print *, "gamdi_grot = ", gamdi_grot
218	print *, "gamdi_h = ", gamdi_h
219
220	IF (.NOT. fxyhypb) THEN
221	IF (ysinus) THEN
222	print *, ' Inigeom, Y = Sinus (Latitude) '
223	! utilisation de f(x, y) avec y = sinus de la latitude
224	CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &
225	rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &
226	xprimm025, rlonp025, xprimp025)
227	ELSE
228	print *, 'Inigeom, Y = Latitude, der. sinusoid .'
229	! utilisation de f(x, y) a tangente sinusoidale, y etant la latit
230
231	pxo = clon * pi / 180.
232	pyo = 2. * clat * pi / 180.
233
234	! determination de transx (pour le zoom) par Newton-Raphson
235
236	itmax = 10
237	eps = .1E-7
238
239	xo1 = 0.
240	DO iter = 1, itmax
241	x1 = xo1
242	f = x1 + alphax * sin(x1-pxo)
243	df = 1. + alphax * cos(x1-pxo)
244	x1 = x1 - f / df
245	xdm = abs(x1-xo1)
246	IF (xdm<=eps) EXIT
247	xo1 = x1
248	END DO
249
250	transx = xo1
251
252	itmay = 10
253	eps = .1E-7
254
255	yo1 = 0.
256	DO iter = 1, itmay
257	y1 = yo1
258	f = y1 + alphay * sin(y1-pyo)
259	df = 1. + alphay * cos(y1-pyo)
260	y1 = y1 - f / df
261	ydm = abs(y1-yo1)
262	IF (ydm<=eps) EXIT
263	yo1 = y1
264	END DO
265
266	transy = yo1
267
268	CALL fxy(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
269	yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
270	rlonp025, xprimp025)
271	END IF
272	ELSE
273	! Utilisation de fxyhyper, f(x, y) à dérivée tangente hyperbolique
274	print *, 'Inigeom, Y = Latitude, dérivée tangente hyperbolique'
275	CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &
276	taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
277	yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
278	rlonp025, xprimp025)
279	END IF
280
281	rlatu(1) = pi / 2.
282	rlatu(jjp1) = -rlatu(1)
283
284	! Calcul aux pôles
285
286	yprimu(1) = 0.
287	yprimu(jjp1) = 0.
288
289	un4rad2 = 0.25 * rad * rad
290
291	! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
292	! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
293	! chaque aire_2d(i, j), ainsi que les quatre élongations
294	! élémentaires cuij et les quatre élongations cvij qui sont
295	! calculées aux mêmes endroits que les aireij.
296
297	coslatm = cos(rlatu1(1))
298	radclatm = 0.5 * rad * coslatm
299
300	aireij1_2d(:iim, 1) = 0.
301	aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
302	aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
303	aireij4_2d(:iim, 1) = 0.
304
305	cuij1(:iim, 1) = 0.
306	cuij2(:iim, 1) = radclatm * xprimp025(:iim)
307	cuij3(:iim, 1) = radclatm * xprimm025(:iim)
308	cuij4(:iim, 1) = 0.
309
310	cvij1(:iim, 1) = 0.
311	cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
312	cvij3(:iim, 1) = cvij2(:iim, 1)
313	cvij4(:iim, 1) = 0.
314
315	do j = 2, jjm
316	coslatm = cos(rlatu1(j))
317	coslatp = cos(rlatu2(j-1))
318	radclatp = 0.5 * rad * coslatp
319	radclatm = 0.5 * rad * coslatm
320	ai14 = un4rad2 * coslatp * yprimu2(j-1)
321	ai23 = un4rad2 * coslatm * yprimu1(j)
322
323	aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
324	aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
325	aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
326	aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
327	cuij1(:iim, j) = radclatp * xprimp025(:iim)
328	cuij2(:iim, j) = radclatm * xprimp025(:iim)
329	cuij3(:iim, j) = radclatm * xprimm025(:iim)
330	cuij4(:iim, j) = radclatp * xprimm025(:iim)
331	cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
332	cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
333	cvij3(:iim, j) = cvij2(:iim, j)
334	cvij4(:iim, j) = cvij1(:iim, j)
335	end do
336
337	coslatp = cos(rlatu2(jjm))
338	radclatp = 0.5 * rad * coslatp
339
340	aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
341	aireij2_2d(:iim, jjp1) = 0.
342	aireij3_2d(:iim, jjp1) = 0.
343	aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)
344
345	cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
346	cuij2(:iim, jjp1) = 0.
347	cuij3(:iim, jjp1) = 0.
348	cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)
349
350	cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
351	cvij2(:iim, jjp1) = 0.
352	cvij3(:iim, jjp1) = 0.
353	cvij4(:iim, jjp1) = cvij1(:iim, jjp1)
354
355	! Périodicité :
356
357	cvij1(iip1, :) = cvij1(1, :)
358	cvij2(iip1, :) = cvij2(1, :)
359	cvij3(iip1, :) = cvij3(1, :)
360	cvij4(iip1, :) = cvij4(1, :)
361
362	cuij1(iip1, :) = cuij1(1, :)
363	cuij2(iip1, :) = cuij2(1, :)
364	cuij3(iip1, :) = cuij3(1, :)
365	cuij4(iip1, :) = cuij4(1, :)
366
367	aireij1_2d(iip1, :) = aireij1_2d(1, :)
368	aireij2_2d(iip1, :) = aireij2_2d(1, :)
369	aireij3_2d(iip1, :) = aireij3_2d(1, :)
370	aireij4_2d(iip1, :) = aireij4_2d(1, :)
371
372	DO j = 1, jjp1
373	DO i = 1, iim
374	aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
375	+ aireij3_2d(i, j) + aireij4_2d(i, j)
376	alpha1_2d(i, j) = aireij1_2d(i, j) / aire_2d(i, j)
377	alpha2_2d(i, j) = aireij2_2d(i, j) / aire_2d(i, j)
378	alpha3_2d(i, j) = aireij3_2d(i, j) / aire_2d(i, j)
379	alpha4_2d(i, j) = aireij4_2d(i, j) / aire_2d(i, j)
380	alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
381	alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
382	alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)
383	alpha3p4_2d(i, j) = alpha3_2d(i, j) + alpha4_2d(i, j)
384	END DO
385
386	aire_2d(iip1, j) = aire_2d(1, j)
387	alpha1_2d(iip1, j) = alpha1_2d(1, j)
388	alpha2_2d(iip1, j) = alpha2_2d(1, j)
389	alpha3_2d(iip1, j) = alpha3_2d(1, j)
390	alpha4_2d(iip1, j) = alpha4_2d(1, j)
391	alpha1p2_2d(iip1, j) = alpha1p2_2d(1, j)
392	alpha1p4_2d(iip1, j) = alpha1p4_2d(1, j)
393	alpha2p3_2d(iip1, j) = alpha2p3_2d(1, j)
394	alpha3p4_2d(iip1, j) = alpha3p4_2d(1, j)
395	END DO
396
397	DO j = 1, jjp1
398	DO i = 1, iim
399	aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
400	aireij4_2d(i + 1, j) + aireij3_2d(i + 1, j)
401	unsaire_2d(i, j) = 1. / aire_2d(i, j)
402	unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
403	unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
404	airesurg_2d(i, j) = aire_2d(i, j) / g
405	END DO
406	aireu_2d(iip1, j) = aireu_2d(1, j)
407	unsaire_2d(iip1, j) = unsaire_2d(1, j)
408	unsair_gam1_2d(iip1, j) = unsair_gam1_2d(1, j)
409	unsair_gam2_2d(iip1, j) = unsair_gam2_2d(1, j)
410	airesurg_2d(iip1, j) = airesurg_2d(1, j)
411	END DO
412
413	DO j = 1, jjm
414	DO i = 1, iim
415	airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
416	aireij1_2d(i, j + 1) + aireij4_2d(i, j + 1)
417	END DO
418	DO i = 1, iim
419	airez = aireij2_2d(i, j) + aireij1_2d(i, j + 1) &
420	+ aireij3_2d(i + 1, j) + aireij4_2d(i + 1, j + 1)
421	unsairez_2d(i, j) = 1. / airez
422	unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
423	fext_2d(i, j) = airez * sin(rlatv(j)) * 2. * omeg
424	END DO
425	airev_2d(iip1, j) = airev_2d(1, j)
426	unsairez_2d(iip1, j) = unsairez_2d(1, j)
427	fext_2d(iip1, j) = fext_2d(1, j)
428	unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
429	END DO
430
431	! Calcul des élongations cu_2d, cv_2d, cvu
432
433	DO j = 1, jjm
434	DO i = 1, iim
435	cv_2d(i, j) = 0.5 * &
436	(cvij2(i, j) + cvij3(i, j) + cvij1(i, j + 1) + cvij4(i, j + 1))
437	cvu(i, j) = 0.5 * (cvij1(i, j) + cvij4(i, j) + cvij2(i, j) &
438	+ cvij3(i, j))
439	cuv(i, j) = 0.5 * (cuij2(i, j) + cuij3(i, j) + cuij1(i, j + 1) &
440	+ cuij4(i, j + 1))
441	unscv2_2d(i, j) = 1. / cv_2d(i, j)**2
442	END DO
443	DO i = 1, iim
444	cuvsurcv_2d(i, j) = airev_2d(i, j) * unscv2_2d(i, j)
445	cvsurcuv_2d(i, j) = 1. / cuvsurcv_2d(i, j)
446	cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
447	cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
448	cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
449	END DO
450	cv_2d(iip1, j) = cv_2d(1, j)
451	cvu(iip1, j) = cvu(1, j)
452	unscv2_2d(iip1, j) = unscv2_2d(1, j)
453	cuv(iip1, j) = cuv(1, j)
454	cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
455	cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
456	cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
457	cuvscvgam2_2d(iip1, j) = cuvscvgam2_2d(1, j)
458	cvscuvgam_2d(iip1, j) = cvscuvgam_2d(1, j)
459	END DO
460
461	DO j = 2, jjm
462	DO i = 1, iim
463	cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i + 1, j) + cuij2(i, j) &
464	+ cuij3(i + 1, j))
465	unscu2_2d(i, j) = 1. / cu_2d(i, j)**2
466	cvusurcu_2d(i, j) = aireu_2d(i, j) * unscu2_2d(i, j)
467	cusurcvu_2d(i, j) = 1. / cvusurcu_2d(i, j)
468	cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
469	cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
470	cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
471	END DO
472	cu_2d(iip1, j) = cu_2d(1, j)
473	unscu2_2d(iip1, j) = unscu2_2d(1, j)
474	cvusurcu_2d(iip1, j) = cvusurcu_2d(1, j)
475	cusurcvu_2d(iip1, j) = cusurcvu_2d(1, j)
476	cvuscugam1_2d(iip1, j) = cvuscugam1_2d(1, j)
477	cvuscugam2_2d(iip1, j) = cvuscugam2_2d(1, j)
478	cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
479	END DO
480
481	! Calcul aux pôles
482
483	cu_2d(:, 1) = 0.
484	unscu2_2d(:, 1) = 0.
485	cvu(:, 1) = 0.
486
487	cu_2d(:, jjp1) = 0.
488	unscu2_2d(:, jjp1) = 0.
489	cvu(:, jjp1) = 0.
490
491	DO j = 1, jjm
492	DO i = 1, iim
493	airvscu2_2d(i, j) = airev_2d(i, j) / (cuv(i, j) * cuv(i, j))
494	aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)
495	END DO
496	airvscu2_2d(iip1, j) = airvscu2_2d(1, j)
497	aivscu2gam_2d(iip1, j) = aivscu2gam_2d(1, j)
498	END DO
499
500	DO j = 2, jjm
501	DO i = 1, iim
502	airuscv2_2d(i, j) = aireu_2d(i, j) / (cvu(i, j) * cvu(i, j))
503	aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)
504	END DO
505	airuscv2_2d(iip1, j) = airuscv2_2d(1, j)
506	aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
507	END DO
508
509	! Calcul des aires aux pôles :
510
511	apoln = sum(aire_2d(:iim, 1))
512	apols = sum(aire_2d(:iim, jjp1))
513	unsapolnga1 = 1. / (apoln**(-gamdi_gdiv))
514	unsapolsga1 = 1. / (apols**(-gamdi_gdiv))
515	unsapolnga2 = 1. / (apoln**(-gamdi_h))
516	unsapolsga2 = 1. / (apols**(-gamdi_h))
517
518	! Changement F. Hourdin calcul conservatif pour fext_2d
519	! constang_2d contient le produit a * cos (latitude) * omega
520
521	DO i = 1, iim
522	constang_2d(i, 1) = 0.
523	END DO
524	DO j = 1, jjm - 1
525	DO i = 1, iim
526	constang_2d(i, j + 1) = rad * omeg * cu_2d(i, j + 1) &
527	* cos(rlatu(j + 1))
528	END DO
529	END DO
530	DO i = 1, iim
531	constang_2d(i, jjp1) = 0.
532	END DO
533
534	! Périodicité en longitude
535
536	DO j = 1, jjm
537	fext_2d(iip1, j) = fext_2d(1, j)
538	END DO
539	DO j = 1, jjp1
540	constang_2d(iip1, j) = constang_2d(1, j)
541	END DO
542
543	call new_unit(unit)
544	open(unit, file="longitude_latitude.txt", status="replace", action="write")
545	write(unit, fmt=) '"longitudes at V points (degrees)"', rlonv 180. / pi
546	write(unit, fmt=) '"latitudes at V points (degrees)"', rlatv 180. / pi
547	write(unit, fmt=) '"longitudes at U points (degrees)"', rlonu 180. / pi
548	write(unit, fmt=) '"latitudes at U points (degrees)"', rlatu 180. / pi
549	close(unit)
550
551	END SUBROUTINE inigeom
552
553	end module comgeom