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