--- trunk/dyn3d/fxysinus.f	2014/03/05 14:57:53	82
+++ trunk/dyn3d/fxysinus.f	2014/04/25 14:58:31	97
@@ -1,111 +1,99 @@
+module fxysinus_m
 
-! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/fxysinus.F,v 1.1.1.1 2004/05/19
-! 12:53:06 lmdzadmin Exp $
-
-SUBROUTINE fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
-    yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, rlonp025, &
-    xprimp025)
-
-
-  USE dimens_m
-  USE paramet_m
-  USE comconst
-  USE nr_util, ONLY: pi
   IMPLICIT NONE
 
-  ! Calcul  des longitudes et des latitudes  pour une fonction f(x,y)
-  ! avec y = Asin( j )  .
-
-  ! Auteur  :  P. Le Van
-
-
-
-  INTEGER i, j
-
-  REAL rlatu(jjp1), yprimu(jjp1), rlatv(jjm), yprimv(jjm), rlatu1(jjm), &
-    yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
-  REAL rlonu(iip1), xprimu(iip1), rlonv(iip1), xprimv(iip1), rlonm025(iip1), &
-    xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)
-
-
-  ! $Header: /home/cvsroot/LMDZ4/libf/grid/fxy_sin.h,v 1.1.1.1 2004/05/19
-  ! 12:53:05 lmdzadmin Exp $
-
-  ! -----------------------------------------------------------------------
-
-  ! ................................................................
-  ! ................  Fonctions in line  ...........................
-  ! ................................................................
-
-  REAL fy, fx, fxprim, fyprim
-  REAL ri, rj
-
-
-  fy(rj) = asin(1.+2.*((1.-rj)/float(jjm)))
-  fyprim(rj) = 1./sqrt((rj-1.)*(jjm+1.-rj))
+  private
+  public fxysinus
 
-  fx(ri) = 2.*pi/float(iim)*(ri-0.5*float(iim)-1.)
-  ! fx    ( ri ) = 2.*pi/FLOAT(iim) * ( ri - 0.5* ( FLOAT(iim) + 1.) )
-  fxprim(ri) = 2.*pi/float(iim)
+contains
 
+  SUBROUTINE fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &
+       yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, rlonp025, &
+       xprimp025)
 
-  ! La valeur de pi est passee par le common/const/ou /const2/ .
-  ! Sinon, il faut la calculer avant d'appeler ces fonctions .
+    ! From LMDZ4/libf/dyn3d/fxysinus.F, version 1.1.1.1, 2004/05/19 12:53:06
+    ! and LMDZ4/libf/grid/fxy_sin.h, v 1.1.1.1, 2004/05/19 12:53:05
 
-  ! ----------------------------------------------------------------
-  ! Fonctions a changer eventuellement, selon x(x) et y(y) choisis .
-  ! -----------------------------------------------------------------
+    ! Author: P. Le Van
 
-  ! .....  ici, on a l'application particuliere suivante   ........
+    ! Calcul des longitudes et des latitudes pour une fonction f(x, y)
+    ! avec y = Arcsin(j).
 
-  ! **************************************
-  ! **     x = 2. * pi/iim *  X         **
-  ! **     y =      pi/jjm *  Y         **
-  ! **************************************
+    USE dimens_m, only: iim, jjm
+    USE nr_util, ONLY: pi
 
-  ! ..................................................................
-  ! ..................................................................
+    INTEGER i, j
 
+    REAL, intent(out):: rlatu(:), yprimu(:) ! (jjp1)
+    REAL, intent(out):: rlatv(:), yprimv(:) ! (jjm)
+    REAL, intent(out):: rlatu1(:) ! (jjm)
+    real, intent(out):: yprimu1(:), rlatu2(:), yprimu2(:) ! (jjm)
+    REAL, intent(out):: rlonu(:), xprimu(:), rlonv(:), xprimv(:) ! (iip1)
+    REAL, intent(out):: rlonm025(:) ! (iip1)
+    real, intent(out):: xprimm025(:), rlonp025(:), xprimp025(:) ! (iip1)
 
+    ! Local:
+    real fxprim
 
-  ! -----------------------------------------------------------------------
+    !-----------------------------------------------------------------------
 
-  ! ......  calcul  des  latitudes  et de y'   .....
+    fxprim = 2. * pi / iim
 
-  DO j = 1, jjm + 1
-    rlatu(j) = fy(float(j))
-    yprimu(j) = fyprim(float(j))
-  END DO
+    ! Calcul des latitudes et de y' 
 
+    forall(j = 1: jjm + 1)
+       rlatu(j) = fy(real(j))
+       yprimu(j) = fyprim(real(j))
+    END forall
 
-  DO j = 1, jjm
+    forall(j = 1: jjm)
+       rlatv(j) = fy(real(j) + 0.5)
+       rlatu1(j) = fy(real(j) + 0.25)
+       rlatu2(j) = fy(real(j) + 0.75)
 
-    rlatv(j) = fy(float(j)+0.5)
-    rlatu1(j) = fy(float(j)+0.25)
-    rlatu2(j) = fy(float(j)+0.75)
+       yprimv(j) = fyprim(real(j) + 0.5)
+       yprimu1(j) = fyprim(real(j) + 0.25)
+       yprimu2(j) = fyprim(real(j) + 0.75)
+    END forall
 
-    yprimv(j) = fyprim(float(j)+0.5)
-    yprimu1(j) = fyprim(float(j)+0.25)
-    yprimu2(j) = fyprim(float(j)+0.75)
+    ! Calcul des longitudes et de x' 
+    forall(i = 1: iim + 1)
+       rlonv(i) = fx(real(i))
+       rlonu(i) = fx(real(i) + 0.5)
+       rlonm025(i) = fx(real(i) - 0.25)
+       rlonp025(i) = fx(real(i) + 0.25)
 
-  END DO
+       xprimv(i) = fxprim
+       xprimu(i) = fxprim
+       xprimm025(i) = fxprim
+       xprimp025(i) = fxprim
+    END forall
 
+  END SUBROUTINE fxysinus
 
-  ! .....  calcul   des  longitudes et de  x'   .....
+  ! Fonctions à changer éventuellement, selon x(x) et y(y) choisis.
+  ! Ici, on a l'application particulière suivante :
 
-  DO i = 1, iim + 1
-    rlonv(i) = fx(float(i))
-    rlonu(i) = fx(float(i)+0.5)
-    rlonm025(i) = fx(float(i)-0.25)
-    rlonp025(i) = fx(float(i)+0.25)
+  ! x = 2 * pi / iim * X
+  ! y = pi / jjm * Y
 
-    xprimv(i) = fxprim(float(i))
-    xprimu(i) = fxprim(float(i)+0.5)
-    xprimm025(i) = fxprim(float(i)-0.25)
-    xprimp025(i) = fxprim(float(i)+0.25)
-  END DO
+  pure REAL function fy(rj)
+    USE dimens_m, only: jjm
+    REAL, intent(in):: rj
+    fy = asin(1. + 2. * ((1. - rj) / jjm))
+  end function fy
 
+  pure real function fx(ri)
+    USE dimens_m, only: iim
+    USE nr_util, ONLY: pi
+    REAL, intent(in):: ri
+    fx = 2.*pi/real(iim)*(ri-0.5*real(iim)-1.)
+  end function fx
 
-  RETURN
-END SUBROUTINE fxysinus
+  pure real function fyprim(rj)
+    USE dimens_m, only: jjm
+    REAL, intent(in):: rj
+    fyprim = 1./sqrt((rj-1.)*(jjm + 1.-rj))
+  end function fyprim
 
+end module fxysinus_m