
      MODULE radtrans
! -------------------------------------------------------------------- !
! Purpose: This module calculates radiation transport
! -------------------------------------------------------------------- !
      use communicate
      IMPLICIT NONE
      REAL(8) :: c, q2
      PUBLIC radtransport
      CONTAINS

! -------------------------------------------------------------------- !
      SUBROUTINE radtransport
! -------------------------------------------------------------------- !
! Purpose:  Solve the radiation transport equation to determine the
!           intensity spectrum.
! -------------------------------------------------------------------- !
      USE grt 
      
      IMPLICIT NONE
      INTEGER(4), PARAMETER :: np = 5000
      INTEGER(4) :: max_theory, max_exp
      INTEGER(4) :: i,j,iagain ,io 
      INTEGER(4) :: ionce    ! change 1->0 for changing hydro mesh
      INTEGER(4) :: unitnum
      INTEGER(4) :: max
      INTEGER(4) :: ilwe,iupe,ilwt,iupt      
! Interpolating
      REAL(8), DIMENSION(np) :: yn1

      REAL(8) :: cfwhm
!     REAL(8) :: c, q2
      REAL(8) , DIMENSION(SG)      :: rlambda,rintsy
      REAL(8) , DIMENSION(SG)      :: xexp,fexp  ! Exp. spectral data for q-squr comp.
      REAL(8) , DIMENSION(SG)      :: emax,emin,optical_dep
      REAL(8) , DIMENSION(SG)      :: eout,riout
      REAL(8) , DIMENSION(MZ-1)    :: dzc,rdzc1
      REAL(8) , DIMENSION(MZ)      :: rdx1
      REAL(8) , DIMENSION(SG,MZ)   :: part1, part2, parts
      REAL(8) , DIMENSION(SG,MZ-1) :: rdzc2
      REAL(8) , DIMENSION(SG,MZ)   :: rdx2, dx2
      REAL(8) , DIMENSION(SG,0:MZ) :: hbdry2   ! Hydro bdry. data for all freq
      REAL(8) , DIMENSION(SG,0:MZ) :: j_y      ! Emis. eval. @ hydro boundries
      REAL(8) , DIMENSION(SG,0:MZ) :: j_slope  ! Slope of line seg. btwn zone centers
      REAL(8) , DIMENSION(SG,MZ)   :: j_slopezb  ! Slope of line seg. btwn zone boundries
      REAL(8) , DIMENSION(SG,MZ)   :: j_yintzb   ! Y-inters. of line seg btwn zone bdrys
      REAL(8) , DIMENSION(SG,MZ)   :: a1         ! Aux. array
      REAL(8) , DIMENSION(SG,MZ)   :: betax0
      REAL(8) , DIMENSION(SG,MZ)   :: betax1
      REAL(8) , DIMENSION(SG,MZ)   :: ttau
      REAL(8) , DIMENSION(SG,MZ)   :: betanegdx  
      REAL(8) , DIMENSION(SG,MZ)   :: exp_alpha  ! Exp[tau%yint]
      REAL(8) , DIMENSION(SG,MZ)   :: exp_ttau
      REAL(8) , DIMENSION(SG,MZ)   :: exp_betanegdx  
      character(len=70) :: envstring
      CHARACTER(len=30) :: ctemp
!     CHARACTER(len=100), EXTERNAL :: dir
      CHARACTER(len=3)  :: cagain
      external fbuff
       
      INTERFACE msexp
         FUNCTION msexp(x)
           IMPLICIT NONE
           REAL(8), DIMENSION(:,:), INTENT(IN) :: x
           REAL(8), DIMENSION(SIZE(x,1),SIZE(x,2)) :: msexp
         END FUNCTION msexp
      END INTERFACE

      INTERFACE qsqrd
         SUBROUTINE qsqrd(x,fexp,max_exp,fthry,c,q2)
           IMPLICIT NONE
           REAL(8), INTENT(IN), DIMENSION(:) :: x  
           REAL(8), INTENT(IN), DIMENSION(:) :: fexp
           INTEGER(4),INTENT(IN) :: max_exp
           REAL(8), INTENT(IN), DIMENSION(:) :: fthry
           REAL(8), INTENT(OUT) :: c
           REAL(8), INTENT(OUT) :: q2
         END SUBROUTINE qsqrd
      END INTERFACE

      DATA ionce  /1/
      DATA iagain /0/
      SAVE ionce, iagain, max_exp, xexp, fexp
      iagain=iagain+1
! -- Reading Exp. Data
      if(ionce.eq.1) then
        envstring = 'WORKDIR'
        unitnum = 36
!       ctemp = 'optdep'//trim(adjustl(get_cmytid()))//"."//cagain
        ctemp = 'expspec.dat'
        call openf(unitnum,ctemp,envstring)                   ! OPEN(unit=36,file=expspec.dat)
        i=1
        do
          read(36,*,end=80) xexp(i),fexp(i)
!         print*,'EXP ',xexp(i),fexp(i)
          i=i+1
        enddo
80      close(36)
        max_exp=i-1
      endif  
!     IF(ionce.EQ.1) THEN
        dx(1:MZ) = hbdry(1:MZ) - hbdry(0:MZ-1)
        dzc(1:MZ-1) = opacity%zc(2:MZ) - opacity%zc(1:MZ-1) 
        rdzc1   = 1.0/dzc  
        rdx1   = 1.0/dx
        rdx2   = SPREAD(rdx1,1,SG)
        dx2    = SPREAD(dx,1,SG)
        hbdry2 = SPREAD(hbdry,1,SG)
        rdzc2  = SPREAD(rdzc1,1,SG)
        ionce = 0
!     ENDIF  NEEDS TO BE MODIFIED

! -- Determine line segments for partial optical depth from trapazoital integration

      a1(:,:)=opacity%a2d(:,:)*dx2(:,:)    ! trap integ.
    
      tau%y(:,MZ)=0.0
!      tau%y(:,MZ-1:0:-1)=a1(:,MZ:1:-1)+tau%y(:,MZ:1:-1)  
      DO i=1,MZ
         tau%y(:,MZ-i)=a1(:,MZ+1-i)+tau%y(:,MZ+1-i)         
      ENDDO         
        
      tau%slope(:,1:MZ)=((tau%y(:,1:MZ)-tau%y(:,0:MZ-1))*rdx2(:,1:MZ))
      tau%yint(:,:)=(tau%y(:,1:MZ)-tau%slope(:,:)*hbdry2(:,1:MZ))
      tau%slope(:,:)=-1.0d0*tau%slope(:,:)
      tau%yint(:,:)=-1.0d0*tau%yint(:,:)

! -- Determine emissivity line segments within hydro zones from zone center data.
      emax(:)= emissivity%a2d(:,1)
      emin(:)= emissivity%a2d(:,MZ)
      j_slope(:,0)  = (emissivity%a2d(:,1) - emin)/(emissivity%zc(1)-hbdry(0))
      j_slope(:,MZ) = (emax - emissivity%a2d(:,MZ))/(hbdry(MZ)-emissivity%zc(MZ))
      j_slope(:,1:MZ-1)=(emissivity%a2d(:,2:MZ)-emissivity%a2d(:,1:MZ-1))*rdzc2(:,1:MZ-1)
      j_y(:,0) = emin
      j_y(:,MZ) = emax
!      j_y(:,1:MZ-1) = emissivity%a2d(:,1:MZ-1) + j_slope(:,1:MZ-1)*&
!     &                (hbdry2(:,1:MZ-1)-emissivity%zc(1:MZ-1))      
      DO i=1,MZ-1 
         j_y(:,i) = emissivity%a2d(:,i) + j_slope(:,i)*(hbdry2(:,i)-emissivity%zc(i))
      ENDDO
      j_slopezb(:,1:MZ) = (j_y(:,1:MZ)-j_y(:,0:MZ-1))*rdx2(:,1:MZ)
      j_yintzb(:,:)  = j_y(:,1:MZ)-j_slopezb(:,:)*hbdry2(:,1:MZ)

! -- Preparing data for piecewise analytical solution
      
      betax0(:,:) = tau%slope(:,:)*hbdry2(:,0:MZ-1)
      betax1(:,:) = tau%slope(:,:)*hbdry2(:,1:MZ)
      betanegdx(:,:) = tau%slope(:,:)*(-dx2(:,:))
      ttau = tau%yint + betax1
      exp_alpha  = msexp(tau%yint)  ! msexp -> EXP with underflow correction 
      exp_betanegdx = msexp(betanegdx)
      exp_ttau = msexp(ttau)

      optical_dep(:)=-1.0d0*tau%yint(:,1)

     WHERE( ABS(tau%slope).GT.1.0e-10 )
!        part1=(j_yintzb*exp_alpha/tau%slope)*(exp_betax1-exp_betax0)
!        part2=(j_slopezb*exp_alpha/tau%slope**2)*(exp_betax1*(betax1-1.0)-exp_betax0*(betax0-1.0))
       part1=(j_yintzb*exp_ttau/tau%slope)*(1.0d0-exp_betanegdx)
       part2=((j_slopezb*exp_ttau/(tau%slope**2))*((betax1-1.0d0)-(exp_betanegdx*(betax0-1.0d0)))) 
     ELSEWHERE
       part1=j_yintzb*exp_alpha*dx2
       part2=0.0
     ENDWHERE
      parts=part1+part2
      rint(:)=SUM((parts),2)

! **************** TESTING *************** !
!        OPEN(40,FILE='/home/manolo/SRMAC/opt_dep.out')
!         OPEN(43,FILE='/home/manolo/SRMAC/tauslope')
!        OPEN(44,FILE='/home/manolo/SRMAC/hbdry2')
!        DO j=1,SG
!            WRITE(40,903) xspec(j), optical_dep(j)
!            write(43,903) xspec(j), (tau%slope(j,i),i=1,mz)
!            write(44,903) xspec(j), (hbdry2(j,i),i=1,mz)
!        enddo
!        close(42)
! **************** TESTING *************** !

         write(cagain,'(i3.3)') iagain

! ---- Writting data for table -----------------------------------------------
         envstring = 'SPECDIR'
         unitnum = 40
         ctemp = 'opd'//trim(adjustl(get_cmytid()))//"."//cagain
         call openf(unitnum,ctemp,envstring)    ! OPEN(unit=40,file=opd..)
         WRITE(40,*) '# ',r8_frm_p(1),' ',r8_frm_p(2),' ',r8_frm_c(1)
         DO j=1,nsize
!         write(40,901)  (12398.0/xspec(j)), optical_dep(j)
!!            write(40,'(3(2x,e16.9))')  xspec(j), opacity%a2d(j,1), emissivity%a2d(j,1)  
         ENDDO    
         close(40)
! -----------------------------------------------------------------------------         
!        write(cagain,'(i3.3)') iagain
!        c_aux = 'spec.out'//cagain
!        c_aux = dir(workdir,c_aux)

         cfwhm=0.0040       ! -- FWHM of line profile(instumental)
         call CONVOL1(xspec,rint,nsize,eout,riout,cfwhm,sg)

         envstring = 'SPECDIR'
         unitnum = 39
         ctemp = 'spec'//trim(adjustl(get_cmytid()))//"."//cagain
         call openf(unitnum,ctemp,envstring)    ! OPEN(unit=39,file=spec.out..)
         rlambda = 0.0d0
         rintsy  = 0.0d0
         io=1
         WRITE(39,*) '# ',r8_frm_p(1),' ',r8_frm_p(2),' ',r8_frm_c(1)
         DO j=nsize,1,-1
            rlambda(io) = 12398.0/eout(j)
            rintsy(io)  = riout(j)
!           write(39,901)  rlambda(io), rintsy(io)
            io=io+1
         ENDDO     
!        close(39)
                  

         max_theory=nsize
         call interpolate(SG,max_theory,rlambda,rintsy,max_exp,xexp,yn1)

!        call qsqrd(xexp,fexp,max_exp,yn1,c,q2)
!        write(1,*) 'c=',c
!        write(1,*) 'q2',q2
!        yn1=yn1*c  ! rescaling results
!        write(39,*) '# q2=',q2,' c=',c
!        # no-rescaling of results see radtrans.f90'
         do i=1,max_exp
!           write(39,901) xexp(i),yn1(i)
!           write(39,901) xexp(i),yn1(i)   !printing yn1 from interpolate
         enddo
         close(39)

902   FORMAT(10(2x,e14.7))
901   FORMAT(1x,e14.7,2x,e14.7)
      END SUBROUTINE radtransport

! *******************************************************************
      END MODULE radtrans


! *******************************************************************
      FUNCTION msexp(x)
! *******************************************************************
! Calcuates e^x and handels underflows
      IMPLICIT NONE
      REAL(8),    DIMENSION(:,:), INTENT(IN)     ::  x
      REAL(8),    DIMENSION(SIZE(x,1),SIZE(x,2)) ::  msexp
      INTEGER(4), DIMENSION(SIZE(x,1),SIZE(x,2)) ::  isize

      isize = -730 
      WHERE(x.le.isize)
        msexp=0.0d0
      ELSEWHERE
        msexp=EXP(x)
      ENDWHERE

      END FUNCTION msexp
! *******************************************************************
! --------------------------------------------------------------- !
    FUNCTION  get_mask(xarr)
! --------------------------------------------------------------- !
    IMPLICIT NONE

    REAL(8), INTENT(IN), DIMENSION(:) :: xarr 
    INTEGER(4), DIMENSION(SIZE(xarr,1)) :: get_mask

    INTEGER(4), PARAMETER :: reg_max = 4
    INTEGER(4) :: i
    TYPE region
       REAL(8) :: xmin
       REAL(8) :: xmax 
    END TYPE region   
    TYPE(region) :: regions(reg_max)

    get_mask=0                       ! Deny all
    regions(1)%xmin=5166
    regions(1)%xmax=5246
    
    regions(2)%xmin=5407
    regions(2)%xmax=5514

    regions(3)%xmin=5975
    regions(3)%xmax=6300

    regions(4)%xmin=6630
    regions(4)%xmax=6780

    DO i=1,reg_max
       WHERE((xarr.GE.regions(i)%xmin).AND.(xarr.LE.regions(i)%xmax))
          get_mask=1 
       ENDWHERE
    ENDDO 
    END FUNCTION get_mask
! --------------------------------------------------------------- !
    SUBROUTINE qsqrd(x,fexp,max_exp,fthry,c,q2)
! --------------------------------------------------------------- !
    IMPLICIT NONE

    REAL(8), INTENT(IN), DIMENSION(:) :: x   ! common x-vals
    REAL(8), INTENT(IN), DIMENSION(:) :: fexp
    INTEGER(4), INTENT(IN) :: max_exp
    REAL(8), INTENT(IN), DIMENSION(:) :: fthry
    REAL(8), INTENT(OUT) :: c
    REAL(8), INTENT(OUT) :: q2

    INTEGER(4) :: i
    REAL(8), DIMENSION(SIZE(x,1)) :: ax, ax1, ax2
    INTEGER(4), DIMENSION(SIZE(x,1)) :: mask 
   
    INTERFACE get_mask
       FUNCTION get_mask(xarr)
         REAL(8), INTENT(IN), DIMENSION(:) :: xarr 
         INTEGER(4), DIMENSION(SIZE(xarr,1)) :: get_mask
       END FUNCTION get_mask
    END INTERFACE

    mask=get_mask(x)
 
!   do i=1,max_exp
!      write(39,*) x(i),fthry(i)*mask(i)
!   enddo
!   c   = (SUM(fexp(1:max_exp)*fthry(1:max_exp)))/&
!   & (SUM(fthry(1:max_exp)*fthry(1:max_exp)))
!   ax(1:max_exp)  = (fexp(1:max_exp)-c*fthry(1:max_exp))
!   ax  = ax(1:max_exp)*ax(1:max_exp)
!   q2  = SUM(ax(1:max_exp))
!   write(1,*) 'max_exp=',max_exp
    ax  = 0.0d0    
    ax1 = 0.0d0
    ax2 = 0.0d0
    WHERE(mask.eq.1)
       ax1 = fexp*fthry
       ax2 = fthry*fthry
    ENDWHERE
       c   = SUM(ax1)/SUM(ax2)
    WHERE(mask.eq.1)
       ax  = fexp-c*fthry
       ax  = ax*ax
    ENDWHERE
       q2  = SUM(ax)
    END SUBROUTINE qsqrd

! ------------------------------------------------------------
     subroutine interpolate(np,max1,x1,y1,max2,x2,yn1)
! ------------------------------------------------------------
     implicit none
     integer, intent(in) :: np
     integer, intent(in) :: max1,max2
     real(8), dimension(np), intent(in) :: x1, y1
     real(8), dimension(np), intent(in) :: x2
     real(8), dimension(np), intent(out) :: yn1
     integer :: i
     integer:: ilower
     real(8) :: func
     external fbuff
!    yn1=0.0d0  Uncomment to crash code !!!!!!
     do i=1,max2
        call locate(x2(i),x1,max1,ilower,np)
        yn1(i)=func(y1,x1,x2(i),ilower)
     enddo
     end subroutine interpolate

! -------------------------------------------------------
      double precision function func(y,x,x0,ilower)
! -------------------------------------------------------
! Region: Profile
! Purpose: Linear Interpolation
!          Using bracket to interpolate at grid point x0
      implicit none
      integer, parameter :: np=5001
      real(8), dimension(np) :: x,y
      real(8) :: slope, x0
      integer :: iupper, ilower

      iupper=ilower+1
      slope=(y(iupper)-y(ilower))/(x(iupper)-x(ilower))
      func=y(iupper)+slope*(x0-x(iupper))
      end function func

