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

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

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

          itmax = 10
          eps = .1E-7

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

          transx = xo1

          itmay = 10
          eps = .1E-7

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

          transy = yo1

          CALL fxy(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
               yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
               rlonp025, xprimp025)
       END IF
    ELSE
       ! Utilisation de fxyhyper, f(x, y) à dérivée tangente hyperbolique
       print *, 'Inigeom, Y = Latitude, dérivée tangente hyperbolique'
       CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &
            taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
            yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &
            rlonp025, xprimp025)
    END IF

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

    ! Calcul aux pôles 

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

    un4rad2 = 0.25 * rad * rad

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

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

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

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

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

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

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

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

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

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

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

    ! Périodicité :

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

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

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

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

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

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

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

    ! Calcul des élongations cu_2d, cv_2d, cvu 

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

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

    ! Calcul aux pôles 

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

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

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

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

    ! Calcul des aires aux pôles :

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

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

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

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