!      
      module gswm09_module
!
! This software is part of the NCAR TIE-GCM.  Use is governed by the 
! Open Source Academic Research License Agreement contained in the file 
! tiegcmlicense.txt.
!      
! 04/13 bae:  Read a single month of GSWM09 boundary perturbations Z, T, U, V, 
!    from X Zhang et al, JGR, 2010 and interpolate to model grid and time.
C   value = amp*cos(nn*pi*(LT-phase)/12), where nn=1D,2S for each UT/glon
!
      use params_module,only: nlon,nlat,nlonp4,nlonp2
      use init_module,only: glon,glat
      use mpi_module,only: lon0,lon1,lat0,lat1
      use addfld_module,only: addfld
      use cons_module,only: pi
      use input_module,only:
     |  gswm09_nlev      ! number of levels from lower boundary to nudge with GSWM09
			 ! if ctmt,gswm09_nlev=0,1, only lbc; if >1 use zonal mean TIDI winds
			 ! and MSIS (or TIEGCM) Tn and %nden/nden (convert to geop z)
      implicit none
      integer,parameter ::    nalt9=13	! for outzav_04080_120f107_90to120km
!
#ifdef SUN
#include <netcdf3.inc>
#else
#include <netcdf.inc>
#endif 
! 
! GSWM09 boundary data at model grid, output by this module,
!   to be used by dt, duv, addiag, etc.

! Easiest assumption for gswm09_nlev=1 only (GSWM09 tides at lbc only)
!  Convert to 97km (use 97.5km at #40=nlb) at UT, for nlat glats and nlonp4 glons
! Interpolation from constant heights 95:2.5:105km for "nudging" on gswm09_nlev p sfcs
!  lbc ~97km is ~94-99km, so interpolate 2.5km in ht for gsmw09_nlev and add
!   zonal means (zm) from files like ~emery/tidi/outzav_04080_120f107_95to105km
      integer,parameter,private :: nmonth= 12, nhour=24, nlat5=37 
!      integer,parameter,private :: nalt2p5=161	! Not needed for partial array
!  dimension of 2 at front is nn=1 diurnal and nn=2 semidiurnal
! dimension at end nalt9 (=5 for 95:2.5:105, =11 for 95:2.5:120) is number of hts to save
!  NOTE:  if gswm09_nlev so large z hts>~105km, must revise dim 5 and zm file

! allocated subdomains:
! Allocated by sub alloc_gswm09 (called by allocdata):
      real,allocatable,dimension(:,:,:) :: ! (lon0:lon1,lat0:lat1,nalt9)
     |  gswm09_z, gswm09_t, gswm09_u, gswm09_v

! Arrays on constant p levels for GSWM09 plus zonal means (zm) (for lbc,dt,duv.F)
      real,allocatable,dimension(:,:,:) :: ! (lon0:lon1,lat0:lat1,gswm09_nlev)
     |  t_zmgswm09, u_zmgswm09, v_zmgswm09
!
! Private module data, read by sub rdgswm09:
! Will be allocated (lon0:lon1,lat0:lat1,nhour) (no month since must be read in)
      real,dimension(nlonp4,nlat,nalt9,2),private :: ! 2=D,S
     |  z_amp, z_ph, t_amp, t_ph, u_amp, u_ph, v_amp, v_ph

      integer,dimension(nlat),private :: jl1
      real,dimension(nlonp4),private :: glonp4

! 6/13 bae:  Define gswm09_alt(nalt9) as 13 alts from 90:2.5:105km
!      real,dimension(nalt9) :: gswm09_alt
!      real,dimension(13) :: gswm09_alt

!
      contains 
!-----------------------------------------------------------------------
      subroutine getgswm09 (itp)
!
! Module driver to read monthly GSWM09 files, and do time interpolations.

      use init_module,only: istep,iyear,iday,secs
!
! Files provided by user via namelist read where have 6 di and 8 sdi tides, 
!  where s0 is migrating tide, nm di tides are w1,w2,e1,e2,e3, and 
!  nm sdi tides are w1,w2,w3,w4,e1,e2,e3,e4 where phases and amps can be revised
      use input_module,only: 
     |  gswm09_di_file,gswm09_sdi_file

!  if gswm09_nlev>0, must read a zonal mean (zm) file for the day and the F107 created
!    from ~emery/msis00/msis2tie.f every 2.5km from 90:2.5:120km (13 alts)
!    and from ~emery/tidi/meankm.m to create zonal mean Tn, %nden/nden, Un, Vn 
!    on glat dres grid of TIEGCM (72 glats -88.75:2.5:88.75) to use in 2.5 or 5deg
!  Sample output: ~emery/tidi/outzav_04080_120f107_90to120km

! interface p levels for Z,DEN and midpoint p levels for TN,UN,VN
!      use fields_module,only: z,	! z(nlev,nlon,nlat,2) usually allocated
!     |  tn,un,vn
!
!
! Driver for obtaining monthly GSWM09 data, called once per timestep from lbc.
!
! Args:
      integer,intent(in) :: itp
! Local:
      integer :: iprint,i
!
      iprint = 0
      if (istep==1) iprint = 1

! Read diurnal and semidiurnal monthly GSWM09 files
        if (istep==1) then
! Read diurnal and semidiurnal monthly GSWM09 files
         if (len_trim(gswm09_di_file)>0) 
     |     call rdgswm09(gswm09_di_file,'diur')
         if (len_trim(gswm09_sdi_file)>0)
     |     call rdgswm09(gswm09_sdi_file,'semi')
        endif		! if (istep==1) then
!
        call mkgswm09(secs,iprint)

      end subroutine getgswm09
!-----------------------------------------------------------------------
      subroutine rdgswm09(monfile,type)
      use input_module,only: mxlen_filename
      implicit none

      integer,external :: nextlu
!
! Args:
      character(len=*),intent(in) :: monfile,type
!
! Local:
      integer :: ncid,istat,lu_typ
      character(len=mxlen_filename) :: dskfile
      real :: adph1,adph2
      character(len=240) :: char240
      real,dimension(8,2) :: tmp	! Read in dres, but save sres or dres
! 01/13 BA Emery: add amplitude multipliers facamp and phase shifts ishift for GSWM09  
! ishift = shift in longitude to shift in phase (negative is earlier)
      integer :: i,il,is,id2,js,nn,k,j,n
! facamp = multiplier for the amplitude of the tide (default=1.)
      real :: facamp,phsft
!
!  Multiply the values of the GSWM 2009 tides by facamp (default=1.) and shift
!   the phases by phsft (+ shift earlier and - shift later, default=0)
!  Set defaults
          facamp = 1.
	  phsft = 0.
!  Set changes from defaults
      select case(trim(type))
        case('diur')    ! diurnal (13 wave numbers +/-6)
          nn = 1
! 05/13 bae:  Shift 2 hour earlier to fit C/NOFS and incr amp 1.5 (negl for day, good for PRE)
! 09/14 bae:  Shift 2 hour earlier to fit C/NOFS Viz but incr amp not nec
	  phsft = 2.
	  facamp = 2.0
        case('semi')   ! semi-diurnal(13 wave numbers +/-6)
          nn = 2
! 05/13 bae:  Shift 2 hour earlier to fit C/NOFS and try 2.0 facamp
! 09/14 bae:  Shift 2 hour earlier to fit C/NOFS Viz but incr amp not nec
	  phsft = 2.
	  facamp = 2.0
        case default
          write(6,"(/,'>>> rdgswm09: unknown type=',a)") type
          call shutdown('rdgswm09')
      end select

! Open monthly file
          lu_typ = nextlu()
        write (6,"(1x,'facamp phsft lu_typ, GSWM09 monfile =',2f6.2,
     |    i3,1x,a)") facamp,phsft,lu_typ,trim(monfile)
          open(lu_typ,file=trim(monfile),status='old')
!  Read monthly file as made in /hao/aim3/emery/gswm/mkpham_gswm09.f
       is = 2
       do i=3,146
        id2 = (i+1)/2
        if (id2*2==i+1) is = is + 1
        do k=1,13	! 90:2.5:120km
         j = 0
         do js=1,36
           j = j + 1
           read (lu_typ,"(4(f6.2,e13.5))") (tmp(n,1),n=1,8)
           if (nlat==72) then
            z_amp(i,j,k,nn) = tmp(2,1)*facamp
            z_ph(i,j,k,nn) = tmp(1,1)-phsft
            t_amp(i,j,k,nn) = tmp(4,1)*facamp
            t_ph(i,j,k,nn) = tmp(3,1)-phsft
            u_amp(i,j,k,nn) = tmp(6,1)*facamp
            u_ph(i,j,k,nn) = tmp(5,1)-phsft
            v_amp(i,j,k,nn) = tmp(8,1)*facamp
            v_ph(i,j,k,nn) = tmp(7,1)-phsft
           endif
           j = j + 1
           read (lu_typ,"(4(f6.2,e13.5))") (tmp(n,2),n=1,8)
! Save in dres or sres
           if (nlat==72) then
            z_amp(i,j,k,nn) = tmp(2,2)*facamp
            z_ph(i,j,k,nn) = tmp(1,2)-phsft
            t_amp(i,j,k,nn) = tmp(4,2)*facamp
            t_ph(i,j,k,nn) = tmp(3,2)-phsft
            u_amp(i,j,k,nn) = tmp(6,2)*facamp
            u_ph(i,j,k,nn) = tmp(5,2)-phsft
            v_amp(i,j,k,nn) = tmp(8,2)*facamp
            v_ph(i,j,k,nn) = tmp(7,2)-phsft
           else
            if (id2*2==i+1) then
            z_amp(is,js,k,nn) = 0.5*(tmp(2,1)+tmp(2,2))*facamp
! Check phase over 0/24 line but do not bother to keep w/in 0-24 range since cosine OK
            adph1 = 0.
            if (tmp(1,1)-tmp(1,2)<-12.) adph1=24.
            adph2 = 0.
            if (tmp(1,1)-tmp(1,2)>12.) adph2=24.
            z_ph(is,js,k,nn) = 0.5*(tmp(1,1)+adph1+tmp(1,2)+adph2)-phsft
            t_amp(is,js,k,nn) = 0.5*(tmp(4,1)+tmp(4,2))*facamp
            adph1 = 0.
            if (tmp(3,1)-tmp(3,2)<-12.) adph1=24.
            adph2 = 0.
            if (tmp(3,1)-tmp(3,2)>12.) adph2=24.
            t_ph(is,js,k,nn) = 0.5*(tmp(3,1)+adph1+tmp(3,2)+adph2)-phsft
            u_amp(is,js,k,nn) = 0.5*(tmp(6,1)+tmp(6,2))*facamp
            adph1 = 0.
            if (tmp(5,1)-tmp(5,2)<-12.) adph1=24.
            adph2 = 0.
            if (tmp(5,1)-tmp(5,2)>12.) adph2=24.
            u_ph(is,js,k,nn) = 0.5*(tmp(5,1)+adph1+tmp(5,2)+adph2)-phsft
            v_amp(is,js,k,nn) = 0.5*(tmp(8,1)+tmp(8,2))*facamp
            adph1 = 0.
            if (tmp(7,1)-tmp(7,2)<-12.) adph1=24.
            adph2 = 0.
            if (tmp(7,1)-tmp(7,2)>12.) adph2=24.
            v_ph(is,js,k,nn) = 0.5*(tmp(7,1)+adph1+tmp(7,2)+adph2)-phsft
            endif
           endif
          enddo
         enddo
! TEMP - write out u_amp,ph for all i at 97.5km (k=4) at equator (nlat/2)
        write (6,"(1x,'97.5km equatorS UT=0, nn glon u_amp,ph =',i2,
     |   3f9.3)") nn,glon(i-2),u_amp(i,nlat/2,4,nn),u_ph(i,nlat/2,4,nn)
       enddo
!        if (nn==2) stop

! Wrap-around points
      do k=1,nalt9
          z_amp(1:2,:,k,nn)           = z_amp(nlonp4-3:nlonp4-2,:,k,nn)
          z_amp(nlonp4-1:nlonp4,:,k,nn) = z_amp(3:4,:,k,nn)
          t_amp(1:2,:,k,nn)           = t_amp(nlonp4-3:nlonp4-2,:,k,nn)
          t_amp(nlonp4-1:nlonp4,:,k,nn) = t_amp(3:4,:,k,nn)
          u_amp(1:2,:,k,nn)           = u_amp(nlonp4-3:nlonp4-2,:,k,nn)
          u_amp(nlonp4-1:nlonp4,:,k,nn) = u_amp(3:4,:,k,nn)
          v_amp(1:2,:,k,nn)           = v_amp(nlonp4-3:nlonp4-2,:,k,nn)
          v_amp(nlonp4-1:nlonp4,:,k,nn) = v_amp(3:4,:,k,nn)
          z_ph(1:2,:,k,nn)             = z_ph(nlonp4-3:nlonp4-2,:,k,nn)
          z_ph(nlonp4-1:nlonp4,:,k,nn) = z_ph(3:4,:,k,nn)
          t_ph(1:2,:,k,nn)             = t_ph(nlonp4-3:nlonp4-2,:,k,nn)
          t_ph(nlonp4-1:nlonp4,:,k,nn) = t_ph(3:4,:,k,nn)
          u_ph(1:2,:,k,nn)             = u_ph(nlonp4-3:nlonp4-2,:,k,nn)
          u_ph(nlonp4-1:nlonp4,:,k,nn) = u_ph(3:4,:,k,nn)
          v_ph(1:2,:,k,nn)             = v_ph(nlonp4-3:nlonp4-2,:,k,nn)
          v_ph(nlonp4-1:nlonp4,:,k,nn) = v_ph(3:4,:,k,nn)
      enddo
!
! Set glonp4 or glon plus wrap-around points
      do il=1,nlon
       glonp4(il+2) = glon(il)
      enddo
       glonp4(1) = glonp4(nlon+1)
       glonp4(2) = glonp4(nlon+2)
       glonp4(nlon+3) = glonp4(3)
       glonp4(nlon+4) = glonp4(4)

      write(6,"('Completed read from GSWM09 monthly data file ',a)") 
     |    trim(monfile)
      write(6,"(/,72('-'))")
      end subroutine rdgswm09
!-----------------------------------------------------------------------
      subroutine mkgswm09 (secs,iprint)
!
! Use data read from gswm09 file to return perturbation in Z
! and T at current model date and time. It is assumed that the gswm09 file
! provides amps and phases for a single month where no interpolation is
! needed since LT can be defined at each glon for UT=secs
!
      use hist_module,only: modeltime
#ifdef MPI
      use mpi_module,only: mp_periodic_f2d
#endif
      implicit none
!
! Args:
      integer,intent(in) :: iprint
      real,intent(in) :: secs
!
! Local:
      integer :: il,jl,k,i0,j0,i1,j1
      real :: uth,slt1
      real :: fgswm09(lon0:lon1,lat0:lat1,4) ! for mp calls
! TEMP for printout 9/14:
      integer :: i,j
      real,dimension(13) :: alt_13

! For addfld calls:
      i0 = lon0 ; i1 = lon1
      j0 = lat0 ; j1 = lat1
! 

! Calculate nalt9 altitudes (90:2.5:120km)
      do k=1,nalt9
       alt_13(k) = 90. + (k-1)*2.5

! tide=amp*cos[n*omega*LT-phase] or amp*cos[n*omega*UT +s*glon -phase
      uth = secs/3600.
      do il=lon0,lon1 
        slt1 = uth + glonp4(il)/15.
        if (slt1.lt.0.) slt1 = slt1 + 24.
        do jl=lat0,lat1
           gswm09_z(il,jl,k) = 
     |       z_amp(il,jl,k,1)*cos(pi*(uth-z_ph(il,jl,k,1))/12.)
     |       + z_amp(il,jl,k,2)*cos(2.*pi*(uth-z_ph(il,jl,k,2))/12.)
           gswm09_t(il,jl,k) = 
     |       t_amp(il,jl,k,1)*cos(pi*(uth-t_ph(il,jl,k,1))/12.)
     |       + t_amp(il,jl,k,2)*cos(2.*pi*(uth-t_ph(il,jl,k,2))/12.)
           gswm09_u(il,jl,k) = 
     |       u_amp(il,jl,k,1)*cos(pi*(uth-u_ph(il,jl,k,1))/12.)
     |       + u_amp(il,jl,k,2)*cos(2.*pi*(uth-u_ph(il,jl,k,2))/12.)
           gswm09_v(il,jl,k) = 
     |       v_amp(il,jl,k,1)*cos(pi*(uth-v_ph(il,jl,k,1))/12.)
     |       + v_amp(il,jl,k,2)*cos(2.*pi*(uth-v_ph(il,jl,k,2))/12.)
        enddo	! jl=lat0,lat1
      enddo	! il=lon0:lon1

! Do mpi periodic points exchange for gswm with f2d(:)
! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4
!
#ifdef MPI
           fgswm09(:,:,1) = gswm09_z(lon0:lon1,lat0:lat1,k)
           fgswm09(:,:,2) = gswm09_t(lon0:lon1,lat0:lat1,k)
           fgswm09(:,:,3) = gswm09_u(lon0:lon1,lat0:lat1,k)
           fgswm09(:,:,4) = gswm09_v(lon0:lon1,lat0:lat1,k)
           call mp_periodic_f2d(fgswm09,lon0,lon1,lat0,lat1,4)
           gswm09_z(lon0:lon1,lat0:lat1,k) = fgswm09(:,:,1)
           gswm09_t(lon0:lon1,lat0:lat1,k) = fgswm09(:,:,2)
           gswm09_u(lon0:lon1,lat0:lat1,k) = fgswm09(:,:,3)
           gswm09_v(lon0:lon1,lat0:lat1,k) = fgswm09(:,:,4)
#else
! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4
      gswm09_z(1:2,:,k)              = gswm09_z(nlonp4-3:nlonp4-2,:,k)
      gswm09_z(nlonp4-1:nlonp4,:,k)  = gswm09_z(3:4,:,k)
      gswm09_t(1:2,:,k)              = gswm09_t(nlonp4-3:nlonp4-2,:,k)
      gswm09_t(nlonp4-1:nlonp4,:,k)  = gswm09_t(3:4,:,k)
      gswm09_u(1:2,:,k)              = gswm09_u(nlonp4-3:nlonp4-2,:,k)
      gswm09_u(nlonp4-1:nlonp4,:,k)  = gswm09_u(3:4,:,k)
      gswm09_v(1:2,:,k)              = gswm09_v(nlonp4-3:nlonp4-2,:,k)
      gswm09_v(nlonp4-1:nlonp4,:,k)  = gswm09_v(3:4,:,k)
#endif
      enddo	! k=1,nalt9

! TEMP 9/14:  Print out values of _lbc for plotting in plot_mon_ut0.m - good for nproc=1
      if (iprint==1) then         
         write (40,"(13f6.1)") (alt_13(i),i=1,nalt9)
         write (40,"(10f8.2)") (glonp4(i),i=lon0,lon1)
         write (40,"(10f8.2)") (glat(i),i=lat0,lat1)
          do i=lon0,lon1
           do j=lat0,lat1
            write (40,"(8e13.5)") (gswm09_z(i,j,k)*0.00001,
     |       gswm09_t(i,j,k),gswm09_u(i,j,k)*0.01,gswm09_v(i,j,k)*0.01,
     |       k=3,4)
           enddo
          enddo
       endif


! TEMP write-out and save altitudes
      do k=1,nalt9
       call addfld('TIDE_Z','TIDE Z','cm',
     |    gswm09_z(i0:i1,j0:j1,k),'lon',i0,i1,'lat',j0,j1,k)
       call addfld('TIDE_T','TIDE TN','K',
     |    gswm09_t(i0:i1,j0:j1,k),'lon',i0,i1,'lat',j0,j1,k)
       call addfld('TIDE_U','TIDE UN','cm/s',
     |    gswm09_u(i0:i1,j0:j1,k),'lon',i0,i1,'lat',j0,j1,k)
       call addfld('TIDE_V','TIDE VN','cm/s',
     |    gswm09_v(i0:i1,j0:j1,k),'lon',i0,i1,'lat',j0,j1,k)
      enddo

      end subroutine mkgswm09
!-----------------------------------------------------------------------
      subroutine alloc_gswm09(lon0,lon1,lat0,lat1)
!
! Args:
      integer,intent(in) :: lon0,lon1,lat0,lat1
!
! Local:
      integer :: istat

! GSWM09 data in altitude to be use-associated to lbc.F:
      allocate(gswm09_z (lon0:lon1,lat0:lat1,nalt9),stat=istat)
      allocate(gswm09_t (lon0:lon1,lat0:lat1,nalt9),stat=istat)
      allocate(gswm09_u (lon0:lon1,lat0:lat1,nalt9),stat=istat)
      allocate(gswm09_v (lon0:lon1,lat0:lat1,nalt9),stat=istat)
      write(6,"('Allocated GSWM09 t,u,v,z ',
     |  '(lon0:lon1,lat0:lat1,nalt9)')")
!
      end subroutine alloc_gswm09
!-----------------------------------------------------------------------
      end module gswm09_module