#define DBG  print*, pid ,__FILE__, __LINE__
#define IAIC_IDX idx
!-------------------------------------------------------------------------------------------------
  module propagat_mod
  use partctl_mod,only:nb8,plist,lat
  !use wamvar_mod
  use mwmvar_mod
  implicit none
!-------------------------------------------------------------------------------------------------
  public :: init_propagat,final_propagat,propagat_preinf !,out_pgt_info
  public :: propagat !,smoth_ee,update_ee,update_e
  private
!-------------------------------------------------------------------------------------------------
  type pgt_info_type
    !integer :: idx=0
    integer :: iwk=1,iwk1=1,jth=1,jth1=1
    integer :: ids(4)
    real(spdp) :: wts(4),wtp(4)
  end type pgt_info_type
  type(pgt_info_type),allocatable :: pgtinf(:,:,:)
!-------------------------------------------------------------------------------------------------
  contains
!-------------------------------------------------------------------------------------------------
  subroutine init_propagat
    if(.not. allocated(pgtinf))allocate(pgtinf(kl,jl,0:np))
  end subroutine init_propagat
!-------------------------------------------------------------------------------------------------
  subroutine final_propagat
    if(allocated(pgtinf))deallocate(pgtinf)
  end subroutine final_propagat
!-------------------------------------------------------------------------------------------------
  subroutine propagat_preinf(IAIC_IDX)
    integer,intent(in) :: IAIC_IDX
    integer :: j,k,i,iyyz,ixx,jyy,ixx1,jyy1,iwk,iwk1,i1,jth,jth1
    real(spdp) :: d0,dddx0,dddy0,duxdx0,duxdy0,duydy0,duydx0,th0
    real(spdp) :: sinth,costh,wk0,ws0,dk0,cg,cgx,cgy
    real(spdp) :: xx,yy,x1,x2,y1,y2,dsidd
    real(spdp) :: ssr1,ssr2,ssrwk,ssrth,wks,fien,fien1,wk1,wk2,dtth0,ths
    real(spdp) :: e11,e12,e13,e14,e1,e21,e22,e23,e24,e2
    real(spdp) :: e32,e33,e34,e3,e41,e42,e43,e44,e4,e31
    real(spdp) :: th1,th2,exxyy,x_d,y_d
    real(spdp) :: dx,dy
    real(spdp) :: wt4(4)
    integer :: idx1,ids(4)
    if(nsp(IAIC_IDX).ne.1) return
    d0=d(IAIC_IDX)
    dddx0 = cnb(IAIC_IDX)%dddx  !dddx(IAIC_IDX)
    dddy0 = cnb(IAIC_IDX)%dddy  !dddy(IAIC_IDX)
    duxdx0=uxx(IAIC_IDX)
    duxdy0=uxy(IAIC_IDX)
    duydx0=uyx(IAIC_IDX)
    duydy0=uyy(IAIC_IDX)
    do j=1,jl
      th0=thet(j)
      sinth=sin(th0)
      costh=cos(th0)
      !sinth=sin_thet(j) !zyd
      !costh=cos_thet(j) !zyd
      do k=1,kl
        wk0=wk(k)
        ws0=zpi*wf(k,IAIC_IDX)
        dk0=d0*wk0
        cg=ccg(k,IAIC_IDX)
        cgx=cg*costh
        cgy=cg*sinth
        !******  1.  "the calculation of wave engery-current spreading"
        !-----------------------------------------------------------------------------------------
        !----------------------------
        ! Points for interp.
        !  (0,dy) ---- (dx,dy)
        !    |            |
        !    |      X     |
        !    |            |
        !  (0,0) ------(dx,0)
        !----------------------------
        !xx=-deltt*(cgx+ux(IAIC_IDX))/rslat(IAIC_IDX)*180./pi
        xx=-deltt*(cgx+ux(IAIC_IDX))/cnb(IAIC_IDX)%rslat*180./pi
        yy=-deltt*(cgy+uy(IAIC_IDX))/rs*180./pi
        if(xx>=0)then
          if(yy>=0)then  ! Part 1
            ids(3)=nb8(IAIC_IDX)%u;ids(4)=nb8(IAIC_IDX)%ur;ids(1)=IAIC_IDX;ids(2)=nb8(IAIC_IDX)%r
          else           ! Part 4
            ids(3)=IAIC_IDX;ids(4)=nb8(IAIC_IDX)%r;ids(1)=nb8(IAIC_IDX)%d;ids(2)=nb8(IAIC_IDX)%dr
            !yy=delty(IAIC_IDX)+yy
            yy=cnb(IAIC_IDX)%delty+yy
          endif
        else
          if(yy>=0)then  ! Part 2
            ids(3)=nb8(IAIC_IDX)%ul;ids(4)=nb8(IAIC_IDX)%u;ids(1)=nb8(IAIC_IDX)%l;ids(2)=IAIC_IDX
            !xx=deltx(IAIC_IDX)+xx
            xx=cnb(IAIC_IDX)%deltx+xx
          else            ! Part 3
            ids(3)=nb8(IAIC_IDX)%l;ids(4)=IAIC_IDX;ids(1)=nb8(IAIC_IDX)%dl;ids(2)=nb8(IAIC_IDX)%d
            !xx=deltx(IAIC_IDX)+xx;yy=delty(IAIC_IDX)+yy
            xx=cnb(IAIC_IDX)%deltx+xx;yy=cnb(IAIC_IDX)%delty+yy
          endif
        endif
        !dx=deltx(IAIC_IDX);dy=delty(IAIC_IDX)
        dx=cnb(IAIC_IDX)%deltx;dy=cnb(IAIC_IDX)%delty
        !-----------------------------------------------------------------------------------------
        !******  2.  "the effect of refraction caused by topography and current"
        if (dk0.lt.40.) then
          dsidd=0.5*g/cosh(dk0)*wk0**2/ws0/cosh(dk0)
        else
          dsidd=0.
        endif
        !ssr1=(dsidd*dddx0+wk0*costh*duxdx0+wk0*sinth*duydx0)*180./pi
        !ssr2=(dsidd*dddy0+wk0*costh*duxdy0+wk0*sinth*duydy0)*180./pi
        !!ssrwk=-(ssr1*costh/rslat(IAIC_IDX)+ssr2*sinth/rs)
        !!ssrth=(ssr1*sinth/rslat(IAIC_IDX)-ssr2*costh/rs)/wk0
        !ssrwk=-(ssr1*costh/cnb(IAIC_IDX)%rslat+ssr2*sinth/rs)
        !ssrth= (ssr1*sinth/cnb(IAIC_IDX)%rslat-ssr2*costh/rs)/wk0
      ssr1=(dsidd*dddx0+wk0*costh*duxdx0+wk0*sinth*duydx0) !*(180.d0/pi)/rslat(idx)
      ssr2=(dsidd*dddy0+wk0*costh*duxdy0+wk0*sinth*duydy0) !*(180.d0/pi)/rs
      ssrwk=-(ssr1*costh+ssr2*sinth)
      ssrwk=-(ssr1*costh+ssr2*sinth)
        !        ssrth=ssrth-cg*costh*tand(90.-(ic-1)*grid)/rs !?yinxq
        !ssrth=ssrth-cg*costh*tand(y(ic))/rs !?yinxq
        !-ssrth=ssrth-cg*costh*tand(lat(plist(IAIC_IDX)%j))/rs !?yinxq
        ssrth=ssrth-cg*costh*tan(lat(plist(IAIC_IDX)%j)*pi/180.)/rs !?yinxq
        wks=wk0-deltt*ssrwk
        if (wks.lt.0.) wks=0.
        if (wks.le.wkmin) then
          iwk=1
          iwk1=1
          fien=0.
          fien1=1.
          wk1=0.
          wk2=wk(iwk1)
        else
          if (wks.lt.wk(kld)) then
            !======= THE FOLLOWING BELONGS TO THE LAGFD-WAM WAVE MODEL
            !            iwk=int(log(wks/wkmin)/log(pwk))+1
            !==========================================================
            !            if (wks.lt.wk(iwk)) iwk=iwk-1
            !            if (wks.gt.wk(iwk+1)) iwk=iwk+1
            !========   The following is modified by Yang Yongzeng =====c
            do iyyz=1,kld
              if(wks>=wk(iyyz) .and. wks<wk(iyyz+1))iwk=iyyz
            enddo
            iwk1=iwk+1
            !=====================================================================c
            wk1=wk(iwk)
            wk2=wk(iwk+1)
            if (iwk.lt.kl) then
              iwk1=iwk+1
              fien=1.
              fien1=1.
            else
              i=iwk-kl+1
              i1=i+1
              fien=wkh(i)
              fien1=wkh(i1)
              iwk=kl
              iwk1=kl
            endif
          else
            wks=wk(kld)
            iwk=kl
            iwk1=kl
            i=kld-kl+1
            i1=i+1
            fien=wkh(i)
            fien1=wkh(i1)
          endif
        endif
        dtth0=deltt*ssrth
        ths=th0-dtth0
        if (ths<-1.0*zpi .or. ths>=2.0*zpi) ths=th0
        if (ths.ge.zpi) ths=ths-zpi
        if (ths.lt. 0.) ths=ths+zpi
        jth=int(ths/deltth)+1;jth1=jth+1
        if(jth.eq.jlp1)then
          jth=jl;jth1=1
        endif
        if (jth1.eq.jlp1) jth1=1
        th1=thet(jth);th2=thet(jth+1)
        !-----------------------------------------------------------------------------------------
!        pgtinf(k,j,IAIC_IDX)%idx=idx1
        pgtinf(k,j,IAIC_IDX)%iwk=iwk
        pgtinf(k,j,IAIC_IDX)%iwk1=iwk1
        pgtinf(k,j,IAIC_IDX)%jth=jth
        pgtinf(k,j,IAIC_IDX)%jth1=jth1
        call inter_wt(0._spdp,dx,0._spdp,dy,xx,yy,wt4)  ! record: idx1,wt4
        pgtinf(k,j,IAIC_IDX)%ids=ids
        pgtinf(k,j,IAIC_IDX)%wts=wt4
        call inter_wt(wk1,wk2,th1,th2,wks,ths,wt4) ! record: iwk,iwk1,jth,jth1,wt4
        wt4(1)=wt4(1)*fien;wt4(2)=wt4(2)*fien1;wt4(3)=wt4(3)*fien;wt4(4)=wt4(4)*fien1
        pgtinf(k,j,IAIC_IDX)%wtp=wt4
        !-----------------------------------------------------------------------------------------
        !          !******  3.  "determing the wave energy at the physical space point (xx,yy)
        !          !******                  and the wave space point (wks,ths)"
        !          !*       3.1 "at physical space point (x1,y1)"
        !          e11=ee(iwk ,jth ,ixx,jyy)*fien
        !          e12=ee(iwk1,jth ,ixx,jyy)*fien1
        !          e13=ee(iwk ,jth1,ixx,jyy)*fien
        !          e14=ee(iwk1,jth1,ixx,jyy)*fien1
        !          call inter(wk1,wk2,th1,th2,e11,e12,e13,e14,wks,ths,e1)
        !          !*       3.2 "at physical space point (x2,y1)"
        !          e21=ee(iwk ,jth ,ixx1,jyy)*fien
        !          e22=ee(iwk1,jth ,ixx1,jyy)*fien1
        !          e23=ee(iwk ,jth1,ixx1,jyy)*fien
        !          e24=ee(iwk1,jth1,ixx1,jyy)*fien1
        !          call inter(wk1,wk2,th1,th2,e21,e22,e23,e24,wks,ths,e2)
        !          !*       3.3 "at physical space point (x1,y2)"
        !          e31=ee(iwk ,jth ,ixx,jyy1)*fien
        !          e32=ee(iwk1,jth ,ixx,jyy1)*fien1
        !          e33=ee(iwk ,jth1,ixx,jyy1)*fien
        !          e34=ee(iwk1,jth1,ixx,jyy1)*fien1
        !          call inter(wk1,wk2,th1,th2,e31,e32,e33,e34,wks,ths,e3)
        !          !*       3.4 "at physical space point (x2,y2)"
        !          e41=ee(iwk ,jth ,ixx1,jyy1)*fien
        !          e42=ee(iwk1,jth ,ixx1,jyy1)*fien1
        !          e43=ee(iwk ,jth1,ixx1,jyy1)*fien
        !          e44=ee(iwk1,jth1,ixx1,jyy1)*fien1
        !          call inter(wk1,wk2,th1,th2,e41,e42,e43,e44,wks,ths,e4)
        !          !*       3.5 "at physical space point (xx,yy)"
        !          call inter(x1,x2,y1,y2,e1,e2,e3,e4,xx,yy,exxyy)
        !          !        e(k,j,IAIC_IDX)=exxyy
        !          e(k,j,IAIC_IDX)=max(exxyy,small)
        !          !        write(6,*) e(k,j,IAIC_IDX)
      enddo
    enddo
  end subroutine propagat_preinf
!-------------------------------------------------------------------------------------------------
  subroutine inter(u1,u2,v1,v2,aa,bb,cc,dd,xt,yt,value)
    real(spdp), intent(in) :: u1,u2,v1,v2,aa,bb,cc,dd,xt,yt
    real(spdp), intent(out) :: value
    real(spdp) :: dox, doy, yvalue, r, q, phia, phib, phic, phid
    dox=u2-u1
    doy=v2-v1
    if(dox.eq.0.0.and.doy.eq.0.0)then
      yvalue=(aa+bb+cc+dd)/4.0
    endif
    if(dox.eq.0.0.and.doy.ne.0.0)then
      r=(yt-v1)/doy
      yvalue=r*(cc+dd)/2.0+(1-r)*(aa+bb)/2.0
    endif
    if(dox.ne.0.0.and.doy.eq.0.0)then
      q=(xt-u1)/dox
      yvalue=q*(bb+dd)/2.0+(1-q)*(aa+cc)/2.0
    endif
    if(dox.ne.0.0.and.doy.ne.0.0)then
      q=(xt-u1)/dox
      r=(yt-v1)/doy
      phia=(1-q)*(1-r)
      phib=q*(1-r)
      phic=(1-q)*r
      phid=q*r
      yvalue=aa*phia
      yvalue=yvalue+bb*phib
      yvalue=yvalue+cc*phic
      yvalue=yvalue+dd*phid
    endif
    value=yvalue
  end subroutine inter
!-------------------------------------------------------------------------------------------------
  subroutine inter_wt(u1,u2,v1,v2,xt,yt,wt4)
    real(spdp), intent(in) :: u1,u2,v1,v2,xt,yt
    real(spdp), intent(out) :: wt4(4)
    real(spdp) :: dox, doy, r, q, phia, phib, phic, phid
    dox=u2-u1
    doy=v2-v1
    if(dox.eq.0.0.and.doy.eq.0.0)then
      !yvalue=(aa+bb+cc+dd)/4.0
      wt4=1.d0/4.d0
    endif
    if(dox.eq.0.0.and.doy.ne.0.0)then
      r=(yt-v1)/doy
      !yvalue=r*(cc+dd)/2.0+(1-r)*(aa+bb)/2.0
      wt4(3:4)=r/2.d0;wt4(1:2)=(1.d0-r)/2.d0
    endif
    if(dox.ne.0.0.and.doy.eq.0.0)then
      q=(xt-u1)/dox
      !yvalue=q*(bb+dd)/2.0+(1-q)*(aa+cc)/2.0
      wt4(2)=q/2.d0;wt4(4)=q/2.d0
      wt4(1)=(1.d0-q)/2.d0;wt4(3)=(1.d0-q)/2.d0
    endif
    if(dox.ne.0.0.and.doy.ne.0.0)then
      q=(xt-u1)/dox
      r=(yt-v1)/doy
      phia=(1-q)*(1-r)
      phib=q*(1-r)
      phic=(1-q)*r
      phid=q*r
      !yvalue=aa*phia
      !yvalue=yvalue+bb*phib
      !yvalue=yvalue+cc*phic
      !yvalue=yvalue+dd*phid
      wt4(1)=phia
      wt4(2)=phib
      wt4(3)=phic
      wt4(4)=phid
    endif
  end subroutine inter_wt
!-------------------------------------------------------------------------------------------------
!  subroutine out_pgt_info2
!    use netcdf_mod
!    integer :: IAIC_IDX
!    real(4),allocatable :: itmp1(:),itmp2d(:,:)
!    integer :: irec,ncid,ii,jj
!    allocate(itmp1(0:np),itmp2d(im,jm))
!    root=0
!    if(pid==root)then
!      call open_nc(ncid,'pgtinfo.nc','c')
!      call dimension_define(ncid,'lon',im,'lon',nf_real)
!      call dimension_define(ncid,'lat',jm,'lat',nf_real)
!      call dimension_define(ncid,'rec',0,'rec',nf_int)
!      call variable_define(ncid,'i1',nf_real,['lon','lat','rec'])
!      call set_attribute(ncid,'units','degrees_north','lat')
!      call set_attribute(ncid,'units','degrees_east','lon')
!      call set_attribute(ncid,'modulo','','lon')
!      call end_define(ncid)
!      call writenc(ncid,'lon',lon)
!      call writenc(ncid,'lat',lat)
!    endif
!    do j=1,jl
!      do k=1,kl        
!      	irec=k+(j-1)*kl
!      	if(pid==root)write(*,*)'out_pgt_info:',irec,k,j
!		    do IAIC_IDX=1,npc
!          !itmp1(IAIC_IDX)=pgtinf(k,j,IAIC_IDX)%wts(1)
!          itmp1(IAIC_IDX)=pgtinf(k,j,IAIC_IDX)%wts(2)
!	      enddo
!	      call irrp_gather(itmp1,itmp2d,root)
!	      if(pid==root)then
!			    do jj=1,jm
!			      do ii=1,im
!			        if(mask(ii,jj)==0)then
!			          itmp2d(ii,jj)=0
!			        elseif(abs(itmp2d(ii,jj))>1000)then
!			        	write(*,*)'bebe',ii,jj,itmp2d(ii,jj)
!			        endif
!			      enddo
!			    enddo
!	      	call writenc(ncid,'i1',itmp2d,irec)
!	      endif
!	    enddo
!	  enddo
!	  if(pid==root)call close_nc(ncid)
!	  deallocate(itmp1,itmp2d)	        
!  end subroutine out_pgt_info2
!!-------------------------------------------------------------------------------------------------
!  subroutine out_pgt_info1
!    use netcdf_mod
!    integer :: IAIC_IDX
!    integer,allocatable :: itmp1(:),itmp2d(:,:)
!    integer :: irec,ncid,ii,jj
!    allocate(itmp1(0:np),itmp2d(im,jm))
!    root=0
!    if(pid==root)then
!      call open_nc(ncid,'pgtinfo.nc','c')
!      call dimension_define(ncid,'lon',im,'lon',nf_real)
!      call dimension_define(ncid,'lat',jm,'lat',nf_real)
!      call dimension_define(ncid,'rec',0,'rec',nf_int)
!      call variable_define(ncid,'i1',nf_int,['lon','lat','rec'])
!      call set_attribute(ncid,'units','degrees_north','lat')
!      call set_attribute(ncid,'units','degrees_east','lon')
!      call set_attribute(ncid,'modulo','','lon')
!      call end_define(ncid)
!      call writenc(ncid,'lon',lon)
!      call writenc(ncid,'lat',lat)
!    endif
!    do j=1,jl
!      do k=1,kl        
!      	irec=k+(j-1)*kl
!      	if(pid==root)write(*,*)'out_pgt_info:',irec,k,j
!		    do IAIC_IDX=1,npc
!          ii=pgtinf(k,j,IAIC_IDX)%ids(1)
!          itmp1(IAIC_IDX)=0;
!          if(ii/=0)then
!          	itmp1(IAIC_IDX)=plist(ii)%i+plist(ii)%j
!	          if(abs(plist(ii)%i)>400)then
!	          	write(*,*)'error: ',plist(ii)%i,pgtinf(k,j,IAIC_IDX)%ids(1),npc
!	          endif
!	        endif
!	      enddo
!	      call irrp_gather(itmp1,itmp2d,root)
!	      if(pid==root)then
!			    do jj=1,jm
!			      do ii=1,im
!			        if(mask(ii,jj)==0)then
!			          itmp2d(ii,jj)=0
!			        elseif(abs(itmp2d(ii,jj))>1000)then
!			        	write(*,*)'bebe',ii,jj,itmp2d(ii,jj)
!			        endif
!			      enddo
!			    enddo
!	      	call writenc(ncid,'i1',itmp2d,irec)
!	      endif
!	    enddo
!	  enddo
!	  if(pid==root)call close_nc(ncid)
!	  deallocate(itmp1,itmp2d)	        
!  end subroutine out_pgt_info1
!  subroutine out_pgt_info
!    use netcdf_mod
!    integer :: IAIC_IDX
!    integer,allocatable :: itmp1(:),itmp2d(:,:)
!    integer :: irec,ncid,ii,jj
!    allocate(itmp1(0:np),itmp2d(im,jm))
!    root=0
!    if(pid==root)then
!      call open_nc(ncid,'pgtinfo.nc','c')
!      call dimension_define(ncid,'lon',im,'lon',nf_real)
!      call dimension_define(ncid,'lat',jm,'lat',nf_real)
!      call dimension_define(ncid,'rec',0,'rec',nf_int)
!      call variable_define(ncid,'i1',nf_int,['lon','lat','rec'])
!      call set_attribute(ncid,'units','degrees_north','lat')
!      call set_attribute(ncid,'units','degrees_east','lon')
!      call set_attribute(ncid,'modulo','','lon')
!      call end_define(ncid)
!      call writenc(ncid,'lon',lon)
!      call writenc(ncid,'lat',lat)
!    endif
!    do j=1,jl
!      do k=1,kl        
!      	irec=k+(j-1)*kl
!      	if(pid==root)write(*,*)'out_pgt_info:',irec,k,j
!		    do IAIC_IDX=1,npc
!        	itmp1(IAIC_IDX)=pgtinf(k,j,IAIC_IDX)%iwk
!	      enddo
!	      call irrp_gather(itmp1,itmp2d,root)
!	      if(pid==root)then
!			    do jj=1,jm
!			      do ii=1,im
!			        if(mask(ii,jj)==0)then
!			          itmp2d(ii,jj)=0
!			        elseif(abs(itmp2d(ii,jj))>1000)then
!			        	write(*,*)'bebe',ii,jj,itmp2d(ii,jj)
!			        endif
!			      enddo
!			    enddo
!	      	call writenc(ncid,'i1',itmp2d,irec)
!	      endif
!	    enddo
!	  enddo
!	  if(pid==root)call close_nc(ncid)
!	  deallocate(itmp1,itmp2d)	        
!  end subroutine out_pgt_info
!-------------------------------------------------------------------------------------------------
  subroutine propagat(IAIC_IDX)
    integer,intent(in) :: IAIC_IDX
    integer :: idx1,iwk,iwk1,jth,jth1,idx0
    real(spdp) :: exxyy
    real(spdp) :: wtp(4),wts(4)
    do j=1,jl
      do k=1,kl
        iwk =pgtinf(k,j,IAIC_IDX)%iwk;iwk1=pgtinf(k,j,IAIC_IDX)%iwk1
        jth =pgtinf(k,j,IAIC_IDX)%jth;jth1=pgtinf(k,j,IAIC_IDX)%jth1
        wtp=pgtinf(k,j,IAIC_IDX)%wtp;wts=pgtinf(k,j,IAIC_IDX)%wts
        exxyy=0.d0
        idx1=pgtinf(k,j,IAIC_IDX)%ids(1)
        exxyy=exxyy+ee(iwk ,jth ,idx1)*wtp(1)*wts(1)+ee(iwk1,jth ,idx1)*wtp(2)*wts(1) &
                   +ee(iwk ,jth1,idx1)*wtp(3)*wts(1)+ee(iwk1,jth1,idx1)*wtp(4)*wts(1)
        !idx1=nb8(idx0)%r
        idx1=pgtinf(k,j,IAIC_IDX)%ids(2)
        exxyy=exxyy+ee(iwk ,jth ,idx1)*wtp(1)*wts(2)+ee(iwk1,jth ,idx1)*wtp(2)*wts(2) &
                   +ee(iwk ,jth1,idx1)*wtp(3)*wts(2)+ee(iwk1,jth1,idx1)*wtp(4)*wts(2)
        !idx1=nb8(idx0)%u
        idx1=pgtinf(k,j,IAIC_IDX)%ids(3)
        exxyy=exxyy+ee(iwk ,jth ,idx1)*wtp(1)*wts(3)+ee(iwk1,jth ,idx1)*wtp(2)*wts(3) &
                   +ee(iwk ,jth1,idx1)*wtp(3)*wts(3)+ee(iwk1,jth1,idx1)*wtp(4)*wts(3)
        !idx1=nb8(idx0)%ur
        idx1=pgtinf(k,j,IAIC_IDX)%ids(4)
        exxyy=exxyy+ee(iwk ,jth ,idx1)*wtp(1)*wts(4)+ee(iwk1,jth ,idx1)*wtp(2)*wts(4) &
                   +ee(iwk ,jth1,idx1)*wtp(3)*wts(4)+ee(iwk1,jth1,idx1)*wtp(4)*wts(4)
        e(k,j,IAIC_IDX)=max(exxyy,small)
      enddo
    enddo
  end subroutine propagat
!-------------------------------------------------------------------------------------------------
  subroutine smoth_ee
    real(spdp),parameter :: a=24.0
    real(spdp) :: n
    integer :: idx
    real(spdp),allocatable :: em(:,:,:)
    allocate(em(kl,jl,0:np));em=0.0
    do idx=1,npc
      em(:,:,idx)=a*ee(:,:,idx);n=a
      if(nb8(idx)%l/=0)then
        em(:,:,idx)=em(:,:,idx)+ee(:,:,nb8(idx)%l);n=n+1
      endif
      if(nb8(idx)%r/=0)then
        em(:,:,idx)=em(:,:,idx)+ee(:,:,nb8(idx)%r);n=n+1
      endif
      if(nb8(idx)%u/=0)then
        em(:,:,idx)=em(:,:,idx)+ee(:,:,nb8(idx)%u);n=n+1
      endif
      if(nb8(idx)%d/=0)then
        em(:,:,idx)=em(:,:,idx)+ee(:,:,nb8(idx)%d);n=n+1
      endif
      em(:,:,idx)=max(small,em(:,:,idx)/n)
    enddo
    ee=em
    deallocate(em)
  end subroutine smoth_ee
!-------------------------------------------------------------------------------------------------
!  subroutine update_ee
!    integer,save :: group_ind_ee
!    integer :: ivar,kk,jj
!    if(group_ind_ee==0)then
!      ivar=1;group_ind_ee=1;kk=kl;jj=jl
!      call irrp_exg_setvar(group_ind_ee,ivar,'ee',ee,kk,jj)
!    endif
!    call irrp_exg_action(group_ind_ee)
!    ee(:,:,0)=small
!  end subroutine update_ee
!-------------------------------------------------------------------------------------------------
!  subroutine update_e
!    integer,save :: group_ind_e
!    integer :: ivar,kk,jj
!    if(group_ind_e==0)then
!      ivar=1;group_ind_e=1;kk=kl;jj=jl
!      call irrp_exg_setvar(group_ind_e,ivar,'e',e,kk,jj)
!    endif
!    call irrp_exg_action(group_ind_e)
!    e(:,:,0)=small
!  end subroutine update_e
!-------------------------------------------------------------------------------------------------
  end module propagat_mod
!-------------------------------------------------------------------------------------------------
