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