module coefficients

   integer :: n_levels
   real, allocatable, dimension(:) :: a, b

   contains

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: read_coeffs
   !
   ! Notes: Obtain table of coefficients for input by this routine from
   !        http://www.ecmwf.int/products/data/technical/model_levels/index.html
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine read_coeffs()
      
      implicit none

      integer :: i, nlvl, istatus

      open(21,file='ecmwf_coeffs',form='formatted',status='old',iostat=istatus)
  
      n_levels = 0

      if (istatus /= 0) then
         write(6,*) 'ERROR: Error opening ecmwf_coeffs' 
         return
      end if

      read(21,*,iostat=istatus) nlvl
      do while (istatus == 0)
         n_levels = n_levels + 1
         read(21,*,iostat=istatus) nlvl
      end do
      
      rewind(21)

      n_levels = n_levels - 1

      allocate(a(0:n_levels))
      allocate(b(0:n_levels))

      write(6,*) ' '
      write(6,*) 'Coefficients for each level:',n_levels
      do i=0,n_levels
         read(21,*,iostat=istatus) nlvl, a(i), b(i)
         write(6,'(i5,5x,f12.6,2x,f12.10)') nlvl, a(i), b(i)
      end do
      write(6,*) ' '

      close(21)
      
   end subroutine read_coeffs


   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: cleanup_coeffs
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine cleanup_coeffs()

      implicit none

      n_levels = 0
      if (allocated(a)) deallocate(a)
      if (allocated(b)) deallocate(b)

   end subroutine cleanup_coeffs

end module coefficients


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! calc_ecmwf_p
!
! The purpose of this program is to compute a 3d pressure field for ECMWF 
!    model-level data sets; the code works in the WPS intermediate file format, 
!    reading a PSFC field from intermediate files, the A and B coefficients 
!    from a text file, ecmwf_coeffs, and writes the pressure data to an 
!    intermediate file.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
program calc_ecmwf_p

   use coefficients
   use date_pack
   use gridinfo_module
   use misc_definitions_module
   use module_debug
   use read_met_module
   use stringutil
   use write_met_module

   implicit none

   ! Local variables
   integer :: i, idiff, n_times, t, istatus, fg_idx
   real :: a_full, b_full
   character (len=19) :: valid_date, temp_date
   character (len=128) :: input_name
   real, allocatable, dimension(:,:) :: psfc
   real, allocatable, dimension(:,:,:) :: tt, qv
   logical :: is_psfc
   type (met_data) :: ecmwf_data, p_data, rh_data


   !
   ! Setup (read namelist and check on time range)
   !
   call get_namelist_params()

   call set_debug_level(WARN)

   ! Compute number of times that we will process
   call geth_idts(end_date(1), start_date(1), idiff)
   call mprintf((idiff < 0),ERROR,'Ending date is earlier than starting date in namelist for domain %i.', i1=1)

   n_times = idiff / interval_seconds

   ! Check that the interval evenly divides the range of times to process
   call mprintf((mod(idiff, interval_seconds) /= 0),WARN, &
                'In namelist, interval_seconds does not evenly divide '// &
                '(end_date - start_date) for domain %i. Only %i time periods '// &
                'will be processed.', i1=1, i2=n_times)

   fg_idx = 1

   input_name = fg_name(fg_idx)

   !
   ! Get coefficients for model level pressures
   !
   call read_coeffs()


   !
   ! Loop over all prefixes listed in namelist for fg_name
   !
   do while (input_name /= '*')

      !
      ! Loop over all times to be processed for this domain
      !
      do t=0,n_times

         ! Get the date string for the current time
         call geth_newdate(valid_date, trim(start_date(1)), t*interval_seconds)
         temp_date = ' '
         write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)

         ! Initialize the module for reading in the met fields
         call read_met_init(trim(input_name), .false., temp_date(1:13), istatus)

         is_psfc = .false.

         if (istatus == 0) then
            call mprintf(.true.,STDOUT,'Reading from %s at time %s', s1=input_name, s2=temp_date(1:13))

            ! Process all fields and levels from the current file; read_next_met_field()
            !   will return a non-zero status when there are no more fields to be read.
            do while (istatus == 0)

               call read_next_met_field(ecmwf_data, istatus)

               if (istatus == 0) then

                  ! Have we found either PSFC or LOGSFP?
                  if ((trim(ecmwf_data%field) == 'PSFC' .and. ecmwf_data%xlvl == 200100.) &
                      .or. trim(ecmwf_data%field) == 'LOGSFP') then
                     p_data = ecmwf_data
                     p_data%field = 'PRESSURE '
                     p_data%desc  = 'Pressure'

                     rh_data = ecmwf_data
                     rh_data%field = 'RH       '
                     rh_data%units = '%'
                     rh_data%desc  = 'Relative humidity'

                     if (.not. allocated(psfc)) then
                        allocate(psfc(ecmwf_data%nx,ecmwf_data%ny))
                     end if
                     call mprintf(.true.,STDOUT,'Found %s field in %s:%s', &
                                  s1=ecmwf_data%field, s2=input_name, s3=temp_date(1:13))

                     is_psfc   = .true.
                     if (trim(ecmwf_data%field) == 'LOGSFP') then
                        psfc(:,:) = exp(ecmwf_data%slab(:,:))
                     else
                        psfc(:,:) = ecmwf_data%slab(:,:)
                     end if

                  ! Have we found temperature?
                  else if (trim(ecmwf_data%field) == 'TT') then

                     if (.not. allocated(tt)) then
                        allocate(tt(ecmwf_data%nx,ecmwf_data%ny,n_levels+1))  ! Extra level is for surface
                     end if

                     if (nint(ecmwf_data%xlvl) >= 1 .and. &
                         nint(ecmwf_data%xlvl) <= n_levels) then 
                        tt(:,:,nint(ecmwf_data%xlvl)) = ecmwf_data%slab
                     else if (nint(ecmwf_data%xlvl) == 200100) then
                        tt(:,:,n_levels+1) = ecmwf_data%slab
                     end if

                  ! Have we found specific humidity?
                  else if (trim(ecmwf_data%field) == 'SPECHUMD') then

                     if (.not. allocated(qv)) then
                        allocate(qv(ecmwf_data%nx,ecmwf_data%ny,n_levels+1))  ! Extra level is for surface
                     end if

                     if (nint(ecmwf_data%xlvl) >= 1 .and. &
                         nint(ecmwf_data%xlvl) <= n_levels) then 
                        qv(:,:,nint(ecmwf_data%xlvl)) = ecmwf_data%slab
                     else if (nint(ecmwf_data%xlvl) == 200100) then
                        qv(:,:,n_levels+1) = ecmwf_data%slab
                     end if

                  end if
                  
                  if (associated(ecmwf_data%slab)) deallocate(ecmwf_data%slab)
   
               end if
   
            end do
   
            call read_met_close()


            ! Now write out, for each level, the pressure field
            if (is_psfc) then

               allocate(p_data%slab(p_data%nx,p_data%ny))
               allocate(rh_data%slab(rh_data%nx,rh_data%ny))

               call write_met_init(trim(get_path(input_name))//'PRES', .false., temp_date(1:13), istatus)

               if (allocated(tt) .and. allocated(qv)) then
                  p_data%xlvl  = 200100.
                  p_data%slab  = psfc
! Surface RH should be computed from surface DEWPT by ungrib
!                  rh_data%xlvl = 200100.
!                  call calc_rh(tt(:,:,n_levels+1), qv(:,:,n_levels+1), psfc, rh_data%slab, rh_data%nx, rh_data%ny)
!                  call write_next_met_field(rh_data, istatus) 
                  call write_next_met_field(p_data, istatus) 
               else
                  call mprintf(.true.,WARN,'Either TT or SPECHUMD not found. No RH will be computed.')
               end if

               do i = 1, n_levels

                  a_full = 0.5 * (a(i-1) + a(i))   ! A and B are dimensioned (0:n_levels)
                  b_full = 0.5 * (b(i-1) + b(i))

                  p_data%xlvl = real(i)
                  p_data%slab = a_full + psfc * b_full

                  if (allocated(tt) .and. allocated(qv)) then
                     rh_data%xlvl = real(i)
                     call calc_rh(tt(:,:,i), qv(:,:,i), p_data%slab, rh_data%slab, rh_data%nx, rh_data%ny)
                     call write_next_met_field(rh_data, istatus) 
                  end if

                  call write_next_met_field(p_data, istatus) 

               end do

               call write_met_close()
  
               deallocate(p_data%slab)
               deallocate(rh_data%slab)

            end if

            if (allocated(tt))   deallocate(tt)
            if (allocated(qv))   deallocate(qv)
            if (allocated(psfc)) deallocate(psfc)

         else
            call mprintf(.true.,ERROR,'Problem opening %s at time %s', s1=input_name, s2=temp_date(1:13))
         end if

      end do 

      fg_idx = fg_idx + 1
      input_name = fg_name(fg_idx)

   end do 

   call cleanup_coeffs()

   stop

end program calc_ecmwf_p


subroutine calc_rh(t, qv, p, rh, nx, ny)

   implicit none

   ! Arguments
   integer, intent(in) :: nx, ny
   real, dimension(nx, ny), intent(in)  :: t, qv, p
   real, dimension(nx, ny), intent(out) :: rh

   ! Constants
   real, parameter :: svp1=0.6112
   real, parameter :: svp2=17.67
   real, parameter :: svp3=29.65
   real, parameter :: svpt0=273.15
   real, parameter :: eps = 0.622

   ! Local variables
   integer :: i, j
   real :: es, e

   do j=1,ny
   do i=1,nx
      es=svp1*10.*exp(svp2*(t(i,j)-svpt0)/(t(i,j)-svp3))
      e=qv(i,j)*p(i,j)/100./(eps+qv(i,j)*(1.-eps)) ! qv is specific humidity
      rh(i,j) = 100.0 * e/es
   end do
   end do

end subroutine calc_rh