!vCIRA2019_0.2.11
subroutine estimate_area_averaged_sst(aakm,sst2d,flon_west,flat_south,fdlon,fdlat,n_lons_max,n_lats_max,&
            n_lons_dat,n_lats_dat,sst,tlon,tlat,izp,sta,ierr_nsta,&
            lu_log_info,lu_log_debug)
        
    !this subroutine calculates area-averaged SST from 5 points
    ! converted from Mark DeMaria's code

    implicit none

    !input ------------------------------------------------
    ! Distance (km) around the center for area-averaged SST
    real :: aakm
    ! point lat 
    real :: tlat
    ! point lon 
    real :: tlon
    ! sst  interpolated for tlat, tlon
    real :: sst
    !max number of lats
    integer :: n_lons_max
    !max number of lons
    integer :: n_lats_max
    !number of lons
    integer :: n_lons_dat
    !number of lats
    integer :: n_lats_dat
    ! lat dist 
    real :: fdlat
    ! lon dist 
    real :: fdlon
    ! left most lon point 
    real :: flon_west
    ! right most lat point 
    real :: flat_south
    ! array sst 
    real :: sst2d(n_lons_max,n_lats_max)
    !  Zonal Periodic flag:
    !       =1 when f1 is periodic in the zonal direction
    !       (normally used only when f1 spans 360 deg
    !       of longitude, starting at 0)
    !       =0 if not periodic
    ! should be izp = 0 for atl, pac  EOHC data 
    integer :: izp 
    ! logicl unit number for log file
    integer :: lu_log_info
    ! logicl unit number for debug file
    integer :: lu_log_debug

    !output -------------------------------------------------------
    ! area averaged sst 
    real :: sta
    ! error getting area averaged NSTA
    ! this error should be used as warning only
    ! if non-zero SST rather than area-averaged  SST is used
    integer :: ierr_nsta
    
    !local ----------------------------------------------------
    real,parameter :: dtr = 3.14159/180.0
    ! avaraging distance along lat in km
    real :: dlataa
    ! avaraging distance along lon in km 
    real :: dlonaa
   
    ! lats/lons of points to use for averaging 
    real :: tlonp
    real :: tlonm
    real :: tlatp
    real :: tlatm
    ! missing value returned by llintsp
    real :: ierr

    ! temp sst
    real :: sstp10
    real :: sstm10
    real :: sstp01
    real :: sstm01

    write(lu_log_debug,*) 'INFO: estimate_area_averaged_sst: ','will now get NSTA for aakm  = ', aakm
    
!   Area averaged SST
    if (aakm .gt. 0.0) then
        !dtr    = 3.14159/180.0
        dlataa = aakm/111.1
        dlonaa = dlataa/cos(abs(dtr*tlat))
        !
        tlonp = tlon + dlonaa
        tlonm = tlon - dlonaa
        tlatp = tlat + dlataa
        tlatm = tlat - dlataa
!
        call llintsp(sst2d,flon_west,flat_south,fdlon,fdlat,&
                n_lons_max,n_lats_max,n_lons_dat,n_lats_dat,sstp10,tlonp,tlat  ,izp,ierr)
        if (ierr .ne. 0) sstp10 = sst
        ierr_nsta = ierr
!
        call llintsp(sst2d,flon_west,flat_south,fdlon,fdlat,&
                n_lons_max,n_lats_max,n_lons_dat,n_lats_dat,sstm10,tlonm,tlat  ,izp,ierr)
        if (ierr .ne. 0) sstm10 = sst
        ierr_nsta = ierr_nsta+ierr
!
        call llintsp(sst2d,flon_west,flat_south,fdlon,fdlat,&
                n_lons_max,n_lats_max,n_lons_dat,n_lats_dat,sstp01,tlon ,tlatp, izp,ierr)
        if (ierr .ne. 0) sstp01 = sst
        ierr_nsta = ierr_nsta+ierr
!
        call llintsp(sst2d,flon_west,flat_south,fdlon,fdlat,&
                n_lons_max,n_lats_max,n_lons_dat,n_lats_dat,sstm01,tlon ,tlatm, izp,ierr)
        if (ierr .ne. 0) sstm01 = sst
        ierr_nsta = ierr_nsta+ierr
!
        sta = 0.2*(sst + sstp10 + sstm10 + sstp01 + sstm01)
        write(lu_log_debug,*)  'INFO: estimate_area_averaged_sst: &
                &NSTA sub estimate sst, nsta, ierr_nsta ', sst, sta, ierr_nsta
        !write(lu_log_debug,*)  'DEBUG: estimate_area_averaged_sst: &
        !        &NSTA sub estimate sst, nsta ', sst, sta, sstp10,sstm10,sstp01,sstm01
    else
        sta = sst
        ierr_nsta = 4
        write(lu_log_debug,*)  'INFO: estimate_area_averaged_sst: &
                &NSTA is equal to SST, aakm = 0, sst, nsta ', sst, sta
    endif

    return

end subroutine estimate_area_averaged_sst