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.

    ! cv_2d est aux points v. cu_2d est aux points
    ! u. 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 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

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

    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

    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))
          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)
       unscv2_2d(iip1, j) = unscv2_2d(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.

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

    ! 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	140	! cv_2d est aux points v. cu_2d est aux points
143			! u. Cf. "inigeom.txt".
144	guez	78
145			USE comconst, ONLY : g, omeg, rad
146			USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
147	guez	139	use dynetat0_m, only: xprimp025, xprimm025, rlatu1, rlatu2, rlatu, rlatv, &
148			yprimu1, yprimu2, rlonu, rlonv
149	guez	78	use jumble, only: new_unit
150			use nr_util, only: pi
151			USE paramet_m, ONLY : iip1, jjp1
152
153	guez	96	! Local:
154	guez	113	INTEGER i, j, unit
155	guez	78	REAL ai14, ai23, airez, un4rad2
156			REAL coslatm, coslatp, radclatm, radclatp
157			REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
158			REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4 ! in m
159			REAL gamdi_gdiv, gamdi_grot, gamdi_h
160			real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &
161			aireij4_2d ! in m2
162
163			!------------------------------------------------------------------
164
165			PRINT *, 'Call sequence information: inigeom'
166
167	guez	125	IF (nitergdiv /= 2) THEN
168			gamdi_gdiv = coefdis / (nitergdiv - 2)
169	guez	78	ELSE
170			gamdi_gdiv = 0.
171			END IF
172	guez	121
173	guez	125	IF (nitergrot /= 2) THEN
174			gamdi_grot = coefdis / (nitergrot - 2)
175	guez	78	ELSE
176			gamdi_grot = 0.
177			END IF
178	guez	121
179	guez	125	IF (niterh /= 2) THEN
180			gamdi_h = coefdis / (niterh - 2)
181	guez	78	ELSE
182			gamdi_h = 0.
183			END IF
184
185			print *, 'gamdi_gdiv = ', gamdi_gdiv
186			print *, "gamdi_grot = ", gamdi_grot
187			print *, "gamdi_h = ", gamdi_h
188
189			un4rad2 = 0.25 * rad * rad
190
191			! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
192			! aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui entourent
193			! chaque aire_2d(i, j), ainsi que les quatre élongations
194			! élémentaires cuij et les quatre élongations cvij qui sont
195			! calculées aux mêmes endroits que les aireij.
196
197			coslatm = cos(rlatu1(1))
198			radclatm = 0.5 * rad * coslatm
199
200			aireij1_2d(:iim, 1) = 0.
201			aireij2_2d(:iim, 1) = un4rad2 * coslatm * xprimp025(:iim) * yprimu1(1)
202			aireij3_2d(:iim, 1) = un4rad2 * coslatm * xprimm025(:iim) * yprimu1(1)
203			aireij4_2d(:iim, 1) = 0.
204
205			cuij1(:iim, 1) = 0.
206			cuij2(:iim, 1) = radclatm * xprimp025(:iim)
207			cuij3(:iim, 1) = radclatm * xprimm025(:iim)
208			cuij4(:iim, 1) = 0.
209
210			cvij1(:iim, 1) = 0.
211			cvij2(:iim, 1) = 0.5 * rad * yprimu1(1)
212			cvij3(:iim, 1) = cvij2(:iim, 1)
213			cvij4(:iim, 1) = 0.
214
215			do j = 2, jjm
216			coslatm = cos(rlatu1(j))
217			coslatp = cos(rlatu2(j-1))
218			radclatp = 0.5 * rad * coslatp
219			radclatm = 0.5 * rad * coslatm
220			ai14 = un4rad2 * coslatp * yprimu2(j-1)
221			ai23 = un4rad2 * coslatm * yprimu1(j)
222
223			aireij1_2d(:iim, j) = ai14 * xprimp025(:iim)
224			aireij2_2d(:iim, j) = ai23 * xprimp025(:iim)
225			aireij3_2d(:iim, j) = ai23 * xprimm025(:iim)
226			aireij4_2d(:iim, j) = ai14 * xprimm025(:iim)
227			cuij1(:iim, j) = radclatp * xprimp025(:iim)
228			cuij2(:iim, j) = radclatm * xprimp025(:iim)
229			cuij3(:iim, j) = radclatm * xprimm025(:iim)
230			cuij4(:iim, j) = radclatp * xprimm025(:iim)
231			cvij1(:iim, j) = 0.5 * rad * yprimu2(j-1)
232			cvij2(:iim, j) = 0.5 * rad * yprimu1(j)
233			cvij3(:iim, j) = cvij2(:iim, j)
234			cvij4(:iim, j) = cvij1(:iim, j)
235			end do
236
237			coslatp = cos(rlatu2(jjm))
238			radclatp = 0.5 * rad * coslatp
239
240			aireij1_2d(:iim, jjp1) = un4rad2 * coslatp * xprimp025(:iim) * yprimu2(jjm)
241			aireij2_2d(:iim, jjp1) = 0.
242			aireij3_2d(:iim, jjp1) = 0.
243			aireij4_2d(:iim, jjp1) = un4rad2 * coslatp * xprimm025(:iim) * yprimu2(jjm)
244
245			cuij1(:iim, jjp1) = radclatp * xprimp025(:iim)
246			cuij2(:iim, jjp1) = 0.
247			cuij3(:iim, jjp1) = 0.
248			cuij4(:iim, jjp1) = radclatp * xprimm025(:iim)
249
250			cvij1(:iim, jjp1) = 0.5 * rad * yprimu2(jjm)
251			cvij2(:iim, jjp1) = 0.
252			cvij3(:iim, jjp1) = 0.
253			cvij4(:iim, jjp1) = cvij1(:iim, jjp1)
254
255			! Périodicité :
256
257			cvij1(iip1, :) = cvij1(1, :)
258			cvij2(iip1, :) = cvij2(1, :)
259			cvij3(iip1, :) = cvij3(1, :)
260			cvij4(iip1, :) = cvij4(1, :)
261
262			cuij1(iip1, :) = cuij1(1, :)
263			cuij2(iip1, :) = cuij2(1, :)
264			cuij3(iip1, :) = cuij3(1, :)
265			cuij4(iip1, :) = cuij4(1, :)
266
267			aireij1_2d(iip1, :) = aireij1_2d(1, :)
268			aireij2_2d(iip1, :) = aireij2_2d(1, :)
269			aireij3_2d(iip1, :) = aireij3_2d(1, :)
270			aireij4_2d(iip1, :) = aireij4_2d(1, :)
271
272			DO j = 1, jjp1
273			DO i = 1, iim
274			aire_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) &
275			+ aireij3_2d(i, j) + aireij4_2d(i, j)
276			alpha1_2d(i, j) = aireij1_2d(i, j) / aire_2d(i, j)
277			alpha2_2d(i, j) = aireij2_2d(i, j) / aire_2d(i, j)
278			alpha3_2d(i, j) = aireij3_2d(i, j) / aire_2d(i, j)
279			alpha4_2d(i, j) = aireij4_2d(i, j) / aire_2d(i, j)
280			alpha1p2_2d(i, j) = alpha1_2d(i, j) + alpha2_2d(i, j)
281			alpha1p4_2d(i, j) = alpha1_2d(i, j) + alpha4_2d(i, j)
282			alpha2p3_2d(i, j) = alpha2_2d(i, j) + alpha3_2d(i, j)
283			alpha3p4_2d(i, j) = alpha3_2d(i, j) + alpha4_2d(i, j)
284			END DO
285
286			aire_2d(iip1, j) = aire_2d(1, j)
287			alpha1_2d(iip1, j) = alpha1_2d(1, j)
288			alpha2_2d(iip1, j) = alpha2_2d(1, j)
289			alpha3_2d(iip1, j) = alpha3_2d(1, j)
290			alpha4_2d(iip1, j) = alpha4_2d(1, j)
291			alpha1p2_2d(iip1, j) = alpha1p2_2d(1, j)
292			alpha1p4_2d(iip1, j) = alpha1p4_2d(1, j)
293			alpha2p3_2d(iip1, j) = alpha2p3_2d(1, j)
294			alpha3p4_2d(iip1, j) = alpha3p4_2d(1, j)
295			END DO
296
297			DO j = 1, jjp1
298			DO i = 1, iim
299			aireu_2d(i, j) = aireij1_2d(i, j) + aireij2_2d(i, j) + &
300			aireij4_2d(i + 1, j) + aireij3_2d(i + 1, j)
301			unsaire_2d(i, j) = 1. / aire_2d(i, j)
302			unsair_gam1_2d(i, j) = unsaire_2d(i, j)**(-gamdi_gdiv)
303			unsair_gam2_2d(i, j) = unsaire_2d(i, j)**(-gamdi_h)
304			airesurg_2d(i, j) = aire_2d(i, j) / g
305			END DO
306			aireu_2d(iip1, j) = aireu_2d(1, j)
307			unsaire_2d(iip1, j) = unsaire_2d(1, j)
308			unsair_gam1_2d(iip1, j) = unsair_gam1_2d(1, j)
309			unsair_gam2_2d(iip1, j) = unsair_gam2_2d(1, j)
310			airesurg_2d(iip1, j) = airesurg_2d(1, j)
311			END DO
312
313			DO j = 1, jjm
314			DO i = 1, iim
315			airev_2d(i, j) = aireij2_2d(i, j) + aireij3_2d(i, j) + &
316			aireij1_2d(i, j + 1) + aireij4_2d(i, j + 1)
317			END DO
318			DO i = 1, iim
319			airez = aireij2_2d(i, j) + aireij1_2d(i, j + 1) &
320			+ aireij3_2d(i + 1, j) + aireij4_2d(i + 1, j + 1)
321			unsairez_2d(i, j) = 1. / airez
322			unsairz_gam_2d(i, j) = unsairez_2d(i, j)**(-gamdi_grot)
323			fext_2d(i, j) = airez * sin(rlatv(j)) * 2. * omeg
324			END DO
325			airev_2d(iip1, j) = airev_2d(1, j)
326			unsairez_2d(iip1, j) = unsairez_2d(1, j)
327			fext_2d(iip1, j) = fext_2d(1, j)
328			unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
329			END DO
330
331	guez	140	! Calcul des élongations cu_2d, cv_2d
332	guez	78
333			DO j = 1, jjm
334			DO i = 1, iim
335			cv_2d(i, j) = 0.5 * &
336			(cvij2(i, j) + cvij3(i, j) + cvij1(i, j + 1) + cvij4(i, j + 1))
337			unscv2_2d(i, j) = 1. / cv_2d(i, j)**2
338			END DO
339			DO i = 1, iim
340			cuvsurcv_2d(i, j) = airev_2d(i, j) * unscv2_2d(i, j)
341			cvsurcuv_2d(i, j) = 1. / cuvsurcv_2d(i, j)
342			cuvscvgam1_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_gdiv)
343			cuvscvgam2_2d(i, j) = cuvsurcv_2d(i, j)**(-gamdi_h)
344			cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
345			END DO
346			cv_2d(iip1, j) = cv_2d(1, j)
347			unscv2_2d(iip1, j) = unscv2_2d(1, j)
348			cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
349			cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
350			cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
351			cuvscvgam2_2d(iip1, j) = cuvscvgam2_2d(1, j)
352			cvscuvgam_2d(iip1, j) = cvscuvgam_2d(1, j)
353			END DO
354
355			DO j = 2, jjm
356			DO i = 1, iim
357			cu_2d(i, j) = 0.5 * (cuij1(i, j) + cuij4(i + 1, j) + cuij2(i, j) &
358			+ cuij3(i + 1, j))
359			unscu2_2d(i, j) = 1. / cu_2d(i, j)**2
360			cvusurcu_2d(i, j) = aireu_2d(i, j) * unscu2_2d(i, j)
361			cusurcvu_2d(i, j) = 1. / cvusurcu_2d(i, j)
362			cvuscugam1_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_gdiv)
363			cvuscugam2_2d(i, j) = cvusurcu_2d(i, j)**(-gamdi_h)
364			cuscvugam_2d(i, j) = cusurcvu_2d(i, j)**(-gamdi_grot)
365			END DO
366			cu_2d(iip1, j) = cu_2d(1, j)
367			unscu2_2d(iip1, j) = unscu2_2d(1, j)
368			cvusurcu_2d(iip1, j) = cvusurcu_2d(1, j)
369			cusurcvu_2d(iip1, j) = cusurcvu_2d(1, j)
370			cvuscugam1_2d(iip1, j) = cvuscugam1_2d(1, j)
371			cvuscugam2_2d(iip1, j) = cvuscugam2_2d(1, j)
372			cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j)
373			END DO
374
375			! Calcul aux pôles
376
377			cu_2d(:, 1) = 0.
378			unscu2_2d(:, 1) = 0.
379
380			cu_2d(:, jjp1) = 0.
381			unscu2_2d(:, jjp1) = 0.
382
383			! Calcul des aires aux pôles :
384
385			apoln = sum(aire_2d(:iim, 1))
386			apols = sum(aire_2d(:iim, jjp1))
387			unsapolnga1 = 1. / (apoln**(-gamdi_gdiv))
388			unsapolsga1 = 1. / (apols**(-gamdi_gdiv))
389			unsapolnga2 = 1. / (apoln**(-gamdi_h))
390			unsapolsga2 = 1. / (apols**(-gamdi_h))
391
392			! Changement F. Hourdin calcul conservatif pour fext_2d
393			! constang_2d contient le produit a * cos (latitude) * omega
394
395			DO i = 1, iim
396			constang_2d(i, 1) = 0.
397			END DO
398			DO j = 1, jjm - 1
399			DO i = 1, iim
400			constang_2d(i, j + 1) = rad * omeg * cu_2d(i, j + 1) &
401			* cos(rlatu(j + 1))
402			END DO
403			END DO
404			DO i = 1, iim
405			constang_2d(i, jjp1) = 0.
406			END DO
407
408			! Périodicité en longitude
409			DO j = 1, jjp1
410			constang_2d(iip1, j) = constang_2d(1, j)
411			END DO
412
413			call new_unit(unit)
414			open(unit, file="longitude_latitude.txt", status="replace", action="write")
415			write(unit, fmt=) '"longitudes at V points (degrees)"', rlonv 180. / pi
416			write(unit, fmt=) '"latitudes at V points (degrees)"', rlatv 180. / pi
417			write(unit, fmt=) '"longitudes at U points (degrees)"', rlonu 180. / pi
418			write(unit, fmt=) '"latitudes at U points (degrees)"', rlatu 180. / pi
419			close(unit)
420
421			END SUBROUTINE inigeom
422
423	guez	3	end module comgeom