!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

 module io_binary

!BOP
! !MODULE: io_binary
! !DESCRIPTION:
! This module provides a binary parallel input/output interface
! for writing fields.
!
! !REVISION HISTORY:
! CVS:$Id: io_binary.F90,v 1.18 2002/12/20 10:21:09 pwjones Exp $
! CVS:$Name: POP_2_0_1 $

! !USES:

   use kinds_mod
   use domain_size
   use domain
   use constants
   use boundary
   use communicate
   use broadcast
   use gather_scatter
   use exit_mod
   use io_types

   implicit none
   private
   save

! !PUBLIC MEMBER FUNCTIONS:

   public :: open_read_binary, &
             open_binary, &
             close_binary, &
             define_field_binary, &
             read_field_binary, &
             write_field_binary

!EOP
!BOC
!-----------------------------------------------------------------------
!
! module variables
!
!-----------------------------------------------------------------------

   character (1), parameter :: &
      attrib_separator=':' ! character separating attrib name and value

!-----------------------------------------------------------------------
!
! binary interfaces for i/o routines
!
!-----------------------------------------------------------------------

   interface write_array
      module procedure write_int_2d, &
                       write_real4_2d, &
                       write_real8_2d, &
                       write_int_3d, &
                       write_real4_3d, &
                       write_real8_3d
   end interface

   interface read_array
     module procedure read_int_2d, &
                      read_real4_2d, &
                      read_real8_2d, &
                      read_int_3d, &
                      read_real4_3d, &
                      read_real8_3d
   end interface

!EOC
!***********************************************************************

 contains

!***********************************************************************
!BOP
! !IROUTINE: open_read_binary
! !INTERFACE:

 subroutine open_read_binary(data_file)

! !DESCRIPTION:
! This routine opens a binary data file (and header file if it
! exists) for reading. It also reads global file attributes.
!
! !REVISION HISTORY:
! same as module

! !INPUT/OUTPUT PARAMETERS:

   type (datafile), intent (inout) :: &
      data_file ! file to be opened

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   character (80) :: &
      path, &! filename to read
      work_line, &! temporary to use for parsing file lines
      comp_line, &! temporary to use for parsing file lines
      data_type, &! temporary to use for data type
      att_name ! attribute name

   integer (i4) :: &
      n, &! loop index
      att_ival, &! value for integer attribute
      hdr_error, &! error for reading header file
      cindx1,cindx2 ! indices for unpacking character strings

   logical (log_kind) :: &
      att_lval ! value for logical attribute

   real (r4) :: &
      att_rval ! value for real attribute

   real (r8) :: &
      att_dval ! value for double attribute

   logical (log_kind) :: &
      header_exists, &! flag to check existence of header file
      lmatch ! flag to use for attribute search

!-----------------------------------------------------------------------
!
! set the readonly flag in the data file descriptor
!
!-----------------------------------------------------------------------

   data_file%readonly = .true.

!-----------------------------------------------------------------------
!
! if unit not assigned already, assign a unit (id) to this file
!
!-----------------------------------------------------------------------

   if (data_file%id(1) == 0) then
     call get_unit(data_file%id(1))
     call get_unit(data_file%id(2))
   endif

!-----------------------------------------------------------------------
!
! check to see whether a header file exists
! if so, open it
! if not, set the unit to -unit so it can be released later
!
!-----------------------------------------------------------------------

   header_exists = .false.

   if (my_task == master_task) then
      path = trim(data_file%full_name)/&
             &/'.hdr'

      inquire(file=path, exist=header_exists)

      if (header_exists) then
         open (unit=data_file%id(2), file=path, status='old')
      else
         write(stdout,*) 'WARNING: Input header file does not exist'
         write(stdout,*) 'for file: ',trim(path)
         write(stdout,*) &
                  'Assuming fields will be read in sequential order'
      endif
   endif

   call broadcast_scalar(header_exists, master_task)

   if (.not. header_exists) then
      data_file%id(2) = -data_file%id(2)
      data_file%current_record = 1
   endif

!-----------------------------------------------------------------------
!
! now open data file
!
!-----------------------------------------------------------------------

   if (my_task < data_file%num_iotasks) then
      path = trim(data_file%full_name)
      open (unit=data_file%id(1),action='read',status='unknown', &
            file=trim(path), form='unformatted',access='direct', &
            recl=data_file%record_length)
   endif

!-----------------------------------------------------------------------
!
! if the header file exists read file attributes from header file
!
!-----------------------------------------------------------------------

   if (header_exists) then

      work_line = char_blank
      hdr_error = 0

      if (my_task == master_task) then
         !*** first read until hit beginning of global attribute
         !*** definition section

         do while (hdr_error == 0 .and. work_line(1:7) /= '&GLOBAL')
            read(data_file%id(2),'(a80)',iostat=hdr_error) comp_line
            work_line = adjustl(comp_line)
         end do
      end if

      !*** now read until hit end of global attribute definition

      att_loop: do while (hdr_error == 0)

         work_line = char_blank
         if (my_task == master_task) then
            read(data_file%id(2),'(a80)',iostat=hdr_error) comp_line
            work_line = adjustl(comp_line)
         endif

         call broadcast_scalar(work_line, master_task)
         call broadcast_scalar(hdr_error, master_task)

         if (work_line(1:1) == '/') exit att_loop

         !*** find location of separator in the string
         !*** to determine name of attribute and extract
         !*** attribute name from beginning of string

         cindx1 = index(work_line,attrib_separator)
         att_name = char_blank
         att_name(1:cindx1-1) = work_line(1:cindx1-1)

         comp_line = work_line
         do n=1,cindx1
            comp_line(n:n) = ' '
         end do
         work_line = adjustl(comp_line)

         !*** now find location of separator in the string
         !*** to determine data type of attribute and extract
         !*** data type from beginning of string

         cindx1 = index(work_line,attrib_separator)
         data_type = char_blank
         data_type(1:cindx1-1) = work_line(1:cindx1-1)

         comp_line = work_line
         do n=1,cindx1
            comp_line(n:n) = ' '
         end do
         work_line = adjustl(comp_line)

         !*** check for standard file attributes

         select case(trim(att_name))
         case ('title','TITLE')
            data_file%title = trim(work_line)
         case ('history','HISTORY')
            data_file%history = trim(work_line)
         case ('conventions','CONVENTIONS')
            data_file%conventions = trim(work_line)
         case default

            !*** if additional attributes exist, add them as
            !*** additional attributes
            !*** note that if they are already defined, the
            !*** call to add_attrib will overwrite value

            select case (trim(data_type))
            case('char','CHAR')
               call add_attrib_file(data_file, trim(att_name), &
                                               trim(work_line))
            case('log','LOG','logical','LOGICAL')
               read(work_line,*) att_lval
               call add_attrib_file(data_file, trim(att_name), att_lval)
            case('int','INT','i4','I4')
               read(work_line,*) att_ival
               call add_attrib_file(data_file, trim(att_name), att_ival)
            case('r4','R4','REAL','real','float','FLOAT')
               read(work_line,*) att_rval
               call add_attrib_file(data_file, trim(att_name), att_rval)
            case('r8','R8','dbl','DBL','double','DOUBLE')
               read(work_line,*) att_dval
               call add_attrib_file(data_file, trim(att_name), att_dval)
            end select

         end select
      end do att_loop

   endif ! header exists

!-----------------------------------------------------------------------
!EOC

 end subroutine open_read_binary

!***********************************************************************
!BOP
! !IROUTINE: open_binary
! !INTERFACE:

 subroutine open_binary(data_file)

! !DESCRIPTION:
! This routine opens a binary data file and header file for writing
! and writes global file attributes to the header file.
!
! !REVISION HISTORY:
! same as module

! !INPUT/OUTPUT PARAMETERS:

   type (datafile), intent (inout) :: &
      data_file ! file to open

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   character (255) :: path

   character (char_len) :: &
      work_line ! temp space for manipulating strings

   character (5), parameter :: &
      hdr_fmt = '(a80)'

   integer (i4) :: &
      n, &! loop index
      cindx1, cindx2, &! indices for character strings
      hdr_unit ! unit number for header file

!-----------------------------------------------------------------------
!
! assign units to file if not assigned already
!
!-----------------------------------------------------------------------

   if (data_file%id(1) == 0) then
      call get_unit(data_file%id(1))
      call get_unit(data_file%id(2))
   endif

!-----------------------------------------------------------------------
!
! open data and header files for writing
!
!-----------------------------------------------------------------------

   !*** open header file from one task only

   if (my_task == master_task) then
      path = trim(data_file%full_name)/&
                                       &/'.hdr'
      open (unit=data_file%id(2), file=path, status='unknown')
      hdr_unit = data_file%id(2)
   end if

   path = trim(data_file%full_name)

   !*** open data file from all io tasks

   if (my_task < data_file%num_iotasks) then
      open(data_file%id(1), file=path, access='direct', &
           form='unformatted', recl=data_file%record_length, &
           status='unknown')
   endif

   !*** initialize record number

   data_file%current_record = 1

!-----------------------------------------------------------------------
!
! write attributes to header file
!
!-----------------------------------------------------------------------

   if (my_task == master_task) then

      work_line = char_blank
      work_line(1:7) = '&GLOBAL'
      write(hdr_unit,hdr_fmt) work_line

      work_line = char_blank
      work_line(1:5) = 'title'
      work_line(6:6) = attrib_separator
      work_line(7:10) = 'char'
      work_line(11:11) = attrib_separator
      write(work_line(12:),*) trim(data_file%title)
      write(hdr_unit,hdr_fmt) work_line

      work_line = char_blank
      work_line(1:7) = 'history'
      work_line(8:8) = attrib_separator
      work_line(9:12) = 'char'
      work_line(13:13) = attrib_separator
      write(work_line(14:),*) trim(data_file%history)
      write(hdr_unit,hdr_fmt) work_line

      work_line = char_blank
      work_line(1:11) = 'conventions'
      work_line(12:12) = attrib_separator
      work_line(11:14) = 'char'
      work_line(15:15) = attrib_separator
      write(work_line(16:),*) trim(data_file%conventions)
      write(hdr_unit,hdr_fmt) work_line

      !*** if additional attributes are defined in the
      !*** file definition, write these as well

      if (associated(data_file%add_attrib_cval)) then
         do n=1,size(data_file%add_attrib_cval)

            work_line = char_blank
            cindx1 = 1
            cindx2 = len_trim(data_file%add_attrib_cname(n))
            work_line(cindx1:cindx2) = trim(data_file%add_attrib_cname(n))
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            cindx2 = cindx1 + 3
            work_line(cindx1:cindx2) = 'char'
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            write(work_line(cindx1:),*) trim(data_file%add_attrib_cval(n))
            write(hdr_unit,hdr_fmt) work_line

         end do
      endif ! cval

      if (associated(data_file%add_attrib_lval)) then
         do n=1,size(data_file%add_attrib_lval)

            work_line = char_blank
            cindx1 = 1
            cindx2 = len_trim(data_file%add_attrib_lname(n))
            work_line(cindx1:cindx2)=trim(data_file%add_attrib_lname(n))
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            cindx2 = cindx1 + 2
            work_line(cindx1:cindx2) = 'log'
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            write(work_line(cindx1:),*) data_file%add_attrib_lval(n)
            write(hdr_unit,hdr_fmt) work_line

         end do
      endif ! lval

      if (associated(data_file%add_attrib_ival)) then
         do n=1,size(data_file%add_attrib_ival)

            work_line = char_blank
            cindx1 = 1
            cindx2 = len_trim(data_file%add_attrib_iname(n))
            work_line(cindx1:cindx2)=trim(data_file%add_attrib_iname(n))
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            cindx2 = cindx1 + 2
            work_line(cindx1:cindx2) = 'int'
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            write(work_line(cindx1:),*) data_file%add_attrib_ival(n)
            write(hdr_unit,hdr_fmt) work_line

         end do
      endif ! ival

      if (associated(data_file%add_attrib_rval)) then
         do n=1,size(data_file%add_attrib_rval)

            work_line = char_blank
            cindx1 = 1
            cindx2 = len_trim(data_file%add_attrib_rname(n))
            work_line(cindx1:cindx2) = &
                                trim(data_file%add_attrib_rname(n))
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            cindx2 = cindx1 + 1
            work_line(cindx1:cindx2) = 'r4'
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            write(work_line(cindx1:),*) data_file%add_attrib_rval(n)
            write(hdr_unit,hdr_fmt) work_line

         end do
      endif ! rval

      if (associated(data_file%add_attrib_dval)) then
         do n=1,size(data_file%add_attrib_dval)

            work_line = char_blank
            cindx1 = 1
            cindx2 = len_trim(data_file%add_attrib_dname(n))
            work_line(cindx1:cindx2)=trim(data_file%add_attrib_dname(n))
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            cindx2 = cindx1 + 1
            work_line(cindx1:cindx2) = 'r8'
            cindx1 = cindx2 + 1
            cindx2 = cindx1
            work_line(cindx1:cindx2) = attrib_separator
            cindx1 = cindx2 + 1
            write(work_line(cindx1:),*) data_file%add_attrib_dval(n)
            write(hdr_unit,hdr_fmt) work_line

         end do
      endif ! dval

   endif ! master task

!-----------------------------------------------------------------------
!EOC

 end subroutine open_binary

!***********************************************************************
!BOP
! !IROUTINE: close_binary
! !INTERFACE:

 subroutine close_binary(data_file)

! !DESCRIPTION:
! This routine closes an open binary data file and the associated
! header file if it exists.
!
! !REVISION HISTORY:
! same as module

! !INPUT/OUTPUT PARAMETERS:

   type (datafile), intent (inout) :: data_file

!EOP
!BOC
!-----------------------------------------------------------------------
!
! close a data file
!
!-----------------------------------------------------------------------

   if (my_task < data_file%num_iotasks) close(data_file%id(1))

   if (my_task == master_task .and. data_file%id(2) > 0) &
      close(data_file%id(2))

!-----------------------------------------------------------------------
!EOC

 end subroutine close_binary

!***********************************************************************
!BOP
! !IROUTINE: define_field_binary
! !INTERFACE:

 subroutine define_field_binary(data_file, io_field)

! !DESCRIPTION:
! This routine defines an io field. When writing a file, this means
! computing the starting record number for the field and then writing
! the field information to the header file. When reading a file,
! the define routine will attempt to fill an io field with information
! from the header file; if a header file does not exist, it will leave
! the io field mostly undefined and will rely on an input record from
! the read routine to locate the field in the file (this latter is
! provided for compatibility with forcing and grid files).
!
! !REVISION HISTORY:
! same as module

! !INPUT/OUTPUT PARAMETERS:

   type (datafile), intent (inout) :: &
      data_file

   type (io_field_desc), intent (inout) :: &
      io_field

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   character (80) :: &
      work_line, &! workspace for manipulating input string
      comp_line, &! comparison string
      att_name, &! attribute name string
      ctype ! attribute data type

   integer (i4) :: &
      n, &! loop index
      att_ival, &! temp value for integer attribute
      hdr_error, &! io error status for header file
      cindx1, &! character string index
      cindx2, &! character string index
      unit ! unit for header file

   logical (log_kind) :: &
      att_lval ! temp value for logical attribute

   real (r4) :: &
      att_rval ! temp value for real attribute

   real (r8) :: &
      att_dval ! temp value for double attribute

!-----------------------------------------------------------------------
!
! for output files, set id to current record number and write
! field attributes to header file.
!
!-----------------------------------------------------------------------

   if (.not. data_file%readonly) then

      !*** define id as current record

      io_field%id = data_file%current_record

      !*** increment record counter for next field

      if (io_field%nfield_dims == 2) then

         data_file%current_record = data_file%current_record + 1

      else if (io_field%nfield_dims == 3) then

         if (associated(io_field%field_i_3d)) then
            n = size(io_field%field_i_3d,dim=3)
         else if (associated(io_field%field_r_3d)) then
            n = size(io_field%field_r_3d,dim=3)
         else if (associated(io_field%field_d_3d)) then
            n = size(io_field%field_d_3d,dim=3)
         endif

         data_file%current_record = data_file%current_record + n

      endif

      !***
      !*** now write attributes to header file
      !***

      unit = data_file%id(2)

      if (my_task == master_task) then

         !*** first write variable name header

         write(unit,'(a1,a)') '&',trim(io_field%short_name)

         !*** write standard attributes to header file

         if (io_field%long_name /= char_blank) &
            write(unit,'(a9,a1,a4,a1,a)') 'long_name', attrib_separator, &
                                          'char', attrib_separator, &
                                          trim(io_field%long_name)

         if (io_field%units /= char_blank) &
            write(unit,'(a5,a1,a4,a1,a)') 'units', attrib_separator, &
                                          'char' , attrib_separator, &
                                          trim(io_field%units)

         if (io_field%grid_loc /= '    ') &
            write(unit,'(a8,a1,a4,a1,a4)') 'grid_loc', attrib_separator, &
                                           'char' , attrib_separator, &
                                           io_field%grid_loc

         if (io_field%missing_value /= undefined) then
            work_line = char_blank
            work_line( 1:13) = 'missing_value'
            work_line(14:14) = attrib_separator
            work_line(15:16) = 'r8'
            work_line(17:17) = attrib_separator
            write(work_line(18:),*) io_field%missing_value
            write(unit,'(a)') trim(work_line)
         endif

         if (any(io_field%valid_range(:) /= undefined)) then
            work_line = char_blank
            work_line( 1:11) = 'valid_range'
            work_line(12:12) = attrib_separator
            work_line(13:14) = 'r8'
            work_line(15:15) = attrib_separator
            write(work_line(16:),*) io_field%valid_range(:)
            write(unit,'(a)') trim(work_line)
         endif

         work_line = char_blank
         work_line(1:2) = 'id'
         work_line(3:3) = attrib_separator
         work_line(4:6) = 'int'
         work_line(7:7) = attrib_separator
         write(work_line(8:),*) io_field%id
         write(unit,'(a)') trim(work_line)

         !**** number of dimensions and dimension info

         work_line = char_blank
         work_line( 1:11) = 'nfield_dims'
         work_line(12:12) = attrib_separator
         work_line(13:15) = 'int'
         work_line(16:16) = attrib_separator
         write(work_line(17:),*) io_field%nfield_dims
         write(unit,'(a)') trim(work_line)

         !**** not sure what to do about dimension info

         ! write field_dim stuff

         !*** write additional attributes to header file

         if (associated(io_field%add_attrib_cval)) then
            do n=1,size(io_field%add_attrib_cval)

               work_line = char_blank
               cindx1 = len_trim(io_field%add_attrib_cname(n))
               cindx2 = cindx1 + 7
               work_line(1:cindx1) = trim(io_field%add_attrib_cname(n))
               work_line(cindx1+1:cindx1+1) = attrib_separator
               work_line(cindx1+2:cindx2-2) = 'char'
               work_line(cindx2-1:cindx2-1) = attrib_separator
               cindx1 = cindx2 + len_trim(io_field%add_attrib_cval(n))
               work_line(cindx2:cindx2) = trim(io_field%add_attrib_cval(n))
               write(unit,'(a)') trim(work_line)
            end do
         endif

         if (associated(io_field%add_attrib_lval)) then
            do n=1,size(io_field%add_attrib_lval)

               work_line = char_blank
               cindx1 = len_trim(io_field%add_attrib_lname(n))
               cindx2 = cindx1 + 6
               work_line(1:cindx1) = trim(io_field%add_attrib_lname(n))
               work_line(cindx1+1:cindx1+1) = attrib_separator
               work_line(cindx1+2:cindx2-2) = 'log'
               work_line(cindx2-1:cindx2-1) = attrib_separator
               write(work_line(cindx2:),*) io_field%add_attrib_lval(n)
               write(unit,'(a)') trim(work_line)
            end do
         endif

         if (associated(io_field%add_attrib_ival)) then
            do n=1,size(io_field%add_attrib_ival)

               work_line = char_blank
               cindx1 = len_trim(io_field%add_attrib_iname(n))
               cindx2 = cindx1 + 6
               work_line(1:cindx1) = trim(io_field%add_attrib_iname(n))
               work_line(cindx1+1:cindx1+1) = attrib_separator
               work_line(cindx1+2:cindx2-2) = 'int'
               work_line(cindx2-1:cindx2-1) = attrib_separator
               write(work_line(cindx2:),*) io_field%add_attrib_ival(n)
               write(unit,'(a)') trim(work_line)
            end do
         endif

         if (associated(io_field%add_attrib_rval)) then
            do n=1,size(io_field%add_attrib_rval)

               work_line = char_blank
               cindx1 = len_trim(io_field%add_attrib_rname(n))
               cindx2 = cindx1 + 5
               work_line(1:cindx1) = trim(io_field%add_attrib_rname(n))
               work_line(cindx1+1:cindx1+1) = attrib_separator
               work_line(cindx1+2:cindx2-2) = 'r4'
               work_line(cindx2-1:cindx2-1) = attrib_separator
               write(work_line(cindx2:),*) io_field%add_attrib_rval(n)
               write(unit,'(a)') trim(work_line)
            end do
         endif

         if (associated(io_field%add_attrib_dval)) then
            do n=1,size(io_field%add_attrib_dval)

               work_line = char_blank
               cindx1 = len_trim(io_field%add_attrib_dname(n))
               cindx2 = cindx1 + 5
               work_line(1:cindx1) = trim(io_field%add_attrib_dname(n))
               work_line(cindx1+1:cindx1+1) = attrib_separator
               work_line(cindx1+2:cindx2-2) = 'r8'
               work_line(cindx2-1:cindx2-1) = attrib_separator
               write(work_line(cindx2:),*) io_field%add_attrib_dval(n)
               write(unit,'(a)') work_line
               write(unit,'(a)') trim(work_line)
            end do
         endif

         !*** end variable attribute section with a delimiter

         write(unit,'(a1)') '/'

      endif ! master task

!-----------------------------------------------------------------------
!
! attempt to define an io field from an input header file
!
!-----------------------------------------------------------------------

   else ! this is an input file

      unit = data_file%id(2)
      if (unit <= 0) then ! no header file, assume fields are defined
                           ! in the order they exist in input file
                           ! set id as current record and increment

         io_field%id = data_file%current_record

         if (associated(io_field%field_i_2d)) then
            data_file%current_record = data_file%current_record + 1
         else if (associated(io_field%field_i_3d)) then
            data_file%current_record = data_file%current_record + &
                                       size(io_field%field_i_3d, dim=3)
         else if (associated(io_field%field_r_2d)) then
            data_file%current_record = data_file%current_record + 1
         else if (associated(io_field%field_r_3d)) then
            data_file%current_record = data_file%current_record + &
                                       size(io_field%field_r_3d, dim=3)
         else if (associated(io_field%field_d_2d)) then
            data_file%current_record = data_file%current_record + 1
         else if (associated(io_field%field_d_3d)) then
            data_file%current_record = data_file%current_record + &
                                       size(io_field%field_d_3d, dim=3)
         else
            call exit_POP(sigAbort, &
             'define: No known binary field descriptor associated')
         end if

      else ! header exists: read all attributes from header file

         if (my_task == master_task) then
            hdr_error = 0
            rewind (unit)
            cindx1 = len_trim(io_field%short_name) + 1
            comp_line(1:1) = '&'
            comp_line(2:cindx1) = trim(io_field%short_name)

            srch_loop: do while (hdr_error == 0) ! look for field in file

               read(unit,'(a80)',iostat=hdr_error) work_line
               work_line = adjustl(work_line)
               if (work_line(1:cindx1) == comp_line(1:cindx1)) &
                  exit srch_loop
            end do srch_loop
         endif ! master_task

         call broadcast_scalar(hdr_error, master_task)
         if (hdr_error /= 0) call exit_POP(sigAbort,&
            'could not find field in binary header file')

         !*** found io_field definition - now extract field attributes

         att_loop: do

            if (my_task == master_task) then
               read(unit,'(a80)',iostat=hdr_error) work_line
               work_line = adjustl(work_line)
            endif

            call broadcast_scalar(hdr_error, master_task)
            call broadcast_scalar(work_line, master_task)

            !*** exit if error reading or reached end of definition

            if (hdr_error /= 0 .or. work_line(1:1) == '/') &
               exit att_loop

            !*** extract attribute name

            cindx1 = index(work_line,attrib_separator)
            att_name = char_blank
            att_name(1:cindx1-1) = work_line(1:cindx1-1)

            !*** strip attribute name from input line

            comp_line = work_line
            do n=1,cindx1
               comp_line(n:n) = ' '
            end do
            work_line = adjustl(comp_line)

            !*** extract attribute data type

            cindx1 = index(work_line,attrib_separator)
            ctype = char_blank
            ctype(1:cindx1-1) = work_line(1:cindx1-1)

            !*** strip data type from input line

            comp_line = work_line
            do n=1,cindx1
               comp_line(n:n) = ' '
            end do
            work_line = adjustl(comp_line)

            !*** if this attribute matches a standard attribute
            !*** set attribute value. otherwise, define an
            !*** additional attribute

            select case (trim(att_name))

            case('long_name','LONG_NAME')
               io_field%long_name = char_blank
               io_field%long_name = trim(work_line)

            case('units','UNITS')
               io_field%units = char_blank
               io_field%units = trim(work_line)

            case('grid_loc','GRID_LOC')
               io_field%grid_loc = char_blank
               io_field%grid_loc = trim(work_line)

            case('missing_value','MISSING_VALUE')
               read(work_line,*) io_field%missing_value

            case('valid_range','VALID_RANGE')
               read(work_line,*) io_field%valid_range(:)

            case('id','ID')
               read(work_line,*) io_field%id

            case('field_dim','FIELD_DIM')
               !**** do not know what to do about dimensions yet

            case('nfield_dims','NFIELD_DIMS')
               read(work_line,*) io_field%nfield_dims

            case default

               !*** add attribute to io field

               select case (ctype)
               case ('char','CHAR','character','CHARACTER')
                  call add_attrib_io_field(io_field, trim(att_name), &
                                                     trim(work_line))

               case ('log','LOG','logical','LOGICAL')
                  read(work_line,*) att_lval
                  call add_attrib_io_field(io_field, trim(att_name), &
                                                     att_lval)

               case ('int','INT','i4','I4','integer','INTEGER')
                  read(work_line,*) att_ival
                  call add_attrib_io_field(io_field, trim(att_name), &
                                                     att_ival)

               case ('r4','R4','real','REAL','float','FLOAT')
                  read(work_line,*) att_rval
                  call add_attrib_io_field(io_field, trim(att_name), &
                                                     att_rval)

               case ('r8','R8','dbl','DBL','double','DOUBLE')
                  read(work_line,*) att_dval
                  call add_attrib_io_field(io_field, trim(att_name), &
                                                     att_dval)

               case default
                  call exit_POP(sigAbort, &
                                'define_io_field: unknown data type')
               end select

            end select ! att_name

         end do att_loop

         if (hdr_error /= 0) call exit_POP(sigAbort, &
            'define_io_field: error reading attribute from header')

      endif ! unit > 0 (header exists)

   endif !readonly

!-----------------------------------------------------------------------
!EOC

 end subroutine define_field_binary

!***********************************************************************
!BOP
! !IROUTINE: write_field_binary
! !INTERFACE:

 subroutine write_field_binary(data_file, io_field)

! !INPUT PARAMETERS:

   type (datafile), intent (in) :: &
      data_file ! file to which data will be written

   type (io_field_desc), intent (in) :: &
      io_field ! field to be written

! !DESCRIPTION:
! This routine writes a binary field to the data file.
!
! !REVISION HISTORY:
! same as module

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   integer (i4) :: &
      record ! starting record for writing array

!-----------------------------------------------------------------------
!
! extract the starting record from the io field
!
!-----------------------------------------------------------------------

   record = io_field%id

!-----------------------------------------------------------------------
!
! write the io field
!
!-----------------------------------------------------------------------

   if (associated(io_field%field_i_2d)) then
      call write_array(data_file,io_field%field_i_2d,record)
   else if (associated(io_field%field_i_3d)) then
      call write_array(data_file,io_field%field_i_3d,record)
   else if (associated(io_field%field_r_2d)) then
      call write_array(data_file,io_field%field_r_2d,record)
   else if (associated(io_field%field_r_3d)) then
      call write_array(data_file,io_field%field_r_3d,record)
   else if (associated(io_field%field_d_2d)) then
      call write_array(data_file,io_field%field_d_2d,record)
   else if (associated(io_field%field_d_3d)) then
      call write_array(data_file,io_field%field_d_3d,record)
   else
      call exit_POP(sigAbort, &
                   'write: No known binary field descriptor associated')
   end if

!-----------------------------------------------------------------------
!EOC

 end subroutine write_field_binary

!***********************************************************************
!BOP
! !IROUTINE: read_field_binary
! !INTERFACE:

 subroutine read_field_binary(data_file, io_field)

! !DESCRIPTION:
! This routine reads a field from a binary input file.
!
! !REVISION HISTORY:
! same as module

! !INPUT/OUTPUT PARAMETERS:

   type (datafile), intent (inout) :: &
      data_file ! file from which to read data

   type (io_field_desc), intent (inout) :: &
      io_field ! data to be read

!EOP
!BOC
!-----------------------------------------------------------------------
!
! if the field location parameters undefined, assume center scalar
!
!-----------------------------------------------------------------------

   if (io_field%field_loc == field_loc_unknown) then
      io_field%field_loc = field_loc_center
      io_field%field_type = field_type_scalar
   endif

!-----------------------------------------------------------------------
!
! read an io field
!
!-----------------------------------------------------------------------

   if (associated(io_field%field_i_2d)) then
      call read_array(data_file,io_field%field_i_2d,io_field%id, &
                                io_field%field_loc, io_field%field_type)
   else if (associated(io_field%field_i_3d)) then
      call read_array(data_file,io_field%field_i_3d,io_field%id, &
                                io_field%field_loc, io_field%field_type)
   else if (associated(io_field%field_r_2d)) then
      call read_array(data_file,io_field%field_r_2d,io_field%id, &
                                io_field%field_loc, io_field%field_type)
   else if (associated(io_field%field_r_3d)) then
      call read_array(data_file,io_field%field_r_3d,io_field%id, &
                                io_field%field_loc, io_field%field_type)
   else if (associated(io_field%field_d_2d)) then
      call read_array(data_file,io_field%field_d_2d,io_field%id, &
                                io_field%field_loc, io_field%field_type)
   else if (associated(io_field%field_d_3d)) then
      call read_array(data_file,io_field%field_d_3d,io_field%id, &
                                io_field%field_loc, io_field%field_type)
   else
      call exit_POP(sigAbort,'read_field: field not associated')
   end if

!-----------------------------------------------------------------------
!EOC

 end subroutine read_field_binary

!***********************************************************************
!BOP
! !IROUTINE: read_int_2d
! !INTERFACE:

 subroutine read_int_2d(data_file, INT2D, start_record, &
                        field_loc, field_type)

! !DESCRIPTION:
! Reads a 2-d horizontal slice of integers from a binary file
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   type (datafile), intent(in) :: &
      data_file ! info about data file

   integer (i4), intent(in) :: &
      start_record, &! starting record of field in file
      field_loc, &! loc of field on horiz grid
      field_type ! type of field (scalar, vector, angle)

! !INPUT/OUTPUT PARAMETERS:

   integer (i4), dimension(:,:,:), intent(inout) :: &
      INT2D ! array to be read

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   integer (i4), dimension(:,:), allocatable :: &
      IOBUFI ! local global-sized buffer

!-----------------------------------------------------------------------
!
! read in global 2-d slice from one processor
!
!-----------------------------------------------------------------------

   if (my_task == master_task) then
      allocate(IOBUFI(nx_global,ny_global))
      read(data_file%id(1),rec=start_record) IOBUFI
   endif

!-----------------------------------------------------------------------
!
! send chunks to processors who own them
!
!-----------------------------------------------------------------------

   call scatter_global(INT2D, IOBUFI, master_task, distrb_clinic, &
                       field_loc, field_type)

   if (my_task == master_task) deallocate(IOBUFI)

!-----------------------------------------------------------------------
!EOC

 end subroutine read_int_2d

!***********************************************************************
!BOP
! !IROUTINE: read_real4_2d
! !INTERFACE:

 subroutine read_real4_2d(data_file, REAL2D, start_record, &
                          field_loc, field_type)

! !DESCRIPTION:
! Reads a 2-d horizontal slice of reals from a binary file
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   type (datafile), intent(in) :: &
      data_file ! info about data file

   integer (i4), intent(in) :: &
      start_record, &! starting record of field in file
      field_loc, &! loc of field on horiz grid
      field_type ! type of field (scalar, vector, angle)

! !INPUT/OUTPUT PARAMETERS:

   real (r4), dimension(:,:,:), intent(inout) :: &
      REAL2D ! array to be read

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   real (r4), dimension(:,:), allocatable :: &
      IOBUFR ! local global-sized buffer

!-----------------------------------------------------------------------
!
! read in global 2-d slice from one processor
!
!-----------------------------------------------------------------------

   if (my_task == master_task) then
      allocate(IOBUFR(nx_global,ny_global))
      read(data_file%id(1),rec=start_record) IOBUFR
   endif

!-----------------------------------------------------------------------
!
! send chunks to processors who own them
!
!-----------------------------------------------------------------------

   call scatter_global(REAL2D, IOBUFR, master_task, distrb_clinic, &
                       field_loc, field_type)

   if (my_task == master_task) deallocate(IOBUFR)

!-----------------------------------------------------------------------
!EOC

 end subroutine read_real4_2d

!***********************************************************************
!BOP
! !IROUTINE: read_real8_2d
! !INTERFACE:

 subroutine read_real8_2d(data_file, DBL2D, start_record, &
                          field_loc, field_type)

! !DESCRIPTION:
! Reads a 2-d horizontal slice of 64-bit reals from a binary file
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   type (datafile), intent(in) :: &
      data_file ! info about data file

   integer (i4), intent(in) :: &
      start_record, &! starting record of field in file
      field_loc, &! loc of field on horiz grid
      field_type ! type of field (scalar, vector, angle)

! !INPUT/OUTPUT PARAMETERS:

   real (r8), dimension(:,:,:), intent(inout) :: &
      DBL2D ! array to be read

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   real (r8), dimension(:,:), allocatable :: &
      IOBUFD ! local global-sized buffer

!-----------------------------------------------------------------------
!
! read in global 2-d slice from one processor
!
!-----------------------------------------------------------------------

   if (my_task == master_task) then
      allocate(IOBUFD(nx_global,ny_global))
      read(data_file%id(1),rec=start_record) IOBUFD
   endif

!-----------------------------------------------------------------------
!
! send chunks to processors who own them
!
!-----------------------------------------------------------------------

   call scatter_global(DBL2D, IOBUFD, master_task, distrb_clinic, &
                       field_loc, field_type)

   if (my_task == master_task) deallocate(IOBUFD)

!-----------------------------------------------------------------------
!EOC

 end subroutine read_real8_2d

!***********************************************************************
!BOP
! !IROUTINE: read_int_3d
! !INTERFACE:

 subroutine read_int_3d(data_file, INT3D, start_record, &
                        field_loc, field_type)

! !DESCRIPTION:
! Reads a 3-d integer array from a binary file
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   type (datafile), intent(in) :: &
      data_file ! info on input data file

   integer (i4), intent(in) :: &
      start_record, &! starting record of field in file
      field_loc, &! loc of field on horiz grid
      field_type ! type of field (scalar, vector, angle)

! !INPUT/OUTPUT PARAMETERS:

   integer (i4), dimension(:,:,:,:), intent(inout) :: &
      INT3D ! array to be read

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   integer (i4) k,n,krecs,klvl,nz

   integer (i4), dimension(:,:), allocatable :: &
      IOBUFI ! global-sized array buffer

!-----------------------------------------------------------------------
!
! determine the number of records each i/o process must read to
! get all the records
!
!-----------------------------------------------------------------------

   nz = size(INT3D, DIM=3)

   if (mod(nz,data_file%num_iotasks) == 0) then
      krecs = nz/data_file%num_iotasks
   else
      krecs = nz/data_file%num_iotasks + 1
   endif

   if (my_task < data_file%num_iotasks) &
      allocate(IOBUFI(nx_global,ny_global))

!-----------------------------------------------------------------------
!
! each i/o process reads a horizontal slab from a record and sends
! the data to processors who own it. read and distribute num_iotasks
! records at a time to keep the messages from getting mixed up.
!
!-----------------------------------------------------------------------

   do k=1,krecs

      if (my_task < data_file%num_iotasks) then
         klvl = (k-1)*data_file%num_iotasks + my_task + 1
         if (klvl <= nz) then
            read(data_file%id(1), rec=start_record+klvl-1) IOBUFI
         endif
      endif

      do n=1,data_file%num_iotasks
         klvl = (k-1)*data_file%num_iotasks + n
         if (klvl <= nz) then
            call scatter_global(INT3D(:,:,klvl,:), IOBUFI, n-1, &
                                distrb_clinic, field_loc, field_type)
         endif
      end do

   end do

   if (my_task < data_file%num_iotasks) deallocate(IOBUFI)

!-----------------------------------------------------------------------
!EOC

 end subroutine read_int_3d

!***********************************************************************
!BOP
! !IROUTINE: read_real4_3d
! !INTERFACE:

 subroutine read_real4_3d(data_file, REAL3D, start_record, &
                          field_loc, field_type)

! !DESCRIPTION:
! Reads a 3-d real array from a binary file
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   type (datafile), intent(in) :: &
      data_file ! info on input data file

   integer (i4), intent(in) :: &
      start_record, &! starting record of field in file
      field_loc, &! loc of field on horiz grid
      field_type ! type of field (scalar, vector, angle)

! !INPUT/OUTPUT PARAMETERS:

   real (r4), dimension(:,:,:,:), intent(inout) :: &
      REAL3D ! array to be read

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   integer (i4) k,n,krecs,klvl,nz

   real (r4), dimension(:,:), allocatable :: &
      IOBUFR ! global-sized array buffer

!-----------------------------------------------------------------------
!
! determine the number of records each i/o process must read to
! get all the records
!
!-----------------------------------------------------------------------

   nz = size(REAL3D, DIM=3)

   if (mod(nz,data_file%num_iotasks) == 0) then
      krecs = nz/data_file%num_iotasks
   else
      krecs = nz/data_file%num_iotasks + 1
   endif

   if (my_task < data_file%num_iotasks) &
      allocate(IOBUFR(nx_global,ny_global))

!-----------------------------------------------------------------------
!
! each i/o process reads a horizontal slab from a record and sends
! the data to processors who own it. read and distribute num_iotasks
! records at a time to keep the messages from getting mixed up.
!
!-----------------------------------------------------------------------

   do k=1,krecs

      if (my_task < data_file%num_iotasks) then
         klvl = (k-1)*data_file%num_iotasks + my_task + 1
         if (klvl <= nz) then
            read(data_file%id(1), rec=start_record+klvl-1) IOBUFR
         endif
      endif

      do n=1,data_file%num_iotasks
         klvl = (k-1)*data_file%num_iotasks + n
         if (klvl <= nz) then
            call scatter_global(REAL3D(:,:,klvl,:), IOBUFR, n-1, &
                                distrb_clinic, field_loc, field_type)
         endif
      end do

   end do

   if (my_task < data_file%num_iotasks) deallocate(IOBUFR)

!-----------------------------------------------------------------------
!EOC

 end subroutine read_real4_3d

!***********************************************************************
!BOP
! !IROUTINE: read_real8_3d
! !INTERFACE:

 subroutine read_real8_3d(data_file, DBL3D, start_record, &
                          field_loc, field_type)

! !DESCRIPTION:
! Reads a 3-d 64-bit real array from a binary file
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   type (datafile), intent(in) :: &
      data_file ! info on input data file

   integer (i4), intent(in) :: &
      start_record, &! starting record of field in file
      field_loc, &! loc of field on horiz grid
      field_type ! type of field (scalar, vector, angle)

! !INPUT/OUTPUT PARAMETERS:

   real (r8), dimension(:,:,:,:), intent(inout) :: &
      DBL3D ! array to be read

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   integer (i4) k,n,krecs,klvl,nz

   real (r8), dimension(:,:), allocatable :: &
      IOBUFD ! global-sized array buffer

!-----------------------------------------------------------------------
!
! determine the number of records each i/o process must read to
! get all the records
!
!-----------------------------------------------------------------------

   nz = size(DBL3D, DIM=3)

   if (mod(nz,data_file%num_iotasks) == 0) then
      krecs = nz/data_file%num_iotasks
   else
      krecs = nz/data_file%num_iotasks + 1
   endif

   if (my_task < data_file%num_iotasks) &
      allocate(IOBUFD(nx_global,ny_global))

!-----------------------------------------------------------------------
!
! each i/o process reads a horizontal slab from a record and sends
! the data to processors who own it. read and distribute num_iotasks
! records at a time to keep the messages from getting mixed up.
!
!-----------------------------------------------------------------------

   do k=1,krecs

      if (my_task < data_file%num_iotasks) then
         klvl = (k-1)*data_file%num_iotasks + my_task + 1
         if (klvl <= nz) then
            read(data_file%id(1), rec=start_record+klvl-1) IOBUFD
         endif
      endif

      do n=1,data_file%num_iotasks
         klvl = (k-1)*data_file%num_iotasks + n
         if (klvl <= nz) then
            call scatter_global(DBL3D(:,:,klvl,:), IOBUFD, n-1, &
                                distrb_clinic, field_loc, field_type)
         endif
      end do

   end do

   if (my_task < data_file%num_iotasks) deallocate(IOBUFD)

!-----------------------------------------------------------------------
!EOC

 end subroutine read_real8_3d

!***********************************************************************
!BOP
! !IROUTINE: write_int_2d
! !INTERFACE:

 subroutine write_int_2d(data_file,INT2D,start_record)

! !DESCRIPTION:
! Writes a 2-d slab of integers to a binary file.
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   type (datafile), intent(in) :: &
      data_file ! file information

   integer (i4), intent(in) :: &
      start_record ! starting record of array in file

   integer (i4), dimension(:,:,:), intent(in) :: &
      INT2D ! array to be written

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   integer (i4), dimension(:,:), allocatable :: &
      IOBUFI ! local global-sized buffer

!-----------------------------------------------------------------------
!
! receive chunks from processors who own them
!
!-----------------------------------------------------------------------

   if (my_task == master_task) allocate(IOBUFI(nx_global,ny_global))

   call gather_global(IOBUFI, INT2D, master_task, distrb_clinic)

!-----------------------------------------------------------------------
!
! write global 2-d slice from one processor
!
!-----------------------------------------------------------------------

   if (my_task == master_task) then
      write(data_file%id(1),rec=start_record) IOBUFI
      deallocate(IOBUFI)
   endif

!-----------------------------------------------------------------------
!EOC
 end subroutine write_int_2d

!***********************************************************************
!BOP
! !IROUTINE: write_real4_2d
! !INTERFACE:

 subroutine write_real4_2d(data_file,REAL2D,start_record)

! !DESCRIPTION:
! Writes a 2-d slab of reals to a binary file.
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   type (datafile), intent(in) :: &
      data_file ! file information

   integer (i4), intent(in) :: &
      start_record ! starting record of array in file

   real (r4), dimension(:,:,:), intent(in) :: &
      REAL2D ! array to be written

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   real (r4), dimension(:,:), allocatable :: &
      IOBUFR ! local global-sized buffer

!-----------------------------------------------------------------------
!
! receive chunks from processors who own them
!
!-----------------------------------------------------------------------

   if (my_task == master_task) allocate(IOBUFR(nx_global,ny_global))
   call gather_global(IOBUFR, REAL2D, master_task, distrb_clinic)

!-----------------------------------------------------------------------
!
! write global 2-d slice from one processor
!
!-----------------------------------------------------------------------

   if (my_task == master_task) then
      write(data_file%id(1),rec=start_record) IOBUFR
      deallocate(IOBUFR)
   endif

!-----------------------------------------------------------------------
!EOC
 end subroutine write_real4_2d

!***********************************************************************
!BOP
! !IROUTINE: write_real8_2d
! !INTERFACE:

 subroutine write_real8_2d(data_file,DBL2D,start_record)

! !DESCRIPTION:
! Writes a 2-d slab of doubles to a binary file.
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   type (datafile), intent(in) :: &
      data_file ! file information

   integer (i4), intent(in) :: &
      start_record ! starting record of array in file

   real (r8), dimension(:,:,:), intent(in) :: &
      DBL2D ! array to be written

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   real (r8), dimension(:,:), allocatable :: &
      IOBUFD ! local global-sized buffer

!-----------------------------------------------------------------------
!
! receive chunks from processors who own them
!
!-----------------------------------------------------------------------

   if (my_task == master_task) allocate(IOBUFD(nx_global,ny_global))
   call gather_global(IOBUFD, DBL2D, master_task, distrb_clinic)

!-----------------------------------------------------------------------
!
! write global 2-d slice from one processor
!
!-----------------------------------------------------------------------

   if (my_task == master_task) then
      write(data_file%id(1),rec=start_record) IOBUFD
      deallocate(IOBUFD)
   endif

!-----------------------------------------------------------------------
!EOC
 end subroutine write_real8_2d

!***********************************************************************
!BOP
! !IROUTINE: write_int_3d
! !INTERFACE:

 subroutine write_int_3d(data_file,INT3D,start_record)

! !DESCRIPTION:
! Writes a 3-d integer array as a series of 2-d slabs to a binary
! output file.
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   type (datafile), intent(in) :: &
      data_file ! info on output data file

   integer (i4), intent(in) :: &
      start_record ! starting position of array in file

   integer (i4), dimension(:,:,:,:), intent(in) :: &
      INT3D ! array to be written

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   integer (i4) :: &
      k, n, & ! dummy counters
      krecs, & ! number of records each iotask must write
      klvl, & ! k index corresponding to current record
      nz ! size of 3rd dimension of 3-d array to be written

   integer (i4), dimension(:,:), allocatable :: &
      IOBUFI ! global-sized buffer array

!-----------------------------------------------------------------------
!
! determine the number of records each i/o process must write to
! get all the records
!
!-----------------------------------------------------------------------

   nz = size(INT3D, DIM=3)

   if (mod(nz,data_file%num_iotasks) == 0) then
      krecs = nz/data_file%num_iotasks
   else
      krecs = nz/data_file%num_iotasks + 1
   endif

   if (my_task < data_file%num_iotasks) &
      allocate(IOBUFI(nx_global,ny_global))

!-----------------------------------------------------------------------
!
! gather and write num_iotasks records at a time to keep the
! messages from getting mixed up.
!
!-----------------------------------------------------------------------

   do k=1,krecs

!-----------------------------------------------------------------------
!
! gather a global slice on each iotask
!
!-----------------------------------------------------------------------

      do n=1,data_file%num_iotasks
         klvl = (k-1)*data_file%num_iotasks + n
         if (klvl <= nz) then
            call gather_global(IOBUFI, INT3D(:,:,klvl,:), &
                               n-1, distrb_clinic)
         endif
      end do

!-----------------------------------------------------------------------
!
! iotasks wait for the gather to be complete and write global
! slabs to the file
!
!-----------------------------------------------------------------------

      if (my_task < data_file%num_iotasks) then
         klvl = (k-1)*data_file%num_iotasks + my_task + 1
         if (klvl <= nz) then
            write(data_file%id(1), rec=start_record+klvl-1) IOBUFI
         endif
      endif

   end do

!-----------------------------------------------------------------------

   if (my_task < data_file%num_iotasks) deallocate(IOBUFI)

!-----------------------------------------------------------------------
!EOC

 end subroutine write_int_3d

!***********************************************************************
!BOP
! !IROUTINE: write_real4_3d
! !INTERFACE:

 subroutine write_real4_3d(data_file,REAL3D,start_record)

! !DESCRIPTION:
! Writes a 3-d real array as a series of 2-d slabs to a binary
! output file.
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   type (datafile), intent(in) :: &
      data_file ! info on output data file

   integer (i4), intent(in) :: &
      start_record ! starting position of array in file

   real (r4), dimension(:,:,:,:), intent(in) :: &
      REAL3D ! array to be written

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   integer (i4) :: &
      k, n, & ! dummy counters
      krecs, & ! number of records each iotask must write
      klvl, & ! k index corresponding to current record
      nz ! size of 3rd dimension of 3-d array to be written

   real (r4), dimension(:,:), allocatable :: &
      IOBUFR ! global-sized buffer array

!-----------------------------------------------------------------------
!
! determine the number of records each i/o process must write to
! get all the records
!
!-----------------------------------------------------------------------

   nz = size(REAL3D, DIM=3)

   if (mod(nz,data_file%num_iotasks) == 0) then
      krecs = nz/data_file%num_iotasks
   else
      krecs = nz/data_file%num_iotasks + 1
   endif

   if (my_task < data_file%num_iotasks) &
      allocate(IOBUFR(nx_global,ny_global))

!-----------------------------------------------------------------------
!
! gather and write num_iotasks records at a time to keep the
! messages from getting mixed up.
!
!-----------------------------------------------------------------------

   do k=1,krecs

!-----------------------------------------------------------------------
!
! gather a global slice on each iotask
!
!-----------------------------------------------------------------------

      do n=1,data_file%num_iotasks
         klvl = (k-1)*data_file%num_iotasks + n
         if (klvl <= nz) then
            call gather_global(IOBUFR, REAL3D(:,:,klvl,:), &
                               n-1, distrb_clinic)
         endif
      end do

!-----------------------------------------------------------------------
!
! iotasks wait for the gather to be complete and write global
! slabs to the file
!
!-----------------------------------------------------------------------

      if (my_task < data_file%num_iotasks) then
         klvl = (k-1)*data_file%num_iotasks + my_task + 1
         if (klvl <= nz) then
            write(data_file%id(1), rec=start_record+klvl-1) IOBUFR
         endif
      endif

   end do

!-----------------------------------------------------------------------

   if (my_task < data_file%num_iotasks) deallocate(IOBUFR)

!-----------------------------------------------------------------------
!EOC

 end subroutine write_real4_3d

!***********************************************************************
!BOP
! !IROUTINE: write_real8_3d
! !INTERFACE:

 subroutine write_real8_3d(data_file,DBL3D,start_record)

! !DESCRIPTION:
! Writes a 3-d 64-bit real array as a series of 2-d slabs to a binary
! output file.
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   type (datafile), intent(in) :: &
      data_file ! info on output data file

   integer (i4), intent(in) :: &
      start_record ! starting position of array in file

   real (r8), dimension(:,:,:,:), intent(in) :: &
      DBL3D ! array to be written

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   integer (i4) :: &
      k, n, & ! dummy counters
      krecs, & ! number of records each iotask must write
      klvl, & ! k index corresponding to current record
      nz ! size of 3rd dimension of 3-d array to be written

   real (r8), dimension(:,:), allocatable :: &
      IOBUFD ! global-sized buffer array

!-----------------------------------------------------------------------
!
! determine the number of records each i/o process must write to
! get all the records
!
!-----------------------------------------------------------------------

   nz = size(DBL3D, DIM=3)

   if (mod(nz,data_file%num_iotasks) == 0) then
      krecs = nz/data_file%num_iotasks
   else
      krecs = nz/data_file%num_iotasks + 1
   endif

   if (my_task < data_file%num_iotasks) &
      allocate(IOBUFD(nx_global,ny_global))

!-----------------------------------------------------------------------
!
! gather and write num_iotasks records at a time to keep the
! messages from getting mixed up.
!
!-----------------------------------------------------------------------

   do k=1,krecs

!-----------------------------------------------------------------------
!
! gather a global slice on each iotask
!
!-----------------------------------------------------------------------

      do n=1,data_file%num_iotasks
         klvl = (k-1)*data_file%num_iotasks + n
         if (klvl <= nz) then
            call gather_global(IOBUFD, DBL3D(:,:,klvl,:), &
                               n-1, distrb_clinic)
         endif
      end do

!-----------------------------------------------------------------------
!
! iotasks wait for the gather to be complete and write global
! slabs to the file
!
!-----------------------------------------------------------------------

      if (my_task < data_file%num_iotasks) then
         klvl = (k-1)*data_file%num_iotasks + my_task + 1
         if (klvl <= nz) then
            write(data_file%id(1), rec=start_record+klvl-1) IOBUFD
         endif
      endif

   end do

!-----------------------------------------------------------------------

   if (my_task < data_file%num_iotasks) deallocate(IOBUFD)

!-----------------------------------------------------------------------
!EOC

 end subroutine write_real8_3d

!***********************************************************************

 end module io_binary

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||


text:
      program ffrtest
      implicit none
      real*8,dimension(0:2000) :: F ! np波函数
      real*8 ::b(0:2000,2)
      real*8 :: FFR4 ! interpolation function
      real*8 :: Y
      integer :: i,j

      open(11,file='fort.7',status='old')
      do i=0,2000
      read(11,*)(b(i,j),j=1,2)
      F(i)=b(i,2)
      end do
      Y=4.5
      write(*,*)FFR4(Y,F,2001)

      end program ffrtest

      FUNCTION FFR4(Y,F,N)
            IMPLICIT REAL*8(A-H,O-Z)
            REAL*8 F(N),P,P1,P2,Q,X,FFR4
            REAL*8 Y
            PARAMETER(X=.16666666666667)
            P=Y
            I=P
            write(*,*)I
            IF(I.LE.0) GO TO 2
            IF(I.GE.N-2) GO TO 4
  1         P=P-I
            P1=P-1.
            P2=P-2.
            Q=P+1.
            FFR4=(-P2*F(I)+Q*F(I+3))*P*P1*X+(P1*F(I+1)-P*F(I+2))*Q*P2*.5
            RETURN
  2         IF(I.LT.0) GO TO 3
            I=1
            GO TO 1
  3         FFR4=F(1)
            RETURN
  4         IF(I.GT.N-2) GO TO 5
            I=N-3
            GO TO 1
  5         FFR4=F(N)
            RETURN
            END function

ext: f

max_stars_repo_name: yazhouLY/samples

max_issues_repo_name: yazhouLY/samples

max_forks_repo_name: yazhouLY/samples


text:
      subroutine testcoupl
      implicit none
c
c     local
c
      integer i,iformat
      character*2 ab(2)
      real*8 ene
      double precision  Zero, One, Two, Three, Four, Half, Rt2
      parameter( Zero = 0.0d0, One = 1.0d0, Two = 2.0d0 )
c
c     include
c
      include 'coupl.inc'
      include 'input.inc'
c
c     Common to lh_readin and printout
c
      double precision  alpha, sin2w, gfermi, alfas
      double precision  mtMS,mbMS,mcMS,mtaMS!MSbar masses
      double precision  Vud,Vus             !CKM matrix elements
      common/values/    alpha,sin2w,gfermi,alfas,
     &                  mtMS,mbMS,mcMS,mtaMS,
     &                  Vud
      open(unit=1,file="couplings_check.txt")
c
c output all info
c
 10   format( 1x,a10,2x,f7.3,' GeV        ',a16,2x,f7.4,' GeV' )
 11   format( 1x,a10,2x,f11.5,2x,f11.5,a3,f11.5,2x,f11.5 )
 12   format( 1x,a10,2x,f6.2,a )
 13   format( 1x,a10,2x,f6.4,a )
 14   format( 1x,2(a10,2x,f10.7,2x,f10.7) )
 15   format( 1x,a10,2x,f9.5,a )
 16   format( 1x,a10,2x,f7.5 )
 17   format( 1x,a10,2x,f8.4 )
 18   format( 1x,a10,2x,f8.4,' GeV' )
 19   format( 1x,a10,2x,f6.4,a13,2x,f6.4 )
 20   format( 1x,a10,2x,f11.5,1x,f11.5 )
 21   format( 1x,a10,2x,f8.4,' GeV',1x,a13,2x,f8.4,' GeV' )
 22   format( 1x,a10,2x,f10.8,a13,2x,f6.4 )
 23   format( 1x,a10,2x,f8.4)
 24   format( 1x,a10,2x,f7.3,' GeV        ',a16,2x,f7.4,' GeV  (calc @ LO)')
 25   format( 1x,a10,2x,f7.3,' GeV        ',a16,2x,f7.4,' GeV')

      write(1,20) 'gwwa ', gwwa
      write(1,20) 'gwwz ', gwwz
      write(1,20) 'gwwh ', gwwh
      write(1,20) 'gzzh ', gzzh
      write(1,20) 'ghhh ', ghhh
      write(1,20) 'gwwhh ', gwwhh
      write(1,20) 'gzzhh ', gzzhh
      write(1,20) 'ghhhh ', ghhhh
      write(1,11) 'gal ',gal(1),' ',gal(2)
      write(1,11) 'gau ',gau(1),' ',gau(2)
      write(1,11) 'gad ',gad(1),' ',gad(2)
      write(1,11) 'gwf ',gwf(1),' ',gwf(2)
      write(1,11) 'gzn ',gzn(1),' ',gzn(2)
      write(1,11) 'gzl ',gzl(1),' ',gzl(2)
      write(1,11) 'gzu ',gzu(1),' ',gzu(2)
      write(1,11) 'gzd ',gzd(1),' ',gzd(2)
      write(1,14) 'gHtop ',ghtop(1),' ',ghtop(2)
      write(1,14) 'gHbot ',ghbot(1),' ',ghbot(2)
      write(1,14) 'gHcha ',ghcha(1),' ',ghcha(2)
      write(1,14) 'gHtau ',ghtau(1),' ',ghtau(2)
      write(1,14) 'gg ',gg(1),' ',gg(2)

      return
      end


text:
  Leona Dance is Users/SummerSong The Bummers older, identical twin sister.
  She plays Capoeira and will be attending UC Berkeley in the fall to study Toxicology.
  Leona is a frequent visitor to The Bakery.
  Leona could be considered a Hippies hippy and would probably live in The Domes if she attended UC Davis. This is because Leona loves to frolic in nature, eat Organic organic vegetables, garden and talk to animals.

Links relating to some of Leonas favorite places/things:
http://www.ncgold.com/Museums_Parks/syrp/ The Yuba River
http://www.briarpatchcoop.com/ Briarpatch Coop
http://www.coolhanduke.com Cool Hand Uke 1
http://www.kvmr.org/personalities/d_scanlan.html Cool Hand Uke 2
http://www.oliveoilsource.com/scripts/company_item.asp?p_com261&ptree Organic Ag Advisors

My sister is obsessed with her cat.  His name is Stinkerton Charlie McGee.   Hes awfully stinky.  Like the Users/MiriamKaufman Wombat. Users/SummerSong The Good Twin

Congrats on being accepted to UC Berkeley!  Whoohoooooo!!!!  What a smartypants!  Users/SummerSong xoxo

20050517 00:08:42 nbsp As I type, youre attempting to teach Miriams eightweek old kitten yoga.   Weirdo! :P Users/SummerSong

20050801 13:08:33 nbsp Its your birthday too, yay! I love twins. Happy Birthday NapoLeona! Users/MiriamKaufman

20050801 13:30:16 nbsp NapoLeona BonaButta! Users/SummerSong

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
Here is the text of the DCR Press Release regarding the endorsement of LEAD for ASUCD President President/VicePresident

DAVIS COLLEGE REPUBLICANS

www.ucdcr.com

Contact: George Andrews  gjandrews@ucdavis.edu

DCR Endorses LEAD Presidential Ticket: Over 5 years of bitter feuds between LEAD and DCR come to a standstill

WEDNESDAY, February 16, 2005   The Davis College Republicans voted unanimously to endorse the LEAD Presidential ticket after hearing from both LEAD and Student Focus slates at this Wednesdays weekly meeting.

We believe that Student Focus has shown over the past year and continues to demonstrate that they are no longer worthy of our staunch support due to the abuse of their veto pen, said DCR Chairman Ryan Clumpner. Despite significant historical differences, the current LEAD slate has demonstrated to share many of the same values as our members, such as fiscal responsibility, increased accountability to the student body and respect for all groups on campus. On the issues that matter most, we are virtually in complete agreement.

The DCRLEAD coalition solidified when the slate made a promise to support both the Pledge bill and the Flag bill during a question and answer session at the DCR meeting. The Pledge bill would make the Pledge of Allegiance a part of ASUCD Senate meetings. The Flag bill would put a U.S and California State flag in the ASUCD Senate room for meetings. Currently, the vast majority of legislative bodies in the United States display the U.S Flag and use the Pledge of Allegiance.

The Pledge and Flag bill have been historically important issues to our membership. Student Focus has yet to even consider this issue since taking control of the Executive office after years of DCRs unwavering support said Clumpner.

DCR intends to throw the full weight of its support behind LEAD, similar to the DCR efforts undertaken during the stunning election of ASUCD Senator Keith Shively in Fall 2004. DCR believes this coalition will be a powerful force in upcoming election and we intend to make our endorsement meaningful, said Clumpner. This fracture in our once concrete Focus voting bloc can result in a major swing in the election. We expect a high turnout in favor of LEAD among our members and look forward to working with LEAD to make ASUCD a respectable and effective student government once again.

30

20050217 09:36:42 nbsp Hey George. Who is DCR endorsing for Senate? Users/RobRoy

20050217 10:15:59 nbsp fox Users/GeorgejAndrews

20050217 18:03:47 nbsp Its funny that after 5 years of voting for the Student Focus lesser evil, Davis College Republicans DCR has finally realized that the Student Focus lesser evil has been useless in providing anyone other than themselves with anything of any value... Users/PaulAmnuaypayoat

20050218 09:15:41 nbsp Paul you squish.  You are just bitter because Focus whooped you when you ran for senate.  I hope this doesnt come back to bite DCR in the butt, but more power to em. Users/ChrisMays

20050218 18:53:34 nbsp Haha Chris, its good that you made it to the Wiki. But seriously, there is much evidence that they have cheated over these years, and tell me, what have they given to you or Davis College Republicans DCR? Most of their slate members are against conservatives, and against fiscal responsability. I might say that voting for Student Focus  communism Users/PaulAmnuaypayoat

20050218 20:45:24 nbsp Funny thing that Student Focus ran and got elected a conservative last year, guess they are all really against conservatives eh?  It is also well known that Focus isnt the only group to do illegal campaigning  cheating, oh well. Users/JennyMays

20050218 23:04:14 nbsp I think the women of Davis College Republicans DCR had a stronger attraction to the LEAD presidential canidates than the males of DCR had to the Focus presidential canidates Users/JamesSchwab

20050218 23:13:17 nbsp Re Jenny: Thats why I said most. Yes, they almost always have a token conservative, but cmon. They have no problem allowing people such as Cari Ham to not only harass conservatives, they reward her with positions of ASUCD Senate influence. They promised the Pledge and the Flag for you guys years ago, something that costs no money, yet have failed to deliver after having all of that control of the government over all of those years. While supporting them at one time in the past may have been the better choice for many of you, it is now time to capitalize on the fact that you can elect your own people, such as Keith Shively. Student Focus is dead, especially after losing the hearts and minds of Davis College Republicans DCR, which is still one of the most influential political groups on campus. I am not saying that LEAD or anyone else is perfect, but at this point any defense of Student Focus is futile. Users/PaulAmnuaypayoat

20050219 23:58:18 nbsp Kalens veto of the 2/3 spending referendum bill is an example of where the LEAD slate is the more conservative, fiscally. Another useful fact: the LEAD slate is running Andrew Wong, who is in ROTC (which means conservative, even though hes not an open conservative). Face it  since the IRS, LEADs dangerous stuff isnt all that dangerous. Users/KenBloom

20050220 00:00:29 nbsp Users/JamesSchwab: The attractiveness of the candidates is the subject of a running debate between the men and women of DCR. They agree to disagree, and so we based our nomination on other factors. Users/KenBloom

I just wanted to say how pleasing it is to see DCR rescind its support for Student Focus, which is the cadre party of student government that serves as the perfect example of how the movie Election and reality collide at UCDavis.  I applaud the criticism at the abuse of the veto, which Kalen Gallagher and Student Focus seems to have repeatedely used to further ruin the state of the organization that we all pay to support.  Repubican or Democrat, Libertarian or Green, no one want their government to waste their money and no rational person could continue to support those that do just that. Users/RobertBaron

20070216 19:35:55 nbsp Oh how the winds of time change the picture.  The majority of current LEAD candidates sadly do not support the ROTC according to DCR.  They didnt like that. Users/GregWebb

20070217 02:09:50 nbsp LOL. Greg is uninformed. Users/JamesSchwab


text:
MODULE double_complex_erf

! Double precision complex error function
! Adapted from the Naval Surface Warfare Center Mathematics Library
! by Alan.Miller @ vic.cmis.csiro.au
! http://www.ozemail.com.au/~milleraj

! Latest revision - 18 September 2000

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)

CONTAINS

FUNCTION dnorm(x, y) RESULT(fn_val)
! Replaces the statement function anorm in the F77 code.

REAL (dp), INTENT(IN) :: x, y
REAL (dp)             :: fn_val

fn_val = MAX(ABS(x), ABS(y))
RETURN
END FUNCTION dnorm

SUBROUTINE dcrec (x, y, u, v)
!-----------------------------------------------------------------------
!             COMPLEX RECIPROCAL U + I*V = 1/(X + I*Y)
!-----------------------------------------------------------------------

REAL (dp), INTENT(IN)   :: x, y
REAL (dp), INTENT(OUT)  :: u, v

REAL (dp) :: d, t

IF (ABS(x) > ABS(y)) GO TO 10
t = x/y
d = y + t*x
u = t/d
v = -1._dp/d
RETURN

10 t = y/x
d = x + t*y
u = 1._dp/d
v = -t/d

RETURN
END SUBROUTINE dcrec

SUBROUTINE cdivid (ar, ai, br, bi, cr, ci)
!-----------------------------------------------------------------------
!     REAL (dp) COMPLEX DIVISION C = A/B AVOIDING OVERFLOW
!-----------------------------------------------------------------------

REAL (dp), INTENT(IN)   :: ar, ai, br, bi
REAL (dp), INTENT(OUT)  :: cr, ci

REAL (dp) :: d, t, u, v

IF (ABS(br) <= ABS(bi)) GO TO 10
t = bi/br
d = br + t*bi
u = (ar + ai*t)/d
v = (ai - ar*t)/d
cr = u
ci = v
RETURN

10 IF (bi == 0._dp) GO TO 20
t = br/bi
d = bi + t*br
u = (ar*t + ai)/d
v = (ai*t - ar)/d
cr = u
ci = v
RETURN

!     DIVISION BY ZERO. C = INFINITY

20 cr = HUGE(1.0_dp)
ci = cr

RETURN
END SUBROUTINE cdivid

SUBROUTINE dcerf (mo, z, w)
!-----------------------------------------------------------------------

!               COMPUTATION OF THE COMPLEX ERROR FUNCTION

!                          -----------------

!                       W = ERF(Z)    IF MO = 0
!                       W = ERFC(Z)   OTHERWISE

!-----------------------------------------------------------------------

INTEGER, INTENT(IN)     :: mo
REAL (dp), INTENT(IN)   :: z(2)
REAL (dp), INTENT(OUT)  :: w(2)

REAL (dp) :: m, n, n2, n4, np1
REAL (dp) :: c2, d, d2, e, eps, r, sn, tol, x, y
REAL (dp) :: a0(2), an(2), b0(2), bn(2)
REAL (dp) :: g0(2), gn(2), h0(2), hn(2)
REAL (dp) :: qf(2), sm(2), sz(2), tm(2), ts(2), w0(2), wn(2)
!------------------------
!     C = 1/SQRT(PI)
!------------------------
REAL (dp), PARAMETER :: c = .56418958354775628694807945156077_dp
!------------------------

!     ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE
!            SMALLEST NUMBER SUCH THAT 1._dp + EPS > 1._dp.

eps = EPSILON(1.0_dp)

!------------------------
x = z(1)
y = z(2)
sn = 1._dp
IF (x >= 0._dp) GO TO 10
x = -x
y = -y
sn = -1._dp

10 r = x*x + y*y
sz(1) = x*x - y*y
sz(2) = 2._dp*x*y

IF (r <= 1._dp) GO TO 20
IF (r >= 144._dp) GO TO 100
IF (ABS(y) > 31.8_dp*x) GO TO 50
IF (ABS(y) > 7.0_dp*x .AND. r < 64._dp) GO TO 50
IF (ABS(y) > 3.2_dp*x .AND. r < 49._dp) GO TO 50
IF (ABS(y) > 2.0_dp*x .AND. r < 36._dp) GO TO 50
IF (ABS(y) > 1.2_dp*x .AND. r < 25._dp) GO TO 50
IF (ABS(y) > 0.9_dp*x .AND. r < 16._dp) GO TO 50
IF (r >= 6.25_dp) GO TO 80
IF (ABS(y) > 0.6_dp*x) GO TO 50
IF (r >= 4.0_dp) GO TO 40

d = x - 2._dp
IF (d*d + y*y < 1._dp) GO TO 40
GO TO 50

!                          TAYLOR SERIES

20 c2 = c + c
tm(1) = c2*x
tm(2) = c2*y
sm(1) = tm(1)
sm(2) = tm(2)
tol = 2._dp*eps
m = 0._dp

21 m = m + 1._dp
d = m + m + 1._dp
ts(1) = tm(1)*sz(1) - tm(2)*sz(2)
ts(2) = tm(1)*sz(2) + tm(2)*sz(1)
tm(1) = -ts(1)/m
tm(2) = -ts(2)/m
ts(1) = tm(1)/d
ts(2) = tm(2)/d
sm(1) = sm(1) + ts(1)
sm(2) = sm(2) + ts(2)
IF (dnorm(ts(1),ts(2)) > tol*dnorm(sm(1),sm(2))) GO TO 21

IF (mo /= 0) GO TO 30
w(1) = sn*sm(1)
w(2) = sn*sm(2)
RETURN

30 IF (sn == 1._dp) GO TO 31
w(1) = 1._dp + sm(1)
w(2) = sm(2)
RETURN

31 w(1) = 0.5_dp + (0.5_dp - sm(1))
w(2) = -sm(2)
RETURN

!                  TAYLOR SERIES AROUND Z0 = 2

40 tm(1) = x
tm(2) = y
CALL erfcm2 (0, tm, w)
IF (mo /= 0) GO TO 41
w(1) = sn*(0.5_dp + (0.5_dp - w(1)))
w(2) = - sn*w(2)
RETURN
41 IF (sn > 0._dp) RETURN
w(1) = 2._dp - w(1)
w(2) = - w(2)
RETURN

!            PADE APPROXIMATION FOR THE TAYLOR SERIES
!                    FOR  (EXP(Z*Z)/Z)*ERF(Z)

50 d = 4._dp
IF (r > 16._dp) d = 16._dp
IF (r > 64._dp) d = 64._dp
d2 = d*d
CALL dcrec (sz(1), sz(2), w(1), w(2))
a0(1) = 1._dp
a0(2) = 0._dp
an(1) = (w(1) + 4._dp/15._dp)*d
an(2) = w(2)*d
b0(1) = 1._dp
b0(2) = 0._dp
bn(1) = (w(1) - 0.4_dp)*d
bn(2) = w(2)*d
CALL cdivid (an(1), an(2), bn(1), bn(2), wn(1), wn(2))
tol = 10._dp*eps
n4 = 0._dp

60 n4 = n4 + 4._dp
e = (n4 + 1._dp)*(n4 + 5._dp)
tm(1) = d*(w(1) - 2._dp/e)
tm(2) = d*w(2)
e = d2*(n4*(n4 + 2.0))/((n4 - 1.0)*(n4 + 3.0)*(n4 + 1.0)**2)

qf(1) = (tm(1)*an(1) - tm(2)*an(2)) + e*a0(1)
qf(2) = (tm(1)*an(2) + tm(2)*an(1)) + e*a0(2)
a0(1) = an(1)
a0(2) = an(2)
an(1) = qf(1)
an(2) = qf(2)
qf(1) = (tm(1)*bn(1) - tm(2)*bn(2)) + e*b0(1)
qf(2) = (tm(1)*bn(2) + tm(2)*bn(1)) + e*b0(2)
b0(1) = bn(1)
b0(2) = bn(2)
bn(1) = qf(1)
bn(2) = qf(2)

w0(1) = wn(1)
w0(2) = wn(2)
CALL cdivid (an(1), an(2), bn(1), bn(2), wn(1), wn(2))
IF (dnorm(wn(1) - w0(1), wn(2) - w0(2)) > tol*dnorm(wn(1), wn(2))) GO TO 60

c2 = c + c
sm(1) = c2*(x*wn(1) - y*wn(2))
sm(2) = c2*(x*wn(2) + y*wn(1))
e = EXP(-sz(1))
qf(1) = e*COS(-sz(2))
qf(2) = e*SIN(-sz(2))
tm(1) = qf(1)*sm(1) - qf(2)*sm(2)
tm(2) = qf(1)*sm(2) + qf(2)*sm(1)

w(1) = sn*tm(1)
w(2) = sn*tm(2)
IF (mo == 0) RETURN
w(1) = 1._dp - w(1)
w(2) = - w(2)
RETURN

!         PADE APPROXIMATION FOR THE ASYMPTOTIC EXPANSION
!                    FOR  Z*EXP(Z*Z)*ERFC(Z)

80 d = 4._dp*r
IF (r < 16._dp) d = 16._dp*r
d2 = d*d
tm(1) = sz(1) + sz(1)
tm(2) = sz(2) + sz(2)
g0(1) = 1._dp
g0(2) = 0._dp
gn(1) = (2._dp + tm(1))/d
gn(2) = tm(2)/d
h0(1) = 1._dp
h0(2) = 0._dp
tm(1) = 3._dp + tm(1)
hn(1) = tm(1)/d
hn(2) = tm(2)/d
CALL cdivid (gn(1), gn(2), hn(1), hn(2), wn(1), wn(2))
np1 = 1._dp
tol = 10._dp*eps

90 n = np1
np1 = n + 1._dp
n2 = n + n
e = (n2*(n2 + 1._dp))/d2
tm(1) = tm(1) + 4._dp
qf(1) = (tm(1)*gn(1) - tm(2)*gn(2))/d - e*g0(1)
qf(2) = (tm(1)*gn(2) + tm(2)*gn(1))/d - e*g0(2)
g0(1) = gn(1)
g0(2) = gn(2)
gn(1) = qf(1)
gn(2) = qf(2)
qf(1) = (tm(1)*hn(1) - tm(2)*hn(2))/d - e*h0(1)
qf(2) = (tm(1)*hn(2) + tm(2)*hn(1))/d - e*h0(2)
h0(1) = hn(1)
h0(2) = hn(2)
hn(1) = qf(1)
hn(2) = qf(2)

w0(1) = wn(1)
w0(2) = wn(2)
CALL cdivid (gn(1), gn(2), hn(1), hn(2), wn(1), wn(2))
IF (dnorm(wn(1) - w0(1), wn(2) - w0(2)) > tol*dnorm(wn(1), wn(2))) GO TO 90

tm(1) = x*hn(1) - y*hn(2)
tm(2) = x*hn(2) + y*hn(1)
CALL cdivid (c*gn(1), c*gn(2), tm(1), tm(2), sm(1), sm(2))
GO TO 130

!                      ASYMPTOTIC EXPANSION

100 CALL dcrec (x, y, tm(1), tm(2))
sm(1) = tm(1)
sm(2) = tm(2)
qf(1) = tm(1)*tm(1) - tm(2)*tm(2)
qf(2) = 2._dp*tm(1)*tm(2)
tol = 2._dp*eps
d = -0.5_dp
110 d = d + 1._dp
ts(1) = tm(1)*qf(1) - tm(2)*qf(2)
ts(2) = tm(1)*qf(2) + tm(2)*qf(1)
tm(1) = -d*ts(1)
tm(2) = -d*ts(2)
sm(1) = sm(1) + tm(1)
sm(2) = sm(2) + tm(2)
IF (dnorm(tm(1), tm(2)) > tol*dnorm(sm(1), sm(2))) GO TO 110
sm(1) = c*sm(1)
sm(2) = c*sm(2)
IF (x < 1.d-2) GO TO 200

!                       TERMINATION

130 e = EXP(-sz(1))
qf(1) = e*COS(-sz(2))
qf(2) = e*SIN(-sz(2))
ts(1) = qf(1)*sm(1) - qf(2)*sm(2)
ts(2) = qf(1)*sm(2) + qf(2)*sm(1)
sm(1) = ts(1)
sm(2) = ts(2)

IF (mo /= 0) GO TO 140
w(1) = sn*(0.5_dp + (0.5_dp - sm(1)))
w(2) = - sn*sm(2)
RETURN

140 IF (sn == 1._dp) GO TO 141
w(1) = 2._dp - sm(1)
w(2) = -sm(2)
RETURN

141 w(1) = sm(1)
w(2) = sm(2)
RETURN

!               MODIFIED ASYMPTOTIC EXPANSION

200 e = EXP(-sz(1))
qf(1) = e*COS(-sz(2))
qf(2) = e*SIN(-sz(2))
w(1) = qf(1)*sm(1) - qf(2)*sm(2)
w(2) = qf(1)*sm(2) + qf(2)*sm(1)
IF (mo == 0) GO TO 210
w(1) = 1._dp + sn*w(1)
w(2) = sn*w(2)
RETURN

210 IF (sn < 0.0) RETURN
w(1) = - w(1)
w(2) = - w(2)

RETURN
END SUBROUTINE dcerf

SUBROUTINE erfcm2 (mo, z, w)
!-----------------------------------------------------------------------
!           CALCULATION OF ERFC(Z) USING THE TAYLOR SERIES
!                          AROUND Z0 = 2
!-----------------------------------------------------------------------

INTEGER, INTENT(IN)     :: mo
REAL (dp), INTENT(IN)   :: z(2)
REAL (dp), INTENT(OUT)  :: w(2)

REAL (dp) :: eps, h(2), t(2), tol, x, y
INTEGER   :: j, n
!------------------------------
!     C = (2/SQRT(PI))*EXP(-4)
!     E = ERFC(2)
!------------------------------
REAL (dp), PARAMETER :: c = .020666985354092053857068941306585476_dp,  &
                        e = .0046777349810472658379307436327470714_dp
!------------------------------
REAL (dp) :: a(63) = (/ .20000000000000000000000000000000000D+01,  &
   .23333333333333333333333333333333333D+01,  .16666666666666666666666666666666667D+01,  &
   .63333333333333333333333333333333333D+00, -.22222222222222222222222222222222222D-01,  &
  -.16349206349206349206349206349206349D+00, -.76984126984126984126984126984126984D-01,  &
  -.24250440917107583774250440917107584D-02,  .12716049382716049382716049382716049D-01,  &
   .50208433541766875100208433541766875D-02, -.25305969750414194858639303083747528D-03,  &
  -.78593217482106370995259884148773038D-03, -.19118154038788959423880058800693721D-03,  &
   .46324144207742091339974937858535742D-04,  .33885549097189308829520469732109944D-04,  &
   .28637897646612243562134629672756034D-05, -.29071891082127275370004560446169188D-05,  &
  -.89674405786490646425523560263096103D-06,  .96069103941908684338469767911200105D-07,  &
   .99432863129093191401848891268744113D-07,  .97610310501460621303387795457283579D-08,  &
  -.65557500375673133822289344530892436D-08, -.18706782059105426900361744016236561D-08,  &
   .20329898993447386223176373714372370D-09,  .16941915827254374668448114614201210D-09,  &
   .10619149520827430973786114446699534D-10, -.10136148256511788733365237088810952D-10,  &
  -.21042890133669970575386166675721692D-11,  .37186985840699828780916522245407087D-12,  &
   .17921843632701679986488128324051002D-12, -.89823991804248069863542565948598397D-16,  &
  -.10533182313660970970232171410372199D-13, -.12340742690978398320850088252659714D-14,  &
   .44315624546581333350568244777175883D-15,  .11584041639989442481950487524296214D-15,  &
  -.10765703619385988116658460442219647D-16, -.70653158723054941879586082239984222D-17,  &
  -.18708903154917138727191793341667090D-18,  .32549879318817103966053527398133297D-18,  &
   .40654116689599228385911733319215613D-19, -.11250074516817311101947327325293424D-19,  &
  -.28923865378584966737386008432031980D-20,  .23653053641701517160704870522922706D-21,  &
   .14665384680061888088099002254334292D-21,  .26971039707314316218154193225264469D-23,  &
  -.58753834789274356433279671015522650D-23, -.59960357240498652932299485494869633D-24,  &
   .18586826578121663981412155416486531D-24,  .38364131854692721721867481914852428D-25,  &
  -.41342210492630142578080062451711039D-26, -.17646283105274988992381528904600860D-26,  &
   .19828685934364181151988692232131608D-28,  .65592252170840353572672782446212733D-28,  &
   .40626551379996340638338449938639730D-29, -.20097984104191034713653294173834095D-29,  &
  -.28104226475997460044096389060743131D-30,  .48705319298749358709127987806547949D-31,  &
   .12664655832830787747161769929972617D-31, -.75168312488894341862391776330113688D-33,  &
  -.45760473722605993842481669806804415D-33, -.56725491529575395930156379514718000D-35,  &
   .13932664042920082608489441616061541D-34,  .10452448992516358449586503951463322D-35 /)
!------------------------------

!     ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE
!            SMALLEST NUMBER SUCH THAT 1._dp + EPS .GT. 1._dp

eps = EPSILON(1.0_dp)

!------------------------------
tol = eps*1.d+12
h(1) = 1._dp + (1._dp - z(1))
h(2) = - z(2)

x = 1._dp
y = 0._dp
w(1) = a(30)
w(2) = 0._dp
DO n = 31,63
  t(1) = x*h(1) - y*h(2)
  t(2) = x*h(2) + y*h(1)
  x = t(1)
  y = t(2)
  t(1) = a(n)*x
  t(2) = a(n)*y
  w(1) = w(1) + t(1)
  w(2) = w(2) + t(2)
  IF (dnorm(t(1), t(2)) <= tol*dnorm(w(1), w(2))) EXIT
END DO

DO j = 1,29
  n = 30 - j
  x = h(1)*w(1) - h(2)*w(2)
  w(2) = h(1)*w(2) + h(2)*w(1)
  w(1) = a(n) + x
END DO
x = h(1)*w(1) - h(2)*w(2)
w(2) = h(1)*w(2) + h(2)*w(1)
w(1) = 1._dp + x

x = c*(h(1)*w(1) - h(2)*w(2))
w(2) = c*(h(1)*w(2) + h(2)*w(1))
w(1) = e + x
IF (mo == 0) RETURN

!                     COMPUTE EXP(Z*Z)*ERFC(Z)

x = z(1)*z(1) - z(2)*z(2)
y = 2._dp*z(1)*z(2)
x = EXP(x)
t(1) = x*COS(y)
t(2) = x*SIN(y)
x = t(1)*w(1) - t(2)*w(2)
y = t(1)*w(2) + t(2)*w(1)
w(1) = x
w(2) = y

RETURN
END SUBROUTINE erfcm2

END MODULE double_complex_erf

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
C    @(#)lcskin.f	20.3 2/13/96
      subroutine lcskin ( s, rr, rf, xf , freq )
      implicit real * 8 (a-h, o-z), integer * 4 (i-n)

      include 'ipfinc/deck25.inc'

      s2 = s * s
      s3 = ( unity - s2 ) * rr
      r2 = freq * valu8 / s3
      rf = rr
      xf = 0.0
      if ( r2 .eq. 0.0 ) go to 270
      qremb = 0.0
      if ( s .lt. tenm6 ) go to 130
      q2 = r2 * s2
      if ( s2 .lt. 0.8 .or. q2 .gt. 64.0 .or. r2 .le. 64.0 ) go to 110
  110 if ( q2 .gt. 64.0 ) qremb = dsqrt ( q2 ) * sqrt2
      xx = dsqrt ( q2 )
      x2 = xx * xx / 64.0
      iback = 2
      if ( x2 .le. unity ) go to 170
      go to 220
  120 a = -berd
      b = -beid
      aremb = gerd
      bremb = geid
  130 xx = dsqrt( r2 )
      x2 = xx * xx / 64.0
      iback = 1
      if ( x2 .le. unity ) go to 170
      go to 220
  140 g = ber
      h = bei
      e = berd
      f = beid
      if ( s .lt. tenm6 ) go to 150
      g = a * ger - b * gei + aremb * ber - bremb *bei
      h = a * gei + b * ger + aremb * bei + bremb * ber
      e = a * gerd - b * geid + aremb * berd - bremb * beid
      f = a * geid + b * gerd + aremb * beid + bremb * berd
  150 e2f2 = e * *2 + f * *2
      s2 = xx * s3 * onehaf / e2f2
      if ( iprsup .ge. 1 ) write (lunit6,160 ) e, f, e2f2, s2
  160 format ( ' "LCSKIN".   E, F, E2F2, S2 =',  4e25.16 )
      rf = ( -h * e + g * f ) * s2
      xf = ( g * e + h * f ) * s2
      go to 270
  170 z = x2
      ber = unity
      bei = 0.0
      berd = 0.0
      beid = onehaf
      gerd = 0.0
      geid = valu9
      ger = - valu10
      gei = 0.0
      ialt = 1
      do 200 k = 1, 14
      if ( ialt .eq. 1 ) go to 180
      ber = ber + fbe(k) * z
      beid = beid + fbed(k) * z
      if ( s .lt. tenm6 ) go to 190
      geid = geid + fked(k) * z
      if ( iback .eq. 2 ) go to 190
      ger = ger + fke(k) * z
      go to 190
  180 bei = bei + fbe(k) * z
      berd = berd + fbed(k) * z
      if ( s .lt. tenm6 ) go to 190
      gerd = gerd + fked(k) * z
      if ( iback .eq. 2 ) go to 190
      gei = gei + fke(k) * z
  190 z = z * x2
  200 ialt = -ialt
      beid = beid * xx
      berd = berd * xx
      if ( s .lt. tenm6 ) go to 210
      xl = dlog( xx * onehaf )
      gerd = -xl * berd - ber / xx + beid * aaa1 + xx * gerd
      geid = -xl * beid - bei / xx - berd * aaa1 + xx * geid
      if ( iback .eq. 2 ) go to 210
      ger = -xl * ber + bei * aaa1 + ger
      gei = -xl * bei - ber * aaa1 + gei
  210 go to (140 ,120 ), iback
  220 x2 = 8.0 / xx
      z = x2
      ber = 0.0
      bei = -valu11
      berd = ber
      beid = bei
      ger = unity / sqrt2
      gei = ger
      gerd = ger
      geid = gei
      ialt = 1
      do 250 k = 1, 6
      thetar = fbe(k+14) * z
      thetai = fbed(k+14) * z
      phir = fke(k+14) * z
      phii = fked(k+14) * z
      ber = ber + thetar
      bei = bei + thetai
      ger = ger + phir
      gei = gei + phii
      if ( ialt .eq. 1 ) go to 230
      berd = berd + thetar
      beid = beid + thetai
      gerd = gerd + phir
      geid = geid + phii
      go to 240
  230 berd = berd - thetar
      beid = beid - thetai
      gerd = gerd - phir
      geid = geid - phii
  240 ialt = -ialt
  250 z = z * x2
      xl = xx * sqrt2
      if ( qremb .lt. 1.0 ) go to 260
      xl = xl - qremb
  260 thetar = -xl + berd
      thetai = -xl + beid
      z = dsqrt ( xx )
      x2 = valu12 / z
      z = valu13 / z * dexp ( thetar )
      fr = z * dcos ( thetai )
      fi = z * dsin ( thetai )
      x2 = x2 * dexp ( ber )
      thetar = x2 * dcos ( bei )
      thetai = x2 * dsin ( bei )
      z = -fr * gerd + fi * geid
      geid = -fr * geid - fi * gerd
      gerd = z
      z = aaa2 * dexp ( -qremb )
      gr = z * dsin ( qremb )
      gi = z * dcos ( qremb )
      berd = thetar * ger - thetai * gei + gerd * gr - geid * gi
      beid = thetar * gei + thetai * ger + gerd * gi + geid * gr
      ger = fr
      gei = fi
      ber = thetar + ger * gr - gei * gi
      bei = thetai + ger * gi + gei * gr
      go to (140 ,120 ), iback
  270 if ( iprsup .ge. 1 ) write (lunit6,280 ) z, thetar, thetai
  280 format ( 36h exit "lcskin".  z, thetar, thetai =,  3e25.16 )
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp