!
!#     #
!#     #   #####     #    #        ####
!#     #     #       #    #       #
!#     #     #       #    #        ####
!#     #     #       #    #            #
!#     #     #       #    #       #    #
! #####      #       #    ######   ####
!
!
!  $Author: pkubota $
!  $Date: 2008/09/23 17:51:54 $
!  $Revision: 1.9 $
!
MODULE Utils
  USE Constants, Only: &
       r8,r4,i4,             &
       i8

  IMPLICIT NONE
  PUBLIC :: tmstmp2
  PUBLIC :: IJtoIBJB
  PUBLIC :: LinearIJtoIBJB
  PUBLIC :: NearestIJtoIBJB
  PUBLIC :: SeaMaskIJtoIBJB
  PUBLIC :: SplineIJtoIBJB
  PUBLIC :: AveBoxIJtoIBJB
  PUBLIC :: FreqBoxIJtoIBJB
  PUBLIC :: lati
  PUBLIC :: vfirec
  PUBLIC :: julday
 
  INTERFACE IJtoIBJB
     MODULE PROCEDURE &
          IJtoIBJB_R, IJtoIBJB_I, &
          IJtoIBJB3_R, IJtoIBJB3_I
  END INTERFACE

  INTERFACE LinearIJtoIBJB
     MODULE PROCEDURE LinearIJtoIBJB_R2D
  END INTERFACE

  INTERFACE NearestIJtoIBJB
     MODULE PROCEDURE &
          NearestIJtoIBJB_I2D, NearestIJtoIBJB_R2D, &
          NearestIJtoIBJB_I3D, NearestIJtoIBJB_R3D
  END INTERFACE

  INTERFACE SeaMaskIJtoIBJB
     MODULE PROCEDURE SeaMaskIJtoIBJB_R2D
  END INTERFACE

  INTERFACE SplineIJtoIBJB
     MODULE PROCEDURE SplineIJtoIBJB_R2D
  END INTERFACE

  INTERFACE AveBoxIJtoIBJB
     MODULE PROCEDURE AveBoxIJtoIBJB_R2D
  END INTERFACE

  INTERFACE FreqBoxIJtoIBJB
     MODULE PROCEDURE FreqBoxIJtoIBJB_I2D, FreqBoxIJtoIBJB_R2D
  END INTERFACE

 
  REAL(KIND=r8), ALLOCATABLE :: lati(:)


CONTAINS
  SUBROUTINE InitUtils(iMax,jMax)
    INTEGER, INTENT(in ) :: iMax
    INTEGER, INTENT(in ) :: jMax
    ALLOCATE(lati(jMax));lati=0.0_r8
  END SUBROUTINE InitUtils
  !
  !***************************************************************************
  !                      (imonth,iday,iyear)
REAL(KIND=r8) FUNCTION julday (imonth,iday,iyear,tod)
  IMPLICIT NONE
  INTEGER, INTENT(IN   ) :: imonth
  INTEGER, INTENT(IN   ) :: iday
  INTEGER, INTENT(IN   ) :: iyear
  REAL(KIND=r8)   , INTENT(IN   ) :: tod
  !
  ! compute the julian day from a normal date
  !
  julday= iday  &
       + MIN(1,MAX(0,imonth-1))*31  &
       + MIN(1,MAX(0,imonth-2))*(28+(1-MIN(1,MOD(iyear,4))))  &
       + MIN(1,MAX(0,imonth-3))*31  &
       + MIN(1,MAX(0,imonth-4))*30  &
       + MIN(1,MAX(0,imonth-5))*31  &
       + MIN(1,MAX(0,imonth-6))*30  &
       + MIN(1,MAX(0,imonth-7))*31  &
       + MIN(1,MAX(0,imonth-8))*31  &
       + MIN(1,MAX(0,imonth-9))*30  &
       + MIN(1,MAX(0,imonth-10))*31  &
       + MIN(1,MAX(0,imonth-11))*30  &
       + MIN(1,MAX(0,imonth-12))*31  &
       + tod/86400.0

END FUNCTION julday

  SUBROUTINE tmstmp2(id, ifday, tod, ihr, iday, mon, iyr)
    !
    !
    !==========================================================================
    !    id(4).......date of current data
    !                id(1)....hour(00/12)
    !                id(2)....month
    !                id(3)....day of month
    !                id(4)....year
    !    ifday.......model forecast day
    !    tod.........todx=tod+swint*f3600, model forecast time of
    !                day in seconds
    !                swint....sw subr. call interval in hours
    !                swint has to be less than or equal to trint
    !                              and mod(trint,swint)=0
    !                f3600=3.6e3
    !    ihr.........hour(00/12)
    !    iday........day of month
    !    mon.........month
    !    iyr.........year
    !    yrl.........length of year in days
    !    monl(12)....length of each month in days
    !==========================================================================
    !

    INTEGER, INTENT(in ) :: id(4)
    INTEGER, INTENT(in ) :: ifday
    REAL(KIND=r8),    INTENT(in ) :: tod
    INTEGER, INTENT(out) :: ihr
    INTEGER, INTENT(out) :: iday
    INTEGER, INTENT(out) :: mon
    INTEGER, INTENT(out) :: iyr

    INTEGER :: kday
    INTEGER :: idaymn
    REAL(KIND=r8)    :: ctim
    REAL(KIND=r8)    :: hrmodl
    INTEGER :: monl(12)

    REAL(KIND=r8), PARAMETER :: yrl =   365.2500
    REAL(KIND=r8), PARAMETER ::  ep = .015625
    DATA MONL/31,28,31,30,31,30,&
         31,31,30,31,30,31/

    ctim=tod+id(1)*3600.0_r8

    IF (ctim >= 86400.e0_r8) THEN
       kday=1
       ctim=ctim-86400.e0_r8
    ELSE
       kday=0
    END IF
    !
    !     adjust time to reduce round off error in divsion
    !
    iday = id(3) + ifday + kday
    hrmodl = (ctim+ep)/3600.0_r8
    ihr = hrmodl
    mon = id(2)
    iyr = id(4)
    DO
       idaymn = monl(mon)
       IF (yrl == 365.25e0_r8 .AND. MOD(iyr,4) == 0 .AND. mon == 2) &
            idaymn=29
       IF (iday <= idaymn) RETURN
       iday = iday - idaymn
       mon = mon + 1
       IF (mon < 13) CYCLE
       mon = 1
       iyr = iyr + 1
    END DO
  END SUBROUTINE tmstmp2




  SUBROUTINE LinearIJtoIBJB_R2D(FieldIn,FieldOut)
    REAL(KIND=r8), INTENT(IN  ) :: FieldIn (: ,:)
    REAL(KIND=r8), INTENT(OUT ) :: FieldOut(:,:)
    INTEGER            :: ib,i
    INTEGER            :: jb,j

    CHARACTER(LEN=*), PARAMETER :: h="**LinearIJtoIBJB**"

       DO jb = 1, SIZE(FieldOut,2)
          DO ib = 1, SIZE(FieldOut,1)
             i=ib
	     j=jb
	     FieldOut(ib,jb)=FieldIn(i,j)
          END DO
       END DO
       
  END SUBROUTINE LinearIJtoIBJB_R2D


  SUBROUTINE NearestIJtoIBJB_I2D(FieldIn,FieldOut)
    INTEGER, INTENT(IN  ) :: FieldIn (:,:)
    INTEGER, INTENT(OUT ) :: FieldOut(:,:)
    INTEGER            :: i
    INTEGER            :: ib
    INTEGER            :: j
    INTEGER            :: jb

    CHARACTER(LEN=*), PARAMETER :: h="**NearestIJtoIBJB**"

       DO jb = 1, SIZE(FieldOut,2)
          DO ib = 1, SIZE(FieldOut,1)
             i=ib
	     j=jb
	     FieldOut(ib,jb)=FieldIn(i,j)
          END DO
       END DO
  END SUBROUTINE NearestIJtoIBJB_I2D



  SUBROUTINE NearestIJtoIBJB_R2D(FieldIn,FieldOut)
    REAL(KIND=r8), INTENT(IN  ) :: FieldIn (:,:)
    REAL(KIND=r8), INTENT(OUT ) :: FieldOut(:,:)

    INTEGER            :: i
    INTEGER            :: ib
    INTEGER            :: j
    INTEGER            :: jb

    CHARACTER(LEN=*), PARAMETER :: h="**NearestIJtoIBJB**"

       DO jb = 1, SIZE(FieldOut,2)
          DO ib = 1, SIZE(FieldOut,1)
             i=ib
	     j=jb
	     FieldOut(ib,jb)=FieldIn(i,j)
          END DO
       END DO
  END SUBROUTINE NearestIJtoIBJB_R2D

  ! 3D version by hmjb
  SUBROUTINE NearestIJtoIBJB_I3D(FieldIn,FieldOut)
    INTEGER, INTENT(IN  ) :: FieldIn (:,:,:)
    INTEGER, INTENT(OUT ) :: FieldOut(:,:,:)
    INTEGER            :: i,k
    INTEGER            :: ib
    INTEGER            :: j
    INTEGER            :: jb

    CHARACTER(LEN=*), PARAMETER :: h="**NearestIJtoIBJB**"

       DO jb = 1, SIZE(FieldOut,2)
          DO ib = 1, SIZE(FieldOut,1)
             i=ib
	     j=jb
	     FieldOut(ib,:,jb)=FieldIn(i,:,j)
          END DO
       END DO
  END SUBROUTINE NearestIJtoIBJB_I3D


  ! 3D version by hmjb
  SUBROUTINE NearestIJtoIBJB_R3D(FieldIn,FieldOut)
    REAL(KIND=r8), INTENT(IN  ) :: FieldIn (: ,:,:)
    REAL(KIND=r8), INTENT(OUT ) :: FieldOut(:,:,:)

    INTEGER            :: i,k
    INTEGER            :: ib
    INTEGER            :: j
    INTEGER            :: jb

    CHARACTER(LEN=*), PARAMETER :: h="**NearestIJtoIBJB**"

    DO jb = 1, SIZE(FieldOut,2)
       DO ib = 1, SIZE(FieldOut,1)
    	  i=ib
     	  j=jb
     	  FieldOut(ib,:,jb)=FieldIn(i,:,j)
       END DO
    END DO
       
  END SUBROUTINE NearestIJtoIBJB_R3D





  SUBROUTINE SeaMaskIJtoIBJB_R2D(FieldIn,FieldOut)
    REAL(KIND=r8)   , INTENT(IN  ) :: FieldIn (:,:)
    REAL(KIND=r8)   , INTENT(OUT ) :: FieldOut(:,:)

    INTEGER            :: i
    INTEGER            :: ib
    INTEGER            :: j
    INTEGER            :: jb

    CHARACTER(LEN=*), PARAMETER :: h="**SeaMaskIJtoIBJB**"

    DO jb = 1, SIZE(FieldOut,2)
       DO ib = 1, SIZE(FieldOut,1)
    	  i=ib
     	  j=jb
     	  FieldOut(ib,jb)=FieldIn(i,j)
       END DO
    END DO

  END SUBROUTINE SeaMaskIJtoIBJB_R2D




  SUBROUTINE SplineIJtoIBJB_R2D(FieldIn,FieldOut)
    REAL(KIND=r8), INTENT(IN  ) :: FieldIn (:,:)
    REAL(KIND=r8), INTENT(OUT ) :: FieldOut(:,:)
    INTEGER            :: i
    INTEGER            :: ib
    INTEGER            :: j
    INTEGER            :: jb

    CHARACTER(LEN=*), PARAMETER :: h="**SplineIJtoIBJB**"

    PRINT *, h

    DO jb = 1, SIZE(FieldOut,2)
       DO ib = 1, SIZE(FieldOut,1)
    	  i=ib
     	  j=jb
     	  FieldOut(ib,jb)=FieldIn(i,j)
       END DO
    END DO
    
  END SUBROUTINE SplineIJtoIBJB_R2D



 SUBROUTINE AveBoxIJtoIBJB_R2D(FieldIn,FieldOut)
    REAL(KIND=r8)   , INTENT(IN  ) :: FieldIn (:,:)
    REAL(KIND=r8)   , INTENT(OUT ) :: FieldOut(:,:)

    INTEGER            :: i
    INTEGER            :: ib
    INTEGER            :: j
    INTEGER            :: jb

    CHARACTER(LEN=*), PARAMETER :: h="**AveBoxIJtoIBJB_R2D**"

    DO jb = 1, SIZE(FieldOut,2)
       DO ib = 1, SIZE(FieldOut,1)
    	  i=ib
     	  j=jb
     	  FieldOut(ib,jb)=FieldIn(i,j)
       END DO
    END DO
  END SUBROUTINE AveBoxIJtoIBJB_R2D


  SUBROUTINE FreqBoxIJtoIBJB_I2D(FieldIn,FieldOut)
    INTEGER(KIND=i8), INTENT(IN  ) :: FieldIn (:,:)
    INTEGER(KIND=i8), INTENT(OUT ) :: FieldOut(:,:)

    INTEGER            :: i
    INTEGER            :: ib
    INTEGER            :: j
    INTEGER            :: jb

    CHARACTER(LEN=*), PARAMETER :: h="**FreqBoxIJtoIBJB**"

    DO jb = 1, SIZE(FieldOut,2)
       DO ib = 1, SIZE(FieldOut,1)
    	  i=ib
     	  j=jb
     	  FieldOut(ib,jb)=FieldIn(i,j)
       END DO
    END DO
  END SUBROUTINE FreqBoxIJtoIBJB_I2D

  SUBROUTINE FreqBoxIJtoIBJB_R2D(FieldIn,FieldOut)
    REAL(KIND=r8)   , INTENT(IN  ) :: FieldIn (:,:)
    REAL(KIND=r8)   , INTENT(OUT ) :: FieldOut(:,:)

    INTEGER            :: i
    INTEGER            :: ib
    INTEGER            :: j
    INTEGER            :: jb

    CHARACTER(LEN=*), PARAMETER :: h="**FreqBoxIJtoIBJB**"

    DO jb = 1, SIZE(FieldOut,2)
       DO ib = 1, SIZE(FieldOut,1)
    	  i=ib
     	  j=jb
     	  FieldOut(ib,jb)=FieldIn(i,j)
       END DO
    END DO
  END SUBROUTINE FreqBoxIJtoIBJB_R2D

 !
 ! maps (i,j) into (ib,jb)
 !
 SUBROUTINE IJtoIBJB_R(FieldIn,FieldOut)
    REAL(KIND=r8)   , INTENT(IN  ) :: FieldIn (:,:)
    REAL(KIND=r8)   , INTENT(OUT ) :: FieldOut(:,:)
    INTEGER               :: i
    INTEGER               :: j
    INTEGER               :: ib
    INTEGER               :: jb
    DO jb = 1, SIZE(FieldOut,2)
       DO ib = 1, SIZE(FieldOut,1)
    	  i=ib
     	  j=jb
     	  FieldOut(ib,jb)=FieldIn(i,j)
       END DO
    END DO
 END SUBROUTINE IJtoIBJB_R

 ! 3D version by hmjb
 SUBROUTINE IJtoIBJB3_R(FieldIn,FieldOut)
    REAL(KIND=r8)   , INTENT(IN  ) :: FieldIn (:,:,:)
    REAL(KIND=r8)   , INTENT(OUT ) :: FieldOut(:,:,:)
    INTEGER               :: i
    INTEGER               :: j
    INTEGER               :: ib
    INTEGER               :: jb
    DO jb = 1, SIZE(FieldOut,2)
       DO ib = 1, SIZE(FieldOut,1)
    	  i=ib
     	  j=jb
     	  FieldOut(ib,:,jb)=FieldIn(i,:,j)
       END DO
    END DO
 END SUBROUTINE IJtoIBJB3_R

 !
 ! maps (i,j) into (ib,jb)
 !
 SUBROUTINE IJtoIBJB_I(FieldIn,FieldOut)
    INTEGER(KIND=i8), INTENT(IN  ) :: FieldIn (:,:)
    INTEGER(KIND=i8), INTENT(OUT ) :: FieldOut(:,:)
    INTEGER               :: i
    INTEGER               :: j
    INTEGER               :: ib
    INTEGER               :: jb
    DO jb = 1, SIZE(FieldOut,2)
       DO ib = 1, SIZE(FieldOut,1)
    	  i=ib
     	  j=jb
     	  FieldOut(ib,jb)=FieldIn(i,j)
       END DO
    END DO
 END SUBROUTINE IJtoIBJB_I
 ! 3D version by hmjb
 SUBROUTINE IJtoIBJB3_I(FieldIn,FieldOut)
    INTEGER(KIND=i8), INTENT(IN  ) :: FieldIn (:,:,:)
    INTEGER(KIND=i8), INTENT(OUT ) :: FieldOut(:,:,:)
    INTEGER               :: i
    INTEGER               :: j
    INTEGER               :: ib
    INTEGER               :: jb
    DO jb = 1, SIZE(FieldOut,2)
       DO ib = 1, SIZE(FieldOut,1)
    	  i=ib
     	  j=jb
     	  FieldOut(ib,:,jb)=FieldIn(i,:,j)
       END DO
    END DO
 END SUBROUTINE IJtoIBJB3_I




!
!------------------------------- VFORMAT ----------------------------------
!
 SUBROUTINE vfirec(iunit,a,n,type)

  INTEGER, INTENT(IN)  :: iunit  !#TO deve ser kind default
  INTEGER, INTENT(IN)  :: n
  REAL(KIND=r8), INTENT(OUT)    :: a(n)
  CHARACTER(len=* ), INTENT(IN) :: type
  !
  ! local
  !
  CHARACTER(len=1 ) :: vc(0:63)
  CHARACTER(len=80) :: line
  CHARACTER(len=1 ) :: cs
  INTEGER           :: ich0
  INTEGER           :: ich9
  INTEGER           :: ichcz
  INTEGER           :: ichca
  INTEGER           :: ichla
  INTEGER           :: ichlz
  INTEGER           :: i
  INTEGER           :: nvalline
  INTEGER           :: nchs
  INTEGER           :: ic
  INTEGER           :: ii
  INTEGER           :: isval
  INTEGER           :: iii
  INTEGER           :: ics
  INTEGER           :: nn
  INTEGER           :: nbits
  INTEGER           :: nc
  REAL(KIND=r8)              :: bias
  REAL(KIND=r8)              :: fact
  REAL(KIND=r8)              :: facti
  REAL(KIND=r8)              :: scfct

  IF (vc(0).ne.'0') CALL vfinit(vc)

  ich0 =ichar('0')
  ich9 =ichar('9')
  ichcz=ichar('Z')
  ichlz=ichar('z')
  ichca=ichar('A')
  ichla=ichar('a')

  READ (iunit,'(2i8,2e20.10)')nn,nbits,bias,fact

  IF (nn.ne.n) THEN
    PRINT*,' Word count mismatch on vfirec record '
    PRINT*,' Words on record - ',nn
    PRINT*,' Words expected  - ',n
    STOP 'vfirec'
  END IF

  nvalline=(78*6)/nbits
  nchs=nbits/6

  DO i=1,n,nvalline
    READ(iunit,'(a78)') line
    ic=0
    DO ii=i,i+nvalline-1
      isval=0
      IF(ii.gt.n) EXIT
      DO iii=1,nchs
         ic=ic+1
         cs=line(ic:ic)
         ics=ichar(cs)
         IF (ics.le.ich9) THEN
            nc=ics-ich0
         ELSE IF (ics.le.ichcz) THEN
            nc=ics-ichca+10
         ELSE
            nc=ics-ichla+36
         END IF
         isval=ior(ishft(nc,6*(nchs-iii)),isval)
      END DO ! loop iii
        a(ii)=isval
    END DO ! loop ii

  END DO ! loop i

  facti=1.0_r8/fact

  IF (type.eq.'LIN') THEN
    DO i=1,n

      a(i)=a(i)*facti-bias

      !print*,'VFM=',i,a(i)
    END DO
  ELSE IF (type.eq.'LOG') THEN
    scfct=2.0_r8**(nbits-1)
    DO i=1,n
        a(i)=sign(1.0_r8,a(i)-scfct)  &
           *(10.0_r8**(abs(20.0_r8*(a(i)/scfct-1.0_r8))-10.0_r8))
    END DO
  END IF
 END SUBROUTINE vfirec
!--------------------------------------------------------
 SUBROUTINE vfinit(vc)
   CHARACTER(len=1), INTENT(OUT  ) :: vc   (*)
   CHARACTER(len=1)                :: vcscr(0:63)
   INTEGER                         :: n

   DATA vcscr/'0','1','2','3','4','5','6','7','8','9'   &
              ,'A','B','C','D','E','F','G','H','I','J'  &
              ,'K','L','M','N','O','P','Q','R','S','T'  &
              ,'U','V','W','X','Y','Z','a','b','c','d'  &
              ,'e','f','g','h','i','j','k','l','m','n'  &
              ,'o','p','q','r','s','t','u','v','w','x'  &
              ,'y','z','{','|'/

  DO n=0,63
      vc(n)=vcscr(n)
  END DO
 END SUBROUTINE vfinit




END MODULE Utils