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	guez	3	module comgeom
2
3			use dimens_m, only: iim, jjm
4
5			implicit none
6
7	guez	79	private iim, jjm
8	guez	3
9	guez	60	real cu_2d(iim + 1, jjm + 1), cv_2d(iim + 1, jjm) ! in m
10	guez	79	real cu((iim + 1) * (jjm + 1)), cv((iim + 1) * jjm) ! in m
11	guez	3	equivalence (cu, cu_2d), (cv, cv_2d)
12
13	guez	60	real unscu2_2d(iim + 1, jjm + 1) ! in m-2
14	guez	79	real unscu2((iim + 1) * (jjm + 1)) ! in m-2
15	guez	3	equivalence (unscu2, unscu2_2d)
16
17	guez	60	real unscv2_2d(iim + 1, jjm) ! in m-2
18	guez	79	real unscv2((iim + 1) * jjm) ! in m-2
19	guez	3	equivalence (unscv2, unscv2_2d)
20
21	guez	79	real aire((iim + 1) * (jjm + 1)), aire_2d(iim + 1, jjm + 1) ! in m2
22			real airesurg_2d(iim + 1, jjm + 1), airesurg((iim + 1) * (jjm + 1))
23	guez	3	equivalence (aire, aire_2d), (airesurg, airesurg_2d)
24
25	guez	57	real aireu_2d(iim + 1, jjm + 1) ! in m2
26	guez	79	real aireu((iim + 1) * (jjm + 1)) ! in m2
27	guez	3	equivalence (aireu, aireu_2d)
28
29	guez	79	real airev((iim + 1) * jjm), airev_2d(iim + 1, jjm) ! in m2
30			real unsaire((iim + 1) * (jjm + 1)), unsaire_2d(iim + 1, jjm + 1) ! in m-2
31	guez	3	equivalence (airev, airev_2d), (unsaire, unsaire_2d)
32
33	guez	60	real apoln, apols ! in m2
34	guez	3
35	guez	46	real unsairez_2d(iim + 1, jjm)
36	guez	79	real unsairez((iim + 1) * jjm)
37	guez	25	equivalence (unsairez, unsairez_2d)
38	guez	3
39	guez	46	real alpha1_2d(iim + 1, jjm + 1)
40	guez	79	real alpha1((iim + 1) * (jjm + 1))
41	guez	25	equivalence (alpha1, alpha1_2d)
42	guez	3
43	guez	46	real alpha2_2d(iim + 1, jjm + 1)
44	guez	79	real alpha2((iim + 1) * (jjm + 1))
45	guez	3	equivalence (alpha2, alpha2_2d)
46
47	guez	46	real alpha3_2d(iim + 1, jjm + 1), alpha4_2d(iim + 1, jjm + 1)
48	guez	79	real alpha3((iim + 1) * (jjm + 1)), alpha4((iim + 1) * (jjm + 1))
49	guez	3	equivalence (alpha3, alpha3_2d), (alpha4, alpha4_2d)
50
51	guez	46	real alpha1p2_2d(iim + 1, jjm + 1)
52	guez	79	real alpha1p2((iim + 1) * (jjm + 1))
53	guez	3	equivalence (alpha1p2, alpha1p2_2d)
54
55	guez	46	real alpha1p4_2d(iim + 1, jjm + 1), alpha2p3_2d(iim + 1, jjm + 1)
56	guez	79	real alpha1p4((iim + 1) * (jjm + 1)), alpha2p3((iim + 1) * (jjm + 1))
57	guez	3	equivalence (alpha1p4, alpha1p4_2d), (alpha2p3, alpha2p3_2d)
58
59	guez	79	real alpha3p4((iim + 1) * (jjm + 1))
60	guez	46	real alpha3p4_2d(iim + 1, jjm + 1)
61	guez	3	equivalence (alpha3p4, alpha3p4_2d)
62
63	guez	46	real fext_2d(iim + 1, jjm), constang_2d(iim + 1, jjm + 1)
64	guez	79	real fext((iim + 1) * jjm), constang((iim + 1) * (jjm + 1))
65	guez	3	equivalence (fext, fext_2d), (constang, constang_2d)
66
67			real rlatu(jjm + 1)
68			! (latitudes of points of the "scalar" and "u" grid, in rad)
69
70			real rlatv(jjm)
71			! (latitudes of points of the "v" grid, in rad, in decreasing order)
72
73			real rlonu(iim + 1) ! longitudes of points of the "u" grid, in rad
74
75			real rlonv(iim + 1)
76			! (longitudes of points of the "scalar" and "v" grid, in rad)
77
78	guez	60	real cuvsurcv_2d(iim + 1, jjm), cvsurcuv_2d(iim + 1, jjm) ! no dimension
79	guez	79	real cuvsurcv((iim + 1) * jjm), cvsurcuv((iim + 1) * jjm) ! no dimension
80	guez	3	equivalence (cuvsurcv, cuvsurcv_2d), (cvsurcuv, cvsurcuv_2d)
81
82	guez	46	real cvusurcu_2d(iim + 1, jjm + 1), cusurcvu_2d(iim + 1, jjm + 1)
83	guez	60	! no dimension
84	guez	79	real cvusurcu((iim + 1) * (jjm + 1)), cusurcvu((iim + 1) * (jjm + 1))
85			! no dimension
86	guez	3	equivalence (cvusurcu, cvusurcu_2d), (cusurcvu, cusurcvu_2d)
87
88	guez	46	real cuvscvgam1_2d(iim + 1, jjm)
89	guez	79	real cuvscvgam1((iim + 1) * jjm)
90	guez	3	equivalence (cuvscvgam1, cuvscvgam1_2d)
91
92	guez	46	real cuvscvgam2_2d(iim + 1, jjm), cvuscugam1_2d(iim + 1, jjm + 1)
93	guez	79	real cuvscvgam2((iim + 1) * jjm), cvuscugam1((iim + 1) * (jjm + 1))
94	guez	3	equivalence (cuvscvgam2, cuvscvgam2_2d), (cvuscugam1, cvuscugam1_2d)
95
96	guez	46	real cvuscugam2_2d(iim + 1, jjm + 1), cvscuvgam_2d(iim + 1, jjm)
97	guez	79	real cvuscugam2((iim + 1) * (jjm + 1)), cvscuvgam((iim + 1) * jjm)
98	guez	3	equivalence (cvuscugam2, cvuscugam2_2d), (cvscuvgam, cvscuvgam_2d)
99
100	guez	79	real cuscvugam((iim + 1) * (jjm + 1))
101	guez	46	real cuscvugam_2d(iim + 1, jjm + 1)
102	guez	3	equivalence (cuscvugam, cuscvugam_2d)
103
104	guez	46	real unsapolnga1, unsapolnga2, unsapolsga1, unsapolsga2
105	guez	3
106	guez	46	real unsair_gam1_2d(iim + 1, jjm + 1), unsair_gam2_2d(iim + 1, jjm + 1)
107	guez	79	real unsair_gam1((iim + 1) * (jjm + 1)), unsair_gam2((iim + 1) * (jjm + 1))
108	guez	3	equivalence (unsair_gam1, unsair_gam1_2d), (unsair_gam2, unsair_gam2_2d)
109
110	guez	46	real unsairz_gam_2d(iim + 1, jjm)
111	guez	79	real unsairz_gam((iim + 1) * jjm)
112	guez	3	equivalence (unsairz_gam, unsairz_gam_2d)
113
114	guez	46	real xprimu(iim + 1), xprimv(iim + 1)
115	guez	3
116			save
117
118	guez	78	contains
119
120			SUBROUTINE inigeom
121
122			! Auteur : P. Le Van
123
124			! Calcul des élongations cuij1, ..., cuij4, cvij1, ..., cvij4 aux mêmes
125			! endroits que les aires aireij1_2d, ..., aireij4_2d.
126
127			! Choix entre une fonction "f(y)" à dérivée sinusoïdale ou à
128			! dérivée tangente hyperbolique. Calcul des coefficients cu_2d,
129			! cv_2d, 1. / cu_2d2, 1. / cv_2d2. Les coefficients cu_2d et cv_2d
130			! permettent de passer des vitesses naturelles aux vitesses
131			! covariantes et contravariantes, ou vice-versa.
132
133			! On a :
134			! u(covariant) = cu_2d * u(naturel), u(contravariant) = u(naturel) / cu_2d
135			! v(covariant) = cv_2d * v(naturel), v(contravariant) = v(naturel) / cv_2d
136
137			! On en tire :
138			! u(covariant) = cu_2d * cu_2d * u(contravariant)
139			! v(covariant) = cv_2d * cv_2d * v(contravariant)
140
141			! On a l'application (x(X), y(Y)) avec - im / 2 + 1 <= X <= im / 2
142			! et - jm / 2 <= Y <= jm / 2
143
144			! x est la longitude du point en radians.
145			! y est la latitude du point en radians.
146			!
147			! On a : cu_2d(i, j) = rad * cos(y) * dx / dX
148			! cv(j) = rad * dy / dY
149			! aire_2d(i, j) = cu_2d(i, j) * cv(j)
150			!
151			! y, dx / dX, dy / dY calculés aux points concernés. cv, bien que
152			! dépendant de j uniquement, sera ici indicé aussi en i pour un
153			! adressage plus facile en ij.
154
155			! xprimu et xprimv sont respectivement les valeurs de dx / dX aux
156			! points u et v. yprimu et yprimv sont respectivement les valeurs
157			! de dy / dY aux points u et v. rlatu et rlatv sont respectivement
158			! les valeurs de la latitude aux points u et v. cvu et cv_2d sont
159			! respectivement les valeurs de cv_2d aux points u et v.
160
161			! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d
162			! aux points u, v, scalaires, et z. Cf. "inigeom.txt".
163
164			USE comconst, ONLY : g, omeg, rad
165			USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
166			use conf_gcm_m, ONLY : fxyhypb, ysinus
167			use fxy_m, only: fxy
168			use fxyhyper_m, only: fxyhyper
169	guez	97	use fxysinus_m, only: fxysinus
170	guez	78	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	guez	96	! Modifiés pxo, pyo, transx, transy
176	guez	78
177	guez	96	! Local:
178	guez	78	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	guez	97	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	guez	78	IF (ysinus) THEN
228	guez	97	print *, 'inigeom: Y = sin(latitude)'
229			! Utilisation de f(x, y) avec y = sinus de la latitude
230	guez	78	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	guez	3	end module comgeom