module sgwf_module
      use params_module,only: nlonp4,nlon,nlat,glon1,dlon
      use init_module,only: glat
      use cons_module,only: dtr,grav
      use cons_module,only: pi
      use init_module,only: iday,uthr
      implicit none

      real, parameter ::
     |      gwsfmax = 1.33353 ,         ! m/s/s
     |      lonpeak = -53.1982,          ! -53 degrees
     |      latpeak = -17.7273,         ! -17 degrees
     |      zpeak   = 1.80e7,       ! in cm, 180 km
     |      parw  = 8.90471e4,   ! in m, 4.5*zonalw=400km (parallel direction)
     |      perw  = 1.41854e5,   ! in m, 4.5*meridw=638km (perpendicular direction)
     |      zw      = 1.918e6,   ! in cm, 4.5*zw=54.6896km
     |      pdir = -89.6312,    ! propagation direction (from east, in degree)
     |      timebeguth = 22.4025,       ! UT time the forcing begins (hour)
     |      timeenduth = 23.0975         ! UT time the forcing ends
      integer, parameter ::
     |      doy = 274                    ! Day of year of occurence

      contains

      subroutine secondary_gwf(z,lev0,lev1,lon0,lon1,lat,gwsfx,gwsfy)

      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) :: z
      real,dimension(lev0:lev1,lon0:lon1),intent(out):: gwsfx,gwsfy

      integer :: k,i,i0,i1,nk,nkm1,nlevs,nlons
      real :: tfac,zfac,fdir,zfacexp
      real :: wrkx,wrky,wrk1,wrk2,facpar,facper
      real, dimension(lon0:lon1):: gwsfx0, gwsfy0
      real, dimension(nlonp4)::glonp4                ! glonp4(3)=-180, glonp4(1:2) and glonp4(nlonp4-1:nlonp4) are periodic grids

      i0 = lon0
      i1 = lon1

      tfac = 0.
      zfac = 0.
      gwsfx(:,:) = 0.
      gwsfy(:,:) = 0.

      do i=3,nlonp4-2
         glonp4(i) = glon1+(i-3)*dlon
      enddo
      glonp4(1)=glonp4(nlonp4-3)
      glonp4(2)=glonp4(nlonp4-2)
      glonp4(nlonp4-1)=glonp4(3)
      glonp4(nlonp4)=glonp4(4)

      if (iday .eq. doy) then
         if ((uthr.gt.timebeguth).and.(uthr.lt.timeenduth)) then
            fdir = pdir*dtr     ! forcing direction, in radian
            
            wrky = (glat(lat)-latpeak)/180.*pi*6.37e6
            do i=i0,i1
               wrkx = (glonp4(i)-lonpeak)/360.*2.*pi*6.37e6
     |              * cos(glat(lat)*dtr)
               wrk1 = wrkx*cos(fdir)+wrky*sin(fdir)
               wrk2 = -wrkx*sin(fdir)+wrky*cos(fdir)
               facpar = exp(-(wrk1/parw)**2./2.)
               facper = exp(-(wrk2/perw)**2./2.)
               gwsfx0(i) = gwsfmax * facpar * facper
     |              * cos(fdir)
               gwsfy0(i) = gwsfmax * facpar * facper
     |               * sin(fdir)
            enddo

            tfac = sin(pi*(uthr-timebeguth)/(timeenduth-timebeguth))**2.
            do i=i0,i1
               do k=lev0,lev1
                  zfacexp = -((z(k,i)-zpeak)/zw)**2./2.
                  if (zfacexp .gt. -10) then
                     zfac = exp(zfacexp)
                  else
                     zfac = 0
                  endif
                  gwsfx(k,i) = gwsfx0(i) * tfac * zfac ! zonal wind tendency
                  gwsfy(k,i) = gwsfy0(i) * tfac * zfac ! meridional wind tendency
               enddo
            enddo
         endif
      endif

      return
      end subroutine secondary_gwf
      end module sgwf_module