      module lbc
!
! Calculate lower boundary conditions for T,U,V,Z.
! Namelist input options for lbc are validated by
! sub inp_lbc in input.F.
!
      use params_module,only: nlonp4,nlat,nlev,dz
      use cons_module,only: pi,atm_amu,gask,grav,freq_semidi,
     |  dt,re,dlamda,tgrad,cs,cor,tn,freq_3m3
      use mpi_module,only: lon0,lon1,lat0,lat1
      use addfld_module,only: addfld
      use init_module,only: istep
      use fields_module,only: 
     |  t_bgrd, z_bgrd,u_bgrd, v_bgrd,                          ! background
     |  u_prt, v_prt,u_prt_tn, v_prt_tn,u_prt_tnm, v_prt_tnm    ! perturbations
      implicit none
!
! Background options (and corresponding namelist input parameters):
!   ECMWF (daily, 2.5 deg only, T,U,V,Z) (ECMWF_NCFILE)
!   NCEP  (daily, 5.0 deg only, T,Z)     (NCEP_NCFILE)
!   Zatmos (TN)                          (ZATMOS_NCFILE)
!   Annual Tide (Z)                      (TIDEANN=0/1)
!   Uvbnd (UN,VN 2.5 or 5.0 deg)         (no namelist input)
!
! Perturbation options (to be added to background):
!   GSWM (mig and nonmig, 2.5 deg)       (GSWM_xx_xxx_NCFILE)
!   GSWM (mig only, 5.0 deg)             (GSWM_MI_xxx_NCFILE)
!   Hough Modes                          (TIDE, TIDE2)
!   Planetary                            (logical planetary)
!   2-day Rossby-Gravity wave (3m3)      (TIDE3M3)
!   Kelvin wave                          (logical Kelvin)
!   
! Rules (enforced by sub inp_lbc in input.F):
!   - Must provide one and only one of background zatmos, ncep, or ecmwf
!     (Addition of annual tide is allowed only with zatmos)
!   - If Zatmos is set, then annual tides must be turned on
!   - GSWM and Hough modes are mutually exclusive
!   - GSWM non-migrating perturbations not allowed at 5 deg resolution
!   - NCEP only available at 5 deg resolution
!   - ECMWF only available at 2.5 deg resolution
!   - Planet, 2-day, Kelvin, and annual tides are not allowed with NCEP or ECMWF
!
! Typical valid configurations (bracketed parameters [] are optional):
!   5.0 deg: Zatmos + Annual + uvbnd + [GSWM (mig) or Hough] + [Planet, 2-day, or Kelvin]
!   5.0 deg: NCEP  + [GSWM (mig) or Hough]
!
!   2.5 deg: Zatmos + Annual + uvbnd + [GSWM (mig and/or non-mig) or Hough] +
!                                      [Planet, 2-day, or Kelvin]
!   2.5 deg: ECMWF + [GSWM (mig and/or non-mig) or Hough]
!
      real :: dday
!
! Total lower boundary conditions returned by this module.
! (dimensioned at full global grid, but defined at subdomains only)
!
      real,dimension(nlonp4,nlat) :: t_lbc, u_lbc, v_lbc, z_lbc ! total LBC
!
! Diurnal and semi-diurnal tidal perturbations using Hough functions: 
      complex,dimension(nlat) :: 
     |  t_di , u_di , v_di , z_di, t_sdi, u_sdi, v_sdi, z_sdi
      complex,parameter :: ci=(0.,1.), expta=1.
      complex :: bnd_sdi(nlonp4), bnd_di(nlonp4)
!
! Annual tides (only zba is used):
      complex,dimension(nlat) :: tba,uba,vba,zba
      complex :: bnda(nlonp4)
!
! For bndry_comp:
      real :: b(nlonp4,2,2),fb(nlonp4,2,nlat)
!
! This t0 is different than the t0 in cons.
      real :: t0(nlev+1) 
!
! Planetary and Kelvin waves:
      real :: zbplanet(nlonp4,nlat)  ! planetary wave (sub bndry_planetary)
      real :: zbkelvin(nlonp4,nlat)  ! kelvin wave (sub bndry_kelvin)
!
! 2-day wave:
      complex,dimension(nlat) :: zb3m3,tb3m3,ub3m3,vb3m3
      complex :: bnd3m3(nlonp4)
!
! Background tn constant for Hough mode background:
      real,parameter :: tbound = 233.7 ! for Hough mode background only
!
! Temporary flag for transition from old to new lbc code:
      character(len=8),parameter :: lbc_code = 'newlbc  '
!     character(len=8),parameter :: lbc_code = 'oldlbc  '
      contains
!-----------------------------------------------------------------------
      subroutine init_lbc 
!
! Called once per run from tgcm.
!
      use cco2gr_module,only: cco2gr,set_cco2_data
      use chemrates_module,only: co2mix
      use cons_module,only: t0cons=>t0
      implicit none
!
! t0 in cons is zeroed out:
      t0cons(:) = 0.
!
! t0 local to this module (different than the t0 in cons) is for
! use in Hough functions:
      t0(:) = 0.
      t0(1) = tbound
      t0(2) = tbound+dz*tgrad
!
      call bndry_diurnal      ! t_di ,u_di ,v_di ,z_di
      call bndry_semidiurnal  ! t_sdi,u_sdi,v_sdi,z_sdi
      call set_cco2_data
      call cco2gr(co2mix)
      end subroutine init_lbc 
!-----------------------------------------------------------------------
      subroutine tuvz_lbc
!
! Update lower boundary subdomains of T,U,V,Z
! This is called every timestep from advance, after reading or calculating
!   ncep, ecmwf, zatmos, or gswm data, and before calling addiag, dt, duv, etc.
!   
! Lower boundary option combinations have been validated by
!   sub inp_lbc in input.F module.
!
      use init_module,only: istep,iday,iyear,iter,secs
      use input_module,only: step, tideann, planetary, kelvin,
     |  gswm_mi_di_ncfile,  ! gswm migrating diurnal data file
     |  gswm_mi_sdi_ncfile, ! gswm migrating semi-diurnal data file
     |  gswm_nm_di_ncfile,  ! gswm non-migrating diurnal data file
     |  gswm_nm_sdi_ncfile, ! gswm non-migrating semi-diurnal data file
     |  ecmwf_ncfile,       ! ECMWF data file
     |  ncep_ncfile,        ! NCEP data file
     |  ncep_reanalysis,    ! NCEP reanalysis data
     |  zatmos_ncfile,      ! ZATMOS data file
     |  tide3m3,            ! 2-day wave amplitude and phase
     |  calendar_advance,
     |  start_day, start_year
      use gswm_module,only: getgswm
      use gswm_module,only: ! (nlonp4,nlat)
     |  gswm_mi_di_z, gswm_mi_sdi_z, gswm_nm_di_z, gswm_nm_sdi_z,
     |  gswm_mi_di_t, gswm_mi_sdi_t, gswm_nm_di_t, gswm_nm_sdi_t,
     |  gswm_mi_di_u, gswm_mi_sdi_u, gswm_nm_di_u, gswm_nm_sdi_u,
     |  gswm_mi_di_v, gswm_mi_sdi_v, gswm_nm_di_v, gswm_nm_sdi_v
      use zatmos_module,only: zatmos_tn, zatmos_bndry
      use ecmwf_module,only: z_ecmwf, t_ecmwf, u_ecmwf, v_ecmwf, ! (nlonp4,nlat)
     |  ecmwf_bndry
      use ncep_module,only: ncep_bndry,
     |  zncep, ! (nlonp4,nlat,npr)
     |  tncep  ! (nlonp4,nlat)
      use ncep_rean_module,only: ncep_rean_bndry,
     |  z_ncep,t_ncep,u_ncep,v_ncep
      use hist_module,only: nstep,modeltime
#ifdef MPI
      use mpi_module,only: mp_periodic_f2d,
     |  mp_bndlons_f2d,mp_bndlats_f2d
#endif
      implicit none
!
! Local:
      integer :: i,j,lat
      real :: rstep
      real :: flbc(lon0-2:lon1+2,lat0-2:lat1+2,7) ! u_prt,v_prt,z_bgrd
      complex :: t_expt_sdi, t_expt_di, uvz_expt_sdi, uvz_expt_di,
     |  expt3m3
!
! Calculate exponentials for Hough functions:
      rstep = float(step)
      t_expt_sdi = cexp(ci*freq_semidi*rstep*iter)
      t_expt_di  = cexp(ci*.5*freq_semidi*rstep*iter)
      uvz_expt_sdi = cexp(ci*freq_semidi*dt*iter)
      uvz_expt_di  = cexp(ci*.5*freq_semidi*dt*iter)
      expt3m3 = cexp(ci*freq_3m3*rstep*iter)
!
! am 8/09: lbc_bgrd save un and vn perturbations at the LB
! from previous timestep
! for istep = 1 see end of the subroutine
      if(istep > 2) then
        u_prt_tnm(lon0:lon1,lat0:lat1)  = u_prt_tn(lon0:lon1,lat0:lat1)
        v_prt_tnm(lon0:lon1,lat0:lat1)  = v_prt_tn(lon0:lon1,lat0:lat1)
      elseif(istep == 2) then ! for istep = 2
        u_prt_tnm(lon0:lon1,lat0:lat1)  = u_prt(lon0:lon1,lat0:lat1)
        v_prt_tnm(lon0:lon1,lat0:lat1)  = v_prt(lon0:lon1,lat0:lat1)
      endif
      if(istep > 1) then
        u_prt_tn(lon0:lon1,lat0:lat1)  = u_prt(lon0:lon1,lat0:lat1)
        v_prt_tn(lon0:lon1,lat0:lat1)  = v_prt(lon0:lon1,lat0:lat1)
      endif
!
! Init subdomain:
      t_lbc(lon0:lon1,lat0:lat1) = 0.
      u_lbc(lon0:lon1,lat0:lat1) = 0.
      v_lbc(lon0:lon1,lat0:lat1) = 0.
      z_lbc(lon0:lon1,lat0:lat1) = 0.
!
      z_bgrd(lon0:lon1,lat0:lat1) = 0.
      t_bgrd(lon0:lon1,lat0:lat1) = 0.
      u_prt(lon0:lon1,lat0:lat1)  = 0.
      v_prt(lon0:lon1,lat0:lat1)  = 0.
!
! Build background for lower boundaries:
      dday = float(iday)+amod(float(iter)*float(step),86400.)/86400.
!
! Zatmos background for T only (redundant in longitude):
! (This is exclusive of NCEP and ECMWF, already verified by sub inp_lbc)
!
      if (len_trim(zatmos_ncfile) > 0) then
        call zatmos_bndry(istep,iday,int(secs))
        do lat=lat0,lat1
          do i=lon0,lon1
            t_lbc(i,lat) = t_lbc(i,lat) + zatmos_tn(lat)
          enddo
        enddo
        if (istep==1) write(6,"('tuvz_lbc: added zatmos_tn to t_lbc')")
      endif
!
! Annual tides are applied to background Z only:
!
      if (tideann > 0) then
        call bndry_annual(dday) ! tba,uba,vba,zba
        do lat=lat0,lat1
          do i=lon0,lon1
            z_lbc(i,lat) = z_lbc(i,lat) + real(zba(lat)*bnda(i)*expta)
          enddo
        enddo
        if (istep==1) write(6,"('tuvz_lbc: added zba to z_lbc')")
      endif
!
! ECMWF daily data (2.5 deg res):
!
      if (len_trim(ecmwf_ncfile) > 0) then 
        if (calendar_advance > 0) then
          call ecmwf_bndry(ecmwf_ncfile,istep,iyear,iday,int(secs))
        else
          call ecmwf_bndry(ecmwf_ncfile,istep,start_year,start_day,0)
        endif
        do lat=lat0,lat1
          do i=lon0,lon1
            t_lbc(i,lat) = t_lbc(i,lat) + t_ecmwf(i,lat)
            u_lbc(i,lat) = u_lbc(i,lat) + u_ecmwf(i,lat)
            v_lbc(i,lat) = v_lbc(i,lat) + v_ecmwf(i,lat)
            z_lbc(i,lat) = z_lbc(i,lat) + z_ecmwf(i,lat)
          enddo
        enddo
        if (istep==1) write(6,"('tuvz_lbc: added ecmwf ',
     |    'to t_lbc,u_lbc,v_lbc,z_lbc')")
      endif
!
! NCEP daily data (5.0 deg res), T,Z only (U,V will be calculated in uvbnd):
!
      if (len_trim(ncep_ncfile) > 0) then
        call ncep_bndry(istep,modeltime,iyear,iday,int(secs)) ! ncep.F
        do lat=lat0,lat1
          do i=lon0,lon1
            t_lbc(i,lat) = t_lbc(i,lat) + tncep(i,lat)
            z_lbc(i,lat) = z_lbc(i,lat) + zncep(i,lat,2)
          enddo
        enddo
        if (istep==1) write(6,"('tuvz_lbc: added ncep to t_lbc,z_lbc')")
      endif
!
! NCEP reanalysis data:
!
      if (len_trim(ncep_reanalysis) > 0) then
        call ncep_rean_bndry(ncep_reanalysis,istep,modeltime,
     |    iyear,iday,int(secs))
        do lat=lat0,lat1
          do i=lon0,lon1
            t_lbc(i,lat) = t_lbc(i,lat) + t_ncep(i,lat)
            u_lbc(i,lat) = u_lbc(i,lat) + u_ncep(i,lat)
            v_lbc(i,lat) = v_lbc(i,lat) + v_ncep(i,lat)
            z_lbc(i,lat) = z_lbc(i,lat) + z_ncep(i,lat)
          enddo
        enddo
        if (istep==1)
     |    write(6,"('tuvz_lbc: added ncep_reanalysis data to ',
     |      't_lbc,u_lbc,v_lbc,z_lbc')")
      endif
!
! End building lower boundary background fields
!
! Add perturbations:
! Add Planetary, 2-day wave, and Kelvin wave perturbations, 
! if requested. (bndry_planetary, bndry_2day and bndry_kelvin 
! were called by init_lbc)
!
! Planetary waves:

      if (planetary > 0) then
        call bndry_planetary  ! zbplanet
        do j=lat0,lat1
          do i=lon0,lon1
            z_lbc(i,j) = z_lbc(i,j)+zbplanet(i,j)
          enddo
        enddo
        if (istep==1) write(6,"('tuvz_lbc: added planetary to z_lbc')")
      endif
!
! Kelvin waves:
      if (kelvin > 0) then
        call bndry_kelvin     ! zbkelvin
        do j=lat0,lat1
          do i=lon0,lon1
            z_lbc(i,j) = z_lbc(i,j)+zbkelvin(i,j)
          enddo
        enddo
        if (istep==1) write(6,"('tuvz_lbc: added kelvin to z_lbc')")
      endif
!
! 2-day wave:
      if (tide3m3(1) /= 0..or.tide3m3(2) /= 0.) then
        call bndry_2day         ! tb3m3,ub3m3,vb3m3,bnd3m3
        do j=lat0,lat1
          do i=lon0,lon1
            t_lbc(i,j) = t_lbc(i,j)+real(tb3m3(j)*bnd3m3(i)*expt3m3)
            z_lbc(i,j) = z_lbc(i,j)+real(zb3m3(j)*bnd3m3(i)*expt3m3)
          enddo
        enddo
        if (istep==1) write(6,"('tuvz_lbc: added 2-day wave to ',
     |    't_lbc, z_lbc')")
      endif
!
! am 8/11/09 lbc_bgrd: save background fields for Z, Tn
!     add Planetary waves,Kelvin waves,2-day wave to the background since there is
!     no explicit neutral wind component like for GSWM and the Hough modes
!     therefore neutral winds have to calculated using the subroutine  uvbnd/uvbgrd
!
      z_bgrd(lon0:lon1,lat0:lat1) = z_lbc(lon0:lon1,lat0:lat1)
      t_bgrd(lon0:lon1,lat0:lat1) = t_lbc(lon0:lon1,lat0:lat1)
!
! Add GSWM or Hough mode perturbations:
! save Un,Vn perturbations in u_prt and v_prt
!
      call getgswm(istep,iyear,iday,secs) ! read gswm data files (istep==1)
!
! GSWM migrating diurnal:
      if (len_trim(gswm_mi_di_ncfile) > 0) then
        do j=lat0,lat1
          do i=lon0,lon1
            t_lbc(i,j) = t_lbc(i,j)+gswm_mi_di_t(i,j)
            z_lbc(i,j) = z_lbc(i,j)+gswm_mi_di_z(i,j)
            u_prt(i,j) = u_prt(i,j)+gswm_mi_di_u(i,j)
            v_prt(i,j) = v_prt(i,j)+gswm_mi_di_v(i,j)
          enddo
        enddo
        if (istep==1) write(6,"('tuvz_lbc: added gswm_mi_di to ',
     |    't_lbc,z_lbc')")
!
! Hough mode diurnal:
      else ! use Hough functions for diurnal tide
        do j=lat0,lat1
          do i=lon0,lon1
            t_lbc(i,j) = t_lbc(i,j)+real(t_di(j)*bnd_di(i)*t_expt_di)
            z_lbc(i,j) = z_lbc(i,j)+real(z_di(j)*bnd_di(i)*uvz_expt_di)
            u_prt(i,j) = u_prt(i,j)+real(u_di(j)*bnd_di(i)*uvz_expt_di)
            v_prt(i,j) = v_prt(i,j)+real(v_di(j)*bnd_di(i)*uvz_expt_di)
          enddo
        enddo
        if (istep==1) then
          write(6,"('tuvz_lbc: Added Hough mode diurnal',
     |      ' tide to t_lbc,z_lbc')")
!         write(6,"(10x,'t_lbc min,max=',2e12.4)") 
!    |      minval(t_lbc(lon0:lon1,lat0:lat1)),
!    |      maxval(t_lbc(lon0:lon1,lat0:lat1))
!         write(6,"(10x,'z_lbc min,max=',2e12.4)") 
!    |      minval(z_lbc(lon0:lon1,lat0:lat1)),
!    |      maxval(z_lbc(lon0:lon1,lat0:lat1))
        endif
      endif
!
! GSWM migrating semi-diurnal:
      if (len_trim(gswm_mi_sdi_ncfile) > 0) then
        do j=lat0,lat1
          do i=lon0,lon1
            t_lbc(i,j) = t_lbc(i,j)+gswm_mi_sdi_t(i,j)
            z_lbc(i,j) = z_lbc(i,j)+gswm_mi_sdi_z(i,j)
            u_prt(i,j) = u_prt(i,j)+gswm_mi_sdi_u(i,j)
            v_prt(i,j) = v_prt(i,j)+gswm_mi_sdi_v(i,j)
          enddo
        enddo
        if (istep==1) write(6,"('tuvz_lbc: added gswm_mi_sdi to ',
     |    't_lbc,z_lbc')")
!
! Hough mode semi-diurnal:
      else ! use Hough functions for semi-diurnal tide
        do j=lat0,lat1
          do i=lon0,lon1
            t_lbc(i,j)=t_lbc(i,j)+real(t_sdi(j)*bnd_sdi(i)*t_expt_sdi)
            z_lbc(i,j)=z_lbc(i,j)+real(z_sdi(j)*bnd_sdi(i)*uvz_expt_sdi)
            u_prt(i,j)=u_prt(i,j)+real(u_sdi(j)*bnd_sdi(i)*uvz_expt_sdi)
            v_prt(i,j)=v_prt(i,j)+real(v_sdi(j)*bnd_sdi(i)*uvz_expt_sdi)
          enddo
        enddo
        if (istep==1) then
          write(6,"('tuvz_lbc: Added Hough mode semi-diurnal',
     |      ' tide to t_lbc,z_lbc')")
!         write(6,"(10x,'t_lbc min,max=',2e12.4)") 
!    |      minval(t_lbc(lon0:lon1,lat0:lat1)),
!    |      maxval(t_lbc(lon0:lon1,lat0:lat1))
!         write(6,"(10x,'z_lbc min,max=',2e12.4)") 
!    |      minval(z_lbc(lon0:lon1,lat0:lat1)),
!    |      maxval(z_lbc(lon0:lon1,lat0:lat1))
        endif
      endif
!
! GSWM non-migrating diurnal:
      if (len_trim(gswm_nm_di_ncfile) > 0) then
        do j=lat0,lat1
          do i=lon0,lon1
            t_lbc(i,j) = t_lbc(i,j)+gswm_nm_di_t(i,j)
            z_lbc(i,j) = z_lbc(i,j)+gswm_nm_di_z(i,j)
            u_prt(i,j) = u_prt(i,j)+gswm_nm_di_u(i,j)
            v_prt(i,j) = v_prt(i,j)+gswm_nm_di_v(i,j)
          enddo
        enddo
        if (istep==1) write(6,"('tuvz_lbc: added gswm_nm_di to ',
     |    't_lbc,z_lbc')")
      endif
!
! GSWM non-migrating semi-diurnal:
      if (len_trim(gswm_nm_sdi_ncfile) > 0) then
        do j=lat0,lat1
          do i=lon0,lon1
            t_lbc(i,j) = t_lbc(i,j)+gswm_nm_sdi_t(i,j)
            z_lbc(i,j) = z_lbc(i,j)+gswm_nm_sdi_z(i,j)
            u_prt(i,j) = u_prt(i,j)+gswm_nm_sdi_u(i,j)
            v_prt(i,j) = v_prt(i,j)+gswm_nm_sdi_v(i,j)
          enddo
        enddo
        if (istep==1) write(6,"('tuvz_lbc: added gswm_nm_sdi to ',
     |    't_lbc,z_lbc')")
      endif
!      
!     write(6,"('End tuvz_lbc: istep=',i5,' z_lbc min,max=',
!    |  2e12.6)") istep,minval(z_lbc(lon0:lon1,lat0:lat1)),
!    |  maxval(z_lbc(lon0:lon1,lat0:lat1))

!     write(6,"('End tuvz_lbc: istep=',i5,' t_lbc min,max=',
!    |  2e12.6)") istep,minval(t_lbc(lon0:lon1,lat0:lat1)),
!    |  maxval(t_lbc(lon0:lon1,lat0:lat1))

! am 9/09 lbc_bgrd: for ecmwf runs the background is in u_lbc, v_lbc and perturbations
! in  u_prt and v_prt, call add_uvbnd is not called since background velocities 
! are set and not calculated, so add them together here
! 11/19/09 btf: ditto for ncep_reanalysis
!
       if(len_trim(ecmwf_ncfile) > 0.or.len_trim(ncep_reanalysis) > 0) 
     |   then 
         u_lbc(lon0:lon1,lat0:lat1) = u_lbc(lon0:lon1,lat0:lat1) + 
     |      u_prt(lon0:lon1,lat0:lat1)
         v_lbc(lon0:lon1,lat0:lat1) = v_lbc(lon0:lon1,lat0:lat1) + 
     |       v_prt(lon0:lon1,lat0:lat1)
       else
! am 8/09 lbc_bgrd: set u_lbc and v_lbc these have only perturbations on them
! and should be equal to u_prt and v_prt
        u_lbc(lon0:lon1,lat0:lat1) = u_prt(lon0:lon1,lat0:lat1)
        v_lbc(lon0:lon1,lat0:lat1) = v_prt(lon0:lon1,lat0:lat1)
      endif
!
! am 8/09 lbc_bgrd update arrays for istep = 1; approximation
      if(istep == 1) then
        u_prt_tn(lon0:lon1,lat0:lat1)  = u_prt(lon0:lon1,lat0:lat1)
        v_prt_tn(lon0:lon1,lat0:lat1)  = v_prt(lon0:lon1,lat0:lat1)
        u_prt_tnm(lon0:lon1,lat0:lat1) = u_prt(lon0:lon1,lat0:lat1)
        v_prt_tnm(lon0:lon1,lat0:lat1) = v_prt(lon0:lon1,lat0:lat1)
      endif
!      
! am 8/09 lbc_bgrd: it's not necessary to do the mpi for v_prt and u_prt
#ifdef MPI
!
! Boundary (halo) latitudes:
      flbc(:,:,1) = z_bgrd(lon0-2:lon1+2,lat0-2:lat1+2)
      call mp_bndlats_f2d(flbc,lon0,lon1,lat0,lat1,1)
      z_bgrd(lon0-2:lon1+2,lat0-2:lat1+2)   = flbc(:,:,1)
!
! Boundary (halo) longitudes:
      flbc(:,:,1) = u_prt(lon0-2:lon1+2,lat0-2:lat1+2)
      flbc(:,:,2) = v_prt(lon0-2:lon1+2,lat0-2:lat1+2)
      flbc(:,:,3) = z_bgrd(lon0-2:lon1+2,lat0-2:lat1+2)
      flbc(:,:,4) = u_prt_tn(lon0-2:lon1+2,lat0-2:lat1+2)
      flbc(:,:,5) = v_prt_tn(lon0-2:lon1+2,lat0-2:lat1+2)
      call mp_bndlons_f2d(flbc,lon0,lon1,lat0,lat1,5)
      u_prt(lon0-2:lon1+2,lat0-2:lat1+2)     = flbc(:,:,1)
      v_prt(lon0-2:lon1+2,lat0-2:lat1+2)     = flbc(:,:,2)
      z_bgrd(lon0-2:lon1+2,lat0-2:lat1+2)    = flbc(:,:,3)
      u_prt_tn(lon0-2:lon1+2,lat0-2:lat1+2)  = flbc(:,:,4)
      v_prt_tn(lon0-2:lon1+2,lat0-2:lat1+2)  = flbc(:,:,5)
!
! Periodic points for lbc:
      flbc(:,:,1) = z_bgrd(lon0-2:lon1+2,lat0-2:lat1+2)
      call mp_periodic_f2d(flbc(lon0:lon1,lat0:lat1,1),
     |       lon0,lon1,lat0,lat1,1)
      z_bgrd(lon0-2:lon1+2,lat0-2:lat1+2) = flbc(:,:,1)
#else
!      do lat=lat0,lat1
!        do i=1,2
!          ulbc(i,lat)        = ulbc(nlon+i,lat)
!          ulbc(nlonp2+i,lat) = ulbc(i+2,lat)
!          vlbc(i,lat)        = vlbc(nlon+i,lat)
!          vlbc(nlonp2+i,lat) = vlbc(i+2,lat)
!          ulbc_nm(i,lat)        = ulbc_nm(nlon+i,lat)
!          ulbc_nm(nlonp2+i,lat) = ulbc_nm(i+2,lat)
!        enddo
!      enddo
#endif
!
! Save to secondary histories:
        call addfld("T_LBC","t_lbc at end of sub tuvz_lbc","deg K",
     |    t_lbc(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)
        call addfld("U_LBC","u_lbc at end of sub tuvz_lbc","cm/s",
     |    u_lbc(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)
        call addfld("V_LBC","v_lbc at end of sub tuvz_lbc","cm/s",
     |    v_lbc(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)
        call addfld("Z_LBC","z_lbc at end of sub tuvz_lbc","cm",
     |    z_lbc(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)
     
        call addfld("V_PRT","v_prt at end of sub tuvz_lbc","cm/s",
     |    v_prt(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)
        call addfld("U_PRT","u_prt at end of sub tuvz_lbc","cm",
     |    u_prt(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)
        call addfld("Z_BGRD","z_bgrd at end of sub tuvz_lbc","cm",
     |    z_bgrd(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)
     
        call addfld("V_PRT_TN","v_prt tn end of sub tuvz_lbc","cm/s",
     | v_prt_tn(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)
        call addfld("U_PRT_TN","u_prt tn end of sub tuvz_lbc","cm",
     | u_prt_tn(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)
     
        call addfld("V_PRT_TNM","v_prt tnm end of sub tuvz_lbc","cm/s",
     | v_prt_tnm(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)
        call addfld("U_PRT_TNM","u_prt tnm end of sub tuvz_lbc","cm",
     | u_prt_tnm(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)
     
      end subroutine tuvz_lbc
!-----------------------------------------------------------------------
      subroutine add_uvbnd
      use input_module,only: zatmos_ncfile
!      use uv_bndry,only: ubnd,vbnd
      integer :: lat,i

! am 8/09 lbc_bgrd: first 2 timesteps old uvbnd is called which has background
! and perturbations on it     
      if(istep < 3) then
        u_lbc = 0.
	v_lbc = 0.
      endif
!      
! am 8/09 lbc_bgrd: background UN,VN at LBC will be calculated in uvbgrd      
!
!      call addfld("U_PRT_ADD","u_lbc before adding ubnd","cm/s",
!     |  u_lbc(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)
!      call addfld("V_PRT_ADD","v_lbc before adding vbnd","cm/s",
!     |  v_lbc(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)
     
! uvbgrd has been called: add ubgrd,vbgrd to the perturbation at the lower boundary 
!   of U,V, if not using ncep_reanalysis or ecmwf. See call uvbnd/uvbgrd and
!   call add_uvbnd in dynamics.
!
      do lat=lat0,lat1
        do i=lon0,lon1
          u_lbc(i,lat) = u_lbc(i,lat)+u_bgrd(i,lat)
          v_lbc(i,lat) = v_lbc(i,lat)+v_bgrd(i,lat)
        enddo
      enddo
      if (istep==1) then
        write(6,"('lbc: added ubgrd,vbgrd to u_lbc,v_lbc')")
!       write(6,"('lbc: u_lbc min,max=',2e12.4)")
!    |    minval(u_lbc(lon0:lon1,lat0:lat1)),
!    |    maxval(u_lbc(lon0:lon1,lat0:lat1))
!       write(6,"('lbc: v_lbc min,max=',2e12.4)")
!    |    minval(v_lbc(lon0:lon1,lat0:lat1)),
!    |    maxval(v_lbc(lon0:lon1,lat0:lat1))

!       write(6,"(11x,'ubnd min,max=',2e12.4)")
!    |    minval(ubnd(lon0:lon1,lat0:lat1)),
!    |    maxval(ubnd(lon0:lon1,lat0:lat1))
!       write(6,"(11x,'vbnd min,max=',2e12.4)")
!    |    minval(vbnd(lon0:lon1,lat0:lat1)),
!    |    maxval(vbnd(lon0:lon1,lat0:lat1))
      endif

!      call addfld("UBGRD_ADD","ubgrd output of sub uvbnd","cm/s",
!     |  u_bgrd(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)
!      call addfld("U_BND_ADD","u_lbc after adding ubnd","cm/s",
!     |  u_lbc(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)
!
!      call addfld("VBGRD_ADD","vbgrd output of sub uvbnd","cm/s",
!     |  v_bgrd(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)
!      call addfld("V_BND_ADD","v_lbc after adding vbnd","cm/s",
!     |  v_lbc(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',lat0,lat1,0)

      end subroutine add_uvbnd
!-----------------------------------------------------------------------
      subroutine bndry_semidiurnal
!
! Lower boundary conditions for semi-diurnal tide, using Hough functions.
! This is called once per run from init, and returns t_sdi, u_sdi, v_sdi, 
!   z_sdi at nlat latitudes. 
!
      use input_module,only: tide
      implicit none
!
! Local:
      integer,parameter :: nalf=19, malf=2
      real :: p(nlat,nalf,malf),hough(nlat,5,malf),cp(nalf/2+1)
      complex :: dzb(nlat)
      real :: b(5,19),rl(5),bhour(5),rlamda,xdot(19),ydot(19),
     |  ptscal,theta,ptjm(2*nlat+1)
      integer :: n,jm,l,lm1,m,mm1,j,ld,i,nm1
!
      complex zee(5),cl(5),expdlm
      data b/
     | 0.969152, 0.0     , 0.216046, 0.0     , 0.093838,
     | 0.0     , 0.909763, 0.0     , 0.342113, 0.0     ,
     |-0.245226, 0.0     , 0.798445, 0.0     , 0.421218,
     | 0.0     ,-0.408934, 0.0     , 0.645517, 0.0     ,
     | 0.024633, 0.0     ,-0.543993, 0.0     , 0.464159,
     | 0.0     , 0.071127, 0.0     ,-0.643189, 0.0     ,
     |-0.001292, 0.0     , 0.139613, 0.0     ,-0.699495,
     | 0.0     ,-0.006673, 0.0     , 0.225090, 0.0     ,
     | 0.000042, 0.0     ,-0.019654, 0.0     , 0.320141,
     | 0.0     , 0.000394, 0.0     ,-0.043345, 0.0     ,
     |-0.000001, 0.0     , 0.001772, 0.0     ,-0.079831,
     | 0.0     ,-0.000016, 0.0     , 0.005401, 0.0     ,
     | 0.0     , 0.0     ,-0.000112, 0.0     , 0.012932,
     | 0.0     , 0.0     , 0.0     ,-0.000476, 0.0     ,
     | 0.0     , 0.0     , 0.000005, 0.0     ,-0.001490,
     | 0.0     , 0.0     , 0.0     , 0.000031, 0.0     ,
     | 0.0     , 0.0     , 0.0     , 0.0     , 0.000129,
     | 0.0     , 0.0     , 0.0     ,-0.000002, 0.0     ,
     | 0.0     , 0.0     , 0.0     , 0.0     ,-0.000009/
      data rl/7.8519E5, 3.6665E5, 2.1098E5, 1.3671E5, 0.9565E5/
      real,external :: sddot ! util.F
!
! Longitudinal structure:
      rlamda = -2.*dlamda
      bnd_sdi(1)=cexp(ci*2.*rlamda)
      expdlm=cexp(ci*2.*dlamda)
      do i=2,nlonp4
        bnd_sdi(i)=bnd_sdi(i-1)*expdlm
      enddo
!
! Zero out if user did not provide amp/phase:
      if (all(tide==0.)) then
        t_sdi = 0.
        u_sdi = 0.
        v_sdi = 0.
        z_sdi = 0.
        return
      endif
!
      bhour = tide(6:10)
      do n=1,5
        zee(n)=tide(n)*cexp(ci*pi*bhour(n)/6.)
        cl(n)=csqrt(cmplx(gask/(atm_amu*grav*rl(n))*
     |    (t0(1)*2./7.+(t0(2)-t0(1))/dz)-.25))-.5*ci
      enddo
      jm=2*nlat+1
!
! Set up hough functions (see sphpac.F)
      do n=2,nalf+1
        nm1 = n-1
        do m=2,malf+1
          mm1=m-1
          call alfk(n,m,cp)
          do j=1,jm
            theta = float(j-1)*pi/float(jm-1)
            call lfpt(n,m,theta,cp,ptscal)
            ptjm(j) = ptscal
          enddo
          do j=1,nlat
            p(j,nm1,mm1) = ptjm(2*(nlat+1-j))
          enddo
        enddo
        do j=1,nlat
          p(j,nm1,2)=sqrt(float(n*(n+1)-6))*p(j,nm1,2)-2.*tn(j)*
     |               p(j,nm1,1)
        enddo
      enddo
!
! util.F: real function sddot(n,x,y)
      do l=1,5
        do ld=1,2
          do j=1,nlat
            xdot(:) = p(j,:,ld)
            ydot(:) = b(l,:)
            hough(j,l,ld)=sddot(19,xdot,ydot)
          enddo
        enddo
      enddo
!
! Define module data:
      do j=1,nlat
        t_sdi(j)=0.
        z_sdi(j)=0.
        dzb(j)=0.
      enddo
      do l=1,5
        do j=1,nlat
          z_sdi(j)=z_sdi(j)+zee(l)*hough(j,l,1)
          dzb(j)=dzb(j)+zee(l)*hough(j,l,2)
          t_sdi(j)=t_sdi(j)+ci*atm_amu*grav/gask*zee(l)*cl(l)*
     |      hough(j,l,1)
        enddo
      enddo
      do j=1,nlat
        u_sdi(j)=freq_semidi*re*(1.-(cor(j)/freq_semidi)**2)
        v_sdi(j)=ci*grav*(dzb(j)-2.*cor(j)/(freq_semidi*cs(j))*
     |    z_sdi(j))/u_sdi(j)
        u_sdi(j)=grav*(cor(j)/freq_semidi*dzb(j)-2./cs(j)*
     |    z_sdi(j))/u_sdi(j)
      enddo

!     write(6,"('bndry_semidiurnal: t_sdi min,max=',2e12.4)")
!    |  minval(real(t_sdi)),maxval(real(t_sdi))
!     write(6,"('bndry_semidiurnal: u_sdi min,max=',2e12.4)")
!    |  minval(real(u_sdi)),maxval(real(u_sdi))
!     write(6,"('bndry_semidiurnal: v_sdi min,max=',2e12.4)")
!    |  minval(real(v_sdi)),maxval(real(v_sdi))
!     write(6,"('bndry_semidiurnal: z_sdi min,max=',2e12.4)")
!    |  minval(real(z_sdi)),maxval(real(z_sdi))

      end subroutine bndry_semidiurnal
!-----------------------------------------------------------------------
      subroutine bndry_diurnal
!
! Lower boundary conditions for diurnal tide, using Hough functions.
! This is called once per run from init, and returns t_di, u_di, v_di, 
!   z_di at nlat latitudes. 
!
      use input_module,only: tide2
!
! Local:
      integer,parameter :: nalf=19, malf=2
      real :: p(nlat,nalf,malf),hough(nlat,5,malf),cp(nalf/2+1)
      complex :: dzb(nlat)
      real :: b(1,19),rl (1),bhour(1),rlamda,xdot(19),ydot(19),
     |  ptscal,theta,ptjm(2*nlat+1),pik
      integer :: l,m,j,n,jm,ld,i
      complex zee(1),cl(1),expdlm
!
      data b/
     | 0.282710,
     | 0.0     ,
     |-0.638229,
     | 0.0     ,
     | 0.620521,
     | 0.0     ,
     |-0.336408,
     | 0.0     ,
     | 0.117021,
     | 0.0     ,
     |-0.028332,
     | 0.0     ,
     | 0.005042,
     | 0.0     ,
     |-0.000686,
     | 0.0     ,
     | 0.000074,
     | 0.0     ,
     |-0.000006/
      data rl/0.6909E5/
      real,external :: sddot ! in util.F
!
! Calculate longitudinal structure
      rlamda = -2.*dlamda
      bnd_di(1)=cexp(ci*rlamda)
      expdlm=cexp(ci*dlamda)
      do i=2,nlonp4
        bnd_di(i)=bnd_di(i-1)*expdlm
      enddo
!
! Zero out if user did not provide amp/phase:
      if (all(tide2==0.)) then
        t_di = 0.
        u_di = 0.
        v_di = 0.
        z_di = 0.
        return
      endif
      bhour(1) = tide2(2)
      pik = 3.14159265358979312 
      do n=1,1
        zee(n)=tide2(n)*cexp(ci*pik*bhour(n)/12.)
        cl(n)=csqrt(cmplx(gask/(atm_amu*grav*rl(n))*
     |  (t0(1)*2./7.+(t0(2)-t0(1))/dz)-.25))-.5*ci
      enddo
      jm=2*nlat+1
!
! Set up hough functions:
!
      do n=1,19
        do m=1,2
          call alfk(n,m,cp)
          do j=1,jm
            theta = float(j-1)*pi/float(jm-1)
            call lfpt(n,m,theta,cp,ptscal)
            ptjm(j)=ptscal
          enddo
          do j=1,nlat
            p(j,n,m) = ptjm(2*(nlat+1-j))
          enddo
        enddo
        do j=1,nlat
          p(j,n,2)=sqrt(float(n*(n+1)-2))*p(j,n,2)-tn(j)*p(j,n,1)
        enddo
      enddo
!
! util.F: real function sddot(n,x,y)
      do l=1,1
        do ld=1,2
          do j=1,nlat
            xdot(:) = p(j,:,ld)
            ydot(:) = b(l,:)
            hough(j,l,ld)=sddot(19,xdot,ydot)
          enddo
        enddo
      enddo
!
! Generate t_di, u_di, v_di, z_di:
      do j=1,nlat
        t_di(j)=0.
        z_di(j)=0.
        dzb(j)=0.
      enddo
      do l=1,1
        do j=1,nlat
          z_di(j)=z_di(j)+zee(l)*hough(j,l,1)
          dzb(j)=dzb(j)+zee(l)*hough(j,l,2)
          t_di(j)=t_di(j)+ci*atm_amu*grav/gask*zee(l)*cl(l)*hough(j,l,1)
        enddo
      enddo
      do j=1,nlat
        u_di(j)=.5*freq_semidi*re*(1.-(cor(j)/(.5*freq_semidi))**2)
        v_di(j)=ci*grav*(dzb(j)-cor(j)/(.5*freq_semidi*cs(j))*z_di(j))/
     |    u_di(j)
        u_di(j)=grav*(cor(j)/(.5*freq_semidi)*dzb(j)-1./cs(j)*z_di(j))/
     |    u_di(j)
      enddo

!     write(6,"('bndry_diurnal: t_di min,max=',2e12.4)")
!    |  minval(real(t_di)),maxval(real(t_di))
!     write(6,"('bndry_diurnal: u_di min,max=',2e12.4)")
!    |  minval(real(u_di)),maxval(real(u_di))
!     write(6,"('bndry_diurnal: v_di min,max=',2e12.4)")
!    |  minval(real(v_di)),maxval(real(v_di))
!     write(6,"('bndry_diurnal: z_di min,max=',2e12.4)")
!    |  minval(real(z_di)),maxval(real(z_di))

      end subroutine bndry_diurnal
!-----------------------------------------------------------------------
      subroutine bndry_annual(dday)
      use input_module,only: tideann
      use cons_module,only: pi,cs,cor,tn,re,grav,gask,
     |  freq_ann,dz,atm_amu
C     ****
C     ****     TIDAL BOUNDARY CONDITION FOR ANNUAL MODE
!
! Args:
      real,intent(in) :: dday
!     integer,intent(in) :: idayin
!
! For 1998 spherepack lib code (sphpac.f)
! (replaces old alfpac.f)
!
      integer,parameter :: nalf=24, malf=2
!
! Local:
      real,parameter :: SCALE=1.
      integer :: n,l,m,mp1,j,jm,ld,i,k
      real :: rm,rt2,factor
      real :: p(nlat,nalf,malf),hough(nlat,0:6,2),cp(nalf/2+1)
      complex dzb(nlat),zzb(nlat)
      real :: b(6,24),rl(0:6),xdot(nalf),ydot(nalf),ptscal,
     |  ptjm(2*nlat+1),theta,w(nlat)
      complex cc(0:6,0:6),cl(0:6),expt
      real,external :: sddot
!
      data ((b(i,j),i = 1,6),j = 1,12)/
     |  -0.882922, 0.000000,-0.345087, 0.000000,-0.202228, 0.000000,
     |   0.000000,-0.930826, 0.000000,-0.301357, 0.000000,-0.152720,
     |  -0.466226, 0.000000, 0.567457, 0.000000, 0.407114, 0.000000,
     |   0.000000,-0.362673, 0.000000, 0.694431, 0.000000, 0.438014,
     |  -0.055436, 0.000000, 0.711847, 0.000000,-0.163050, 0.000000,
     |   0.000000,-0.044983, 0.000000, 0.625545, 0.000000,-0.325772,
     |  -0.002909, 0.000000, 0.225723, 0.000000,-0.749160, 0.000000,
     |   0.000000,-0.002773, 0.000000, 0.186467, 0.000000,-0.723674,
     |  -0.000086, 0.000000, 0.034940, 0.000000,-0.435919, 0.000000,
     |   0.000000,-0.000103, 0.000000, 0.029425, 0.000000,-0.379254,
     |  -0.000002, 0.000000, 0.003267, 0.000000,-0.122687, 0.000000,
     |   0.000000,-0.000003, 0.000000, 0.002928, 0.000000,-0.104008/
      data ((b(i,J),i = 1,6),j = 13,24)/
     |   0.0     , 0.000000, 0.000206, 0.000000,-0.021267, 0.000000,
     |   0.0     , 0.0     , 0.000000, 0.000202, 0.000000,-0.018228,
     |   0.0     , 0.0     , 0.000009, 0.000000,-0.002540, 0.000000,
     |   0.0     , 0.0     , 0.000000, 0.000010, 0.000000,-0.002252,
     |   0.0     , 0.0     , 0.0     , 0.000000,-0.000223, 0.000000,
     |   0.0     , 0.0     , 0.0     , 0.0     , 0.000000,-0.000208,
     |   0.0     , 0.0     , 0.0     , 0.0     ,-0.000015, 0.000000,
     |   0.0     , 0.0     , 0.0     , 0.0     , 0.000000,-0.000015,
     |   0.0     , 0.0     , 0.0     , 0.0     ,-0.000001, 0.000000,
     |   0.0     , 0.0     , 0.0     , 0.0     , 0.000000,-0.000001,
     |   0.0     , 0.0     , 0.0     , 0.0     , 0.0     , 0.000000,
     |   0.0     , 0.0     , 0.0     , 0.0     , 0.0     , 0.0     /
      data (rl(n),n=1,6)/
     |  -10.8409E5,-7.0243E5,-2.4874E5,-1.9696E5,-1.0694E5,
     |  -0.9119E5/
!
! Annual boundary coefficients for lower boundary at z = -17.
!
      data ((cc(k,n),n=0,6),k=0,3)/
     |   ( 0.421595E+02, 0.000000E+00),(-0.226810E+00, 0.000000E+00),
     |   ( 0.415692E+00, 0.000000E+00),( 0.547413E-01, 0.000000E+00),
     |   (-0.312575E-01, 0.000000E+00),( 0.280035E-01, 0.000000E+00),
     |   (-0.185304E-01, 0.000000E+00),( 0.300941E-01,-0.734916E-03),
     |   ( 0.100061E+01,-0.168070E+00),(-0.212277E+00, 0.133133E+00),
     |   ( 0.226725E-01,-0.553615E-01),(-0.426632E-01, 0.941169E-03),
     |   (-0.109421E-02,-0.877230E-02),(-0.576746E-02,-0.113662E-01),
     |   ( 0.215525E-01, 0.269012E-01),( 0.212126E-01, 0.377161E-01),
     |   (-0.329359E-01, 0.112051E-01),(-0.117842E-01, 0.145664E-01),
     |   (-0.813882E-02, 0.424029E-01),(-0.705682E-02, 0.100935E-01),
     |   (-0.589241E-02, 0.224780E-01),(-0.956909E-03, 0.169400E-02),
     |   ( 0.283082E-02,-0.579849E-02),( 0.141201E-02,-0.351546E-02),
     |   ( 0.288114E-02, 0.600793E-04),(-0.896708E-04, 0.917448E-04),
     |   ( 0.146016E-02, 0.152191E-02),( 0.911123E-04,-0.141497E-03)/
      data ((cc(k,n),n=0,6),k=4,6)/
     |   (-2.200025E-02,-0.159748E-02),( 0.122976E-03,-0.180022E-03),
     |   ( 0.218866E-02, 0.217849E-02),( 0.361138E-03, 0.264706E-03),
     |   ( 0.803136E-04, 0.595681E-04),(-0.711489E-04, 0.138662E-04),
     |   (-0.205981E-04, 0.237149E-05),( 0.211351E-04, 0.146409E-04),
     |   (-0.195730E-02,-0.810895E-03),(-0.442840E-04, 0.292211E-04),
     |   (-0.123450E-04,-0.538276E-04),(-0.253900E-04,-0.621475E-04),
     |   ( 0.767875E-04,-0.422957E-05),( 0.206095E-04, 0.178031E-05),
     |   (-0.191415E-03, 0.722780E-04),(-0.275580E-04,-0.427594E-04),
     |   ( 0.132229E-03,-0.182378E-03),( 0.341214E-04, 0.355103E-05),
     |   ( 0.326348E-04, 0.543373E-04),(-0.109405E-04, 0.189210E-04),
     |   ( 0.489606E-04, 0.396533E-05)/
!
! 8/20/98: This routine now called at every time step from advnce if
!   the run is advancing calendar day. The single argument is current
!   decimal model day.
! 
      rt2 = sqrt(2.)
      if (tideann==0) then
!
! Zero boundary condition, except for geopotential:
!
        do j = 1,nlat
          zba(j) = cc(0,0)*1./rt2*1.e5
          tba(j) = 0.
          uba(j) = 0.
          vba(j) = 0.
        enddo
        bnda = 1. ! whole array op
        return
      endif
      jm = 2*nlat+1
!
! Height variation
!
      do n = 1,6
        cl(n) = -csqrt(cmplx(gask/(atm_amu*grav*rl(n))*
     |    (t0(1)*2./7.+(t0(2)-t0(1))/dz)-.25))-.5*ci
      enddo
      cl(0) = 0.
!
! Set up legendre polynomials
!
! Using new (1998) spherepack (sphpac.f):
! (nalf=24, see cbndrya.h)
!
      do n=1,nalf
        do m=0,1
          mp1=m+1
          call alfk(n,m,cp)
          do j=1,jm
            theta = float(j-1)*pi/float(jm-1)
            call lfpt(n,m,theta,cp,ptscal)
            ptjm(j) = ptscal
          enddo
          do j=1,nlat
            p(j,n,mp1) = ptjm(2*(nlat+1-j))
          enddo
        enddo
        m = 0
        rm = float(m)
        do j=1,nlat
          p(j,n,2)=sqrt(float(n*(n+1)-m*(m+1)))*p(j,n,2)-rm*tn(j)*
     |      p(j,n,1)
        enddo
      enddo
!
! Now evaluate hough functions
!
      do l=1,6
        do ld=1,2
          do j=1,nlat
            xdot(:) = p(j,:,ld)
            ydot(:) = b(l,:)
            hough(j,l,ld)=sddot(24,xdot,ydot)
          enddo
        enddo
      enddo
!
! Hough function of order zero
!
      do j=1,nlat
        hough(j,0,1) = 1./rt2
        hough(j,0,2) = 0.
      enddo
!
! Generate zba, tba, uba, vba
!
      do j=1,nlat
        zba(j) = 0.
        tba(j) = 0.
        uba(j) = 0.
        vba(j) = 0.
      enddo
!
! Summation over frequency, k, calculation of phase factor
!
      do k = 0,6
        expt = cexp(ci*float(k)*(dday-1.)*86400.*freq_ann)
!
! Summation over order, n
!
        do n = 0,6
          factor = scale
          if (k.eq.0.and.n.eq.0) factor = 1.
          do j = 1,nlat
            zzb(j) = cc(k,n)*hough(j,n,1)*expt*1.e5*factor
            dzb(j) = cc(k,n)*hough(j,n,2)*expt*1.e5*factor
            w(j) = grav/(re*((float(k)*freq_ann)**2-cor(j)**2))
            zba(j) = zba(j)+zzb(j)
            tba(j) = tba(j)+ci*atm_amu*grav/gask*cl(n)*zzb(j)
            uba(j) = uba(j)+w(j)*(cor(j)*dzb(j)-rm*float(k)*freq_ann/
     |        cs(j)*zzb(j))
            vba(j) = vba(j)+ci*w(j)*(float(k)*freq_ann*dzb(j)-rm*cor(j)/
     |        cs(j)*zzb(j))
          enddo
        enddo
      enddo
!
! Longitudinal structure
!
      bnda = 1. ! whole array op
      end subroutine bndry_annual
!-----------------------------------------------------------------------
      subroutine bndry_2day
      use input_module,only: tide3m3
      use cons_module,only: pi,cs,cor,tn,t0,re,dz,dlamda,
     |  gask,freq_3m3,grav,atm_amu
C     ****
C     ****     Tidal boundary condition for the two day mixed
C     ****     Rossby-Gravity wave.  This is denoted by indices (3,-3)
C     ****     and is asymmetric.  The charactistic depth is 10.5km.
C     ****     See Forbes, Jeffrey M. Tidal and Planetary Waves. In:
C     ****     Upper Mesosphere and Lower Thermosphere:  A Review of
C     ****     Experiment and Theory.  Geophysical Monograph 87,
C     ****     Copyright 1995 by American Geophysical Union.
C     ****
C     ****     Local variables
C     ****
      integer,parameter :: nalf=19, malf=2
      real :: p(nlat,nalf,malf),hough(nlat,5,malf),
     |  cp(nalf/2+1),xdot(10),ydot(10),ptjm(2*nlat+1),ptscal,
     |  theta,w(nlat)
      complex :: dzb(nlat)
      complex zee,cl,expdlm
!
! Where:
!
! B(10) = Coefficients of normalised associated Legendre
!         polynomials needed to build Hough(3,-3).
!         i.e.: P(l,3), where l = 4,6,8,10,12 (higher
!         terms are too small to merit consideration)
! 
! RL = equivalent depth for wave
!    = 10.5km = 0.105E7 cm
!
! BHOUR = phase of wave measured in hours with constraint:
!         0 .le. BHOUR .lt. p
!         ( p = period of wave in hours = 2.07412*24)
!
! TIDE = amplitude of wave in cm
!
! ZEE = complex amplitude incorporating phase (BHOUR) and
!       amplitude (TIDE)
!
! EXPDLM = exp(2*pi*i*dlamda), where dlamda is
!          longitudinal grid spacing (radians)
!
      real :: b(10) = (/
     |  0.0     ,
     | -0.997138,
     |  0.0     ,
     |  0.075553,
     |  0.0     ,
     | -0.002905,
     |  0.0     ,
     |  0.000069,
     |  0.0     ,
     | -0.000001/)
      real,parameter :: rl=0.105e7
      real :: bhour,tide,rlamda
      integer :: jm,n,nm2,m,mm2,j,ld,i
      real,external :: sddot
      integer :: itide3m3
!
      itide3m3 = 0
      if (tide3m3(1) /= 0..or.tide3m3(2) /= 0.) itide3m3 = 1
!
!  Two-day tide is present
!
!  zee = complex amplitude of wave at time zero.
!      = amp * exp(2*pi*i*BHOUR/P2DAY)
!
!  cl = vertical wave length
!     = ((K*H + dH/ds)/hn - 1/4)**(1/2) -i/2
!
!   where:
!    K = 2/7
!    H = scale height (cm)
!    hn = equivalent depth for mode = 1.05E7 cm
!    s = vertical pressure coordinate
!  
! Init:
      do j=1,nlat
        zb3m3(j)=0.
        tb3m3(j)=0.
        ub3m3(j)=0.
        vb3m3(j)=0.
      enddo
      if (itide3m3.eq.0) return
!
      tide = tide3m3(1)
      bhour = tide3m3(2)
      zee = tide*cexp(ci*freq_3m3*bhour*60.*60.)
      cl = csqrt(cmplx(gask/(atm_amu*grav*rl)*
     |  (t0(1)*2./7.+(t0(2)-t0(1))/dz)-.25))-.5*ci
      jm=2*nlat+1
!
! Set up Hough functions
!
! Using new (1998) spherepack (sphpac.f):
      do n=3,12
        nm2 = n-2
        do m=3,4
          mm2=m-2
          call alfk(n,m,cp)
          do j=1,jm
            theta = float(j-1)*pi/float(jm-1)
            call lfpt(n,m,theta,cp,ptscal)
            ptjm(j) = ptscal
          enddo
          do j=1,nlat
            p(j,nm2,mm2) = ptjm(2*(nlat+1-j))
          enddo
        enddo
        m = 3
        do j=1,nlat
          p(j,nm2,2)=sqrt(float(n*(n+1)-m*(m+1)))*p(j,nm2,2)-
     |               float(m)*tn(j)*p(j,nm2,1)
        enddo
      enddo
!
      do ld = 1,2
        do j = 1,nlat
          xdot(:) = p(j,1:10,ld)
          ydot(:) = b(:)
          hough(j,1,ld) = sddot(10,xdot,ydot)
!         hough(j,1,ld) = sddot(10,p(j,1,ld),nlat,b,1)
        enddo
      enddo
!
! Generate zb3m3, ub3m3, vb3m3 AND tb3m3
!
      do j = 1,nlat
        zb3m3(j) = zee*hough(j,1,1)
        dzb(j) = zee*hough(j,1,2)
        tb3m3(j) = ci*atm_amu*grav/gask*zee*cl*hough(j,1,1)
      enddo
      do j = 1,nlat
        ub3m3(j) = freq_3m3*re*(1.-(cor(j)/freq_3m3)**2)
        vb3m3(j) = ci*grav*(dzb(j)-float(m)*cor(j)/(freq_3m3*cs(j))
     |             *zb3m3(j))/ub3m3(j)
        ub3m3(j) = grav*(cor(j)/freq_3m3*dzb(j)-float(m)/cs(j)*
     |             zb3m3(j))/ub3m3(j)
      enddo
!
! Calculate longitudinal structure
!
      rlamda = -2.*dlamda
      bnd3m3(1)=cexp(ci*float(m)*rlamda)
      expdlm=cexp(ci*float(m)*dlamda)
      do i=2,nlonp4
        bnd3m3(i)=bnd3m3(i-1)*expdlm
      enddo
      end subroutine bndry_2day
!-----------------------------------------------------------------------
      subroutine bndry_planetary
      use init_module,only: iter,glat
      use cons_module,only: pi,dlamda,dt
!
! Calculate contribution to ZB from planetary waves
! 7/31/09 btf: This version obtained from Hanli.
!
!      real,parameter :: amplan=0.0, t0plan=6.0E+5
      real,parameter :: amplan=5.00E4, t0plan=6.0E+5
!     real,parameter :: amplan=11.0E4, t0plan=2.5E+5
!     data amplan,t0plan/ 5.0E4,2.5E+5/
!     data amplan,t0plan/ 11.0E4,2.5E+5/
!
!     real,parameter :: frqplan=1.212e-5  ! 6 day
      real,parameter :: frqplan=3.356e-5   ! 52 hour
!
! Local:
      real :: time,fac1
      integer :: j,i,istartlat,iendlat,iterstart,iterend
      real,parameter :: startlat=32.5
      real,parameter :: endlat=82.5
      real,parameter :: lonphs=-140./180.*3.14159
!
!      iterstart = (150.*86400.+.1)/dt
!      iterstart = 0.0
!      iterstart = (15.*86400.+.1)/dt
!      iterend = (25.*86400.+.1)/dt     ! caseB
      iterstart = (30.*86400.+.1)/dt
      iterend = (40.*86400.+.1)/dt

      istartlat = 0
      iendlat = 0
      do j=1,nlat-1
        if (startlat >= glat(j) .and. startlat < glat(j+1)) istartlat=j
        if (endlat > glat(j) .and. endlat <= glat(j+1)) iendlat=j
      enddo
      if (istartlat==0 .or. iendlat==0) then
        write(6,"(/,'>>> bndry_planetary: could not find index to',
     |    ' startlat=',f8.2,': glat=',/,(6f8.2))") startlat,glat
        call shutdown('startlat')
      endif
      zbplanet = 0. ! init whole-array
!
! 10/19/10 btf and Jia Yue: Always define time (sub bndry_planetary)
      time = (iter-iterstart)*dt
      if (iter < iterstart .or. iter > iterend) then
         fac1 = exp(-(time/t0plan)**2)
      else
         fac1 = 1.
      endif
!
! 10/8/04 btf: To implement planetary waves, uncomment below code.
!
      do j = istartlat,iendlat
        do i = 1,nlonp4
! **** PLANETARY WAVE 1
          zbplanet(i,j) = amplan*(sin((glat(j)-30.)*pi/(60.-
     |           (glat(j)-60.)/3.)))**2*
     |           fac1*cos(-pi+(i-3)*dlamda-lonphs)   ! i=3 corresponds to -pi

!         zbplanet(i,j) = amplan*(sin((glat(j)-30.)*pi/(60.-
!    |           (glat(j)-60.)/3.)))**2*
!    |           (1.-exp(-time/t0plan))*sin(1.*(-pi+(i-3)*dlamda))

!         zbplanet(i,j) = amplan*(sin((glat(j)-30.)*pi/(60.-
!    |           (glat(j)-60.)/3.)))**2*
!    |           sin(1.*(-pi+(i-3)*dlamda))
!         zbplanet(i,j) = 0.
!
! **** PLANETARY WAVE 2
!         zbplanet(i,j) = amplan*(sin((glat(j)-30.)*pi/(60.-
!    |           (glat(j)-60.)/3.)))**2*
!    |           (1.-exp(-time/t0plan))*sin(2.*(-pi+(i-3)*dlamda))
!         zbplanet(i,j) = amplan*(sin((glat(j)-30.)*pi/(60.-
!    |           (glat(j)-60.)/3.)))**2*
!    |           sin(2.*(-pi+(i-3)*dlamda))
! **** PW 3
!          zbplanet(i,j) = amplan*(sin((glat(j)-30.)*pi/(60.-
!     |           (glat(j)-60.)/3.)))**2*
!     |           fac1*sin(-pi+3.*(i-3)*dlamda+frqplan*time)


        enddo
      enddo
!
! Top north latitude is zero:
      zbplanet(:,nlat) = 0.
      zbplanet(:,1) = 0.

!     write(6,"('bndry_planet: zbplanet min,max=',2e12.4)")
!    |  minval(zbplanet),maxval(zbplanet)

      end subroutine bndry_planetary
!-----------------------------------------------------------------------
      subroutine bndry_kelvin
!
! Calculate contribution to ZB from ultra-fast Kelvin Waves
!
      use init_module,only: iter,glat
      use cons_module,only: pi,dt
      use hist_module,only: nstep,modeltime ! for print only
!
! Local:
! 10/7/04 btf: amplan=1.0e5 resulted in NaN's in most fields after ~1/2 hour.
! 10/7/04 btf: amplan=1.0e4 resulted in NaN's in most fields after ~14 hours.
!  2/1/05 btf: changed t0plan from 2.5e4 to 2.5e5 for timegcm1.2
!
!     real,parameter :: amplan=1.0E5, t0plan=2.5E+4
!     real,parameter :: amplan=1.0E4, t0plan=2.5E+4 ! timegcm1.1
      real,parameter :: amplan=1.0E4, t0plan=2.5E+5 ! timegcm1.2
      real :: omg3d,wvx1,time
      integer :: j,i
!
      omg3d = 2.*pi/86400./3.    ! radian frequency of 3 days (rad/sec)
      wvx1 = 2.*pi/360.          ! wavenumber 1 (rad/degree)
      zbkelvin = 0.              ! init whole-array
      time = iter*dt

      do j = 2,nlat-1
        do i = 1,nlonp4
          zbkelvin(i,j) = zbkelvin(i,j)+
     |           amplan*exp(-(glat(j)/30.)**2)*
     |           sin(omg3d*time-wvx1*(i-3)*5.)*
     |           (1.-exp(-time/t0plan))
        enddo
!
!       write(6,"('bndry_kelvin: mtime=',3i4,' nstep=',i4,' dt=',f8.2,
!    |    ' iter=',i5,' time=',e12.4,' j=',i3,' zbkelvin(:,j)=',/,
!    |    (6e12.4))") modeltime(1:3),nstep,dt,iter,time,j,zbkelvin(:,j)
!
      enddo
      zbkelvin(:,1) = 0.    ! highest south latitude
      zbkelvin(:,nlat) = 0. ! highest north latitude

!     write(6,"('bndry_kelvin: zbkelvin min,max=',2e12.4)")
!    |  minval(zbkelvin),maxval(zbkelvin)

      end subroutine bndry_kelvin
!-----------------------------------------------------------------------
      subroutine bndry_comp(barm,lev0,lev1,lon0,lon1,lat)
      use cons_module,only: rmass_o3
      use solgar_module,only: xoxlb
!
! This is called from subdomain latitude loop in dynamics.
! Note fb is calculated at longitude subdomains only.
!
! CALCULATE MATRICES B(nlonp4,2,2) AND VECTORS FB(nlonp4,2)
!   REPRESENTING BOUNDARY CONDITION IN COMP, WHERE PSI1
!   AND PSI2 ARE EVALUATED.
!
! CURRENT BOUNDARY CONDITION IS:
!
!   0.5*(PSI1(-1/2)+PSI1(1/2)) = PSO2LB
!   0.5*(PSI2(-1/2)+PSI2(1/2)) = XOXLB*RMOX(REAL)/MBAR
!   WHERE:
!     XOXLB = OX NUMBER DENSITY MIXING RATIO AT LOWER
!       BOUNDARY
!     PSO2LB= O2 MASS MIXING RATIO AT LOWER BOUNDARY
!
!   THIS GIVES:
!     PSI1(-1/2) = B(1,1)*PSI1(1/2)+B(1,2)*PSI2(1/2)+FB(1)
!     PSI2(-1/2) = B(2,1)*PSI1(1/2)+B(2,2)*PSI2(1/2)+FB(2)
!     WHERE:
!       B(1,1) = B(2,2) = -1.
!       B(1,2) = B(2,1) =  0.
!       FB(1)  = 2.*PSO2LB
!       FB(2)  = 2.*XOXLB*RMOX(TRUE)/MBAR
!
! Args:
      integer,intent(in) ::
     |  lev0,lev1,           ! first and last level indices, this task
     |  lon0,lon1,           ! first and last longitude indices, this task
     |  lat                  ! latitude index
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  barm                 ! mbar (from addiag)
!
! Local:
      real,parameter :: pso2lb=0.24
      integer :: i,k,kk,nmsk
      real :: rmo3a(nlonp4),rmo3b(nlonp4)
!
! SET MATRIX B
!
      do i = 1,nlonp4
        b(i,1,1) = -1.
        b(i,2,2) = -1.
        b(i,1,2) =  0.
        b(i,2,1) =  0.
      enddo
!
! rmo3a, rmo3b = RMOX(TRUE),  K = 1/2, 3/2
!
      rmo3a(:) = rmass_o3
      rmo3b(:) = rmass_o3
!
! rmo3a = RMOX(TRUE), K=0  (EXTRAPOLATION)
! Set fb(1) AND fb(2). Note fb is defined at subdomains only.
!
      do i = lon0,lon1
        rmo3a(i) = 1.5*rmo3a(i)-0.5*rmo3b(i)
        fb(i,1,lat) = 2.*pso2lb
        fb(i,2,lat) = 2.*xoxlb(lat)*rmo3a(i)/barm(1,i)

!       write(6,"('bndry_comp: i=',i3,' lat=',i3,' xoxlb(lat)=',e12.4,
!    |    ' rmo3a(i)=',e12.4,' barm(1,i)=',e12.4,' fb(i,2)=',
!    |    e12.4)") i,lat,xoxlb(lat),rmo3a(i),barm(1,i),
!    |    fb(i,2)

      enddo

!     write(6,"('bndry_comp: lat=',i3,' fb(:,1,lat)=',/,(6e12.4))")
!    |  lat,fb(:,1)
!     write(6,"('bndry_comp: lat=',i3,' fb(:,2,lat)=',/,(6e12.4))")
!    |  lat,fb(:,2)

      end subroutine bndry_comp
!-----------------------------------------------------------------------
      end module lbc
