MODULE process_domain_module
   
   USE input_module
   USE output_module
   USE module_model_basics
   USE module_arrays
   
   integer                             :: next_file
   logical                             :: run_out_of_files
   logical                             :: file_is_open
   integer                             :: vertical       !!! needed for v5dcreation
 

   CONTAINS

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: process_domain
   ! Purpose: Process the input data
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   SUBROUTINE process_domain ()
   
      USE date_pack
      USE gridinfo_module
      USE misc_definitions_module
      USE module_debug
      USE v5d_module
   
      IMPLICIT NONE

      ! ------ some vis5d stuff, contain MAXVARS,MAXTIMES,MAXROWS,MAXCOLUMNS,MAXLEVELS
      include 'v5df.h'
   
      ! Local variables
      integer                                :: idiff, n_times, i
      character (len=19)                     :: valid_date, temp_date, datestr
      character (len=128)                    :: cname, stagger, cunits, cdesc
      integer                                :: istatus, reclength
      integer                                :: ret                         !!! v5d stuff
      integer, dimension(MAXLEVELS)          :: v5d_nl
      real, dimension(MAXLEVELS)             :: v5d_levels
      integer, dimension(MAXTIMES)           :: timestamp, datestamp
      integer                                :: projection
      real, dimension(100)                   :: proj_args
   

      ! Compute number of times that we will process
      CALL geth_idts(end_date, start_date, idiff)
      CALL mprintf((idiff < 0),ERROR,'Ending date is earlier than starting date in namelist.')
   
   
      ! Check that the interval evenly divides the range of times to process
      n_times = idiff / interval_seconds
      CALL mprintf((mod(idiff, interval_seconds) /= 0),WARN, &
                   'In namelist, interval_seconds does not evenly divide '// &
                   '(end_date - start_date). Only %i time periods '// &
                   'will be processed.', i2=n_times)
   
      ! 
      ! DO TIME_INDEPENDANT PROCESSING
      ! 

      ! Initialize the input module to read input data
      CALL input_init(1, istatus)
      CALL mprintf((istatus /= 0),ERROR, 'input_init(): Error opening input.')
   
      ! Read global attributes from the input file 
      CALL read_global_attrs ()
      IF ( iprogram .le. 1 ) n_times = 0


      IF ( output_type == 'grads' ) THEN    !!! We are making grads type files

        ! Open .dat and .ctl files
        DO cunit=10,100
          INQUIRE(unit=cunit, opened=is_used)
          if (.not. is_used) exit
        END DO
        ctlfile = trim(output_root_name)//'.ctl'
        OPEN(cunit,file=ctlfile)

        DO dunit=10,100
          INQUIRE(unit=dunit, opened=is_used)
          IF (.not. is_used) EXIT
        END DO

        IF      ( west_east_dim .gt. 2 .AND. south_north_dim .gt. 2 ) THEN
          reclength = (west_east_dim)*(south_north_dim)
        ELSE IF ( west_east_dim .gt. 2 .AND. south_north_dim .eq. 2 ) THEN
          reclength = west_east_dim
        ELSE IF ( west_east_dim .eq. 2 .AND. south_north_dim .gt. 2 ) THEN
          reclength = south_north_dim
        ELSE IF ( west_east_dim .eq. 2 .AND. south_north_dim .eq. 2 ) THEN
          reclength = 1
        END IF
#ifdef RECL4
        reclength = reclength*4
#endif
        datfile = trim(output_root_name)//'.dat'
!        OPEN(dunit,file=datfile,form='unformatted',access="direct", recl=reclength)
        OPEN(dunit,file=datfile)

      ELSE IF ( output_type == 'v5d' ) THEN    !!! We are making grads type files

        v5dfile = trim(output_root_name)//'.v5d'
  
      END IF


      CALL arw_get_next_time(handle, datestr, istatus)

!!!
!!!
!!! MARS: we go to next time to init, because first step of PHTOT and PTOT are not filled
!!! MARS: caca caca caca. il vaut mieux le geop init pour calculer les hauteurs
!!!
!CALL arw_get_next_time(handle, datestr, istatus)



      vertical = 0      !!! needed for v5d creation (equally spaced levels in generic units)
      ! If interp_method /= 0, then we need to know if this can be done, for now only doing it for wrfout data
      IF ( iprogram .ge. 6 .AND. interp_method /= 0 ) THEN
         CALL get_interp_info (datestr)
      ELSE
        vertical_type = 'n'   !! In case user set this for other programs
        number_of_zlevs = bottom_top_dim
      END IF



#ifdef V5D
      IF ( output_type == 'v5d' ) THEN   !!! Create the vis5d file
        v5d_nl = IMISSING
        CALL v5d_varnames (v5d_nl)
        CALL v5d_times (n_times+1, timestamp, datestamp)
        CALL v5d_proj (datestr, projection, proj_args)
        v5d_levels = MISSING
        IF (vertical == 0 ) THEN
          v5d_levels(1) = 0.0   !!! Same basics set up for sigma level data
          v5d_levels(2) = 1.0
        ELSE
          DO i=1,number_of_zlevs
            v5d_levels(i) = interp_levels(i)
          END DO
        END IF
        ret = v5dCreate( v5dfile, n_times+1, numvars,                     &
                         !!west_east_dim, south_north_dim, v5d_nl,        &
                         south_north_dim, west_east_dim, v5d_nl,        &
                         varnames, timestamp, datestamp, 1, projection, & 
                         proj_args, vertical, v5d_levels )
        IF ( ret == 0 ) THEN 
          CALL mprintf(.true.,ERROR, 'v5dCreate')
        ELSE
          CALL mprintf(.true.,STDOUT, '   ')
          CALL mprintf(.true.,STDOUT, '   SUCCESS: v5d file has been Created')
        END IF

      END IF 
#endif

      CALL input_close()       !!! Will be opened again by next process if needed

      !
      ! BEGIN TIME-DEPENDANT PROCESSING
        next_file = 1
        ivars     = 0
      
      ! Loop over all times to be processed
      could_not_find = plot_these_fields      !!! Keep a list of requested fields
      run_out_of_files = .FALSE.
      file_is_open = .FALSE.
      rec = 0
      all_files : DO time=0,n_times
   
         IF ( iprogram == 1 ) then
           temp_date = "0000-00-00_00:00:00" 
         ELSE 
           CALL geth_newdate(valid_date, trim(start_date), time*interval_seconds)
           temp_date = ' '
           write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)
         ENDIF
         IF ( time == 0 ) tdef_date = temp_date
   
         CALL get_fields(temp_date)
   
         IF (run_out_of_files) EXIT all_files
   
      ENDDO all_files  ! Loop over n_times

      CALL mprintf(.true.,STDOUT, ' ')
      CALL mprintf(.true.,STDOUT, 'DONE Processing Data')

      IF ( output_type == 'grads' ) THEN    !!! We are making grads type files
        ! Create the .ctl file
        CALL mprintf(.true.,STDOUT, ' ')
        CALL mprintf(.true.,STDOUT, 'CREATING .ctl file')
        CALL create_ctl ()
      END IF  


#ifdef V5D
      IF ( output_type == 'v5d' ) THEN   !!! Create the vis5d file
        ret = v5dclose()
      END IF
#endif

   
   END SUBROUTINE process_domain


   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: get_interp_info
   ! Purpose: Get extra info if we need to interpolate to pressure / height    
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   SUBROUTINE get_interp_info (valid_date)

      IMPLICIT NONE

      ! Arguments
      character (len=19)                 :: valid_date
    
      ! Local variables
      integer                             :: istatus, wrftype
      real, pointer, dimension(:,:,:)     :: real_array
      character (len=3)                   :: memorder
      character (len=128)                 :: cname, stagger
      character (len=128), dimension(3)   :: dimnames
      integer, dimension(3)               :: domain_start, domain_end
      integer, dimension(2)               :: loc_of_min_z
      real, allocatable, dimension(:,:,:) :: ph_tmp, phb_tmp
      integer                             :: found


      found = 0

      wrftype = 104
      memorder = "XYZ"
      domain_start = 1
      domain_end(1) = west_east_dim
      domain_end(2) = south_north_dim

      IF (iprogram == 6 .AND. vertical_type == 'p') THEN
        domain_end(3) = bottom_top_dim
        stagger = '-'
        istatus = 0
        CALL read_spec_field(domain_start, domain_end, 'QVAPOR', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        found = found + istatus

        memorder = "XY"
        domain_end(3) = 1
        istatus = 0
        CALL read_spec_field(domain_start, domain_end, 'MU', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        found = found + istatus

        istatus = 0
        CALL read_spec_field(domain_start, domain_end, 'MUB', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        found = found + istatus

        memorder = "0"
        domain_end(1) = 1
        domain_end(2) = 1
        istatus = 0
        CALL read_spec_field(domain_start, domain_end, 'P_TOP', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        found = found + istatus

        memorder = "Z"
        domain_end(1) = bottom_top_dim
        istatus = 0
        CALL read_spec_field(domain_start, domain_end, 'ZNU', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        found = found + istatus

        domain_end(1) = bottom_top_dim + 1
        stagger = 'Z'
        istatus = 0
        CALL read_spec_field(domain_start, domain_end, 'ZNW', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        found = found + istatus

        IF ( found == 0 ) THEN
          vertical = 3    !!! (unequally spaced levels in mb - v5d)
          CALL mprintf(.true.,STDOUT, '   Interpolating to PRESSURE levels ')
          CALL mprintf(.true.,STDOUT, '   ')
        ELSE
          CALL mprintf(.true.,STDOUT, '   WARNING: Asked to interpolate to PRESSURE, but we do not have enough information ' )
          CALL mprintf(.true.,STDOUT, '            Will output data on MODEL LEVELS ' )
          vertical_type = 'n'
        ENDIF

      ENDIF

      IF (iprogram == 8 .AND. vertical_type == 'p') THEN
        domain_end(3) = bottom_top_dim
        stagger = '-'

        istatus = 0
        CALL read_spec_field(domain_start, domain_end, 'P', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        found = found + istatus

        istatus = 0
        CALL read_spec_field(domain_start, domain_end, 'PB', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        found = found + istatus

        IF ( found == 0 ) THEN
          vertical = 3    !!! (unequally spaced levels in mb - v5d)
          CALL mprintf(.true.,STDOUT, '   Interpolating to PRESSURE levels ')
          CALL mprintf(.true.,STDOUT, '   ')
        ELSE
          CALL mprintf(.true.,STDOUT, '   WARNING: Asked to interpolate to PRESSURE, but we do not have enough information ' )
          CALL mprintf(.true.,STDOUT, '            Will output data on MODEL LEVELS ' )
          vertical_type = 'n'
        ENDIF

      ENDIF


      IF (vertical_type == 'z' .AND. interp_method == 1 ) THEN

!!!
!!! MARS MARS: caca caca caca: corriger ce probleme de staggering ?
!!!
        domain_end(3) = bottom_top_dim + 1
        stagger = 'Z'

!         domain_end(3) = bottom_top_dim
!         stagger = ''
!!!
!!! MARS MARS: caca caca caca
!!!

        istatus = 0
        CALL read_spec_field(domain_start, domain_end, 'PH', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        found = found + istatus

        istatus = 0
        CALL read_spec_field(domain_start, domain_end, 'PHB', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        found = found + istatus

!!!print *, 'get_interp_info'
!
!        istatus = 0
!        CALL read_spec_field(domain_start, domain_end, 'PHTOT', wrftype, &
!                             memorder, stagger, dimnames, real_array, valid_date, istatus)
!        found = found + istatus

        IF ( found == 0 ) THEN
          vertical = 2    !!! (unequally spaced levels in km - v5d)
          CALL mprintf(.true.,STDOUT, '   Interpolating to USER SPECIFIED HEIGHT levels ')
          CALL mprintf(.true.,STDOUT, '   ')
        ELSE
          CALL mprintf(.true.,STDOUT, '   WARNING: Asked to interpolate&
& to USER SPECIFIED HEIGHT, but we do not have enough information ' )
          CALL mprintf(.true.,STDOUT, '            Will output data on MODEL LEVELS ' )
          vertical_type = 'n'
        ENDIF

      ENDIF

      IF (vertical_type == 'z' .AND. interp_method == -1 ) THEN
        domain_end(3) = bottom_top_dim + 1
        stagger = 'Z'

        istatus = 0
        CALL read_spec_field(domain_start, domain_end, 'PH', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        found = found + istatus
        IF ( istatus == 0 ) THEN
          allocate(ph_tmp(domain_end(1),domain_end(2),domain_end(3)-1))
          ph_tmp = 0.5*(real_array(:,:,1:domain_end(3)-1)+real_array(:,:,2:domain_end(3)))
        ENDIF

        istatus = 0
        CALL read_spec_field(domain_start, domain_end, 'PHB', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        found = found + istatus
        IF ( istatus == 0 ) THEN
          allocate(phb_tmp(domain_end(1),domain_end(2),domain_end(3)-1))
          phb_tmp = 0.5*(real_array(:,:,1:domain_end(3)-1)+real_array(:,:,2:domain_end(3)))
        ENDIF

        IF ( found == 0 ) THEN
          vertical = 2    !!! (unequally spaced levels in km - v5d)
!!          ph_tmp = ( (ph_tmp+phb_tmp) / 9.81)/1000.        !!! convert to height in km

          ph_tmp = ( (ph_tmp+phb_tmp) / 3.72)/1000.        !!! convert to height in km

          !!! Generate nice heights. 
          !!! Get location of lowest height on model grid
          !!! Use column above this point as our heights to interpolate to
          !!! Adjust lowerst and heighest heights a bit
          number_of_zlevs = bottom_top_dim
          loc_of_min_z = minloc(ph_tmp(:,:,1))
          interp_levels(1:number_of_zlevs) =                                   &
                    ph_tmp(loc_of_min_z(1),loc_of_min_z(2),1:number_of_zlevs)

	  !!--- .002 problematic on Mars with low pressure
          !!interp_levels(1) = interp_levels(1) + 0.002
          !!interp_levels(1) = MAX(interp_levels(1), interp_levels(2)/2.0)  !! no neg values
          !!interp_levels(number_of_zlevs) = interp_levels(number_of_zlevs) - 0.002
          vertical = 2
          CALL mprintf(.true.,STDOUT, '   Interpolating to GENERATED HEIGHT levels ')
          CALL mprintf(.true.,STDOUT, '   ')
        ELSE
          CALL mprintf(.true.,STDOUT, '   WARNING: Asked to interpolate to GENERATED &
&HEIGHT, but we do not have enough information ' )
          CALL mprintf(.true.,STDOUT, '            Will output data on MODEL LEVELS ' )
          vertical_type = 'n'
        ENDIF

        IF (ALLOCATED(ph_tmp))  DEALLOCATE(ph_tmp)
        IF (ALLOCATED(phb_tmp)) DEALLOCATE(phb_tmp)
      ENDIF
    

   END SUBROUTINE get_interp_info

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: get_interp_array
   ! Purpose: Get array we will use when interpolting 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   SUBROUTINE get_interp_array (valid_date)

      USE module_pressure

      IMPLICIT NONE

      ! Arguments
      character (len=19)                 :: valid_date
    
      ! Local variables
      integer                             :: istatus, wrftype
      real, pointer, dimension(:,:,:)     :: real_array
      character (len=3)                   :: memorder
      character (len=128)                 :: cname, stagger
      character (len=128), dimension(3)   :: dimnames
      integer, dimension(3)               :: domain_start, domain_end
      integer, dimension(2)               :: loc_of_min_z
      real, allocatable, dimension(:,:,:) :: p_tmp, pb_tmp, ph_tmp, phb_tmp
      real, pointer, dimension(:,:,:)     :: SCR3


      wrftype = 104
      memorder = "XYZ"
      domain_start = 1
      domain_end(1) = west_east_dim
      domain_end(2) = south_north_dim

      IF (iprogram == 6 .AND. vertical_type == 'p') THEN
        domain_end(3) = bottom_top_dim
        stagger = '-'
        CALL read_spec_field(domain_start, domain_end, 'QVAPOR', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        allocate(QV(west_east_dim,south_north_dim,bottom_top_dim))
        QV = real_array

        memorder = "XY"
        domain_end(3) = 1
        stagger = '-'
        CALL read_spec_field(domain_start, domain_end, 'MU', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        allocate(MU(west_east_dim,south_north_dim))
        MU = real_array(:,:,1)
        CALL read_spec_field(domain_start, domain_end, 'MUB', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        allocate(MUB(west_east_dim,south_north_dim))
        MUB = real_array(:,:,1)

        memorder = "0"
        domain_end(1) = 1
        domain_end(2) = 1
        CALL read_spec_field(domain_start, domain_end, 'P_TOP', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        PTOP = real_array(1,1,1)

        memorder = "Z"
        domain_end(1) = bottom_top_dim
        CALL read_spec_field(domain_start, domain_end, 'ZNU', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        allocate(ZNU(bottom_top_dim))
        ZNU = real_array(:,1,1)

        domain_end(1) = bottom_top_dim + 1
        stagger = 'Z'
        CALL read_spec_field(domain_start, domain_end, 'ZNW', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        allocate(ZNW(bottom_top_dim+1))
        ZNW = real_array(:,1,1)

       
        CALL pressure(SCR3)
        ALLOCATE(vert_array(west_east_dim,south_north_dim,bottom_top_dim))
        vert_array = 0.01 * SCR3      !! pressure array in HPa

        IF (ALLOCATED(QV))  DEALLOCATE(QV)
        IF (ALLOCATED(MU))  DEALLOCATE(MU)
        IF (ALLOCATED(MUB)) DEALLOCATE(MUB)
        IF (ALLOCATED(ZNU)) DEALLOCATE(ZNU)
        IF (ALLOCATED(ZNW)) DEALLOCATE(ZNW)
      ENDIF



      IF (iprogram == 8 .AND. vertical_type == 'p') THEN
        domain_end(3) = bottom_top_dim
        stagger = '-'
        CALL read_spec_field(domain_start, domain_end, 'P', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        allocate(p_tmp(west_east_dim,south_north_dim,bottom_top_dim))
        p_tmp = real_array
        CALL read_spec_field(domain_start, domain_end, 'PB', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        allocate(pb_tmp(west_east_dim,south_north_dim,bottom_top_dim))
        pb_tmp = real_array

        IF (ALLOCATED(vert_array)) DEALLOCATE(vert_array)
        ALLOCATE(vert_array(west_east_dim,south_north_dim,bottom_top_dim))
        vert_array = (p_tmp+pb_tmp)/100.      !! pressure array in HPa

        IF (ALLOCATED(p_tmp))   DEALLOCATE(p_tmp)
        IF (ALLOCATED(pb_tmp))  DEALLOCATE(pb_tmp)
      ENDIF



      IF (vertical_type == 'z' ) THEN

!!!
!!! MARS MARS: caca caca caca: corriger ce probleme de staggering ?
!!!
        domain_end(3) = bottom_top_dim + 1
        stagger = 'Z'

!         domain_end(3) = bottom_top_dim
!         stagger = ''
!!!
!!! MARS MARS: caca caca caca
!!!

        CALL read_spec_field(domain_start, domain_end, 'PH', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        allocate(ph_tmp(west_east_dim,south_north_dim,bottom_top_dim))
        ph_tmp = 0.5*(real_array(:,:,1:domain_end(3)-1)+real_array(:,:,2:domain_end(3)))
        CALL read_spec_field(domain_start, domain_end, 'PHB', wrftype, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
        allocate(phb_tmp(west_east_dim,south_north_dim,bottom_top_dim))
        phb_tmp = 0.5*(real_array(:,:,1:domain_end(3)-1)+real_array(:,:,2:domain_end(3)))


!        CALL read_spec_field(domain_start, domain_end, 'PHTOT', wrftype, &
!                             memorder, stagger, dimnames, real_array, valid_date, istatus)
!
!!print *, 'get_interp_array'
!!!!
!!!! MARS MARS: caca caca caca
!!!!
!        allocate(ph_tmp(west_east_dim,south_north_dim,bottom_top_dim))
!        !!ph_tmp = 0.5*0.5*(real_array(:,:,1:domain_end(3)-1)+real_array(:,:,2:domain_end(3)))
!ph_tmp = 0.5*real_array
!
!        allocate(phb_tmp(west_east_dim,south_north_dim,bottom_top_dim))
!        !!phb_tmp = 0.5*0.5*(real_array(:,:,1:domain_end(3)-1)+real_array(:,:,2:domain_end(3)))
!phb_tmp = 0.5*real_array
!!!!
!!!! MARS MARS: caca caca caca
!!!!

        IF (ALLOCATED(vert_array)) DEALLOCATE(vert_array)
        ALLOCATE(vert_array(west_east_dim,south_north_dim,bottom_top_dim))
!!        vert_array = ( (ph_tmp+phb_tmp) / 9.81)/1000.        !! height array in km
        vert_array = ( (ph_tmp+phb_tmp) / 3.72)/1000.        !! height array in km

        IF (ALLOCATED(ph_tmp))  DEALLOCATE(ph_tmp)
        IF (ALLOCATED(phb_tmp)) DEALLOCATE(phb_tmp)
      ENDIF


   END SUBROUTINE get_interp_array

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: get_fields
   ! Purpose: Read all fields in input file and process required output fields
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   SUBROUTINE get_fields (valid_date)

      USE module_get_file_names
      USE module_diagnostics
      USE module_interp

      IMPLICIT NONE

      ! Local variables
      integer                            :: istatus, i, j, k, ii, jj, kk
      real, pointer, dimension(:,:,:)    :: real_array
      character (len=3)                  :: memorder
      character (len=19)                 :: datestr, valid_date
      character (len=128)                :: cname, stagger, cunits, cdesc
      character (len=128), dimension(3)  :: dimnames
      integer, dimension(3)              :: domain_start, domain_end
      real, pointer, dimension(:,:,:)    :: data_out 
      real, allocatable, dimension(:,:,:):: SCR
      integer                            :: nxout, nyout, nzout
      character (len=20)                 :: dummy
      integer                            :: is_there


      ! Initialize the input module to read static input data for this domain
      CALL mprintf(.true.,STDOUT, ' ')
      CALL mprintf(.true.,STDOUT, ' Processing  time --- %s', s1=trim(valid_date))
      get_right_time : DO
        IF ( next_file > number_of_input_files ) THEN
           CALL mprintf(.true.,STDOUT, '   --- WE HAVE RUN OUT OF INPUT FILES ---')
           run_out_of_files = .TRUE.
           RETURN
        END IF
        IF ( .not. file_is_open) THEN
          CALL input_init(next_file, istatus)
          CALL mprintf((istatus /= 0),ERROR, 'input_init(): Error opening input.')
          file_is_open = .TRUE.
        END IF
        CALL arw_get_next_time(handle, datestr, istatus)
        IF ( istatus /= 0 ) THEN     ! might be in a next file
           CALL input_close()
           file_is_open = .FALSE.
           CALL mprintf(.true.,STDOUT, '   Date not in this file - see if there are more files ')
           next_file = next_file + 1
           CYCLE get_right_time
        END IF
        IF ( TRIM(datestr) .LT. TRIM(valid_date) ) THEN
           CALL mprintf(.true.,STDOUT, '   Not looking for %s - skipping to next time' , s1=datestr )
           CYCLE get_right_time
        ELSEIF ( TRIM(datestr) .EQ. TRIM(valid_date) ) THEN
           CALL mprintf(.true.,STDOUT, '   Found the right date - continue ' )
           EXIT get_right_time
        ELSEIF ( TRIM(datestr) .GT. TRIM(valid_date) ) THEN
           CALL mprintf(.true.,STDOUT, '   Found  %s  before  %s', s1=trim(datestr),s2=trim(valid_date))
           run_out_of_files = .TRUE.
           RETURN
        ENDIF
      ENDDO get_right_time

      CALL clobber_arrays
#ifdef IO_GRIB1
      IF (io_form_input == GRIB1) THEN
        REWIND (13)
      END IF
#endif
   

      !! IF we are interpolting we need the pressure/height array to interpolate to
      IF ( vertical_type /= 'n' ) THEN
         CALL get_interp_array (valid_date)
      ENDIF

   
      ! Read fields using the input module; we know that there are no more
      !   fields to be read when read_next_field() returns a non-zero status.
      istatus = 0
      process_all_fields : DO WHILE (istatus == 0)  
        CALL read_next_field(domain_start, domain_end, cname, cunits, cdesc, &
                             memorder, stagger, dimnames, real_array, valid_date, istatus)
#ifdef IO_GRIB1
        IF (io_form_input == GRIB1) THEN
          IF (istatus == -5) THEN
            istatus = 0
            CYCLE process_all_fields
          END IF
        END IF
#endif
        IF (istatus == 0) THEN
          IF (cname(len_trim(cname)-1:len_trim(cname)) == "_U" ) CYCLE process_all_fields
          IF (cname(len_trim(cname)-1:len_trim(cname)) == "_V" ) CYCLE process_all_fields

            CALL keep_arrays(cname, real_array)  !!! Keep some fields around

!!
!! DECIDE IF WE WANT THE FIELD AND IF YES, INTERPOLATE TO CORRECT GRID IF NEEEDED
!! assume we got the grid info up at the top somewhere
!!

            dummy = ","//trim(cname)//","
            is_there = INDEX(plot_these_fields,trim(dummy))
            IF ( ( INDEX(plot,'all') /= 0 .OR. is_there /= 0) .AND. domain_end(2) /= 1) THEN
              IF (ALLOCATED(SCR)) DEALLOCATE(SCR)
              ALLOCATE(SCR(domain_end(1),domain_end(2),domain_end(3)))
              SCR = real_array
              CALL interp( SCR, domain_end(1), domain_end(2), domain_end(3), &
                           data_out, nxout, nyout, nzout, &
                           vert_array, interp_levels, number_of_zlevs) 
                           
            ! Write the fields we want out to the .dat file, also keeps a list of what is written out
              CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
            ENDIF

        ENDIF    !! end "istatus==0"

      END DO process_all_fields



      IF (  INDEX(plot,'list') /= 0 .OR. INDEX(plot,'file') /= 0 ) THEN
!! Do we have any DIAGNOSTICS to process?
        CALL process_diagnostics ()
      END IF


!! We are DONE for this time




   !! print a list of the requested fields we could not find
   IF ( len_trim(could_not_find) > 1 ) THEN
     CALL mprintf(.true.,STDOUT, '   ')
     CALL mprintf(.true.,STDOUT, '   WARNING: The following requested fields could not be found ' )

     could_not_find = could_not_find(2:len_trim(could_not_find))
     is_there = INDEX(could_not_find,",")
#ifdef V5D
     IF ( output_type == 'v5d' .AND. is_there > 1 ) THEN   !!! Warning for possible bad data file
       CALL mprintf(.true.,STDOUT, '   WARNING: This will create a problem for VIS5D, please ' )
       CALL mprintf(.true.,STDOUT, '            remove these fields from your input list and ' )
       CALL mprintf(.true.,STDOUT, '            run the code again.                          ' )
       CALL mprintf(.true.,STDOUT, '                                                         ' )
     END IF
#endif
     DO WHILE ( is_there > 1 ) 
       PRINT*,"           ", could_not_find(1:is_there-1)
       could_not_find = could_not_find(is_there+1:len_trim(could_not_find))
       is_there = INDEX(could_not_find,",")
     END DO
   
   CALL mprintf(.true.,STDOUT, '   ')
   ENDIF
    

   END SUBROUTINE get_fields

END MODULE process_domain_module
