Sources/dyn3d/comgeom.f

module comgeom

  use dimens_m, only: iim, jjm

  implicit none

  private iim, jjm

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

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

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

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

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

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

  real apoln, apols ! in m2

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  real unsapolnga1, unsapolnga2, unsapolsga1, unsapolsga2                

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

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

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

  save

contains

  SUBROUTINE inigeom

    ! Auteur : P. Le Van

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

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

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

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

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

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

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

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

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

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

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

    PRINT *, 'Call sequence information: inigeom'

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

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

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

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

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

          itmax = 10
          eps = .1E-7

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

          transx = xo1

          itmay = 10
          eps = .1E-7

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

          transy = yo1

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

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

    ! Calcul aux pôles 

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

    un4rad2 = 0.25 * rad * rad

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

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

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

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

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

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

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

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

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

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

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

    ! Périodicité :

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

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

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

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

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

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

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

    ! Calcul des élongations cu_2d, cv_2d, cvu 

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

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

    ! Calcul aux pôles 

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

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

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

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

    ! Calcul des aires aux pôles :

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

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

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

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

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

  END SUBROUTINE inigeom

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