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 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)

  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.

    ! Fonction "f(y)" à 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.

    ! 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 dynetat0_m, only: xprimp025, xprimm025, rlatu1, rlatu2, rlatu, rlatv, &
         yprimu1, yprimu2, rlonu, rlonv
    use jumble, only: new_unit
    use nr_util, only: pi
    USE paramet_m, ONLY : iip1, jjp1

    ! Local:
    INTEGER i, j, unit
    REAL cvu(iip1, jjp1), cuv(iip1, jjm)
    REAL ai14, ai23, airez, un4rad2
    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 gamdi_gdiv, gamdi_grot, gamdi_h
    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 / (nitergdiv - 2)
    ELSE
       gamdi_gdiv = 0.
    END IF

    IF (nitergrot /= 2) THEN
       gamdi_grot = coefdis / (nitergrot - 2)
    ELSE
       gamdi_grot = 0.
    END IF

    IF (niterh /= 2) THEN
       gamdi_h = coefdis / (niterh - 2)
    ELSE
       gamdi_h = 0.
    END IF

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

    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	guez	60	real cuvsurcv_2d(iim + 1, jjm), cvsurcuv_2d(iim + 1, jjm) ! no dimension
68	guez	79	real cuvsurcv((iim + 1) * jjm), cvsurcuv((iim + 1) * jjm) ! no dimension
69	guez	3	equivalence (cuvsurcv, cuvsurcv_2d), (cvsurcuv, cvsurcuv_2d)
70
71	guez	46	real cvusurcu_2d(iim + 1, jjm + 1), cusurcvu_2d(iim + 1, jjm + 1)
72	guez	60	! no dimension
73	guez	79	real cvusurcu((iim + 1) * (jjm + 1)), cusurcvu((iim + 1) * (jjm + 1))
74			! no dimension
75	guez	3	equivalence (cvusurcu, cvusurcu_2d), (cusurcvu, cusurcvu_2d)
76
77	guez	46	real cuvscvgam1_2d(iim + 1, jjm)
78	guez	79	real cuvscvgam1((iim + 1) * jjm)
79	guez	3	equivalence (cuvscvgam1, cuvscvgam1_2d)
80
81	guez	46	real cuvscvgam2_2d(iim + 1, jjm), cvuscugam1_2d(iim + 1, jjm + 1)
82	guez	79	real cuvscvgam2((iim + 1) * jjm), cvuscugam1((iim + 1) * (jjm + 1))
83	guez	3	equivalence (cuvscvgam2, cuvscvgam2_2d), (cvuscugam1, cvuscugam1_2d)
84
85	guez	46	real cvuscugam2_2d(iim + 1, jjm + 1), cvscuvgam_2d(iim + 1, jjm)
86	guez	79	real cvuscugam2((iim + 1) * (jjm + 1)), cvscuvgam((iim + 1) * jjm)
87	guez	3	equivalence (cvuscugam2, cvuscugam2_2d), (cvscuvgam, cvscuvgam_2d)
88
89	guez	79	real cuscvugam((iim + 1) * (jjm + 1))
90	guez	46	real cuscvugam_2d(iim + 1, jjm + 1)
91	guez	3	equivalence (cuscvugam, cuscvugam_2d)
92
93	guez	46	real unsapolnga1, unsapolnga2, unsapolsga1, unsapolsga2
94	guez	3
95	guez	46	real unsair_gam1_2d(iim + 1, jjm + 1), unsair_gam2_2d(iim + 1, jjm + 1)
96	guez	79	real unsair_gam1((iim + 1) * (jjm + 1)), unsair_gam2((iim + 1) * (jjm + 1))
97	guez	3	equivalence (unsair_gam1, unsair_gam1_2d), (unsair_gam2, unsair_gam2_2d)
98
99	guez	46	real unsairz_gam_2d(iim + 1, jjm)
100	guez	79	real unsairz_gam((iim + 1) * jjm)
101	guez	3	equivalence (unsairz_gam, unsairz_gam_2d)
102
103			save
104
105	guez	78	contains
106
107			SUBROUTINE inigeom
108
109			! Auteur : P. Le Van
110
111			! Calcul des élongations cuij1, ..., cuij4, cvij1, ..., cvij4 aux mêmes
112			! endroits que les aires aireij1_2d, ..., aireij4_2d.
113
114	guez	121	! Fonction "f(y)" à dérivée tangente hyperbolique. Calcul des
115			! coefficients cu_2d, cv_2d, 1. / cu_2d2, 1. / cv_2d2. Les
116			! coefficients cu_2d et cv_2d permettent de passer des vitesses
117			! naturelles aux vitesses covariantes et contravariantes, ou
118			! vice-versa.
119	guez	78
120			! On a :
121			! u(covariant) = cu_2d * u(naturel), u(contravariant) = u(naturel) / cu_2d
122			! v(covariant) = cv_2d * v(naturel), v(contravariant) = v(naturel) / cv_2d
123
124			! On en tire :
125			! u(covariant) = cu_2d * cu_2d * u(contravariant)
126			! v(covariant) = cv_2d * cv_2d * v(contravariant)
127
128			! On a l'application (x(X), y(Y)) avec - im / 2 + 1 <= X <= im / 2
129			! et - jm / 2 <= Y <= jm / 2
130
131			! x est la longitude du point en radians.
132			! y est la latitude du point en radians.
133			!
134			! On a : cu_2d(i, j) = rad * cos(y) * dx / dX
135			! cv(j) = rad * dy / dY
136			! aire_2d(i, j) = cu_2d(i, j) * cv(j)
137			!
138			! y, dx / dX, dy / dY calculés aux points concernés. cv, bien que
139			! dépendant de j uniquement, sera ici indicé aussi en i pour un
140			! adressage plus facile en ij.
141
142	guez	139	! cvu et cv_2d sont respectivement les valeurs de
143			! cv_2d aux points u et v.
144	guez	78
145			! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d
146			! aux points u, v, scalaires, et z. Cf. "inigeom.txt".
147
148			USE comconst, ONLY : g, omeg, rad
149			USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
150	guez	139	use dynetat0_m, only: xprimp025, xprimm025, rlatu1, rlatu2, rlatu, rlatv, &
151			yprimu1, yprimu2, rlonu, rlonv
152	guez	78	use jumble, only: new_unit
153			use nr_util, only: pi
154			USE paramet_m, ONLY : iip1, jjp1
155
156	guez	96	! Local:
157	guez	113	INTEGER i, j, unit
158	guez	78	REAL cvu(iip1, jjp1), cuv(iip1, jjm)
159			REAL ai14, ai23, airez, un4rad2
160			REAL coslatm, coslatp, radclatm, radclatp
161			REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
162			REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4 ! in m
163			REAL gamdi_gdiv, gamdi_grot, gamdi_h
164			real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &
165			aireij4_2d ! in m2
166			real airuscv2_2d(iim + 1, jjm)
167			real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)
168			real aivscu2gam_2d(iim + 1, jjm)
169
170			!------------------------------------------------------------------
171
172			PRINT *, 'Call sequence information: inigeom'
173
174	guez	125	IF (nitergdiv /= 2) THEN
175			gamdi_gdiv = coefdis / (nitergdiv - 2)
176	guez	78	ELSE
177			gamdi_gdiv = 0.
178			END IF
179	guez	121
180	guez	125	IF (nitergrot /= 2) THEN
181			gamdi_grot = coefdis / (nitergrot - 2)
182	guez	78	ELSE
183			gamdi_grot = 0.
184			END IF
185	guez	121
186	guez	125	IF (niterh /= 2) THEN
187			gamdi_h = coefdis / (niterh - 2)
188	guez	78	ELSE
189			gamdi_h = 0.
190			END IF
191
192			print *, 'gamdi_gdiv = ', gamdi_gdiv
193			print *, "gamdi_grot = ", gamdi_grot
194			print *, "gamdi_h = ", gamdi_h
195
196			un4rad2 = 0.25 * rad * rad
197
198			! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
199			! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
200			! chaque aire_2d(i, j), ainsi que les quatre élongations
201			! élémentaires cuij et les quatre élongations cvij qui sont
202			! calculées aux mêmes endroits que les aireij.
203
204			coslatm = cos(rlatu1(1))
205			radclatm = 0.5 * rad * coslatm
206
207			aireij1_2d(:iim, 1) = 0.
208			aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
209			aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
210			aireij4_2d(:iim, 1) = 0.
211
212			cuij1(:iim, 1) = 0.
213			cuij2(:iim, 1) = radclatm * xprimp025(:iim)
214			cuij3(:iim, 1) = radclatm * xprimm025(:iim)
215			cuij4(:iim, 1) = 0.
216
217			cvij1(:iim, 1) = 0.
218			cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
219			cvij3(:iim, 1) = cvij2(:iim, 1)
220			cvij4(:iim, 1) = 0.
221
222			do j = 2, jjm
223			coslatm = cos(rlatu1(j))
224			coslatp = cos(rlatu2(j-1))
225			radclatp = 0.5 * rad * coslatp
226			radclatm = 0.5 * rad * coslatm
227			ai14 = un4rad2 * coslatp * yprimu2(j-1)
228			ai23 = un4rad2 * coslatm * yprimu1(j)
229
230			aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
231			aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
232			aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
233			aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
234			cuij1(:iim, j) = radclatp * xprimp025(:iim)
235			cuij2(:iim, j) = radclatm * xprimp025(:iim)
236			cuij3(:iim, j) = radclatm * xprimm025(:iim)
237			cuij4(:iim, j) = radclatp * xprimm025(:iim)
238			cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
239			cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
240			cvij3(:iim, j) = cvij2(:iim, j)
241			cvij4(:iim, j) = cvij1(:iim, j)
242			end do
243
244			coslatp = cos(rlatu2(jjm))
245			radclatp = 0.5 * rad * coslatp
246
247			aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
248			aireij2_2d(:iim, jjp1) = 0.
249			aireij3_2d(:iim, jjp1) = 0.
250			aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)
251
252			cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
253			cuij2(:iim, jjp1) = 0.
254			cuij3(:iim, jjp1) = 0.
255			cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)
256
257			cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
258			cvij2(:iim, jjp1) = 0.
259			cvij3(:iim, jjp1) = 0.
260			cvij4(:iim, jjp1) = cvij1(:iim, jjp1)
261
262			! Périodicité :
263
264			cvij1(iip1, :) = cvij1(1, :)
265			cvij2(iip1, :) = cvij2(1, :)
266			cvij3(iip1, :) = cvij3(1, :)
267			cvij4(iip1, :) = cvij4(1, :)
268
269			cuij1(iip1, :) = cuij1(1, :)
270			cuij2(iip1, :) = cuij2(1, :)
271			cuij3(iip1, :) = cuij3(1, :)
272			cuij4(iip1, :) = cuij4(1, :)
273
274			aireij1_2d(iip1, :) = aireij1_2d(1, :)
275			aireij2_2d(iip1, :) = aireij2_2d(1, :)
276			aireij3_2d(iip1, :) = aireij3_2d(1, :)
277			aireij4_2d(iip1, :) = aireij4_2d(1, :)
278
279			DO j = 1, jjp1
280			DO i = 1, iim
281			aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
282			+ aireij3_2d(i, j) + aireij4_2d(i, j)
283			alpha1_2d(i, j) = aireij1_2d(i, j) / aire_2d(i, j)
284			alpha2_2d(i, j) = aireij2_2d(i, j) / aire_2d(i, j)
285			alpha3_2d(i, j) = aireij3_2d(i, j) / aire_2d(i, j)
286			alpha4_2d(i, j) = aireij4_2d(i, j) / aire_2d(i, j)
287			alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
288			alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
289			alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)
290			alpha3p4_2d(i, j) = alpha3_2d(i, j) + alpha4_2d(i, j)
291			END DO
292
293			aire_2d(iip1, j) = aire_2d(1, j)
294			alpha1_2d(iip1, j) = alpha1_2d(1, j)
295			alpha2_2d(iip1, j) = alpha2_2d(1, j)
296			alpha3_2d(iip1, j) = alpha3_2d(1, j)
297			alpha4_2d(iip1, j) = alpha4_2d(1, j)
298			alpha1p2_2d(iip1, j) = alpha1p2_2d(1, j)
299			alpha1p4_2d(iip1, j) = alpha1p4_2d(1, j)
300			alpha2p3_2d(iip1, j) = alpha2p3_2d(1, j)
301			alpha3p4_2d(iip1, j) = alpha3p4_2d(1, j)
302			END DO
303
304			DO j = 1, jjp1
305			DO i = 1, iim
306			aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
307			aireij4_2d(i + 1, j) + aireij3_2d(i + 1, j)
308			unsaire_2d(i, j) = 1. / aire_2d(i, j)
309			unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
310			unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
311			airesurg_2d(i, j) = aire_2d(i, j) / g
312			END DO
313			aireu_2d(iip1, j) = aireu_2d(1, j)
314			unsaire_2d(iip1, j) = unsaire_2d(1, j)
315			unsair_gam1_2d(iip1, j) = unsair_gam1_2d(1, j)
316			unsair_gam2_2d(iip1, j) = unsair_gam2_2d(1, j)
317			airesurg_2d(iip1, j) = airesurg_2d(1, j)
318			END DO
319
320			DO j = 1, jjm
321			DO i = 1, iim
322			airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
323			aireij1_2d(i, j + 1) + aireij4_2d(i, j + 1)
324			END DO
325			DO i = 1, iim
326			airez = aireij2_2d(i, j) + aireij1_2d(i, j + 1) &
327			+ aireij3_2d(i + 1, j) + aireij4_2d(i + 1, j + 1)
328			unsairez_2d(i, j) = 1. / airez
329			unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
330			fext_2d(i, j) = airez * sin(rlatv(j)) * 2. * omeg
331			END DO
332			airev_2d(iip1, j) = airev_2d(1, j)
333			unsairez_2d(iip1, j) = unsairez_2d(1, j)
334			fext_2d(iip1, j) = fext_2d(1, j)
335			unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
336			END DO
337
338			! Calcul des élongations cu_2d, cv_2d, cvu
339
340			DO j = 1, jjm
341			DO i = 1, iim
342			cv_2d(i, j) = 0.5 * &
343			(cvij2(i, j) + cvij3(i, j) + cvij1(i, j + 1) + cvij4(i, j + 1))
344			cvu(i, j) = 0.5 * (cvij1(i, j) + cvij4(i, j) + cvij2(i, j) &
345			+ cvij3(i, j))
346			cuv(i, j) = 0.5 * (cuij2(i, j) + cuij3(i, j) + cuij1(i, j + 1) &
347			+ cuij4(i, j + 1))
348			unscv2_2d(i, j) = 1. / cv_2d(i, j)**2
349			END DO
350			DO i = 1, iim
351			cuvsurcv_2d(i, j) = airev_2d(i, j) * unscv2_2d(i, j)
352			cvsurcuv_2d(i, j) = 1. / cuvsurcv_2d(i, j)
353			cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
354			cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
355			cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
356			END DO
357			cv_2d(iip1, j) = cv_2d(1, j)
358			cvu(iip1, j) = cvu(1, j)
359			unscv2_2d(iip1, j) = unscv2_2d(1, j)
360			cuv(iip1, j) = cuv(1, j)
361			cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
362			cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
363			cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
364			cuvscvgam2_2d(iip1, j) = cuvscvgam2_2d(1, j)
365			cvscuvgam_2d(iip1, j) = cvscuvgam_2d(1, j)
366			END DO
367
368			DO j = 2, jjm
369			DO i = 1, iim
370			cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i + 1, j) + cuij2(i, j) &
371			+ cuij3(i + 1, j))
372			unscu2_2d(i, j) = 1. / cu_2d(i, j)**2
373			cvusurcu_2d(i, j) = aireu_2d(i, j) * unscu2_2d(i, j)
374			cusurcvu_2d(i, j) = 1. / cvusurcu_2d(i, j)
375			cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
376			cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
377			cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
378			END DO
379			cu_2d(iip1, j) = cu_2d(1, j)
380			unscu2_2d(iip1, j) = unscu2_2d(1, j)
381			cvusurcu_2d(iip1, j) = cvusurcu_2d(1, j)
382			cusurcvu_2d(iip1, j) = cusurcvu_2d(1, j)
383			cvuscugam1_2d(iip1, j) = cvuscugam1_2d(1, j)
384			cvuscugam2_2d(iip1, j) = cvuscugam2_2d(1, j)
385			cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
386			END DO
387
388			! Calcul aux pôles
389
390			cu_2d(:, 1) = 0.
391			unscu2_2d(:, 1) = 0.
392			cvu(:, 1) = 0.
393
394			cu_2d(:, jjp1) = 0.
395			unscu2_2d(:, jjp1) = 0.
396			cvu(:, jjp1) = 0.
397
398			DO j = 1, jjm
399			DO i = 1, iim
400			airvscu2_2d(i, j) = airev_2d(i, j) / (cuv(i, j) * cuv(i, j))
401			aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)
402			END DO
403			airvscu2_2d(iip1, j) = airvscu2_2d(1, j)
404			aivscu2gam_2d(iip1, j) = aivscu2gam_2d(1, j)
405			END DO
406
407			DO j = 2, jjm
408			DO i = 1, iim
409			airuscv2_2d(i, j) = aireu_2d(i, j) / (cvu(i, j) * cvu(i, j))
410			aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)
411			END DO
412			airuscv2_2d(iip1, j) = airuscv2_2d(1, j)
413			aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)
414			END DO
415
416			! Calcul des aires aux pôles :
417
418			apoln = sum(aire_2d(:iim, 1))
419			apols = sum(aire_2d(:iim, jjp1))
420			unsapolnga1 = 1. / (apoln**(-gamdi_gdiv))
421			unsapolsga1 = 1. / (apols**(-gamdi_gdiv))
422			unsapolnga2 = 1. / (apoln**(-gamdi_h))
423			unsapolsga2 = 1. / (apols**(-gamdi_h))
424
425			! Changement F. Hourdin calcul conservatif pour fext_2d
426			! constang_2d contient le produit a * cos (latitude) * omega
427
428			DO i = 1, iim
429			constang_2d(i, 1) = 0.
430			END DO
431			DO j = 1, jjm - 1
432			DO i = 1, iim
433			constang_2d(i, j + 1) = rad * omeg * cu_2d(i, j + 1) &
434			* cos(rlatu(j + 1))
435			END DO
436			END DO
437			DO i = 1, iim
438			constang_2d(i, jjp1) = 0.
439			END DO
440
441			! Périodicité en longitude
442			DO j = 1, jjp1
443			constang_2d(iip1, j) = constang_2d(1, j)
444			END DO
445
446			call new_unit(unit)
447			open(unit, file="longitude_latitude.txt", status="replace", action="write")
448			write(unit, fmt=) '"longitudes at V points (degrees)"', rlonv 180. / pi
449			write(unit, fmt=) '"latitudes at V points (degrees)"', rlatv 180. / pi
450			write(unit, fmt=) '"longitudes at U points (degrees)"', rlonu 180. / pi
451			write(unit, fmt=) '"latitudes at U points (degrees)"', rlatu 180. / pi
452			close(unit)
453
454			END SUBROUTINE inigeom
455
456	guez	3	end module comgeom