    module surface
! SURFACE - support code for reading the MRF surface analysis files
! and for computing means of the relevant quantities.
! SYNOPSIS
!    use surface
!
! DEPENDENCIES
!    fixio - in current directory
!    resolution, precision, io, ferror, command_line_support - in support dir
!
! MODIFICATION HISTORY
! 08/11/05 - Changed name of module to "surface" and gave a matching filename,
! 07/21/04 - First implementation by Eric to compute and write an averaged
!            surface analysis file for the MRF code.
! 07/26/04 - Istvan made a major reduction (simplification) of the code by
!            eliminating several subroutines and variables, and
!            shortening the parameter list of a few subroutines. The
!            eliminated structures were mainly associated with the
!            mpi implementation of the original MRF code. The most extensive
!            changes were done in fixi.f90. 
!
    use fixio
    use resolution
    use precision
    use io
    use ferror
    use command_line_support,only:program_exit
    implicit none
!
!  On Lahey Intel Linux, input and output binary data require byte-swapping.
!
    integer,save,private::nread=INSWAP
    integer,save,private::nwrite=OUTSWAP
!
    type sfcanl
       real(D)::smc(LONF,LATG,LSOIL)
       real(D)::stc(LONF,LATG,LSOIL)
       real(D)::tsea(LONF,LATG)
       real(D)::sheleg(LONF,LATG)
       real(D)::tg3(LONF,LATG)
       real(D)::zorl(LONF,LATG)
       real(D)::cv(LONF,LATG)
       real(D)::cvb(LONF,LATG)
       real(D)::cvt(LONF,LATG)
       real(D)::alvsf(LONF,LATG)
       real(D)::alvwf(LONF,LATG)
       real(D)::alnsf(LONF,LATG)
       real(D)::alnwf(LONF,LATG)
       real(D)::slmsk(LONF,LATG)
       real(D)::vfrac(LONF,LATG)
       real(D)::canopy(LONF,LATG)
       real(D)::f10m(LONF,LATG)
       real(D)::vtype(LONF,LATG)
       real(D)::stype(LONF,LATG)
       real(D)::facsf(LONF,LATG)
       real(D)::facwf(LONF,LATG)
       real(D)::uustar(LONF,LATG)
       real(D)::ffmm(LONF,LATG)
       real(D)::ffhh(LONF,LATG)
    end type sfcanl
!
!  We assume no parallel processing (1 cpu) 
!
!  No explicit routines are needed for assignment, as the default Fortran
!  operator suffices for sfcanl.
!
    interface operator(+)
       module procedure add_sfcanl
    end interface
    interface operator(/)
       module procedure div_sfcanl_n
    end interface
    contains
!------------------------------
    subroutine read_sfcanl(fn,sfc,fhour,idate)
    implicit none
    character(*),intent(in)::fn
    type(sfcanl),intent(out)::sfc
    real(R),intent(out)::fhour
    integer,intent(out)::idate(4)  ! hh,mm,dd,yyyy
!
    integer::ioerr
!
    open(unit=nread,file=fn,form='unformatted',position='rewind', &
       status='old',iostat=ioerr)
    if(ioerr.ne.0) then
       write(IERR,900) trim(fn),nread
       call program_exit(1)
    endif
900 format('read_sfcanl: cannot open ',a/ &
       'On Lahey, setting the environment variable "FORT90L=-Wl,-T',i0, &
       '"is required if the input files are byte-swapped')
!
    call para_fixio_r(LONF,LONF2,LATG,LATG2,sfc%tsea,sfc%smc,sfc%sheleg, &
       sfc%stc,sfc%tg3, &
       sfc%zorl,sfc%cv,sfc%cvb,sfc%cvt,sfc%alvsf,sfc%alvwf,sfc%alnsf,sfc%alnwf,&
       sfc%slmsk,sfc%vfrac,sfc%canopy,sfc%f10m,sfc%vtype,sfc%stype,sfc%facsf, &
       sfc%facwf,sfc%uustar,sfc%ffmm,sfc%ffhh,nread, &
       fhour,idate)
    close(nread)
    return
    end subroutine read_sfcanl
!------------------------------
    subroutine write_sfcanl(fn,sfc,fhour,cdate)
    implicit none
    character(*),intent(in)::fn
    type(sfcanl),intent(in)::sfc
    real(R),intent(in)::fhour
    character(*),intent(in)::cdate  ! yyyyddmmhh format
!
    integer::ioerr,idate(4)
!
    open(unit=nwrite,file=fn,form='unformatted',position='rewind',iostat=ioerr)
    if(ioerr.ne.0) then
       write(IERR,900) trim(fn)
       call program_exit(1)
    endif
900 format('write_sfcanl: cannot create ',a/a)
!
!    nlats=LATG2
!
!  The input date string has the format yyyymmddhh, but the date
!  format in the MRF sigma files is hh mm dd yyyy.
!
    if(len_trim(cdate).eq.10) then
       read(cdate,'(i4,3i2)') idate(4),idate(2),idate(3),idate(1)
    else
       idate=0
    endif

    call para_fixio_w(LONF,LONF2,LATG,LATG2,sfc%tsea,sfc%smc,sfc%sheleg, &
       sfc%stc,sfc%tg3, &
       sfc%zorl,sfc%cv,sfc%cvb,sfc%cvt,sfc%alvsf,sfc%alvwf,sfc%alnsf,sfc%alnwf,&
       sfc%slmsk,sfc%vfrac,sfc%canopy,sfc%f10m,sfc%vtype,sfc%stype,sfc%facsf, &
       sfc%facwf,sfc%uustar,sfc%ffmm,sfc%ffhh,nwrite, &
       fhour,idate)
    close(nwrite)
    return
    end subroutine write_sfcanl
!------------------------------
    elemental function add_sfcanl(x,y) result(z)
    implicit none
    type(sfcanl),intent(in)::x,y
    type(sfcanl)::z
!
!  Add all the array components of sfcanl pairwise
!
    z%smc=x%smc+y%smc
    z%stc=x%stc+y%stc
    z%tsea=x%tsea+y%tsea
    z%sheleg=x%sheleg+y%sheleg
    z%tg3=x%tg3+y%tg3
    z%zorl=x%zorl+y%zorl
    z%cv=x%cv+y%cv
    z%cvb=x%cvb+y%cvb
    z%cvt=x%cvt+y%cvt
    z%alvsf=x%alvsf+y%alvsf
    z%alvwf=x%alvwf+y%alvwf
    z%alnsf=x%alnsf+y%alnsf
    z%alnwf=x%alnwf+y%alnwf
    z%slmsk=x%slmsk+y%slmsk
    z%vfrac=x%vfrac+y%vfrac
    z%canopy=x%canopy+y%canopy
    z%f10m=x%f10m+y%f10m
    z%vtype=x%vtype+y%vtype
    z%stype=x%stype+y%stype
    z%facsf=x%facsf+y%facsf
    z%facwf=x%facwf+y%facwf
    z%uustar=x%uustar+y%uustar
    z%ffmm=x%ffmm+y%ffmm
    z%ffhh=x%ffhh+y%ffhh
    return
    end function add_sfcanl
!------------------------------
    elemental function div_sfcanl_n(x,n) result(z)
    implicit none
    type(sfcanl),intent(in)::x
    integer,intent(in)::n
    type(sfcanl)::z
!
!  Divide each component array element by n
!
    z%smc=x%smc/n
    z%stc=x%stc/n
    z%tsea=x%tsea/n
    z%sheleg=x%sheleg/n
    z%tg3=x%tg3/n
    z%zorl=x%zorl/n
    z%cv=x%cv/n
    z%cvb=x%cvb/n
    z%cvt=x%cvt/n
    z%alvsf=x%alvsf/n
    z%alvwf=x%alvwf/n
    z%alnsf=x%alnsf/n
    z%alnwf=x%alnwf/n
    z%slmsk=x%slmsk/n
    z%vfrac=x%vfrac/n
    z%canopy=x%canopy/n
    z%f10m=x%f10m/n
    z%vtype=x%vtype/n
    z%stype=x%stype/n
    z%facsf=x%facsf/n
    z%facwf=x%facwf/n
    z%uustar=x%uustar/n
    z%ffmm=x%ffmm/n
    z%ffhh=x%ffhh/n
    return
    end function div_sfcanl_n
!------------------------------
    end module surface
