!--Derived type for Date/time handle
 use datedef

 Implicit none

!Dimension defines
 Integer :: Nx, Ny, Nsig
 Integer,parameter :: Np=4

 Real  ::   plev(Np)
 Data    plev/850.,700.,500.,300./

 Real, allocatable,dimension(:) :: Rx
 Real :: Rx0, Dx
 Real, allocatable, dimension(:) :: Ry

! Variables on sigma levels
 Real, allocatable,dimension(:) ::  sig,sigf,sp1d
 Real, allocatable,dimension(:,:,:) :: u,v,w,t,q,qc,h
 Real, allocatable,dimension(:,:) :: ht, ps, rain, tgb, swt, rno

! Variables on Pressure levels and sea level, ptp: potential temp
 Real, allocatable,dimension(:,:,:) :: up,vp,wp,tp,hp
 Real, allocatable,dimension(:,:) :: Slp1, Slp2
 Real :: vartemp

! Date and time and data file name
 Type(datetime) :: s_date,c_date,e_date,temp_date
 Character(len=20) :: input_dir,inputfile,output_dir
 integer, parameter :: fu=111, fv=112, fw=113, ft=114,fh=115,fslp=116  ! File ID for output
 integer :: lrec, orec1, orec2

!--Local varialbes
!--Variables of header
 integer mdate0,ibltyp,icup,ipptls,iboudy,il,jl,kl
 real   dxsp,ptsp,clat,clon,plat,plon
 character*6 proj
 real    dto,dtb,dtr,dtc, ptop
 integer iotyp,idate
 
 integer :: i,j,k,i1,j1

!Namelist
 namelist /Datetimes/ s_date,e_date
 namelist /Domains/ Nx,Rx0,Dx,Ny,Nsig, input_dir,output_dir


 open(1,file="info.regcm")
 read(1,datetimes)
 read(1,domains)
 close(1) 

write(*,datetimes)
write(*,domains)

write(*,*) "Allocating fields..."
allocate(Rx(Nx)); allocate(Ry(Ny))

allocate(sigf(Nsig+1));allocate(sp1d(Nsig+1));allocate(sig(Nsig))
allocate(U(Nx,Ny,Nsig));allocate(V(Nx,Ny,Nsig));allocate(W(Nx,Ny,Nsig));allocate(T(Nx,Ny,Nsig))
allocate(q(Nx,Ny,Nsig));allocate(qc(Nx,Ny,Nsig));allocate(h(Nx,Ny,Nsig));
allocate(ht(Nx,Ny));allocate(ps(Nx,Ny));allocate(rain(Nx,Ny));allocate(tgb(Nx,Ny));allocate(swt(Nx,Ny));allocate(rno(Nx,Ny));

allocate(Up(Nx,Ny,Np));allocate(Vp(Nx,Ny,Np));allocate(Wp(Nx,Ny,Np));allocate(Tp(Nx,Ny,Np));allocate(Hp(Nx,Ny,Np));
allocate(Slp1(Nx,Ny));allocate(Slp2(Nx,Ny));

!--Initialized dimensions
do i=1,Nx 
  Rx(i)=Rx0 + (i-1)*Dx
enddo
open(1,file="info.ry")
do j=1,Ny
  read(1,*)Ry(j)
enddo
close(1) 

write(*,*) "Reading header..."

!File list
! Read header information
 open(1,file=trim(input_dir)//'/OUT_HEAD',form='unformatted',access='direct',recl=Ny*Nx*4)
 read(1,rec=1) mdate0,ibltyp,icup,ipptls,iboudy                      &
        , il,jl,kl,sp1d,dxsp,ptsp,clat,clon,plat,plon,proj         &
        , dto,dtb,dtr,dtc,iotyp
 if(il-2.ne.Ny.or.jl-2.ne.Nx.or.kl.ne.Nsig) then
    write(*,*) 'domain parameters are not consistent'
    write(*,*) 'il,jl,kl,Ny,Nx,Nsig',il,jl,kl,Ny,Nx,Nsig
    stop
 endif
 read(1,rec=2) ((ht(i,j),i=1,Nx),j=1,Ny)

 do k=1,Nsig+1
    sigf(k)=sp1d(Nsig+2-k)
 enddo

 do k=1,Nsig
    sig(k) = 0.5*(sigf(k)+sigf(k+1))
    write(*,*)k,sig(k)
 enddo

 ptop=ptsp*10.
 close(1)
!--end read header

write(*,*) "Prepare reading data.."

!output files
open(fu,file=trim(output_dir)//"/u.dat",form='unformatted',recl=Ny*Nx*4,access='direct')
open(fv,file=trim(output_dir)//"/v.dat",form='unformatted',recl=Ny*Nx*4,access='direct')
open(fw,file=trim(output_dir)//"/w.dat",form='unformatted',recl=Ny*Nx*4,access='direct')
open(ft,file=trim(output_dir)//"/t.dat",form='unformatted',recl=Ny*Nx*4,access='direct')
open(fh,file=trim(output_dir)//"/h.dat",form='unformatted',recl=Ny*Nx*4,access='direct')
open(fslp,file=trim(output_dir)//"/slp.dat",form='unformatted',recl=Ny*Nx*4,access='direct')
orec1=1   ! For 3D fields
orec2=1   ! For 2D fields 


! Loop through time and read data on sigma level
 c_date = s_date
 do while(date_comp(c_date,e_date) /=1 )
!---open data for the first day of month
    if ((c_date%d==1).and.(c_date%h==0)) then
       write(inputfile,'(a,i4.4,i2.2,i2.2,i2.2)')"ATM.",c_date%y,c_date%m,c_date%d,c_date%h
       write(*,*)"Open data file month ",inputfile
       if(iotyp.eq.1) then
          open(55,file=trim(input_dir)//"/"//inputfile,form='unformatted',recl=Ny*Nx*4,access='direct')
          lrec=0
       else if(iotyp.eq.2) then
          open(55,file=trim(input_dir)//"/"//inputfile,form='unformatted')
       endif
    endif

    write(*,*)c_date

!---Read data on sigma levels
         if(iotyp.eq.1) then
            do k=Nsig,1,-1
               lrec=lrec+1
               read(55,rec=lrec) ((u(i,j,k),i=1,Nx),j=1,Ny)
            enddo
            do k=Nsig,1,-1
               lrec=lrec+1
               read(55,rec=lrec) ((v(i,j,k),i=1,Nx),j=1,Ny)
            enddo
            do k=Nsig,1,-1
               lrec=lrec+1
               read(55,rec=lrec) ((w(i,j,k),i=1,Nx),j=1,Ny)
            enddo
            do k=Nsig,1,-1
               lrec=lrec+1
               read(55,rec=lrec) ((t(i,j,k),i=1,Nx),j=1,Ny)
            enddo
            do k=Nsig,1,-1
               lrec=lrec+1
               read(55,rec=lrec) ((q(i,j,k),i=1,Nx),j=1,Ny)
            enddo
            do k=Nsig,1,-1
               lrec=lrec+1
               read(55,rec=lrec) ((qc(i,j,k),i=1,Nx),j=1,Ny)
            enddo
            lrec=lrec+1
            read(55,rec=lrec) ((ps(i,j),i=1,Nx),j=1,Ny)
            lrec=lrec+1
            read(55,rec=lrec) ((rain(i,j),i=1,Nx),j=1,Ny)
            lrec=lrec+1
            read(55,rec=lrec) ((tgb(i,j),i=1,Nx),j=1,Ny)
            lrec=lrec+1
            read(55,rec=lrec) ((swt(i,j),i=1,Nx),j=1,Ny)
            lrec=lrec+1
            read(55,rec=lrec) ((rno(i,j),i=1,Nx),j=1,Ny)
         else if(iotyp.eq.2) then
            read(55) idate
            print*,' IDATE = ',idate
            do k=Nsig,1,-1
               read(55) ((u(i,j,k),i=1,Nx),j=1,Ny)
            enddo
            do k=Nsig,1,-1
               read(55) ((v(i,j,k),i=1,Nx),j=1,Ny)
            enddo
            do k=Nsig,1,-1
               read(55) ((w(i,j,k),i=1,Nx),j=1,Ny)
            enddo
            do k=Nsig,1,-1
               read(55) ((t(i,j,k),i=1,Nx),j=1,Ny)
            enddo
            do k=Nsig,1,-1
               read(55) ((q(i,j,k),i=1,Nx),j=1,Ny)
            enddo
            do k=Nsig,1,-1
               read(55) ((qc(i,j,k),i=1,Nx),j=1,Ny)
            enddo
            read(55) ((ps(i,j),i=1,Nx),j=1,Ny)
            read(55) ((rain(i,j),i=1,Nx),j=1,Ny)
            read(55) ((tgb(i,j),i=1,Nx),j=1,Ny)
            read(55) ((swt(i,j),i=1,Nx),j=1,Ny)
            read(55) ((rno(i,j),i=1,Nx),j=1,Ny)
         endif

!---Convert to pressure levels and mean sea level
    ! -- Calc. Height of sigma levels
    call HTSIG(T,h,PS,ht,sig,Nx,Ny,Nsig)
    ! -- Calc. Sea-Level Pressure using
    !      1. ERRICO's solution described in height
    !      2. a simple formulae
    !      3. MM5 method
      CALL SLPRES(H,T,PS,HT,TGB,SLP1,SLP2,SIG,Nx,Ny,Nsig)

    ! to interpolate H,U,V,W,T
    !     1. For Heights of Pressure surface
      CALL HEIGHT(HP,H,T,PS,HT,SIG,Nx,Ny,Nsig,PLEV,NP)
    !     2. For Zonal and Meridional Winds
      CALL INTLIN(UP,U,PS,SIG,Nx,Ny,Nsig,PLEV,NP)
      CALL INTLIN(VP,V,PS,SIG,Nx,Ny,Nsig,PLEV,NP)
      CALL INTLIN(WP,W,PS,SIG,Nx,Ny,Nsig,PLEV,NP)
    !     3. For Temperatures
      CALL INTLOG(TP,T,PS,SIG,Nx,Ny,Nsig,PLEV,NP)
    !     4. Potential temperature

    if ((c_date%d==days_of_month(c_date%y,c_date%m)).and.(c_date%h==18)) then
       write(*,*)"Close data file for month ",c_date%m
       close(55)
    endif

    !write data on plevels to grads files
    do k=1,Np
       write(fu,rec=orec1)up(:,:,k)
       write(fv,rec=orec1)vp(:,:,k)
       write(fw,rec=orec1)wp(:,:,k)
       write(ft,rec=orec1)tp(:,:,k)
       write(fh,rec=orec1)hp(:,:,k)
       orec1=orec1+1
    enddo

    write(fslp,rec=orec2)slp1
    orec2=orec2+1
    write(fslp,rec=orec2)slp2
    orec2=orec2+1

    call add_hrs(c_date,6)

 enddo

 close(11)


 Stop 
 end