trunk/dyn3d/comgeom.f

module comgeom

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

  implicit none

  private iim, jjm, ip1jmp1, ip1jm

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

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

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

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

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

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

  real apoln, apols ! in m2

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  real unsapolnga1, unsapolnga2, unsapolsga1, unsapolsga2                

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

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

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

  save

contains

  SUBROUTINE inigeom

    ! Auteur : P. Le Van

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

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

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

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

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

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

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

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

    USE comconst, ONLY : g, omeg, rad
    USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
    use conf_gcm_m, ONLY : fxyhypb, ysinus
    USE dimens_m, ONLY : iim, jjm
    use fxy_m, only: fxy
    use fxyhyper_m, only: fxyhyper
    use jumble, only: new_unit
    use nr_util, only: pi
    USE paramet_m, ONLY : iip1, jjp1
    USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &
         grossismy, pxo, pyo, taux, tauy, transx, transy
    ! Modifies pxo, pyo, transx, transy

    ! Variables locales

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

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

    PRINT *, 'Call sequence information: inigeom'

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

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

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

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

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

          itmax = 10
          eps = .1E-7

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

          transx = xo1

          itmay = 10
          eps = .1E-7

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

          transy = yo1

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

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

    ! Calcul aux pôles 

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

    un4rad2 = 0.25 * rad * rad

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

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

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

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

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

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

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

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

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

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

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

    ! Périodicité :

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

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

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

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

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

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

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

    ! Calcul des élongations cu_2d, cv_2d, cvu 

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

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

    ! Calcul aux pôles 

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

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

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

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

    ! Calcul des aires aux pôles :

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

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

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

    ! Périodicité en longitude

    DO j = 1, jjm
       fext_2d(iip1, j) = fext_2d(1, j)
    END DO
    DO j = 1, jjp1
       constang_2d(iip1, j) = constang_2d(1, j)
    END DO

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

  END SUBROUTINE inigeom

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