! This module replaces grid/fxy*.h

MODULE lmdz_fxy
  IMPLICIT NONE; PRIVATE
  PUBLIC fxy_new, fxy_reg, fxy_sin, fxy_prim
CONTAINS
  SUBROUTINE fxy_new(ripx, fx, fxprim, fy, fyprim, ri, rj, bigy)
    INCLUDE "dimensions.h"
    REAL :: ripx, fx, fxprim, fy, fyprim, ri, rj, bigy

    !....stretching in x...
    ripx(ri) = (ri - 1.0) * 2. * pi / REAL(iim)
    fx  (ri) = ripx(ri) + transx + &
            alphax * SIN(ripx(ri) + transx - pxo) - pi
    fxprim(ri) = 2. * pi / REAL(iim) * &
            (1. + alphax * COS(ripx(ri) + transx - pxo))

    !....stretching in y...
    bigy(rj) = 2. * (REAL(jjp1) - rj) * pi / jjm
    fy(rj) = (bigy(rj) + transy + &
            alphay * SIN(bigy(rj) + transy - pyo)) / 2. - pi / 2.
    fyprim(rj) = (pi / jjm) * (1. + &
            alphay * COS(bigy(rj) + transy - pyo))

    ! fy(rj)= pyo-pisjjm*(rj-transy)+coefalpha*SIN(depisjm*(rj-
    ! *  transy ))
    !   fyprim(rj)= pisjjm-pisjjm*coefy2* COS(depisjm*(rj-transy))
    !--------------------------------------------------------------
  END SUBROUTINE fxy_new

  SUBROUTINE fxy_reg(fy, fx, fxprim, fyprim, ri, rj)
    INCLUDE "dimensions.h"
    REAL :: fy, fx, fxprim, fyprim, ri, rj

    fy    (rj) = pi / REAL(jjm) * (0.5 * REAL(jjm) + 1. - rj)
    fyprim(rj) = pi / REAL(jjm)

    ! fy(rj)=ASIN(1.+2.*((1.-rj)/REAL(jjm)))
    ! fyprim(rj)=1./SQRT((rj-1.)*(jjm+1.-rj))

    fx    (ri) = 2. * pi / REAL(iim) * (ri - 0.5 * REAL(iim) - 1.)
    ! fx    ( ri ) = 2.*pi/REAL(iim) * ( ri - 0.5* ( REAL(iim) + 1.) )
    fxprim(ri) = 2. * pi / REAL(iim)
    !
    !    La valeur de pi est passee par le common/const/ou /const2/ .
    !    Sinon, il faut la calculer avant d'appeler ces fonctions .
    !
    !   ----------------------------------------------------------------
    ! Fonctions a changer eventuellement, selon x(x) et y(y) choisis .
    !   -----------------------------------------------------------------
    !
    !    .....  ici, on a l'application particuliere suivante   ........
    !
    !            **************************************
    !            **     x = 2. * pi/iim *  X         **
    !            **     y =      pi/jjm *  Y         **
    !            **************************************
    !
    !   ..................................................................
    !   ..................................................................
    !
    !
    !
    !-----------------------------------------------------------------------

  END SUBROUTINE fxy_reg

  SUBROUTINE fxy_sin(fy, fx, fxprim, fyprim, ri, rj)
    REAL :: fy, fx, fxprim, fyprim, ri, rj

    fy(rj) = ASIN(1. + 2. * ((1. - rj) / REAL(jjm)))
    fyprim(rj) = 1. / SQRT((rj - 1.) * (jjm + 1. - rj))

    fx    (ri) = 2. * pi / REAL(iim) * (ri - 0.5 * REAL(iim) - 1.)
    ! fx    ( ri ) = 2.*pi/REAL(iim) * ( ri - 0.5* ( REAL(iim) + 1.) )
    fxprim(ri) = 2. * pi / REAL(iim)
    !
    !
    !    La valeur de pi est passee par le common/const/ou /const2/ .
    !    Sinon, il faut la calculer avant d'appeler ces fonctions .
    !
    !   ----------------------------------------------------------------
    ! Fonctions a changer eventuellement, selon x(x) et y(y) choisis .
    !   -----------------------------------------------------------------
    !
    !    .....  ici, on a l'application particuliere suivante   ........
    !
    !            **************************************
    !            **     x = 2. * pi/iim *  X         **
    !            **     y =      pi/jjm *  Y         **
    !            **************************************
    !
    !   ..................................................................
    !   ..................................................................
    !
    !
    !
    !-----------------------------------------------------------------------
  END SUBROUTINE fxy_sin

  SUBROUTINE fxy_prim(fy, fx, fxprim, fyprim, ri, rj)
    REAL :: fy, fx, fxprim, fyprim, ri, rj

    fy    (rj) = pi / REAL(jjm) * (0.5 * REAL(jjm) + 1. - rj)
    fyprim(rj) = pi / REAL(jjm)

    ! fy(rj)=ASIN(1.+2.*((1.-rj)/REAL(jjm)))
    ! fyprim(rj)=1./SQRT((rj-1.)*(jjm+1.-rj))

    fx    (ri) = 2. * pi / REAL(iim) * (ri - 0.5 * REAL(iim) - 1.)
    ! fx    ( ri ) = 2.*pi/REAL(iim) * ( ri - 0.5* ( REAL(iim) + 1.) )
    fxprim(ri) = 2. * pi / REAL(iim)
    !
    !
    !    La valeur de pi est passee par le common/const/ou /const2/ .
    !    Sinon, il faut la calculer avant d'appeler ces fonctions .
    !
    !   ----------------------------------------------------------------
    ! Fonctions a changer eventuellement, selon x(x) et y(y) choisis .
    !   -----------------------------------------------------------------
    !
    !    .....  ici, on a l'application particuliere suivante   ........
    !
    !            **************************************
    !            **     x = 2. * pi/iim *  X         **
    !            **     y =      pi/jjm *  Y         **
    !            **************************************
    !
    !   ..................................................................
    !   ..................................................................
    !
    !
    !
    !-----------------------------------------------------------------------
  END SUBROUTINE fxy_prim
END MODULE lmdz_fxy