!
!  $Author: tomita $
!  $Date: 2007/08/01 20:09:58 $
!  $Revision: 1.1.1.1 $
!
PROGRAM SoilMoistureWeeklyNCEP
   USE InputParameters, ONLY: r4, r8,i4, &
                              nferr, nfclm, nfout, nfctl,nfvgm, &
                              Imax, Jmax, Idim, Jdim,kdim,Undef, &
                              nLats,DirPreIn, DirPreOut, DirModelIn, &
                              VarName, VarNameOut,VarNameOutSN, DirMain, GrADS, Linear, &
                              InitInputParameters,Date,VarNameVeg

   USE LinearInterpolation, ONLY: gLatsL=>LatOut, &
       InitLinearInterpolation, DoLinearInterpolation

   USE AreaInterpolation, ONLY: gLatsA=>gLats, &
       InitAreaInterpolation, DoAreaInterpolation

   use MOD_NETCDF
   ! First Point of Initial Data is near South Pole and near Greenwhich
   ! First Point of Output  Data is near North Pole and near Greenwhich

   IMPLICIT NONE


   INTEGER :: Month, LRec, ios,k,i,j,LRecOut


   CHARACTER (LEN=12) :: TimeGrADS='  Z         '


   REAL (KIND=r4), PARAMETER :: porosity=0.45_r4

   CHARACTER (LEN=3), DIMENSION(12) :: MonthChar = &
             (/ 'JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', &
                'JUL', 'AUG', 'SEP', 'OCT', 'NOV', 'DEC' /)

   REAL (KIND=r8), DIMENSION (:,:)  , ALLOCATABLE :: VarField
   REAL (KIND=r8), DIMENSION (:,:)  , ALLOCATABLE :: VarFieldOutput

   INTEGER       , DIMENSION (:,:)  , ALLOCATABLE :: VegetationMask
   REAL (KIND=r4), DIMENSION (:)    , ALLOCATABLE :: Latitude

   REAL (KIND=r4), DIMENSION (:,:,:), ALLOCATABLE :: SoilMoistureWeeklyOut
   REAL (KIND=r4), DIMENSION (:,:,:), ALLOCATABLE :: SoilMoistureWeekly
   REAL (KIND=r4), DIMENSION (:,:)  , ALLOCATABLE :: SnownWeekly
   REAL (KIND=r4), DIMENSION (:,:)  , ALLOCATABLE :: SnownWeeklyOut

   INTEGER :: nv
   INTEGER :: nfinp=5    ! Standard Read In
   INTEGER :: nfprt=6    ! Standard Print Out
   INTEGER :: nfsno=10   ! To Read Unformatted Weekly SoilMoisture Data
   CHARACTER (LEN=16)   :: NCEPName='gdas1'
   CHARACTER (LEN=1024) :: DGDInp,GDASInp
   CHARACTER (LEN=2)    :: UTC
   CHARACTER (LEN=10)   :: DateLabel ='YYYYMMDDHH' ! Date: yyyymmddhh or DateLabel='        hh'
                                                   ! If Year (yyyy), Month (mm) and Day (dd) Are Unknown
   CHARACTER (len = 526) :: FILE_NAME



   CALL InitInputParameters ()

   DateLabel=Date
   DGDInp=TRIM(DirMain)//'pre/datain/ '
   UTC='00'
   UTC=DateLabel(9:10)
   GDASInp=TRIM(NCEPName)//'.T'//UTC//'Z.sfcanl.netcdf.'//DateLabel ! Input NCEP Analysis File Name

   FILE_NAME=TRIM(DGDInp)//TRIM(GDASInp)

   ! Close netcdf files

   CALL Init_Netcdf_Parameter(FILE_NAME)

   Idim=xSize
   Jdim=ySize

   IF (Linear) THEN
      CALL InitLinearInterpolation ()
   ELSE
      CALL InitAreaInterpolation ()
   END IF

   ALLOCATE (Latitude (Jmax))
   ALLOCATE (VegetationMask          (Imax,Jmax))
   ALLOCATE (SoilMoistureWeekly      (Imax,Jmax,4))
   ALLOCATE (SnownWeekly             (Imax,Jmax))
   ALLOCATE (SnownWeeklyOut          (Imax,Jmax))
   ALLOCATE (SoilMoistureWeeklyOut   (Imax,Jmax,kdim))
   ALLOCATE (VarFieldOutput(Imax,Jmax))
   ALLOCATE (VarField(Idim,Jdim))
   nv=36
   DO k=1,4

      ! Read In Input SoilMoistureWeekly
     
      CALL READ_Netcdf_Field(VarField,nfield=nVarNetcdf(nv))
      nv=nv+1
      ! Interpolate Input Regular Grid Deep Soil Temperature To Gaussian Grid on Output

      IF (Linear) THEN
         CALL DoLinearInterpolation (VarField(:,:), VarFieldOutput(:,:))
      ELSE
         CALL DoAreaInterpolation   (VarField(:,:), VarFieldOutput(:,:))
      END IF
      WRITE (UNIT=*, FMT='(A20,2G12.5)') TRIM(VAR_NAME(nv-1,1)),MINVAL(VarFieldOutput),MAXVAL(VarFieldOutput)

      DO j=1,jMax
         DO i=1,iMax
            SoilMoistureWeekly(i,j,k)=REAL(VarFieldOutput(i,j),r4)
            IF (SoilMoistureWeekly(i,j,k) <=0.0_r4)THEN
                SoilMoistureWeekly(i,j,k) = 0.9_r4
            END IF
         END DO
      END DO
   END DO   


   nv=27

   ! Read In Input Snown
     
   CALL READ_Netcdf_Field(VarField,nfield=nVarNetcdf(nv))
      nv=nv+1
      ! Interpolate Input Regular Grid Deep Soil Temperature To Gaussian Grid on Output

   IF (Linear) THEN
      CALL DoLinearInterpolation (VarField(:,:), VarFieldOutput(:,:))
   ELSE
      CALL DoAreaInterpolation   (VarField(:,:), VarFieldOutput(:,:))
   END IF
   WRITE (UNIT=*, FMT='(A20,2G12.5)') TRIM(VAR_NAME(nv-1,1)),MINVAL(VarFieldOutput),MAXVAL(VarFieldOutput)

   DO j=1,jMax
      DO i=1,iMax
         SnownWeekly(i,j)=REAL(VarFieldOutput(i,j),r4)
         SnownWeeklyOut(i,j)=SnownWeekly(i,j)*1000
      END DO
   END DO

   ! Close netcdf files

   CALL Finalize_Netcdf_Parameter()

   ! Read  Adjusted Interpolated VegetationMask
   INQUIRE (IOLENGTH=LRecOut) VegetationMask
   OPEN (UNIT=nfvgm, FILE=TRIM(DirMain)//DirModelIn//TRIM(VarNameVeg)//nLats, &
         FORM='UNFORMATTED', ACCESS='DIRECT', RECL=LRecOut, ACTION='READ', &
         STATUS='OLD', IOSTAT=ios)
   IF (ios /= 0) THEN
      WRITE (UNIT=nferr, FMT='(3A,I4)') &
            ' ** (Error) ** Open file ', &
              TRIM(DirMain)//DirModelIn//TRIM(VarNameVeg)//nLats, &
            ' returned IOStat = ', ios
      STOP  ' ** (Error) **'
   END IF
   READ (UNIT=nfvgm, REC=1) VegetationMask
   CLOSE (UNIT=nfvgm)
   
   
   
   DO k=1,kdim
      DO j=1,jMax
         DO i=1,Imax
           !  volumetric soil moisture 
           IF(k >=1 .and. k<=2)THEN
            SoilMoistureWeeklyOut(i,j,k)= MIN(( (SoilMoistureWeekly(i,j,1))),1.0)
           ELSE IF(k >2 .and. k<=4)THEN
            SoilMoistureWeeklyOut(i,j,k)= MIN(( (SoilMoistureWeekly(i,j,2))),1.0)
           ELSE IF(k >4 .and. k<=6)THEN
            SoilMoistureWeeklyOut(i,j,k)= MIN(( (SoilMoistureWeekly(i,j,3))),1.0)
           ELSE IF(k >6 .and. k<=8)THEN
            SoilMoistureWeeklyOut(i,j,k)= MIN(( (SoilMoistureWeekly(i,j,4))),1.0)
           ELSE
            SoilMoistureWeeklyOut(i,j,k)= MIN(0.25*((SoilMoistureWeekly(i,j,4)) + (SoilMoistureWeekly(i,j,3)) +   &
                                                    (SoilMoistureWeekly(i,j,2)) + (SoilMoistureWeekly(i,j,1))),1.0)
           END IF

           if(VegetationMask(i,j) >= 1_i4 )THEN 
              SoilMoistureWeeklyOut(i,j,k)=MAX(SoilMoistureWeeklyOut(i,j,k)/porosity,0.0000001)
              if(SoilMoistureWeeklyOut(i,j,k) == 0.0) SoilMoistureWeeklyOut(i,j,k)=1.0000
              if(SoilMoistureWeeklyOut(i,j,k) >= 1.0) SoilMoistureWeeklyOut(i,j,k)=1.0000

           ELSE
              SoilMoistureWeeklyOut(i,j,k)=UNDEF
           endif 

         END DO
      END DO   
   END DO

   
   ! Generates Soil Moisture Input File for the AGCM Model

   INQUIRE (IOLENGTH=LRecOut) SoilMoistureWeeklyOut(1:Imax,1:Jmax,1)
   OPEN (FILE=TRIM(DirMain)//DirModelIn//TRIM(VarNameOut)//'.'//TRIM(Date(1:8))//nLats, &
         UNIT=nfout, FORM='UNFORMATTED', ACCESS='DIRECT', RECL=LRecOut, &
         ACTION='WRITE', STATUS='REPLACE', IOSTAT=ios)
   IF (ios /= 0) THEN
      WRITE (UNIT=nferr, FMT='(3A,I4)') &
            ' ** (Error) ** Open file ', &
              TRIM(DirMain)//DirModelIn//TRIM(VarNameOut)//nLats, &
            ' returned IOStat = ', ios
      STOP  ' ** (Error) **'
   END IF
   DO k=1,kdim
      WRITE (UNIT=nfout, REC=k) SoilMoistureWeeklyOut(1:Imax,1:Jmax,kdim+1-k)
   END DO
   CLOSE (UNIT=nfout)

   ! Generates Snow deep Input File for the AGCM Model

   INQUIRE (IOLENGTH=LRecOut) SnownWeeklyOut(1:Imax,1:Jmax)
   OPEN (FILE=TRIM(DirMain)//DirModelIn//TRIM(VarNameOutSN)//'.'//TRIM(Date(1:8))//nLats, &
         UNIT=nfout, FORM='UNFORMATTED', ACCESS='DIRECT', RECL=LRecOut, &
         ACTION='WRITE', STATUS='REPLACE', IOSTAT=ios)
   IF (ios /= 0) THEN
      WRITE (UNIT=nferr, FMT='(3A,I4)') &
            ' ** (Error) ** Open file ', &
              TRIM(DirMain)//DirModelIn//TRIM(VarNameOutSN)//nLats, &
            ' returned IOStat = ', ios
      STOP  ' ** (Error) **'
   END IF
   WRITE (UNIT=nfout, REC=1) SnownWeeklyOut
   CLOSE (UNIT=nfout)

   IF (GrADS) THEN
      TimeGrADS(1:2)       =Date(9:10)! hour
      TimeGrADS(4:5)       =Date(7:8) ! day
      TimeGrADS(9:12)      =Date(1:4) ! year
      READ (Date(5:6), FMT='(I2)') Month
      TimeGrADS(6:8)=MonthChar(Month)

      OPEN (UNIT=nfctl, FILE=TRIM(DirMain)//TRIM(DirPreOut)//TRIM(VarNameOut)//'.'//Date(1:8)//'.ctl', &
            FORM='FORMATTED', ACCESS='SEQUENTIAL', ACTION='WRITE', &
            STATUS='REPLACE', IOSTAT=ios)
      IF (ios /= 0) THEN
         WRITE (UNIT=nferr, FMT='(3A,I4)') &
               ' ** (Error) ** Open file ', &
                 TRIM(DirMain)//TRIM(DirPreOut)//TRIM(VarName)//'.'//Date(1:8)//'.ctl', &
               ' returned IOStat = ', ios
         STOP  ' ** (Error) **'
      END IF
      WRITE (UNIT=nfctl, FMT='(A)') 'DSET '// &
             TRIM(DirMain)//DirModelIn//TRIM(VarNameOut)//'.'//TRIM(Date(1:8))//nLats
      WRITE (UNIT=nfctl, FMT='(A)') '*'
      WRITE (UNIT=nfctl, FMT='(A)') 'OPTIONS YREV BIG_ENDIAN'
      WRITE (UNIT=nfctl, FMT='(A)') '*'
      WRITE (UNIT=nfctl, FMT='(A)') 'UNDEF -999.0'
      WRITE (UNIT=nfctl, FMT='(A)') '*'
      WRITE (UNIT=nfctl, FMT='(A)') 'TITLE SoilMoistureWeekly on a Gaussian Grid'
      WRITE (UNIT=nfctl, FMT='(A)') '*'
      WRITE (UNIT=nfctl, FMT='(A,I5,A,F8.3,F15.10)') &
                          'XDEF ',Imax,' LINEAR ',0.0_r8,360.0_r8/REAL(Imax,r8)
      WRITE (UNIT=nfctl, FMT='(A,I5,A)') 'YDEF ',Jmax,' LEVELS '
      IF (Linear) THEN
         Latitude(1:Jmax:1) = gLatsL(Jmax:1:-1)
         WRITE (UNIT=nfctl, FMT='(8F10.5)')Latitude
      ELSE
         Latitude(1:Jmax:1) = gLatsA(Jmax:1:-1)
         WRITE (UNIT=nfctl, FMT='(8F10.5)') Latitude
      END IF
      WRITE (UNIT=nfctl, FMT='(A)') 'ZDEF 8 LEVELS 1 2 3 4 5 6 7 8'
      WRITE (UNIT=nfctl, FMT='(A)') 'TDEF 1 LINEAR '//TimeGrADS//' 1MO'
      WRITE (UNIT=nfctl, FMT='(A)') 'VARS 1'
      WRITE (UNIT=nfctl, FMT='(A)') 'SoilMoisture 8 99 SoilMoistureWeekly [kg/m3]'
      WRITE (UNIT=nfctl, FMT='(A)') 'ENDVARS'
      CLOSE (UNIT=nfctl)
   END IF

PRINT *, "*** SoilMoistureWeeklyNCEP ENDS NORMALLY ***"

   IF (GrADS) THEN
   TimeGrADS(1:2)=Date(9:10)
   TimeGrADS(4:5)=Date(7:8)
   TimeGrADS(9:12)=Date(1:4)
   READ (Date(5:6), FMT='(I2)') Month
   TimeGrADS(6:8)=MonthChar(Month)

      ! Write GrADS Control File
      OPEN (FILE=TRIM(DirMain)//DirPreOut//TRIM(VarNameOutSN)//'.'//TRIM(Date(1:8))//nLats//'.ctl', &
            UNIT=nfctl, FORM='FORMATTED', ACCESS='SEQUENTIAL', &
            ACTION='WRITE', STATUS='REPLACE', IOSTAT=ios)
      IF (ios /= 0) THEN
         WRITE (UNIT=nferr, FMT='(3A,I4)') &
               ' ** (Error) ** Open file ', &
                 TRIM(DirMain)//DirPreOut//TRIM(VarNameOutSN)//'.'//TRIM(Date(1:8))//nLats//'.ctl', &
               ' returned IOStat = ', ios
         STOP  ' ** (Error) **'
      END IF
      WRITE (UNIT=nfctl, FMT='(A)') 'DSET '// &
             TRIM(DirMain)//DirModelIn//TRIM(VarNameOutSN)//'.'//TRIM(Date(1:8))//nLats
      WRITE (UNIT=nfctl, FMT='(A)') '*'
      WRITE (UNIT=nfctl, FMT='(A)') 'OPTIONS YREV BIG_ENDIAN'
      WRITE (UNIT=nfctl, FMT='(A)') '*'
      WRITE (UNIT=nfctl, FMT='(A)') 'UNDEF -999.0'
      WRITE (UNIT=nfctl, FMT='(A)') '*'
      WRITE (UNIT=nfctl, FMT='(A)') 'TITLE SNOWWeekly on a Gaussian Grid'
      WRITE (UNIT=nfctl, FMT='(A)') '*'
      WRITE (UNIT=nfctl, FMT='(A,I5,A,F8.3,F15.10)') &
                          'XDEF ',Imax,' LINEAR ',0.0_r8,360.0_r8/REAL(Imax,r8)
      WRITE (UNIT=nfctl, FMT='(A,I5,A)') 'YDEF ',Jmax,' LEVELS '
      IF (Linear) THEN
         Latitude(1:Jmax:1) = gLatsL(Jmax:1:-1)
         WRITE (UNIT=nfctl, FMT='(8F10.5)') Latitude
      ELSE
         Latitude(1:Jmax:1) = gLatsA(Jmax:1:-1)
         WRITE (UNIT=nfctl, FMT='(8F10.5)') Latitude
      END IF
      WRITE (UNIT=nfctl, FMT='(A)') 'ZDEF 1 LEVELS 1000'
      WRITE (UNIT=nfctl, FMT='(A)') 'TDEF 1 LINEAR '//TimeGrADS//' 1MO'
      WRITE (UNIT=nfctl, FMT='(A)') 'VARS 1'
      WRITE (UNIT=nfctl, FMT='(A)') 'snow 0 99 SNOWWeekly [kg/m3]'
      WRITE (UNIT=nfctl, FMT='(A)') 'ENDVARS'
      CLOSE (UNIT=nfctl)
   END IF

PRINT *, "*** SNOWWeekly ENDS NORMALLY ***"

CONTAINS

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

  INTEGER, INTENT(IN)  :: iunit  !#TO deve ser kind default
  INTEGER, INTENT(IN)  :: n
  REAL(KIND=r4), INTENT(OUT)    :: aa(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
  REAL(KIND=r8)     :: a(n)

  vc='0'
  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
  aa=REAL(a,kind=r4)
 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 PROGRAM SoilMoistureWeeklyNCEP
