diff --git a/kap/make/makefile_base b/kap/make/makefile_base index 834a638f2..684458ca5 100644 --- a/kap/make/makefile_base +++ b/kap/make/makefile_base @@ -17,9 +17,9 @@ SRCS = \ kap_def.f90 \ op_def.f90 \ op_load.f \ - op_load_master.f \ + op_load_master.f90 \ op_common.f \ - op_ev.f \ + op_ev.f90 \ op_osc.f \ op_radacc.f \ op_eval.f \ @@ -100,7 +100,7 @@ else $(COMPILE_NO_CHECKS) $(FCfree) $< endif -op_load.o op_common.o op_ev.o op_osc.o op_radacc.o op_load_master.o : %.o : %.f +op_load.o op_common.o op_osc.o op_radacc.o : %.o : %.f ifneq ($(QUIET),) @echo COMPILE_LEGACY_NOCHECKS $< @$(COMPILE_LEGACY_NOCHECKS) $(FCfixed) $< diff --git a/kap/private/op_common.f b/kap/private/op_common.f index c181dd564..afd8d3ec4 100644 --- a/kap/private/op_common.f +++ b/kap/private/op_common.f @@ -2,8 +2,8 @@ module op_common use math_lib use op_def - - contains + + contains subroutine xindex(flt, ilab, xi, ih, i3, ierr) implicit none integer, intent(in) :: i3 @@ -30,26 +30,26 @@ subroutine xindex(flt, ilab, xi, ih, i3, ierr) do i = 1, 4 ih(i) = ih2 + i - 2 ilab(i) = i3*ih(i) - enddo + end do xi = 2.*(x-ih2) - 1 return end subroutine xindex c********************************************************************** - subroutine jrange(ih, jhmin, jhmax, i3) + subroutine jrange(ih, jhmin, jhmax, i3) implicit none integer, intent(in) :: ih(4), i3 - integer, intent(out) :: jhmin, jhmax + integer, intent(out) :: jhmin, jhmax integer :: i - + jhmin = 0 jhmax = 1000 do i = 1, 4 jhmin = max(jhmin, js(ih(i)*i3)/i3) jhmax = min(jhmax, je(ih(i)*i3)/i3) - enddo - + end do + return end subroutine jrange c********************************************************************** @@ -66,11 +66,11 @@ subroutine findne(ilab, fa, nel, nkz, jhmin, jhmax, ih, real, intent(inout) :: flrho c local variables integer :: i, j, n, jh, jm, itt, jne, jnn - real :: flrmin, flrmax, flr(4,4), uyi(4), efa(4, 7:118), + real :: flrmin, flrmax, flr(4,4), uyi(4), efa(4, 7:118), : flrh(4, 7:118), u(4), flnei(4), y, zeta, efa_temp -c declare variables in common block, by default: real (a-h, o-z), integer (i-n) +c declare variables in common block, by default: real (a-h, o-z), integer (i-n) ! integer :: ite1, ite2, ite3, jn1, jn2, jne3, ntot, nc, nf, int, -! : ne1, ne2, np, kp1, kp2, kp3, npp, mx, nx +! : ne1, ne2, np, kp1, kp2, kp3, npp, mx, nx ! real :: umin, umax, epatom, oplnck, fion, yy1, yy2, yx ! common /atomdata/ ite1,ite2,ite3,jn1(91),jn2(91),jne3,umin,umax,ntot, ! + nc,nf,int(17),epatom(17,91,25),oplnck(17,91,25),ne1(17,91,25), @@ -82,20 +82,20 @@ subroutine findne(ilab, fa, nel, nkz, jhmin, jhmax, ih, c efa(i,jh)=sum_n epa(i,jh,n)*fa(n) c flrh(i,jh)=log10(rho(i,jh)) c -c Get efa +c Get efa do i = 1, 4 itt = (ilab(i)-ite1)/2 + 1 do jne = jn1(itt), jn2(itt), i3 jnn = (jne-jn1(itt))/2 + 1 jh = jne/i3 efa_temp = 0. - do n = 1, nel + do n = 1, nel efa_temp = efa_temp + epatom(nkz(n), itt, jnn)*fa(n) - enddo !n - efa(i, jh) = efa_temp - enddo !jne - enddo !i -c + end do !n + efa(i, jh) = efa_temp + end do !jne + end do !i +c c Get range for efa.gt.0 do i = 1, 4 do jh = jhmin, jhmax @@ -103,18 +103,18 @@ subroutine findne(ilab, fa, nel, nkz, jhmin, jhmax, ih, jm = jh - 1 goto 3 endif - enddo + end do goto 4 3 jhmax = MIN(jhmax, jm) 4 continue - enddo -c + end do +c c Get flrh do jh=jhmin,jhmax do i=1,4 flrh(i, jh) = flmu + 0.25*i3*jh - log10(dble(efa(i, jh))) end do - end do + end do c c Find flrmin and flrmax flrmin = -1000 @@ -122,8 +122,8 @@ subroutine findne(ilab, fa, nel, nkz, jhmin, jhmax, ih, do i = 1, 4 flrmin = max(flrmin, flrh(i,jhmin)) flrmax = min(flrmax, flrh(i,jhmax)) - enddo -c + end do +c c Check range of flrho if(flrho .lt. flrmin .or. flrho .gt. flrmax)then ierr = 101 @@ -136,34 +136,34 @@ subroutine findne(ilab, fa, nel, nkz, jhmin, jhmax, ih, jm = jh - 1 goto 5 endif - enddo + end do print*,' Interpolations in j for flne' print*,' Not found, i=',i stop 5 jm=max(jm,jhmin+1) jm=min(jm,jhmax-2) -c +c do i = 1, 4 do j = 1, 4 u(j) = flrh(i, jm+j-2) flr(i,j) = flrh(i, jm+j-2) - enddo + end do call solve(u, flrho, zeta, uyi(i), ierr) if (ierr /= 0) return y = jm + 0.5*(zeta+1) flnei(i) = .25*i3*y - enddo + end do c c Interpolations in i flne = fint(flnei, xi) uy = fint(uyi, xi) -c Get epa +c Get epa epa = exp10(dble(flne + flmu - flrho)) -c +c return -c +c 601 format(' For flt=',1p,e11.3,', flrho=',e11.3,' is out of range'/ - + ' Allowed range for flrho is ',e11.3,' to ',e11.3) + + ' Allowed range for flrho is ',e11.3,' to ',e11.3) end subroutine findne c*********************************************************************** @@ -173,7 +173,7 @@ subroutine yindex(jhmin, jhmax, flne, jh, i3, eta) real, intent(in) :: flne integer, intent(out) :: jh(4) real, intent(out) :: eta -c local variables +c local variables integer :: j, k real :: y c @@ -183,9 +183,9 @@ subroutine yindex(jhmin, jhmax, flne, jh, i3, eta) j = min(j,jhmax-3) do k = 1, 4 jh(k) = j + k - 2 - enddo + end do eta = 2.*(y-j)-1 -c +c return end subroutine yindex @@ -194,20 +194,20 @@ subroutine findux(flr, xi, eta, ux) implicit none real, intent(in) :: flr(4, 4), xi, eta real, intent(out) :: ux -c local variables - integer :: i, j +c local variables + integer :: i, j real :: uxj(4), u(4) c do j = 1, 4 do i = 1, 4 u(i) = flr(i, j) - enddo + end do uxj(j) = fintp(u, xi) - enddo + end do ux = fint(uxj, eta) c return - end subroutine findux + end subroutine findux C************************************** function fint(u,r) @@ -216,7 +216,7 @@ function fint(u,r) c If P(R) = u(1) u(2) u(3) u(4) c for R = -3 -1 1 3 c then a cubic fit is: - P(R)=( + P(R)=( + 27*(u(3)+u(2))-3*(u(1)+u(4)) +R*( + 27*(u(3)-u(2))-(u(4)-u(1)) +R*( + -3*(u(2)+u(3))+3*(u(4)+u(1)) +R*( @@ -233,7 +233,7 @@ function fintp(u,r) c If P(R) = u(1) u(2) u(3) u(4) c for R = -3 -1 1 3 c then a cubic fit to the derivative is: - PP(R)=( + PP(R)=( + 27*(u(3)-u(2))-(u(4)-u(1)) +2.*R*( + -3*(u(2)+u(3))+3*(u(4)+u(1)) +3.*R*( + -3*(u(3)-u(2))+(u(4)-u(1)) )))/48. @@ -243,14 +243,14 @@ function fintp(u,r) return end function fintp C -c*********************************************************************** +c*********************************************************************** subroutine scatt(ih, jh, rion, uf, f, umesh, semesh, dscat, ntot, epa, ierr) use op_load, only: BRCKR integer, intent(inout) :: ierr real :: umesh(:), semesh(:) ! (nptot) real :: f(:,:,:) ! (nptot,4,4) dimension rion(28, 4, 4),uf(0:100), - + fscat(0:100),p(nptot),rr(28),ih(4),jh(4) + + fscat(0:100),p(nptot),rr(28),ih(4),jh(4) integer i,j,k,n c HH: always use meshtype q='m' ite3=2 @@ -264,16 +264,16 @@ subroutine scatt(ih, jh, rion, uf, f, umesh, semesh, dscat, ntot, epa, ierr) fne = real(exp10(ITE3*dble(jh(j))/4d0)) do k = 1, ntot p(k) = f(k, i, j) - enddo + end do do m = 1, 28 rr(m) = rion(m, i, j) - enddo + end do CALL BRCKR(FT,FNE,RR,28,UF,100,FSCAT,ierr) if (ierr /= 0) return do n = 0, 100 u = uf(n) fscat(n) = cscat*(fscat(n)-1) - enddo + end do do n = 2, ntot-1 u=umesh(n) se=semesh(n) @@ -281,27 +281,27 @@ subroutine scatt(ih, jh, rion, uf, f, umesh, semesh, dscat, ntot, epa, ierr) ua=umin+dscat*m ub=ua+dscat p(n)=p(n)+((ub-u)*fscat(m)+(u-ua)*fscat(m+1))/(dscat*se) - enddo + end do u=umesh(ntot) se=semesh(ntot) p(ntot)=p(ntot)+fscat(100)/se p(1)=p(1)+fscat(1)/(1.-.5*umin) do k=1,ntot f(k,i,j)=p(k) - enddo - enddo - enddo + end do + end do + end do C return end subroutine scatt -c*********************************************************************** +c*********************************************************************** subroutine screen1(ih,jh,rion,umesh,ntot,epa,f) use op_load, only: screen2 real :: umesh(:) ! (nptot) real, pointer :: f(:,:,:) ! (nptot,4,4) - dimension uf(0:100), - + fscat(0:100), ih(4), jh(4) - real, target :: rion(28, 4, 4) + dimension uf(0:100), + + fscat(0:100), ih(4), jh(4) + real, target :: rion(28, 4, 4) integer i, j, k, m real, pointer :: p(:), rr(:) c @@ -315,16 +315,16 @@ subroutine screen1(ih,jh,rion,umesh,ntot,epa,f) fne = real(exp10(ITE3*dble(jh(j))/4d0)) ! do k=1,ntot p => f(1:ntot,i,j) -! enddo +! end do ! do m=1,28 rr => rion(1:28,i,j) -! enddo +! end do call screen2(ft,fne,rr,epa,ntot,umin,umax,umesh,p) ! do k=1,ntot ! f(k,i,j)=p(k) -! enddo - enddo - enddo +! end do + end do + end do C return end subroutine screen1 diff --git a/kap/private/op_def.f90 b/kap/private/op_def.f90 index 9a869e080..516f30186 100644 --- a/kap/private/op_def.f90 +++ b/kap/private/op_def.f90 @@ -24,59 +24,61 @@ ! *********************************************************************** -MODULE op_def - IMPLICIT NONE +module op_def + implicit none integer, parameter :: nptot = 10000 integer, parameter :: nrad = 17 integer, parameter :: ipe = 17 - integer,dimension(140:320),parameter :: JS=(/14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 15, 18, 19, 22, 23,& - 26, 27, 30, 31, 34, 34, 34, 34, 36, 36, 36, 36, 36, 36, 38, 38,& - 38, 38, 38, 39, 40, 40, 40, 40, 40, 41, 42, 42, 42, 42, 42, 43,& - 44, 44, 44, 44, 44, 45, 46, 46, 46, 46, 46, 47, 48, 48, 48, 48,& - 48, 49, 50, 50, 50, 50, 52, 52, 52, 52, 52, 52, 54, 54, 54, 54,& - 54, 54, 56, 56, 56, 56, 56, 56, 58, 58, 58, 58, 58, 58, 60, 60,& - 60, 60, 60, 60, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 65,& - 66, 66, 66, 66, 66, 67, 68, 68, 68, 68, 68, 69, 70, 70, 70, 70,& - 70, 71, 72, 72, 72, 72, 72, 73, 74, 74, 74, 74, 76, 76, 76, 76,& - 76, 76, 78, 78, 78, 78, 78, 78, 80, 80, 80, 80, 80, 80, 82, 82,& - 82, 82, 82, 82, 84, 84, 84, 84, 84, 84, 86, 86, 86, 86, 86, 87,& - 88, 88, 88, 88, 88/) + integer, dimension(140:320), parameter :: JS=(/ & + 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 15, 18, 19, 22, 23,& + 26, 27, 30, 31, 34, 34, 34, 34, 36, 36, 36, 36, 36, 36, 38, 38,& + 38, 38, 38, 39, 40, 40, 40, 40, 40, 41, 42, 42, 42, 42, 42, 43,& + 44, 44, 44, 44, 44, 45, 46, 46, 46, 46, 46, 47, 48, 48, 48, 48,& + 48, 49, 50, 50, 50, 50, 52, 52, 52, 52, 52, 52, 54, 54, 54, 54,& + 54, 54, 56, 56, 56, 56, 56, 56, 58, 58, 58, 58, 58, 58, 60, 60,& + 60, 60, 60, 60, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 65,& + 66, 66, 66, 66, 66, 67, 68, 68, 68, 68, 68, 69, 70, 70, 70, 70,& + 70, 71, 72, 72, 72, 72, 72, 73, 74, 74, 74, 74, 76, 76, 76, 76,& + 76, 76, 78, 78, 78, 78, 78, 78, 80, 80, 80, 80, 80, 80, 82, 82,& + 82, 82, 82, 82, 84, 84, 84, 84, 84, 84, 86, 86, 86, 86, 86, 87,& + 88, 88, 88, 88, 88/) - integer,dimension(140:320),parameter :: JE=(/& - 52, 56, 56, 58, 58, 60, 60, 60, 60, 62, 62, 64, 64, 66, 66, 68,& - 68, 70, 70, 72, 72, 74, 74, 74, 74, 76, 76, 78, 78, 80, 80, 80,& - 80, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 84, 86, 86, 86,& - 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 90, 88, 88, 88, 88, 92,& - 92, 92, 92, 94, 94, 94, 94, 94, 94, 89, 90, 90, 90, 90, 90, 91,& - 92, 92, 92, 92, 92, 94, 94, 90, 90, 92, 92, 92, 92, 92, 92, 93,& - 94, 94, 94, 94, 94, 94, 94, 96, 96, 96, 96, 96, 96, 98, 98, 98,& - 98, 98, 98, 99,100,100,100,100,100,100,100,102,102,102,102,102,& - 102,103,104,104,104,104,104,104,104,106,106,106,106,106,106,107,& - 108,108,108,108,108,108,108,110,110,110,110,110,110,111,112,112,& - 112,112,112,112,112,114,114,114,114,114,114,115,116,116,116,116,& - 116,116,116,118,118/) + integer,dimension(140:320),parameter :: JE=(/ & + 52, 56, 56, 58, 58, 60, 60, 60, 60, 62, 62, 64, 64, 66, 66, 68,& + 68, 70, 70, 72, 72, 74, 74, 74, 74, 76, 76, 78, 78, 80, 80, 80,& + 80, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 84, 86, 86, 86,& + 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 90, 88, 88, 88, 88, 92,& + 92, 92, 92, 94, 94, 94, 94, 94, 94, 89, 90, 90, 90, 90, 90, 91,& + 92, 92, 92, 92, 92, 94, 94, 90, 90, 92, 92, 92, 92, 92, 92, 93,& + 94, 94, 94, 94, 94, 94, 94, 96, 96, 96, 96, 96, 96, 98, 98, 98,& + 98, 98, 98, 99,100,100,100,100,100,100,100,102,102,102,102,102,& + 102,103,104,104,104,104,104,104,104,106,106,106,106,106,106,107,& + 108,108,108,108,108,108,108,110,110,110,110,110,110,111,112,112,& + 112,112,112,112,112,114,114,114,114,114,114,115,116,116,116,116,& + 116,116,116,118,118/) -! - INTEGER,DIMENSION(17),parameter :: kz=(/1,2,6,7,8,10,11,12,13,14,16,18,20,24,25,26,28/) - character(len=2), dimension(17),parameter :: name=(/'H ','He','C ','N ','O ','Ne','Na',& - 'Mg','Al','Si','S ','Ar','Ca','Cr','Mn',& - 'Fe','Ni'/) - REAL,DIMENSION(17),PARAMETER :: AMASS=(/1.0080,4.0026,12.0111,14.0067,15.9994,20.179,& - 22.9898,24.305,26.9815,28.086,32.06,39.948, & - 40.08,51.996,54.9380,55.847,58.71/) + integer, dimension(17), parameter :: kz=(/1,2,6,7,8,10,11,12,13,14,16,18,20,24,25,26,28/) + character(len=2), dimension(17), parameter :: name=(/ & + 'H ','He','C ','N ','O ','Ne','Na',& + 'Mg','Al','Si','S ','Ar','Ca','Cr','Mn',& + 'Fe','Ni'/) + real, dimension(17), parameter :: AMASS=(/ & + 1.0080d0, 4.0026d0, 12.0111d0, 14.0067d0, 15.9994d0, 20.179d0, & + 22.9898d0, 24.305d0, 26.9815d0, 28.086d0, 32.06d0, 39.948d0, & + 40.08d0, 51.996d0, 54.9380d0, 55.847d0, 58.71d0/) - integer,save :: ite1, ite2, ite3, jne3 , ntotp, nc, nf - integer,dimension(91),save :: jn1, jn2 - integer,dimension(17),save :: int - real,save :: umin, umax - real,dimension(17,91,25),save :: epatom, oplnck - integer,dimension(17,91,25),save :: ne1p, ne2p,np,kp1,kp2,kp3,npp - real,dimension(-1:28,28,91,25),save :: fionp - real,allocatable,DIMENSION(:),save :: yy2,yx - INTEGER,allocatable,DIMENSION(:),save :: nx + integer, save :: ite1, ite2, ite3, jne3 , ntotp, nc, nf + integer, dimension(91), save :: jn1, jn2 + integer, dimension(17), save :: int + real, save :: umin, umax + real, dimension(17, 91, 25), save :: epatom, oplnck + integer, dimension(17, 91, 25), save :: ne1p, ne2p, np, kp1, kp2, kp3, npp + real, dimension(-1:28, 28, 91, 25), save :: fionp + real, allocatable, dimension(:), save :: yy2, yx + integer, allocatable, dimension(:), save :: nx integer, parameter :: op_cache_version = 1 -END module op_def +end module op_def diff --git a/kap/private/op_ev.f b/kap/private/op_ev.f90 similarity index 54% rename from kap/private/op_ev.f rename to kap/private/op_ev.f90 index 72ff12fb3..ccb106f37 100644 --- a/kap/private/op_ev.f +++ b/kap/private/op_ev.f90 @@ -23,132 +23,128 @@ ! ! *********************************************************************** - module op_ev - use math_lib - use op_def - use const_def - use kap_def, only: kap_test_partials, kap_test_partials_val, kap_test_partials_dval_dx - - contains -c*********************************************************************** - subroutine abund(nel, izz, fa, flmu, fmu1, nkz) - implicit none - integer, intent(in) :: nel, izz(ipe) - real, intent(in) :: fa(ipe) - real, intent(out) :: flmu, fmu1(nrad) - integer, intent(out) :: nkz(ipe) -c local variables +module op_ev + use math_lib + use op_def + use const_def + use kap_def, only: kap_test_partials, kap_test_partials_val, kap_test_partials_dval_dx + + implicit none + +contains + + subroutine abund(nel, izz, fa, flmu, fmu1, nkz) + integer, intent(in) :: nel, izz(ipe) + real, intent(in) :: fa(ipe) + real, intent(out) :: flmu, fmu1(nrad) + integer, intent(out) :: nkz(ipe) + integer :: k, k1, k2, m real :: amamu(ipe), fmu, a1, c1, fmu0 -c -c Get k1,get amamu(k) + + ! Get k1,get amamu(k) do k = 1, nel do m = 1, ipe if(izz(k).eq.kz(m))then amamu(k) = amass(m) nkz(k) = m - goto 1 - endif - enddo - print*,' k=',k,', izz(k)=',izz(k) - print*,' kz(m) not found' + goto 1 + end if + end do + write(*,*) ' k=',k,', izz(k)=',izz(k) + write(*,*) ' kz(m) not found' stop 1 continue - enddo -c -c Mean atomic weight = fmu - fmu = 0. + end do + + ! Mean atomic weight = fmu + fmu = 0d0 do k = 1, nel fmu = fmu + fa(k)*amamu(k) - enddo -c + end do + do k2 = 1, nel a1 = fa(k2) - c1 = 1./(1.-a1) + c1 = 1d0/(1d0-a1) fmu0 = c1*(fmu - fa(k2)*amamu(k2)) - fmu1(k2) = a1*(amamu(k2)-fmu0)*1.660531e-24 !dmu/dlog xi - enddo -c - fmu = fmu*1.660531e-24 ! Convert to cgs + fmu1(k2) = a1*(amamu(k2)-fmu0)*1.660531d-24 !dmu/dlog xi + end do + + fmu = fmu*1.660531d-24 ! Convert to cgs flmu = log10(dble(fmu)) -c + return - end subroutine abund -c********************************************************************** + end subroutine abund - subroutine rd(nel, nkz, izz, ilab, jh, n_tot, ff, rr, i3, umesh, fac) - implicit none - integer, intent(in) :: nel, nkz(ipe), izz(ipe), ilab(4), - > jh(4), n_tot, i3 + + subroutine rd(nel, nkz, izz, ilab, jh, n_tot, ff, rr, i3, umesh, fac) + integer, intent(in) :: nel, nkz(ipe), izz(ipe), ilab(4), jh(4), n_tot, i3 real, intent(in) :: umesh(:) ! (nptot) real(dp), intent(in) :: fac(nel) real, intent(out) :: ff(:,:,:,:) ! (nptot, ipe, 4, 4) real, intent(out) :: rr(28, ipe, 4, 4) -c local variables + ! local variables integer :: i, j, k, l, m, n, itt, jnn, izp, ne1, ne2, ne, ib, ia real :: fion(-1:28), yb, ya, d -c declare variables in common block (instead of by default: real (a-h, o-z), integer (i-n)) -! integer :: ite1, ite2, ite3, jn1, jn2, jne3, ntot, nc, nf, int, -! : ne1p, ne2p, np, kp1, kp2, kp3, npp, mx, nx -! real :: umin, umax, epatom, oplnck, fionp, yy1, yy2, yx -! common /atomdata/ ite1,ite2,ite3,jn1(91),jn2(91),jne3,umin,umax,ntot, -! + nc,nf,int(17),epatom(17,91,25),oplnck(17,91,25),ne1p(17,91,25), -! + ne2p(17,91,25),fionp(-1:28,28,91,25),np(17,91,25),kp1(17,91,25), -! + kp2(17,91,25),kp3(17,91,25),npp(17,91,25),mx(33417000), -! + yy1(33417000),yy2(120000000),nx(19305000),yx(19305000) -! save /atomdata/ -c -c i=temperature inex -c j=density index -c k=frequency index -c n=element index -c Get: -c mono opacity cross-section ff(k,n,i,j) -c modified cross-section for selected element, ta(k,i,j) -c -c Initialisations - rr=0. -c -c Start loop on i (temperature index) + ! declare variables in common block (instead of by default: real (a-h, o-z), integer (i-n)) + ! integer :: ite1, ite2, ite3, jn1, jn2, jne3, ntot, nc, nf, int, + ! : ne1p, ne2p, np, kp1, kp2, kp3, npp, mx, nx + ! real :: umin, umax, epatom, oplnck, fionp, yy1, yy2, yx + ! common /atomdata/ ite1,ite2,ite3,jn1(91),jn2(91),jne3,umin,umax,ntot, + ! + nc,nf,int(17),epatom(17,91,25),oplnck(17,91,25),ne1p(17,91,25), + ! + ne2p(17,91,25),fionp(-1:28,28,91,25),np(17,91,25),kp1(17,91,25), + ! + kp2(17,91,25),kp3(17,91,25),npp(17,91,25),mx(33417000), + ! + yy1(33417000),yy2(120000000),nx(19305000),yx(19305000) + ! save /atomdata/ + + ! i=temperature inex + ! j=density index + ! k=frequency index + ! n=element index + ! Get: + ! mono opacity cross-section ff(k,n,i,j) + ! modified cross-section for selected element, ta(k,i,j) + ! + ! Initialisations + rr = 0d0 + + ! Start loop on i (temperature index) do i = 1, 4 itt = (ilab(i) - ite1)/2 + 1 do j = 1, 4 jnn = (jh(j)*i3 - jn1(itt))/2 + 1 -c Read mono opacities + ! Read mono opacities do n = 1, nel izp = izz(n) ne1 = ne1p(nkz(n), itt, jnn) ne2 = ne2p(nkz(n), itt, jnn) do ne = ne1, ne2 fion(ne) = fionp(ne, nkz(n), itt, jnn) - enddo + end do do ne = ne1, min(ne2, izp-2) rr(izp-1-ne, n, i, j) = fion(ne) - enddo -! call zetbarp(izp, ne1, ne2, fion, zet, i3) -! zetal(n, i, j) = zet + end do + ! call zetbarp(izp, ne1, ne2, fion, zet, i3) + ! zetal(n, i, j) = zet do k = 1, n_tot ff(k, n, i, j) = yy2(k+kp2(nkz(n), itt, jnn)) - enddo + end do if (fac(n) /= 1d0) then do k = 1, size(ff,dim=1) ff(k,n,i,j) = fac(n)*ff(k,n,i,j) end do end if - enddo !n - enddo !j - enddo !i -c - + end do !n + end do !j + end do !i return -c - end subroutine rd -c*********************************************************************** - subroutine ross(flmu, fmu1, dv, ntot, rs, rossl) - implicit none + end subroutine rd + + + subroutine ross(flmu, fmu1, dv, ntot, rs, rossl) integer, intent(in) :: ntot real, intent(in) :: flmu, dv, fmu1(nrad) real, intent(in) :: rs(:,:,:) ! (nptot, 4, 4) @@ -156,89 +152,89 @@ subroutine ross(flmu, fmu1, dv, ntot, rs, rossl) integer :: i, j, n, k2 real(dp) :: drs, dd, oross real :: fmu, tt, ss, dd2, drsp(nrad), exp10_flmu - real :: log10_bohr_radius_sqr = -16.55280 -c -c oross=cross-section in a.u. -c rossl=log10(ROSS in cgs) + real :: log10_bohr_radius_sqr = -16.55280d0 + + ! oross=cross-section in a.u. + ! rossl=log10(ROSS in cgs) exp10_flmu = real(exp10(dble(flmu))) do i = 1, 4 do j = 1, 4 drs = 0.d0 do n = 1, ntot - dd = 1.d0/rs(n, i, j) + dd = 1d0/rs(n, i, j) dd2 = dd*dd drs = drs + dd - enddo - oross = 1.d0/(drs*dv) + end do + oross = 1d0/(drs*dv) rossl(i, j) = log10(dble(oross)) + log10_bohr_radius_sqr - flmu !log10(fmu) - enddo !j - enddo !i -c + end do !j + end do !i + return - end subroutine ross -c*********************************************************************** - subroutine mix(ntot, nel, fa, ff, rs, rr, rion) - implicit none + end subroutine ross + + + subroutine mix(ntot, nel, fa, ff, rs, rr, rion) integer, intent(in) :: ntot, nel real, intent(in) :: ff(:,:,:,:) ! (nptot, ipe, 4, 4) real, intent(in) :: fa(ipe), rr(28, 17, 4, 4) real, intent(out) :: rs(:,:,:) ! (nptot, 4, 4) real, intent(out) :: rion(28, 4, 4) -c local variables + + ! local variables integer :: i, j, k, n, m real :: rs_temp, rion_temp -c + do j = 1, 4 do i = 1, 4 do n = 1, ntot rs_temp = ff(n,1,i,j)*fa(1) do k = 2, nel rs_temp = rs_temp + ff(n,k,i,j)*fa(k) - enddo + end do rs(n, i, j) = rs_temp !rs(n, i, j) = dot_product(ff(n,1:nel,i,j),fa(1:nel)) - enddo + end do do m = 1, 28 rion_temp = rr(m, 1, i, j)*fa(1) do k = 2, nel rion_temp = rion_temp + rr(m,k,i,j)*fa(k) - enddo + end do rion(m,i,j) = rion_temp !rion(m,i,j) = dot_product(rr(m,1:nel,i,j),fa(1:nel)) - enddo - enddo - enddo -c + end do + end do + end do + return - end subroutine mix -C*********************************************************************** - subroutine interp(nel, rossl, xi, eta, g, i3, ux, uy, gx, gy) + end subroutine mix + + + subroutine interp(nel, rossl, xi, eta, g, i3, ux, uy, gx, gy) use op_common, only: fint, fintp - implicit none + integer, intent(in) :: nel, i3 real, intent(in) :: ux, uy, xi, eta, rossl(4, 4) real, intent(out) :: gx, gy, g -c local variables + ! local variables integer :: i, j, l, k2 real :: V(4), U(4), vyi(4) -C - DO I = 1, 4 - DO J = 1, 4 - U(J) = rossl(I,J) - ENDDO - V(I) = FINT(U,ETA) + + do i = 1, 4 + do j = 1, 4 + U(J) = rossl(i,j) + end do + V(I) = fint(U,eta) vyi(i) = fintp(u,eta) - ENDDO - g=FINT(V,XI) - gy=fint(vyi,xi) - gx=fintp(v,xi) -c + end do + g = fint(V, xi) + gy = fint(vyi, xi) + gx = fintp(v, xi) + gy=gy/uy - gx=(80./real(I3))*(gx-gy*ux) + gx=(80d0/real(I3))*(gx-gy*ux) - RETURN - end subroutine interp - + return + end subroutine interp -c*********************************************************************** - end module op_ev +end module op_ev diff --git a/kap/private/op_eval.f b/kap/private/op_eval.f index 9585525a5..8b94c174d 100644 --- a/kap/private/op_eval.f +++ b/kap/private/op_eval.f @@ -127,8 +127,8 @@ subroutine eval_op_radacc( ierr = 5 return endif - enddo - enddo + end do + end do c outer: do i = 1, nel inner: do n = 1, ipe @@ -142,11 +142,11 @@ subroutine eval_op_radacc( endif cycle outer endif - enddo inner + end do inner write(6,*)'OP - CHEM. ELEMENT CANNOT BE INCLUDED: Z = ', izzp(i) ierr = 8 return - enddo outer + end do outer c c Calculate mean atomic weight (flmu) and c array kzz indicating elements for which to calculate g_rad @@ -247,7 +247,7 @@ subroutine eval_op_radacc( const = 13.30295d0 + log10(flux) ! = -log10(c) - log10(amu) + log10(flux) do k2 = 1, kk grl1(k2) = const + flmu - log10(dble(am1(k2))) + f(k2) + g ! log g_rad - enddo !k2 + end do !k2 c g1 = g ! log kappa @@ -321,12 +321,12 @@ subroutine eval_op_ev( endif cycle outer endif - enddo inner + end do inner write(6,*)'OP - CHEM. ELEMENT CANNOT BE INCLUDED: Z = ', + izzp(i) ierr=8 return - enddo outer + end do outer c Calculate mean atomic weight (flmu) call abund(nel, izz, fa, flmu, fmu1, nkz) @@ -480,12 +480,12 @@ subroutine eval_alt_op( endif cycle outer endif - enddo inner + end do inner write(6,*)'OP - CHEM. ELEMENT CANNOT BE INCLUDED: Z = ', + izzp(i) ierr=8 return - enddo outer + end do outer c Calculate mean atomic weight (flmu) call abund(nel, izz, fa, flmu, nkz) diff --git a/kap/private/op_eval_mombarg.f90 b/kap/private/op_eval_mombarg.f90 index 1d74e04b3..3dce72eea 100644 --- a/kap/private/op_eval_mombarg.f90 +++ b/kap/private/op_eval_mombarg.f90 @@ -74,7 +74,7 @@ subroutine compute_grad(k, fk, logT_face, logRho_face,l, r,lkap_ross_cell, lgrad epa_mix_cell = 0d0 do i=imin,imax epa_mix_cell(i) = dot_product(fk,epatom(1:nel,i)) - enddo + end do amu_mix_cell = dot_product(fk,amamu) @@ -136,9 +136,9 @@ subroutine compute_grad(k, fk, logT_face, logRho_face,l, r,lkap_ross_cell, lgrad !',', logRho_grid(ii,jj) missing_point(ii,jj) = 0 endif - enddo - enddo - enddo + end do + end do + end do if (SUM(missing_point) > 0) then retry = .true. if (SUM(missing_point(2:4,1:3)) == 0) then @@ -179,7 +179,7 @@ subroutine compute_grad(k, fk, logT_face, logRho_face,l, r,lkap_ross_cell, lgrad return endif tries = tries + 1 - enddo + end do !write(*,*) 'Points for interpolation selected', k !!! Compute the monochromatic cross-section for the local mixture. @@ -190,21 +190,21 @@ subroutine compute_grad(k, fk, logT_face, logRho_face,l, r,lkap_ross_cell, lgrad ik = i_grid(ii,jj)!+ 1648*(ke-1) do m=1,nptot sig_mix_cell(ii,jj,m) = dot_product(fk,sig(:,ik,m)) - enddo - enddo - enddo + end do + end do + end do !!! Compute the Rossland mean cross-section by integrating over variable v (mesh equally spaced in v). sig_Ross = 0d0 do jj=1,4 do ii=1,4 - sig_int(1:nptot-1) = (1/sig_mix_cell(ii,jj,1:nptot-1) + 1/sig_mix_cell(ii,jj,2:nptot))/2. * dv + sig_int(1:nptot-1) = (1d0/sig_mix_cell(ii,jj,1:nptot-1) + 1d0/sig_mix_cell(ii,jj,2:nptot))*0.5d0 * dv do m=1, nptot-1 !sig_Ross(ii,jj) = sig_Ross(ii,jj) + (1/sig_mix_cell(ii,jj,m) + 1/sig_mix_cell(ii,jj,m+1))/2. * dv sig_Ross(ii,jj) = sig_Ross(ii,jj) + sig_int(m) - enddo - enddo - enddo + end do + end do + end do lkap_Ross = log10_bohr_radius_sqr - log_amu_mix_cell - log10(sig_Ross) @@ -212,8 +212,8 @@ subroutine compute_grad(k, fk, logT_face, logRho_face,l, r,lkap_ross_cell, lgrad do jj = 1, 4 do ii = 1, 4 call rossl_interpolator% add_point(logT_grid(ii,jj), logRho_grid(ii,jj), lkap_Ross(ii,jj)) - enddo - enddo + end do + end do lkap_ross_cell = rossl_interpolator% evaluate(logT_face,logRho_face) @@ -228,15 +228,15 @@ subroutine compute_grad(k, fk, logT_face, logRho_face,l, r,lkap_ross_cell, lgrad gam = 0 do m=1, nptot-1 gam = gam + ((eumesh(ke,ik,m)/ sig_mix_cell(ii,jj,m)) +(eumesh(ke,ik,m+1)/ sig_mix_cell(ii,jj,m+1))) - enddo - gamma_k(ii,jj,ke) = gam/2. * dv - enddo - enddo - enddo + end do + gamma_k(ii,jj,ke) = 0.5d0*gam*dv + end do + end do + end do deallocate(sig_mix_cell,sig_int) - where (gamma_k < 0.) + where (gamma_k < 0d0) gamma_k = 10d-30 endwhere lgamm = log10(gamma_k) @@ -246,21 +246,21 @@ subroutine compute_grad(k, fk, logT_face, logRho_face,l, r,lkap_ross_cell, lgrad do jj = 1, 4 do ii = 1, 4 call gaml_interpolators(ke)% add_point(logT_grid(ii,jj), logRho_grid(ii,jj), lgamm(ii,jj,ke)) - enddo - enddo + end do + end do lgamm_cell(ke) = gaml_interpolators(ke)% evaluate(logT_face, logRho_face) !cell - enddo - lgamm_cell(1) = -30. - lgamm_cell(2) = -30. + end do + lgamm_cell(1) = -30d0 + lgamm_cell(2) = -30d0 !write(*,*) 'lgamm_cell computed'!, ke, lgamm_cell(ke) ! Compute log g_rad. flux = l / (4d0*pi*r*r) sf = lkap_ross_cell - log_c + log10(flux) + log10(amu_mix_cell) lgrad_cell = sf + lgamm_cell - log10(amamu) - lgrad_cell(1) = -30. - lgrad_cell(2) = -30. + lgrad_cell(1) = -30d0 + lgrad_cell(2) = -30d0 !write(*,*) 'log kappa + log g_rad Fe cell', k, lkap_ross_cell, lgrad_cell(16) !write(*,*) 'compute_grad done', k @@ -318,7 +318,7 @@ subroutine compute_gamma_grid(ngp, fk_all,lgamm_pcg, lkap_face_pcg, logT_pcg, lo epa_mix_cell = 0d0 do i=imin,imax epa_mix_cell(i) = dot_product(fk,epatom(1:nel,i)) - enddo + end do amu_mix_cell = dot_product(fk,amamu) @@ -333,19 +333,19 @@ subroutine compute_gamma_grid(ngp, fk_all,lgamm_pcg, lkap_face_pcg, logT_pcg, lo do i=1,1648 do m=1,nptot sig_mix_cell(i,m) = dot_product(fk,sig(:,i,m)) - enddo - enddo + end do + end do !$OMP END PARALLEL DO !!! Compute the Rossland mean cross-section by integrating over variable v (mesh equally spaced in v). sig_Ross = 0d0 !$OMP PARALLEL DO PRIVATE(i,m,sig_int) SCHEDULE(guided) do i=1,1648 - sig_int(1:nptot-1) = (1/sig_mix_cell(i,1:nptot-1) + 1/sig_mix_cell(i,2:nptot))/2. * dv !inv_sig_mix_cell(ii,jj,1:nptot-1) + inv_sig_mix_cell(ii,jj,2:nptot) - do m=1, nptot-1 - sig_Ross(i) = sig_Ross(i) + sig_int(m) - enddo - enddo + sig_int(1:nptot-1) = (1d0/sig_mix_cell(i,1:nptot-1) + 1d0/sig_mix_cell(i,2:nptot))*0.5d0 * dv !inv_sig_mix_cell(ii,jj,1:nptot-1) + inv_sig_mix_cell(ii,jj,2:nptot) + do m=1, nptot-1 + sig_Ross(i) = sig_Ross(i) + sig_int(m) + end do + end do !$OMP END PARALLEL DO @@ -359,13 +359,13 @@ subroutine compute_gamma_grid(ngp, fk_all,lgamm_pcg, lkap_face_pcg, logT_pcg, lo gam = 0d0 do m=1, nptot-1 gam = gam + ((eumesh(ke,i,m)/ sig_mix_cell(i,m))+(eumesh(ke,i,m+1)/ sig_mix_cell(i,m+1))) - enddo - gamma_k(ke,i) = gam/2. * dv - enddo - enddo + end do + gamma_k(ke,i) = 0.5d0*gam * dv + end do + end do !$OMP END PARALLEL DO - where (gamma_k < 0.) + where (gamma_k < 0d0) gamma_k = 10d-30 endwhere @@ -480,9 +480,9 @@ subroutine compute_grad_fast(k,fk, logT_face, logRho_face, l, r,lgrad_cell, ierr missing_point(ii,jj) = 0 endif - enddo - enddo - enddo + end do + end do + end do if (SUM(missing_point) > 0) then retry = .true. @@ -536,9 +536,9 @@ subroutine compute_grad_fast(k,fk, logT_face, logRho_face, l, r,lgrad_cell, ierr i_grid(ii,jj) = i endif - enddo - enddo - enddo + end do + end do + end do retry = .false. else @@ -548,25 +548,25 @@ subroutine compute_grad_fast(k,fk, logT_face, logRho_face, l, r,lgrad_cell, ierr endif endif tries = tries + 1 - enddo + end do do jj=1,4 do ii=1,4 ik = i_grid(ii,jj) do ke=3,nel lgamm(ii,jj,ke) = lgamm_pcg(ke,ik) - enddo + end do lkap_Ross(ii,jj) = lkap_face_pcg(ik) - enddo - enddo + end do + end do call rossl_interpolator% initialize() do jj = 1, 4 do ii = 1, 4 call rossl_interpolator% add_point(logT_grid(ii,jj), logRho_grid(ii,jj), lkap_Ross(ii,jj)) - enddo - enddo + end do + end do lkap_ross_cell = rossl_interpolator% evaluate(logT_face,logRho_face) @@ -575,19 +575,19 @@ subroutine compute_grad_fast(k,fk, logT_face, logRho_face, l, r,lgrad_cell, ierr do jj = 1, 4 do ii = 1, 4 call gaml_interpolators(ke)% add_point(logT_grid(ii,jj), logRho_grid(ii,jj), lgamm(ii,jj,ke)) - enddo - enddo + end do + end do lgamm_cell(ke) = gaml_interpolators(ke)% evaluate(logT_face, logRho_face) !cell - enddo + end do lgamm_cell(1) = -30d0 lgamm_cell(2) = -30d0 flux = l / (4d0*pi*r*r) sf = lkap_ross_cell - log_c + log10(flux) + log10(amu_mix_cell) lgrad_cell = sf + lgamm_cell - log10(amamu) - lgrad_cell(1) = -30. - lgrad_cell(2) = -30. + lgrad_cell(1) = -30d0 + lgrad_cell(2) = -30d0 end subroutine compute_grad_fast @@ -648,7 +648,7 @@ subroutine compute_kappa(k,fk, logT_cntr, logRho_cntr, lkap_ross_cell, dlnkap_ra epa_mix_cell = 0d0 do i=imin,imax epa_mix_cell(i) = dot_product(fk,epatom(1:nel,i)) - enddo + end do amu_mix_cell = dot_product(fk,amamu) @@ -709,9 +709,9 @@ subroutine compute_kappa(k,fk, logT_cntr, logRho_cntr, lkap_ross_cell, dlnkap_ra missing_point(ii,jj) = 0 endif - enddo - enddo - enddo + end do + end do + end do if (SUM(missing_point) > 0) then retry = .true. @@ -752,7 +752,7 @@ subroutine compute_kappa(k,fk, logT_cntr, logRho_cntr, lkap_ross_cell, dlnkap_ra return endif tries = tries + 1 - enddo + end do !!! Compute the monochromatic cross-section for the local mixture. @@ -765,20 +765,20 @@ subroutine compute_kappa(k,fk, logT_cntr, logRho_cntr, lkap_ross_cell, dlnkap_ra ik = i_grid(ii,jj) do m=1,nptot sig_mix_cell(ii,jj,m) = dot_product(fk,sig(:,ik,m)) - enddo - enddo - enddo + end do + end do + end do !!! Compute the Rossland mean cross-section by integrating over variable v (mesh equally spaced in v). sig_Ross = 0d0 do jj=1,4 do ii=1,4 - sig_int = (1/sig_mix_cell(ii,jj,1:nptot-1) + 1/sig_mix_cell(ii,jj,2:nptot))/2. * dv !(1:nptot-1) inv_sig_mix_cell(ii,jj,1:nptot-1) + inv_sig_mix_cell(ii,jj,2:nptot) + sig_int = (1/sig_mix_cell(ii,jj,1:nptot-1) + 1/sig_mix_cell(ii,jj,2:nptot))*0.5d0 * dv !(1:nptot-1) inv_sig_mix_cell(ii,jj,1:nptot-1) + inv_sig_mix_cell(ii,jj,2:nptot) do m=1, nptot-1 sig_Ross(ii,jj) = sig_Ross(ii,jj) + sig_int(m) - enddo - enddo - enddo + end do + end do + end do lkap_Ross = log10_bohr_radius_sqr - log_amu_mix_cell - log10(sig_Ross) @@ -787,8 +787,8 @@ subroutine compute_kappa(k,fk, logT_cntr, logRho_cntr, lkap_ross_cell, dlnkap_ra do jj = 1, 4 do ii = 1, 4 call rossl_interpolator% add_point(logT_grid(ii,jj), logRho_grid(ii,jj), lkap_Ross(ii,jj)) - enddo - enddo + end do + end do lkap_ross_cell = rossl_interpolator% evaluate(logT_cntr,logRho_cntr) dlnkap_rad_dlnT = rossl_interpolator% evaluate_deriv(logT_cntr,logRho_cntr, .true., .false.) @@ -847,7 +847,7 @@ subroutine compute_kappa_grid(fk,lkap_ross_pcg, ierr,ite,jne,epatom,amamu,sig) epa_mix_cell = 0d0 do i=imin,imax epa_mix_cell(i) = dot_product(fk,epatom(1:nel,i)) - enddo + end do amu_mix_cell = dot_product(fk,amamu) @@ -862,19 +862,19 @@ subroutine compute_kappa_grid(fk,lkap_ross_pcg, ierr,ite,jne,epatom,amamu,sig) do i=1,1648 do m=1,nptot sig_mix_cell(i,m) = dot_product(fk,sig(:,i,m)) - enddo - enddo + end do + end do !$OMP END PARALLEL DO !!! Compute the Rossland mean cross-section by integrating over variable v (mesh equally spaced in v). sig_Ross = 0 !$OMP PARALLEL DO PRIVATE(i,m,sig_int) SCHEDULE(guided) do i=1,1648 - sig_int(1:nptot-1) = (1/sig_mix_cell(i,1:nptot-1) + 1/sig_mix_cell(i,2:nptot))/2. * dv !inv_sig_mix_cell(ii,jj,1:nptot-1) + inv_sig_mix_cell(ii,jj,2:nptot) + sig_int(1:nptot-1) = (1/sig_mix_cell(i,1:nptot-1) + 1/sig_mix_cell(i,2:nptot))*0.5d0 * dv !inv_sig_mix_cell(ii,jj,1:nptot-1) + inv_sig_mix_cell(ii,jj,2:nptot) do m=1, nptot-1 sig_Ross(i) = sig_Ross(i) + sig_int(m) - enddo - enddo + end do + end do !$OMP END PARALLEL DO deallocate(sig_mix_cell,sig_int) @@ -937,7 +937,7 @@ subroutine compute_kappa_fast(k,fk, logT_cntr, logRho_cntr,lkap_ross_cell, dlnka epa_mix_cell = 0d0 do i=imin,imax epa_mix_cell(i) = dot_product(fk,epatom(:,i)) - enddo + end do amu_mix_cell = dot_product(fk,amamu) @@ -999,9 +999,9 @@ subroutine compute_kappa_fast(k,fk, logT_cntr, logRho_cntr,lkap_ross_cell, dlnka missing_point(ii,jj) = 0 endif - enddo - enddo - enddo + end do + end do + end do if (SUM(missing_point) > 0) then retry = .true. @@ -1055,9 +1055,9 @@ subroutine compute_kappa_fast(k,fk, logT_cntr, logRho_cntr,lkap_ross_cell, dlnka logRho_grid(ii,jj) = logRho(i) i_grid(ii,jj) = i endif - enddo - enddo - enddo + end do + end do + end do retry = .false. else write(*,*) 'Cannot find points for interpolation compute_kappa_fast', k, ite_min, jne_min, logT_min, logRho_min,logT_cntr, logRho_cntr, missing_point, ii_min, jj_min, offset1, offset2,imin, imax @@ -1067,21 +1067,21 @@ subroutine compute_kappa_fast(k,fk, logT_cntr, logRho_cntr,lkap_ross_cell, dlnka endif endif tries = tries + 1 - enddo + end do do jj=1,4 do ii=1,4 ik = i_grid(ii,jj) lkap_Ross(ii,jj) = lkap_ross_pcg(ik) - enddo - enddo + end do + end do call rossl_interpolator% initialize() do jj = 1, 4 do ii = 1, 4 call rossl_interpolator% add_point(logT_grid(ii,jj), logRho_grid(ii,jj), lkap_Ross(ii,jj)) - enddo - enddo + end do + end do lkap_ross_cell = rossl_interpolator% evaluate(logT_cntr,logRho_cntr) dlnkap_rad_dlnT = rossl_interpolator% evaluate_deriv(logT_cntr, logRho_cntr, .true., .false.) @@ -1111,8 +1111,8 @@ subroutine interpolate_kappa(k,logT_cntr, logRho_cntr,lkap_ross_cell, dlnkap_rad do jj = 1, 4 do ii = 1, 4 call rossl_interpolator% add_point(logT_grid(ii,jj), logRho_grid(ii,jj), lkap_grid(ii,jj)) - enddo - enddo + end do + end do lkap_ross_cell = rossl_interpolator% evaluate(logT_cntr,logRho_cntr) dlnkap_rad_dlnT = rossl_interpolator% evaluate_deriv(logT_cntr,logRho_cntr, .true., .false.) diff --git a/kap/private/op_load.f b/kap/private/op_load.f index ea274dfa5..1fa43b110 100644 --- a/kap/private/op_load.f +++ b/kap/private/op_load.f @@ -22,7 +22,7 @@ ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ! ! *********************************************************************** - + c FORTRAN 90 module for calculation of radiative accelerations, c based on the Opacity Project (OP) code "OPserver". c See CHANGES_HU for changes made to the original code. @@ -33,7 +33,7 @@ module op_load use math_lib use op_def logical :: have_loaded_op = .false. - + contains C****************************************************************** @@ -41,10 +41,10 @@ subroutine op_dload(path, cache_filename, ierr) implicit none character (len=*), intent(in) :: path, cache_filename integer, intent(out) :: ierr - - - - + + + + integer,parameter :: ipz=28 real :: am,amm,delp,dpack integer :: ios,it,ite11,ite22,ite33,itt,itte1,itte2,itte3,izz,jne,ite @@ -56,7 +56,7 @@ subroutine op_dload(path, cache_filename, ierr) real :: dv,dv1 integer :: cache_version - + common /mesh/ ntotv,dv,dv1,umesh,semesh ! common /atomdata/ ! common/atomdata/ ite1,ite2,ite3,jn1(91),jn2(91),jne3,umin,umax,ntot, @@ -64,8 +64,8 @@ subroutine op_dload(path, cache_filename, ierr) ! + ne2(17,91,25),fion(-1:28,28,91,25),np(17,91,25),kp1(17,91,25), ! + kp2(17,91,25),kp3(17,91,25),npp(17,91,25),mx(33417000), ! + yy1(33417000),yy2(120000000),nx(19305000),yx(19305000) -! - +! + integer,dimension(ipe) :: ifl,iflp character num(0:9)*1,zlab(ipe)*3,tlab*6,zlabp(ipe)*3 DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ @@ -78,7 +78,7 @@ subroutine op_dload(path, cache_filename, ierr) real :: y_temp(nptot) integer :: nx_index, left_n, right_n, n_index real :: left_val, right_val, cross_section, slope - + if(allocated(yy2) .eqv. .false.) then ! yy2 actually needs 29,563 x 10,000 length ALLOCATE(yy2(30000*10000),nx(19305000),yx(19305000),stat=ierr) @@ -88,16 +88,16 @@ subroutine op_dload(path, cache_filename, ierr) yx=0.0 ! write(*,*) "ierr",ierr end if - - + + ierr=0 if (have_loaded_op) return - + !$omp critical (critial_do_op_dload) - + if (have_loaded_op) goto 1001 - + !path = '../OP4STARS_1.3' !call getenv("oppath", path) !if (len(trim(path)) == 0) then @@ -142,7 +142,7 @@ subroutine op_dload(path, cache_filename, ierr) zlab(n)='m'//num(kz(n)/10)//num(kz(n)-10*(kz(n)/10)) iflp(n)=70+n zlabp(n)='a'//num(kz(n)/10)//num(kz(n)-10*(kz(n)/10)) - enddo + end do write(*,*) 'loading OP mono data...' @@ -171,7 +171,7 @@ subroutine op_dload(path, cache_filename, ierr) NTOTP=NF IF(NTOTP.GT.nptot)then write(6,6002)ntotp,nptot - ierr=2 + ierr=2 goto 1001 endif INT(1)=1 @@ -183,7 +183,7 @@ subroutine op_dload(path, cache_filename, ierr) c ITE1=MAX(ITE1,ITTE1) c ITE2=MIN(ITE2,ITTE2) C -c READ MESH FILES +c READ MESH FILES OPEN(1,FILE=trim(path)//'/'//ZLAB(1)//'.mesh',status='old', + form='unformatted',iostat=ios) if (ios /= 0) then @@ -191,13 +191,13 @@ subroutine op_dload(path, cache_filename, ierr) ierr = -1 goto 1001 end if - READ(1)DV,NTOTV,(UMESH(N),N=1,NTOTV) + READ(1)DV,NTOTV,(UMESH(N),N=1,NTOTV) umin=umesh(1) umax=umesh(ntotv) DV1=DV - CLOSE(1) -C -C GET MESH FOR SCREEN + CLOSE(1) +C +C GET MESH FOR SCREEN CALL IMESH(UMESH,NTOTV) C C SUBSEQUENT FILES @@ -285,7 +285,7 @@ subroutine op_dload(path, cache_filename, ierr) goto 1001 end if endif - enddo + end do C READ HEADINGS NN=1 READ(IFL(1))IZZ,ITE,AM,UM,UX,NCCC,NFFF,DelP,JNE1,JNE2,JNE3 @@ -301,7 +301,7 @@ subroutine op_dload(path, cache_filename, ierr) endif JNE1=MAX(JNE1,JNE11) JNE2=MIN(JNE2,JNE22) - enddo + end do itt=(it-ite1)/2+1 jn1(itt)=jne1 jn2(itt)=jne2 @@ -328,14 +328,14 @@ subroutine op_dload(path, cache_filename, ierr) left_n = nx_temp(nx_index-1) right_n = nx_temp(nx_index) slope = (right_val - left_val)/float(right_n - left_n) - + do n_index = left_n, right_n cross_section = left_val + (n_index-left_n)*slope yy2(ncount2 + n_index) = cross_section - enddo + end do yy2(ncount2 + left_n) = left_val yy2(ncount2 + right_n) = right_val - enddo + end do kp2(n, itt, jnn) = ncount2 ncount2 = ncount2 + ntotp else @@ -351,9 +351,9 @@ subroutine op_dload(path, cache_filename, ierr) ncount3=ncount3+npp(n,itt,jnn) endif endif - enddo - enddo -c + end do + end do +c c write(6,610)it c write(6,*)'ncount1 = ',ncount1 c write(6,*)'ncount2 = ',ncount2 @@ -366,25 +366,25 @@ subroutine op_dload(path, cache_filename, ierr) close(iflp(n)) 150 CONTINUE c - enddo - + end do + write(*,*) 'done loading OP mono data' have_loaded_op = .true. - + !write(*,*)'ncount1 = ',ncount1 !write(6,*)'ncount2 = ',ncount2 !write(6,*)'ncount3 = ',ncount3 ios = 0 - open(1, file=trim(cache_filename), iostat=ios, + open(1, file=trim(cache_filename), iostat=ios, > action='write', form='unformatted') if (ios == 0) then write(*,*) 'write ' // trim(cache_filename) write(1) op_cache_version, ntotv,dv,dv1,umesh, > ite1,ite2,ite3,jn1,jn2,jne3,umin,umax,ntotp,nc,nf,int,epatom,oplnck, ne1p, - > ne2p,fionp,np,kp1,kp2,kp3,npp,yy2,nx,yx + > ne2p,fionp,np,kp1,kp2,kp3,npp,yy2,nx,yx close(1) end if - + 1001 continue C pre-calculate semesh @@ -428,20 +428,20 @@ end subroutine op_dload c*********************************************************************** SUBROUTINE IMESH(UMESH,NTOT) -C +C DIMENSION UMESH(nptot) COMMON/CIMESH/U(100),AA(nptot),BB(nptot),IN(nptot),ITOT,NN save /cimesh/ - + UMIN=UMESH(1) UMAX=UMESH(NTOT) -c +c II=100 A=(II*UMIN-UMAX)/REAL(II-1) B=(UMAX-UMIN)/REAL(II-1) DO I=1,II U(I)=A+B*I - ENDDO + ENDDO c ib=2 ub=u(ib) @@ -458,21 +458,21 @@ SUBROUTINE IMESH(UMESH,NTOT) nn=n-1 ibb=ib-1 goto 1 - endif - endif + endif + endif in(n)=ib aa(n)=(ub-umesh(n))/d bb(n)=(umesh(n)-ua)/d - enddo + end do c 1 ib=ibb do n=nn+1,ntot ib=ib+1 in(n)=ib u(ib)=umesh(n) - enddo + end do itot=ib -c +c return end SUBROUTINE IMESH @@ -484,8 +484,8 @@ subroutine msh(dv, ntot, umesh, semesh, uf, dscat) real, intent(out) :: umesh(:), semesh(:) ! (nptot) integer :: i, k, ntotv real :: dvp, dv1, umin, umax, umeshp(nptot), semeshp(nptot) - common /mesh/ ntotv, dvp, dv1, umeshp, semeshp - save /mesh/ + common /mesh/ ntotv, dvp, dv1, umeshp, semeshp + save /mesh/ c ntot = ntotv dv = dvp @@ -502,13 +502,13 @@ subroutine msh(dv, ntot, umesh, semesh, uf, dscat) dscat = (umax - umin)*0.01 do i = 0, 100 uf(i) = umin + i*dscat - enddo + end do c return c end subroutine msh - + subroutine solve(u,v,z,uz,ierr) integer, intent(inout) :: ierr dimension u(4) @@ -516,13 +516,13 @@ subroutine solve(u,v,z,uz,ierr) c If P(R) = u(1) u(2) u(3) u(4) c for R = -3 -1 1 3 c then a cubic fit is: - P(R)=( + P(R)=( + 27*(u(3)+u(2))-3*(u(1)+u(4)) +R*( + 27*(u(3)-u(2))-(u(4)-u(1)) +R*( + -3*(u(2)+u(3))+3*(u(4)+u(1)) +R*( + -3*(u(3)-u(2))+(u(4)-u(1)) ))))/48. c First derivative is: - PP(R)=( + PP(R)=( + 27*(u(3)-u(2))-(u(4)-u(1))+ 2*R*( + -3*(u(2)+u(3))+3*(u(4)+u(1)) +3*R*( + -3*(u(3)-u(2))+(u(4)-u(1)) )))/48. @@ -537,19 +537,19 @@ subroutine solve(u,v,z,uz,ierr) d=(v-p(z))/uz z=z+d if(abs(d).lt.1.e-4)return - enddo -c + end do +c ! print*,' Not converged after 10 iterations in SOLVE' ! print*,' v=',v ! DO N=1,4 ! PRINT*,' N, U(N)=',N,U(N) -! ENDDO +! ENDDO ierr = 10 return ! stop -c +c end subroutine solve -c*********************************************************************** +c*********************************************************************** SUBROUTINE BRCKR(T,FNE,RION,NION,U,NFREQ,SF, ierr) integer, intent(inout) :: ierr @@ -674,7 +674,7 @@ SUBROUTINE FDETA(X,ETA, ierr) C integer n,k c -! ierr = 0 +! ierr = 0 a=x*0.88622693 c IF(X.LT.1)THEN @@ -694,7 +694,7 @@ SUBROUTINE FDETA(X,ETA, ierr) 11 ETA=LOG(dble(U)) c ELSE - if(a.lt.2)then + if(a.lt.2)then E=LOG(dble(X)) else e=pow(1.5d0*a,2d0/3d0) @@ -760,20 +760,20 @@ subroutine screen2(ft,fne,rion,epa,ntot,umin,umax,umesh,p) data twopi/6.283185/ COMMON/CIMESH/U(100),AA(nptot),BB(nptot),IN(nptot),ITOT,NN save /cimesh/ -c +c rydt=ft/157894. aune=1.48185e-25*fne -c -c get alp2=1/(Debye)**2 +c +c get alp2=1/(Debye)**2 b=0 do i=1,ipz b=b+rion(i)*i**2 - enddo + end do alp2=(5.8804e-19)*fne*b/(epa*ft) if(alp2/ft.lt.5e-8)return !!!!!!!!!!! -c +c c=1.7337*aune/sqrt(rydt) -c +c do i=1,itot w=u(i)*rydt f(i)=0. @@ -790,22 +790,22 @@ subroutine screen2(ft,fne,rion,epa,ntot,umin,umax,umesh,p) q=(1./x2-1./x1+LOG(dble(x1/x2)))* + (fkp*(1.-exp(dble(-twopi*k/fkp))))/(fk*(1.-exp(dble(-twopi*k/fk)))) ff=ff+wt(j)*q - enddo + end do f(i)=f(i)+crz*ff 1 continue - enddo + end do c p(1)=f(1) do n=2,nn w=umesh(n)*rydt p(n)=p(n)+(aa(n)*f(in(n)-1)+bb(n)*f(in(n)))/(w*w*w) - enddo + end do do n=nn+1,ntot w=umesh(n)*rydt p(n)=p(n)+f(in(n))/(w*w*w) - enddo + end do c return - end subroutine screen2 + end subroutine screen2 end module op_load diff --git a/kap/private/op_load_master.f b/kap/private/op_load_master.f deleted file mode 100644 index c41405cb3..000000000 --- a/kap/private/op_load_master.f +++ /dev/null @@ -1,58 +0,0 @@ - module op_load_master - - use const_def, only: dp - - - implicit none - logical :: loaded_op_master = .false. - contains - - - - subroutine load_op_master(emesh_data_for_op_mono_path, iz,ite,jne,epatom,amamu,sig,eumesh,ierr) - - character (len=*), intent(in) :: emesh_data_for_op_mono_path - - integer, intent(inout) :: ierr - integer , pointer, intent(out) :: iz(:),ite(:),jne(:) - real(dp), pointer, intent(out) :: sig(:,:,:) - real(dp), pointer, intent(out):: epatom(:,:),amamu(:),eumesh(:,:,:) - - integer :: n, m, ke, ik - CHARACTER(LEN=72) :: FMT - integer :: nel, nptot, np - parameter(nel = 17, nptot = 10000, np=1648) !number of elements and number of u-mesh points. - real(dp), allocatable :: amamu_f(:,:) - integer, allocatable :: iz_f(:,:) - - - if (loaded_op_master) return - - allocate(iz_f(nel,np),iz(nel),ite(np),jne(np),stat=ierr) - allocate(sig(nel,np,nptot),stat=ierr) - allocate(epatom(nel,np),amamu_f(nel,np),amamu(nel),eumesh(nel,np,nptot),stat=ierr) - - FMT = '(i2,1x,i3,1x,i3,1x,F14.10,1x,F14.10,10000(1x,E12.6E3),10000(1x,E13.6E3))' - - write(*,*) 'Opening file...' - open(1,file = emesh_data_for_op_mono_path,form = 'formatted', action ='read') - write(*,*) 'Loading OP mono data...' - - do ke =1, nel - do n =1,np - read(1,FMT)iz_f(ke,n),ite(n),jne(n),epatom(ke,n),amamu_f(ke,n),(sig(ke,n,m), m=1,nptot),(eumesh(ke,n,m), m=1,nptot) - enddo - enddo - - close(1) - - do ke=1,nel - amamu(ke) = amamu_f(ke,1) - iz(ke) = iz_f(ke,1) - enddo - - write(*,*) 'OP mono data loaded.' - ierr = 0 - loaded_op_master = .true. - end subroutine load_op_master - end module op_load_master diff --git a/kap/private/op_load_master.f90 b/kap/private/op_load_master.f90 new file mode 100644 index 000000000..7cf061c8d --- /dev/null +++ b/kap/private/op_load_master.f90 @@ -0,0 +1,59 @@ +module op_load_master + + use const_def, only: dp + + implicit none + + logical :: loaded_op_master = .false. + + contains + + subroutine load_op_master(emesh_data_for_op_mono_path, iz,ite,jne,epatom,amamu,sig,eumesh,ierr) + + character (len=*), intent(in) :: emesh_data_for_op_mono_path + + integer, intent(inout) :: ierr + integer , pointer, intent(out) :: iz(:),ite(:),jne(:) + real(dp), pointer, intent(out) :: sig(:,:,:) + real(dp), pointer, intent(out) :: epatom(:,:),amamu(:),eumesh(:,:,:) + + integer :: n, m, ke, ik + CHARACTER(LEN=72) :: fmt + integer, parameter :: nel = 17 ! number of elements + integer, parameter :: nptot = 10000 ! number of u-mesh points + integer, parameter :: np = 1648 + real(dp), allocatable :: amamu_f(:,:) + integer, allocatable :: iz_f(:,:) + + if (loaded_op_master) return + + allocate(iz_f(nel,np),iz(nel),ite(np),jne(np),stat=ierr) + allocate(sig(nel,np,nptot),stat=ierr) + allocate(epatom(nel,np),amamu_f(nel,np),amamu(nel),eumesh(nel,np,nptot),stat=ierr) + + fmt = '(i2,1x,i3,1x,i3,1x,F14.10,1x,F14.10,10000(1x,E12.6E3),10000(1x,E13.6E3))' + + write(*,*) 'Opening file...' + open(1,file = emesh_data_for_op_mono_path,form = 'formatted', action ='read') + write(*,*) 'Loading OP mono data...' + + do ke =1, nel + do n =1,np + read(1,fmt)iz_f(ke,n),ite(n),jne(n),epatom(ke,n),amamu_f(ke,n),(sig(ke,n,m), m=1,nptot),(eumesh(ke,n,m), m=1,nptot) + end do + end do + + close(1) + + do ke=1,nel + amamu(ke) = amamu_f(ke,1) + iz(ke) = iz_f(ke,1) + end do + + write(*,*) 'OP mono data loaded.' + ierr = 0 + loaded_op_master = .true. + + end subroutine load_op_master + +end module op_load_master diff --git a/kap/private/op_osc.f b/kap/private/op_osc.f index 295f5325d..0c676fc0d 100644 --- a/kap/private/op_osc.f +++ b/kap/private/op_osc.f @@ -22,7 +22,7 @@ ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ! ! *********************************************************************** - + module op_osc use math_lib use op_def @@ -37,34 +37,34 @@ subroutine abund(nel, izz, fa, flmu, nkz) real, intent(in) :: fa(ipe) real, intent(out) :: flmu integer, intent(out) :: nkz(ipe) -c local variables +c local variables integer :: k, k1, k2, m real :: amamu(ipe), fmu c c Get k1,get amamu(k) - do k = 1, nel + do k = 1, nel do m = 1, ipe if(izz(k).eq.kz(m))then amamu(k) = amass(m) nkz(k) = m goto 1 endif - enddo + end do print*,' k=',k,', izz(k)=',izz(k) print*,' kz(m) not found' stop - 1 continue - enddo -c + 1 continue + end do +c c Mean atomic weight = fmu fmu = 0. do k = 1, nel fmu = fmu + fa(k)*amamu(k) - enddo + end do c fmu = fmu*1.660531e-24 ! Convert to cgs flmu = log10(dble(fmu)) -c +c return end subroutine abund c********************************************************************** @@ -86,7 +86,7 @@ subroutine xindex(flt, ilab, xi, ih, i3, ierr) ierr = 102 return endif -c +c x = 40.*flt/real(i3) ih2 = x ih2 = max(ih2, 140/i3+2) @@ -94,25 +94,25 @@ subroutine xindex(flt, ilab, xi, ih, i3, ierr) do i = 0, 5 ih(i) = ih2 + i - 2 ilab(i) = i3*ih(i) - enddo + end do xi = 2.*(x-ih2) - 1 c return end subroutine xindex c********************************************************************** - subroutine jrange(ih, jhmin, jhmax, i3) + subroutine jrange(ih, jhmin, jhmax, i3) implicit none integer, intent(in) :: ih(0:5), i3 - integer, intent(out) :: jhmin, jhmax + integer, intent(out) :: jhmin, jhmax integer :: i -c +c jhmin = 0 jhmax = 1000 do i = 0, 5 jhmin = max(jhmin, js(ih(i)*i3)/i3) jhmax = min(jhmax, je(ih(i)*i3)/i3) - enddo -c + end do +c return end subroutine jrange c********************************************************************** @@ -120,7 +120,7 @@ subroutine findne(ilab, fa, nel, nkz, jhmin, jhmax, ih, + flrho, flt, xi, flne, flmu, flr, epa, uy, i3, ierr) use op_load, only : solve implicit none - integer, intent(in) :: ilab(0:5), nel, nkz(ipe), jhmin, + integer, intent(in) :: ilab(0:5), nel, nkz(ipe), jhmin, > ih(0:5), i3 integer, intent(inout) :: jhmax integer, intent(out) :: ierr @@ -129,11 +129,11 @@ subroutine findne(ilab, fa, nel, nkz, jhmin, jhmax, ih, real, intent(inout) :: flrho c local variables integer :: i, j, n, jh, jm, itt, jne, jnn - real :: flrmin, flrmax, flr(4,4), uyi(4), efa(0:5, 7:118), + real :: flrmin, flrmax, flr(4,4), uyi(4), efa(0:5, 7:118), : flrh(0:5, 7:118), u(4), flnei(4), y, zeta, efa_temp -c declare variables in common block, by default: real (a-h, o-z), integer (i-n) +c declare variables in common block, by default: real (a-h, o-z), integer (i-n) ! integer :: ite1, ite2, ite3, jn1, jn2, jne3, ntot, nc, nf, int, -! : ne1, ne2, np, kp1, kp2, kp3, npp, mx, nx +! : ne1, ne2, np, kp1, kp2, kp3, npp, mx, nx ! real :: umin, umax, epatom, oplnck, fion, yy1, yy2, yx ! common /atomdata/ ite1,ite2,ite3,jn1(91),jn2(91),jne3,umin,umax,ntot, ! + nc,nf,int(17),epatom(17,91,25),oplnck(17,91,25),ne1(17,91,25), @@ -145,20 +145,20 @@ subroutine findne(ilab, fa, nel, nkz, jhmin, jhmax, ih, c efa(i,jh)=sum_n epa(i,jh,n)*fa(n) c flrh(i,jh)=log10(rho(i,jh)) c -c Get efa +c Get efa do i = 0, 5 itt = (ilab(i)-ite1)/2 + 1 do jne = jn1(itt), jn2(itt), i3 jnn = (jne-jn1(itt))/2 + 1 jh = jne/i3 efa_temp = 0. - do n = 1, nel + do n = 1, nel efa_temp = efa_temp + epatom(nkz(n), itt, jnn)*fa(n) - enddo !n - efa(i, jh) = efa_temp - enddo !jne - enddo !i -c + end do !n + efa(i, jh) = efa_temp + end do !jne + end do !i +c c Get range for efa.gt.0 do i = 0, 5 do jh = jhmin, jhmax @@ -166,18 +166,18 @@ subroutine findne(ilab, fa, nel, nkz, jhmin, jhmax, ih, jm = jh - 1 goto 3 endif - enddo + end do goto 4 3 jhmax = MIN(jhmax, jm) 4 continue - enddo -c + end do +c c Get flrh do jh = jhmin,jhmax do i = 0,5 flrh(i, jh) = flmu + 0.25*i3*jh - log10(dble(efa(i,jh))) - enddo - enddo + end do + end do c c Find flrmin and flrmax flrmin = -1000 @@ -185,8 +185,8 @@ subroutine findne(ilab, fa, nel, nkz, jhmin, jhmax, ih, do i = 0, 5 flrmin = max(flrmin, flrh(i,jhmin)) flrmax = min(flrmax, flrh(i,jhmax)) - enddo -c + end do +c c Check range of flrho if(flrho .lt. flrmin .or. flrho .gt. flrmax)then ierr = 101 @@ -199,34 +199,34 @@ subroutine findne(ilab, fa, nel, nkz, jhmin, jhmax, ih, jm = jh - 1 goto 5 endif - enddo + end do print*,' Interpolations in j for flne' print*,' Not found, i=',i stop 5 jm=max(jm,jhmin+1) jm=min(jm,jhmax-2) -c +c do i = 1, 4 do j = 1, 4 u(j) = flrh(i, jm+j-2) flr(i,j) = flrh(i, jm+j-2) - enddo + end do call solve(u, flrho, zeta, uyi(i), ierr) if (ierr /= 0) return y = jm + 0.5*(zeta+1) flnei(i) = .25*i3*y - enddo + end do c c Interpolations in i flne = fint(flnei, xi) uy = fint(uyi, xi) -c Get epa - epa = exp10(dble(flne + flmu - flrho)) -c +c Get epa + epa = exp10(dble(flne + flmu - flrho)) +c return -c +c 601 format(' For flt=',1p,e11.3,', flrho=',e11.3,' is out of range'/ - + ' Allowed range for flrho is ',e11.3,' to ',e11.3) + + ' Allowed range for flrho is ',e11.3,' to ',e11.3) end subroutine findne c*********************************************************************** subroutine yindex(jhmin, jhmax, flne, jh, i3, eta) @@ -235,7 +235,7 @@ subroutine yindex(jhmin, jhmax, flne, jh, i3, eta) real, intent(in) :: flne integer, intent(out) :: jh(0:5) real, intent(out) :: eta -c local variables +c local variables integer :: j, k real :: y c @@ -245,9 +245,9 @@ subroutine yindex(jhmin, jhmax, flne, jh, i3, eta) j = min(j,jhmax-3) do k = 0, 5 jh(k) = j + k - 2 - enddo + end do eta = 2.*(y-j)-1 -c +c return end subroutine yindex c*********************************************************************** @@ -255,41 +255,41 @@ subroutine findux(flr, xi, eta, ux) implicit none real, intent(in) :: flr(4, 4), xi, eta real, intent(out) :: ux -c local variables - integer :: i, j +c local variables + integer :: i, j real :: uxj(4), u(4) c do j = 1, 4 do i = 1, 4 u(i) = flr(i, j) - enddo + end do uxj(j) = fintp(u, xi) - enddo + end do ux = fint(uxj, eta) c return - end subroutine findux + end subroutine findux c********************************************************************** subroutine rd(nel, nkz, izz, ilab, jh, n_tot, ff, rr, i3, umesh, fac) implicit none - integer, intent(in) :: nel, nkz(ipe), izz(ipe), ilab(0:5), + integer, intent(in) :: nel, nkz(ipe), izz(ipe), ilab(0:5), > jh(0:5), n_tot, i3 real, intent(in) :: umesh(nptot) real(dp), intent(in) :: fac(nel) real, intent(out) :: ff(:,:,0:,0:) ! (nptot, ipe, 6, 6) real, intent(out) :: rr(28, ipe, 0:5, 0:5) -c local variables +c local variables integer :: i, j, k, l, m, n, itt, jnn, izp, ne1, ne2, ne, ib, ia real :: fion(-1:28), yb, ya, d -c declare variables in common block (instead of by default: real (a-h, o-z), integer (i-n)) +c declare variables in common block (instead of by default: real (a-h, o-z), integer (i-n)) ! integer :: ite1, ite2, ite3, jn1, jn2, jne3, ntot, nc, nf, int, -! : ne1p, ne2p, np, kp1, kp2, kp3, npp, mx, nx -! real :: umin, umax, epatom, oplnck, fionp, yy1, yy2, yx +! : ne1p, ne2p, np, kp1, kp2, kp3, npp, mx, nx +! real :: umin, umax, epatom, oplnck, fionp, yy1, yy2, yx ! common /atomdata/ ite1,ite2,ite3,jn1(91),jn2(91),jne3,umin,umax,ntot, ! + nc,nf,int(17),epatom(17,91,25),oplnck(17,91,25),ne1p(17,91,25), ! + ne2p(17,91,25),fionp(-1:28,28,91,25),np(17,91,25),kp1(17,91,25), ! + kp2(17,91,25),kp3(17,91,25),npp(17,91,25),mx(33417000), -! + yy1(33417000),yy2(120000000),nx(19305000),yx(19305000) +! + yy1(33417000),yy2(120000000),nx(19305000),yx(19305000) ! save /atomdata/ c c i=temperature inex @@ -300,45 +300,45 @@ subroutine rd(nel, nkz, izz, ilab, jh, n_tot, ff, rr, i3, umesh, fac) c mono opacity cross-section ff(k,n,i,j) c modified cross-section for selected element, ta(k,i,j) c -c Initialisations +c Initialisations rr=0. ff=0. -c +c c Start loop on i (temperature index) do i = 0, 5 - itt = (ilab(i) - ite1)/2 + 1 + itt = (ilab(i) - ite1)/2 + 1 do j = 0, 5 - jnn = (jh(j)*i3 - jn1(itt))/2 + 1 -c Read mono opacities + jnn = (jh(j)*i3 - jn1(itt))/2 + 1 +c Read mono opacities do n = 1, nel izp = izz(n) ne1 = ne1p(nkz(n), itt, jnn) ne2 = ne2p(nkz(n), itt, jnn) do ne = ne1, ne2 fion(ne) = fionp(ne, nkz(n), itt, jnn) - enddo + end do do ne = ne1, min(ne2, izp-2) rr(izp-1-ne, n, i, j) = fion(ne) - enddo - + end do + do k = 1, n_tot ff(k, n, i, j) = yy2(k+kp2(nkz(n), itt, jnn)) - enddo - + end do + if (fac(n) /= 1d0) then do k = 1, size(ff,dim=1) ff(k,n,i,j) = fac(n)*ff(k,n,i,j) end do end if - enddo !n - enddo !j - enddo !i + end do !n + end do !j + end do !i c return c end subroutine rd c*********************************************************************** - subroutine ross(flmu, dv, ntot,rs, rossl) + subroutine ross(flmu, dv, ntot,rs, rossl) implicit none integer, intent(in) :: ntot real, intent(in) :: flmu, dv, rs(nptot, 0:5, 0:5) @@ -355,11 +355,11 @@ subroutine ross(flmu, dv, ntot,rs, rossl) do n = 1, ntot dd = 1.d0/rs(n, i, j) drs = drs + dd - enddo + end do oross = 1.d0/(drs*dv) rossl(i, j) = log10(oross) - 16.55280d0 - flmu !log10(fmu) - enddo !j - enddo !i + end do !j + end do !i c return end subroutine ross @@ -369,7 +369,7 @@ subroutine mix(ntot, nel, fa, ff, rs, rr, rion) integer, intent(in) :: ntot, nel real, intent(in) :: ff(nptot, ipe, 0:5, 0:5), fa(ipe), rr(28, 17, 0:5, 0:5) real, intent(out) :: rs(nptot, 0:5, 0:5), rion(28, 0:5, 0:5) -c local variables +c local variables integer :: i, j, k, n, m real :: rs_temp, rion_temp c @@ -379,24 +379,24 @@ subroutine mix(ntot, nel, fa, ff, rs, rr, rion) !rs_temp = ff(n,1,i,j)*fa(1) !do k = 2, nel ! rs_temp = rs_temp + ff(n,k,i,j)*fa(k) - !enddo - !rs(n,i,j) = rs_temp + !end do + !rs(n,i,j) = rs_temp rs(n, i, j) = dot_product(ff(n,1:nel,i,j),fa(1:nel)) - enddo + end do do m = 1, 28 !rion_temp = rr(m, 1, i, j)*fa(1) !do k = 2, nel ! rion_temp = rion_temp + rr(m,k,i,j)*fa(k) - !enddo + !end do !rion(m,i,j) = rion_temp rion(m,i,j) = dot_product(rr(m,1:nel,i,j),fa(1:nel)) - enddo - enddo - enddo + end do + end do + end do c return end subroutine mix -C*********************************************************************** +C*********************************************************************** subroutine interp(nel, rossl, xi, eta, g, i3, ux, uy, gx, gy) implicit none integer, intent(in) :: nel, i3 @@ -406,11 +406,11 @@ subroutine interp(nel, rossl, xi, eta, g, i3, ux, uy, gx, gy) integer :: i, j, l real :: V(4), U(4), vyi(4) real :: x3(3), fx3!, fxxy(0:5, 0:5), fyxy(0:5, 0:5) -c pointers and targets +c pointers and targets real, target :: fx(0:5, 0:5), fy(0:5, 0:5), fxy(0:5, 0:5), fyx(0:5, 0:5), : fxx(0:5, 0:5), fyy(0:5, 0:5), rossl(0:5, 0:5) - real, pointer :: f3(:), fin(:, :), finx(:, :), finy(:, :) -c + real, pointer :: f3(:), fin(:, :), finx(:, :), finy(:, :) +c c interpolation of g (=rosseland mean opacity) c Use refined techniques of bi-cubic spline interpolation (Seaton 1993): ! call deriv(rossl, fx, fy, fxy) @@ -418,87 +418,87 @@ subroutine interp(nel, rossl, xi, eta, g, i3, ux, uy, gx, gy) ! gy = 0.5*gy/uy ! gx = (80./real(i3))*(0.5*gx-gy*ux) c Alternatively, use interpolation techniques by M.-A. Dupret to ensure smooothness required -c for pulsation studies: +c for pulsation studies: do i = 0, 3 x3(1) = i x3(2) = i+1 - x3(3) = i+2 + x3(3) = i+2 do j = 1, 4 f3 => rossl(i:i+2, j) call deriv3(f3, x3, fx3) fx(i+1, j) = fx3 - enddo - enddo + end do + end do do i = 0, 3 x3(1) = i x3(2) = i+1 - x3(3) = i+2 - do j = 1, 4 + x3(3) = i+2 + do j = 1, 4 f3 => rossl(j, i:i+2) call deriv3(f3, x3, fx3) fy(j, i+1) = fx3 - enddo - enddo + end do + end do do i = 1, 2 x3(1) = i x3(2) = i+1 - x3(3) = i+2 + x3(3) = i+2 do j = 2, 3 f3 => fx(i:i+2, j) call deriv3(f3, x3, fx3) fxx(i+1, j) = fx3 - enddo - enddo + end do + end do do i = 1, 2 x3(1) = i x3(2) = i+1 - x3(3) = i+2 + x3(3) = i+2 do j = 2, 3 f3 => fx(j, i:i+2) call deriv3(f3, x3, fx3) fxy(j, i+1) = fx3 - enddo - enddo + end do + end do do i = 1, 2 x3(1) = i x3(2) = i+1 - x3(3) = i+2 + x3(3) = i+2 do j = 2, 3 f3 => fy(i:i+2, j) call deriv3(f3, x3, fx3) fyx(i+1, j) = fx3 - enddo - enddo + end do + end do do i = 1, 2 x3(1) = i x3(2) = i+1 - x3(3) = i+2 + x3(3) = i+2 do j = 2, 3 f3 => fy(j, i:i+2) call deriv3(f3, x3, fx3) fyy(j, i+1) = fx3 - enddo - enddo -! call deriv(rossl, fx, fy, fxy) -! call deriv(fx, fxx, fxy, fxxy) -! call deriv(fy, fyx, fyy, fyxy) + end do + end do +! call deriv(rossl, fx, fy, fxy) +! call deriv(fx, fxx, fxy, fxxy) +! call deriv(fy, fyx, fyy, fyxy) fin => rossl(2:3, 2:3) finx => fx(2:3, 2:3) finy => fy(2:3, 2:3) - call interp3(fin, finx, finy, xi, eta, g) + call interp3(fin, finx, finy, xi, eta, g) fin => fx(2:3, 2:3) finx => fxx(2:3, 2:3) - finy => fxy(2:3, 2:3) - call interp3(fin, finx, finy, xi, eta, gx) + finy => fxy(2:3, 2:3) + call interp3(fin, finx, finy, xi, eta, gx) fin => fy(2:3, 2:3) finx => fyx(2:3, 2:3) - finy => fyy(2:3, 2:3) - call interp3(fin, finx, finy, xi, eta, gy) + finy => fyy(2:3, 2:3) + call interp3(fin, finx, finy, xi, eta, gy) gy = 0.5*gy/uy gx = (80./real(i3))*(0.5*gx-gy*ux) c RETURN - end subroutine interp + end subroutine interp C************************************** function fint(u,r) dimension u(4) @@ -506,7 +506,7 @@ function fint(u,r) c If P(R) = u(1) u(2) u(3) u(4) c for R = -3 -1 1 3 c then a cubic fit is: - P(R)=( + P(R)=( + 27*(u(3)+u(2))-3*(u(1)+u(4)) +R*( + 27*(u(3)-u(2))-(u(4)-u(1)) +R*( + -3*(u(2)+u(3))+3*(u(4)+u(1)) +R*( @@ -523,7 +523,7 @@ function fintp(u,r) c If P(R) = u(1) u(2) u(3) u(4) c for R = -3 -1 1 3 c then a cubic fit to the derivative is: - PP(R)=( + PP(R)=( + 27*(u(3)-u(2))-(u(4)-u(1)) +2.*R*( + -3*(u(2)+u(3))+3*(u(4)+u(1)) +3.*R*( + -3*(u(3)-u(2))+(u(4)-u(1)) )))/48. @@ -537,53 +537,53 @@ end function fintp SUBROUTINE DERIV(f, fx, fy, fxy) C real, intent(in) :: f(0:5, 0:5) - real, intent(out) :: fx(0:5, 0:5), fy(0:5, 0:5), fxy(0:5, 0:5) + real, intent(out) :: fx(0:5, 0:5), fy(0:5, 0:5), fxy(0:5, 0:5) real :: C(6) C C GET FX - DO 70 J = 0, 5 + DO J = 0, 5 L=0 - DO 50 I = 0, 5 + DO I = 0, 5 L=L+1 C(L)=F(I,J) - 50 CONTINUE + END DO CALL GET(C,L) L=0 - DO 60 I = 0, 5 + DO I = 0, 5 L = L + 1 FX(I, J) = C(L) - 60 CONTINUE - 70 CONTINUE + END DO + END DO C C GET FY - DO 100 I = 0, 5 + DO I = 0, 5 L=0 - DO 80 J = 0, 5 + DO J = 0, 5 L = L + 1 C(L) = F(I, J) - 80 CONTINUE + END DO CALL GET(C,L) L=0 - DO 90 J = 0, 5 + DO J = 0, 5 L = L + 1 FY(I,J) = C(L) - 90 CONTINUE - 100 CONTINUE + END DO + END DO C C GET FXY - DO 130 I = 0, 5 + DO I = 0, 5 L = 0 - DO 110 J = 0, 5 + DO J = 0, 5 L = L + 1 C(L) = FX(I, J) - 110 CONTINUE + END DO CALL GET(C,L) L=0 - DO 120 J = 0, 5 + DO J = 0, 5 L = L + 1 FXY(I,J) = C(L) - 120 CONTINUE - 130 CONTINUE + END DO + END DO c RETURN C @@ -622,20 +622,20 @@ SUBROUTINE GET(F,N) D(1)=-.5 T(1)=.5*(-F(1)+F(2)-FP1) C - DO 10 J=2,N-1 + DO J=2,N-1 D(J)=-1./(4.+D(J-1)) T(J)=-D(J)*(F(J-1)-2.*F(J)+F(J+1)-T(J-1)) - 10 CONTINUE + END DO C D(N)=(FPN+F(N-1)-F(N)-T(N-1))/(2.+D(N-1)) C - DO 20 J=N-1,1,-1 + DO J=N-1,1,-1 D(J)=D(J)*D(J+1)+T(J) - 20 CONTINUE + END DO C - DO 30 J=2,N-1 + DO J=2,N-1 F(J)=-F(J)+F(J+1)-2.*D(J)-D(J+1) - 30 CONTINUE + END DO F(1)=FP1 F(N)=FPN C @@ -666,13 +666,13 @@ subroutine INTERP2(f, fx, fy, fxy, xi, eta, g, gx, gy) C C Y = (eta + 5.)/2. - X = (xi + 5.)/2. + X = (xi + 5.)/2. ! i = floor(x) ! j = floor(y) I = X + 1.E-5 IF(ABS(X-I).LE.1.E-5) X = I J = Y + 1.E-5 - IF(ABS(Y-J).LE.1.E-5) Y = J + IF(ABS(Y-J).LE.1.E-5) Y = J C C INTERPOLATE C @@ -726,36 +726,36 @@ subroutine INTERP2(f, fx, fy, fxy, xi, eta, g, gx, gy) C END SUBROUTINE INTERP2 c************************************************************* -c This subroutine estimates the partial derivative of a function +c This subroutine estimates the partial derivative of a function c as described in the PhD thesis by M.-A. Dupret. subroutine deriv3(f, x, fx) real, intent(in) :: f(3), x(3) real, intent(out) :: fx real :: a, a1, a2, b, x1, x2 -c +c x1 = (f(3)- f(2))/(x(3)-x(2)) x2 = (f(2)-f(1))/(x(2)-x(1)) b = abs(x1/x2) a = (2.*x(2) - x(1) - x(3))/((x(2)-x(1))*(x(2)-x(3))) a1 = (x(2)-x(3))/((x(1)-x(3))*(x(1)-x(2))) a2 = (x(2)-x(1))/((x(3)-x(2))*(x(3)-x(1))) - if (b .ge. 0.2 .and. b .le. 5.) then + if (b .ge. 0.2 .and. b .le. 5.) then fx = a*f(2) + a1*f(1) + a2*f(3) else if (abs(x1) f(1:ntot,i,j) -! enddo +! end do ! do m=1,28 rr => rion(1:28,i,j) -! enddo +! end do call screen2(ft,fne,rr,epa,ntot,umin,umax,umesh,p) ! do k=1,ntot ! f(k,i,j)=p(k) -! enddo - enddo - enddo +! end do + end do + end do C return end subroutine screen1 diff --git a/kap/private/op_radacc.f b/kap/private/op_radacc.f index f102b2add..8d538b210 100644 --- a/kap/private/op_radacc.f +++ b/kap/private/op_radacc.f @@ -22,14 +22,14 @@ ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ! ! *********************************************************************** - + module op_radacc - + use math_lib use op_def use const_def ! use op_load, only: nptot, ipe - + logical, parameter :: dbg = .false. @@ -42,24 +42,24 @@ subroutine abund(nel, izz, kk, iz1, fa, kzz, flmu, am1, fmu1, nkz) real, intent(in) :: fa(ipe) integer, intent(out) :: kzz(nrad), nkz(ipe) real, intent(out) :: flmu, am1(nrad), fmu1(nrad) -c local variables +c local variables integer :: k, k1, k2, m real :: fmu, a1, c1, fmu0, amamu(ipe) c c Get k1,get amamu(k) - do k = 1, nel + do k = 1, nel do m = 1, ipe if(izz(k).eq.kz(m))then amamu(k) = amass(m) nkz(k) = m goto 1 endif - enddo + end do print*,' k=',k,', izz(k)=',izz(k) print*,' kz(m) not found' stop - 1 continue - enddo + 1 continue + end do c k1=-1 do k2 = 1, kk @@ -67,48 +67,48 @@ subroutine abund(nel, izz, kk, iz1, fa, kzz, flmu, am1, fmu1, nkz) if(izz(k).eq.iz1(k2))then k1 = k endif - enddo + end do kzz(k2) = k1 - am1(k2) = amamu(k1) - enddo -c + am1(k2) = amamu(k1) + end do +c c Mean atomic weight = fmu - fmu = 0. + fmu = 0d0 do k = 1, nel fmu = fmu + fa(k)*amamu(k) - enddo + end do c do k2 = 1, kk a1 = fa(kzz(k2)) - c1 = 1./(1.-a1) + c1 = 1d0/(1d0-a1) fmu0 = c1*(fmu - fa(kzz(k2))*amamu(kzz(k2))) fmu1(k2) = a1*(am1(k2)-fmu0)*1.660531d-24 !dmu/dlog xi - enddo -c + end do +c fmu = fmu*1.660531d-24 ! Convert to cgs flmu = log10(dble(fmu)) -c +c return end subroutine abund c********************************************************************** - subroutine rd(i3, kk, kzz, nel, nkz, izz, ilab, jh, ntot, umesh, + subroutine rd(i3, kk, kzz, nel, nkz, izz, ilab, jh, ntot, umesh, : semesh, ff, rr, ta, fac) implicit none - integer, intent(in) :: i3, kk, kzz(nrad), nel, nkz(ipe), + integer, intent(in) :: i3, kk, kzz(nrad), nel, nkz(ipe), : izz(ipe), ilab(4), jh(4), ntot real, intent(in) :: umesh(:), semesh(:) ! (nptot) real(dp), intent(in) :: fac(nel) real, intent(out) :: ff(:,:,:,:) ! (nptot, ipe, 4, 4) real, intent(out) :: ta(:,:,:,:) ! (nptot, nrad, 4, 4) real, intent(out) :: rr(28, ipe, 4, 4) -c local variables +c local variables integer :: i, j, k, k2, l, n, m, itt, jnn, izp, ne1, ne2, ne, ib, ia real :: ya, yb, d, se, u, fion(-1:28) !, ff_temp(nptot), ta_temp(nptot) -c declare variables in common block, by default: real (a-h, o-z), integer (i-n) +c declare variables in common block, by default: real (a-h, o-z), integer (i-n) ! integer :: ite1, ite2, ite3, jn1, jn2, jne3, ntotp, nc, nf, int, -! : ne1p, ne2p, np, kp1, kp2, kp3, npp, mx, nx +! : ne1p, ne2p, np, kp1, kp2, kp3, npp, mx, nx ! real :: umin, umax, epatom, oplnck, fionp, yy1, yy2, yx ! common /atomdata/ ite1,ite2,ite3,jn1(91),jn2(91),jne3,umin,umax,ntotp, ! + nc,nf,int(17),epatom(17,91,25),oplnck(17,91,25),ne1p(17,91,25), @@ -116,7 +116,7 @@ subroutine rd(i3, kk, kzz, nel, nkz, izz, ilab, jh, ntot, umesh, ! + kp2(17,91,25),kp3(17,91,25),npp(17,91,25),mx(33417000), ! + yy1(33417000),yy2(120000000),nx(19305000),yx(19305000) ! save /atomdata/ - + include 'formats' c @@ -129,7 +129,7 @@ subroutine rd(i3, kk, kzz, nel, nkz, izz, ilab, jh, ntot, umesh, c modified cross-section for selected element, ta(k,i,j) c zet(i,j) for diffusion coefficient c -c Initialisations +c Initialisations if (dbg) then write(*,2) 'size(ff,dim=1)', size(ff,dim=1) write(*,2) 'size(ff,dim=2)', size(ff,dim=2) @@ -147,21 +147,21 @@ subroutine rd(i3, kk, kzz, nel, nkz, izz, ilab, jh, ntot, umesh, write(*,*) 'rd: ff = 0' end if - ff = 0. + ff = 0d0 if (dbg) write(*,*) 'rd: rr = 0' - rr = 0. + rr = 0d0 if (dbg) write(*,*) 'rd: ta = 0' - ta = 0. -c + ta = 0d0 +c c Start loop on i (temperature index) do i = 1, 4 if (dbg) write(*,2) 'rd: i', i - itt = (ilab(i) - ite1)/2 + 1 -c Read mono opacities + itt = (ilab(i) - ite1)/2 + 1 +c Read mono opacities do j = 1, 4 if (dbg) write(*,2) 'rd: j', j - jnn = (jh(j)*i3 - jn1(itt))/2 + 1 - do n = 1, nel + jnn = (jh(j)*i3 - jn1(itt))/2 + 1 + do n = 1, nel if (dbg) write(*,2) 'rd: n', n izp = izz(n) ne1 = ne1p(nkz(n), itt, jnn) @@ -169,11 +169,11 @@ subroutine rd(i3, kk, kzz, nel, nkz, izz, ilab, jh, ntot, umesh, do ne = ne1, ne2 fion(ne) = fionp(ne, nkz(n), itt, jnn) if (ne .le. min(ne2, izp-2)) rr(izp-1-ne, n, i, j) = fion(ne) - enddo - + end do + do k = 1, ntot ff(k, n, i, j) = yy2(k+kp2(nkz(n), itt, jnn)) - enddo + end do if (fac(n) /= 1d0) then do k = 1, size(ff,dim=1) @@ -181,13 +181,13 @@ subroutine rd(i3, kk, kzz, nel, nkz, izz, ilab, jh, ntot, umesh, end do end if - do k2 = 1, kk + do k2 = 1, kk if (dbg) write(*,2) 'rd: k2', k2 if (kzz(k2)== n) then ib = 1 yb = yx(1+kp3(nkz(kzz(k2)), itt, jnn)) u = umesh(ib) - se = semesh(ib) + se = semesh(ib) ta(1, k2, i, j) = se*ff(1, n, i, j) - yb do m = 2, npp(nkz(kzz(k2)), itt, jnn) ia = ib @@ -199,24 +199,24 @@ subroutine rd(i3, kk, kzz, nel, nkz, izz, ilab, jh, ntot, umesh, u = umesh(l) se = semesh(l) ta(l, k2, i, j) = se*ff(l, n, i, j) -(ya + (l-ia)*d) - enddo + end do u = umesh(ib) se = semesh(ib) ta(ib, k2, i, j) = se*ff(ib, n, i, j) - yb - enddo - goto 101 ! get out of k2-loop - endif - enddo !k2 - 101 continue - enddo !n - enddo !j - enddo !i + end do + goto 101 ! get out of k2-loop + endif + end do !k2 + 101 continue + end do !n + end do !j + end do !i c return c end subroutine rd - + c*********************************************************************** subroutine mix(kk, kzz, ntot, nel, fa, ff, rr, rs, rion) implicit none @@ -225,7 +225,7 @@ subroutine mix(kk, kzz, ntot, nel, fa, ff, rr, rs, rion) real, intent(in) :: ff(:,:,:,:) ! (nptot, ipe, 4, 4) real, intent(out) :: rs(:,:,:) ! (nptot, 4, 4) real, intent(out) :: rion(28, 4, 4) -c local variables +c local variables integer :: i, j, k, n, m, k2 c @@ -233,18 +233,18 @@ subroutine mix(kk, kzz, ntot, nel, fa, ff, rr, rs, rion) do j = 1, 4 do n = 1, ntot rs(n, i, j) = dot_product(ff(n,1:nel,i,j),fa(1:nel)) - enddo - + end do + do m = 1, 28 rion(m,i,j) = dot_product(rr(m,1:nel,i,j),fa(1:nel)) - enddo - enddo - enddo + end do + end do + end do c return end subroutine mix c*********************************************************************** - subroutine ross(kk, flmu, fmu1, dv, ntot, rs, rossl, gaml, ta) + subroutine ross(kk, flmu, fmu1, dv, ntot, rs, rossl, gaml, ta) implicit none integer, intent(in) :: kk, ntot real, intent(in) :: rs(:,:,:) ! (nptot, 4, 4) @@ -252,7 +252,7 @@ subroutine ross(kk, flmu, fmu1, dv, ntot, rs, rossl, gaml, ta) real, intent(in) :: flmu, dv, fmu1(nrad) real, intent(out) :: rossl(4, 4), gaml(4, 4, nrad) -c local variables +c local variables integer :: k2, i, j, n real(dp) :: drs, dd, dgm(nrad) real :: fmu, oross, tt, exp10_flmu @@ -262,29 +262,29 @@ subroutine ross(kk, flmu, fmu1, dv, ntot, rs, rossl, gaml, ta) exp10_flmu = real(exp10(dble(flmu))) do i = 1, 4 do j = 1, 4 - drs = 0.d0 - dgm(:) = 0.d0 + drs = 0d0 + dgm(:) = 0d0 do n = 1, ntot !10000 - dd = 1.d0/rs(n, i, j) - drs = drs + dd + dd = 1d0/rs(n, i, j) + drs = drs + dd do k2 = 1, kk - tt = ta(n, k2, i, j) + tt = ta(n, k2, i, j) dgm(k2) = dgm(k2) + tt*dd - enddo - enddo - oross = 1./(drs*dv) - rossl(i, j) = log10(dble(oross)) - 16.55280 - flmu + end do + end do + oross = 1d0/(drs*dv) + rossl(i, j) = log10(dble(oross)) - 16.55280d0 - flmu do k2 = 1, kk if(dgm(k2).gt.0) then - dgm(k2) = dgm(k2)*dv + dgm(k2) = dgm(k2)*dv gaml(i, j, k2) = log10(dgm(k2)) else - gaml(i, j, k2) = -30. - endif - enddo - enddo !j - enddo !i -c + gaml(i, j, k2) = -30d0 + endif + end do + end do !j + end do !i +c return end subroutine ross c*********************************************************************** @@ -298,13 +298,13 @@ subroutine interp(nel, kk, rossl, gaml, xi, eta, g, i3, f) real :: v(4), u(4), vyi(4), xi c interpolation of g (=rosseland mean opacity) - DO i = 1, 4 - DO j = 1, 4 + do i = 1, 4 + do j = 1, 4 u(j) = rossl(i, j) - ENDDO + end do v(i) = fint(u, eta) vyi(i) = fintp(u, eta) - ENDDO + end do g = fint(v, xi) do k2 = 1, kk @@ -313,17 +313,17 @@ subroutine interp(nel, kk, rossl, gaml, xi, eta, g, i3, f) do i = 1, 4 do j = 1, 4 ! HH: This gives irregularities, perhaps it is preferable to assign nonnegative values of -! neighbouring interpolation points (in subroutine "ross") ? +! neighbouring interpolation points (in subroutine "ross") ? u(j) = gaml(i, j, k2) - enddo + end do v(i) = fint(u, eta) vyi(i) = fintp(u, eta) - enddo - + end do + f(k2) = fint(v, xi) - - 1 continue - enddo !k2 + + 1 continue + end do !k2 c return end subroutine interp diff --git a/star/test_suite/1.4M_ms_op_mono/inlist_1.4M_ms_op_mono b/star/test_suite/1.4M_ms_op_mono/inlist_1.4M_ms_op_mono index 2c59bbafe..341b0e3bb 100644 --- a/star/test_suite/1.4M_ms_op_mono/inlist_1.4M_ms_op_mono +++ b/star/test_suite/1.4M_ms_op_mono/inlist_1.4M_ms_op_mono @@ -44,7 +44,7 @@ use_gold2_tolerances = .true. ! limit max_model_number as part of test_suite - max_model_number = 10 + max_model_number = 30 ! we expect to be stopped by max_model_number before this xa_central_lower_limit_species(1) = 'h1' diff --git a/star/test_suite/do1_test_source b/star/test_suite/do1_test_source index d0b2ef9dd..79d1c6450 100755 --- a/star/test_suite/do1_test_source +++ b/star/test_suite/do1_test_source @@ -17,7 +17,7 @@ do_one pisn "termination code: Star is unbound" "final.mod" skip do_one 15M_dynamo "all values are within tolerances" "final.mod" auto do_one 1.3M_ms_high_Z "stop because log_surface_luminosity >= log_L_upper_limit" "final.mod" auto -do_one 1.4M_ms_op_mono "successfully used OP_mono opacities" skip skip +do_one 1.4M_ms_op_mono "successfully used OP_mono opacities" "final.mod" auto do_one 1.5M_with_diffusion "stop because have dropped below central lower limit for h1" "final.mod" auto do_one 5M_cepheid_blue_loop "crossed blue edge to start 3rd crossing" "final.mod" auto do_one 7M_prems_to_AGB "stop because log_surface_luminosity >= log_L_upper_limit" "final.mod" auto @@ -69,7 +69,7 @@ do_one ns_h "stop because time >= max_age_in_seconds" "final.mod" auto do_one ns_he "stop because power_nuc_burn >= power_nuc_burn_upper_limit" "final.mod" auto do_one other_physics_hooks "stop because have dropped below central lower limit for h1" "final.mod" auto do_one R_CrB_star "stop because Teff <= Teff_lower_limit" "final.mod" auto -do_one radiative_levitation "found expected effects of radiative levitation" skip skip +do_one radiative_levitation "found expected effects of radiative levitation" "final.mod" auto do_one relax_composition_j_entropy "termination code: xa_central_lower_limit" "final.mod" auto do_one rsp_BEP "good match for period" "final.mod" skip