diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py index 52174cc8..afdb2aa6 100644 --- a/python/ctem/ctem/driver.py +++ b/python/ctem/ctem/driver.py @@ -59,6 +59,11 @@ def __init__(self, param_file="ctem.in", isnew=True): 'ejc': 'surface_ejc.dat', 'pos': 'surface_pos.dat', 'time': 'surface_time.dat', + 'stack': 'surface_stacknum.dat', + 'rego': 'surface_rego.dat', + 'melt': 'surface_melt.dat', + 'comp': 'surface_comp.dat', + 'age' : 'surface_age.dat', 'ocum': 'ocumulative.dat', 'odist': 'odistribution.dat', 'pdist': 'pdistribution.dat', @@ -309,6 +314,8 @@ def process_output(self): if (self.user['savetruelist'].upper() == 'T'): self.redirect_outputs(['tcum'], 'dist') self.redirect_outputs(['impmass'], 'misc') + if (self.user['saverego'].upper() == 'T') : + self.redirect_outputs(['stack','rego','age','melt','comp'], 'rego') diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index c3c779c3..50ac5760 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -89,7 +89,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt ! doregotrack & age simulation test real(DP) :: melt, clock, age, thick real(SP),dimension(user%gridsize, user%gridsize) :: agetop - real(SP),dimension(60) :: agetot + real(SP),dimension(MAXAGEBINS) :: agetot type(regolisttype),pointer :: current => null() real(DP) :: age_resolution diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index f29b7600..6da3fb25 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -63,6 +63,7 @@ module module_globals real(DP),parameter :: PF = 5.0e9 ! The shock pressure exceeding the hungoit elastic limit of common geological materials in Pa. (5 GPa) real(DP),parameter :: RAD_GP = 20.0_DP ! The maximum radial position of producing impact glass spherules within a transient crater (unit of crater radii, crater%rad) + type regodatatype real(SP),dimension(MAXAGEBINS) :: age real(DP) :: thickness @@ -275,7 +276,7 @@ module module_globals character(*),parameter :: TIMEFILE = 'surface_time.dat' character(*),parameter :: EJCOVFILE = 'surface_ejc.dat' character(*),parameter :: DEMFILE = 'surface_dem.dat' -character(*),parameter :: REGOFILE = 'surface_regotop.dat' +character(*),parameter :: REGOFILE = 'surface_rego.dat' character(*),parameter :: MELTFILE = 'surface_melt.dat' character(*),parameter :: COMPFILE = 'surface_comp.dat' character(*),parameter :: STACKNUMFILE = 'surface_stacknum.dat' diff --git a/src/io/io_read_regotrack.f90 b/src/io/io_read_regotrack.f90 index ed93300c..3acdbb14 100644 --- a/src/io/io_read_regotrack.f90 +++ b/src/io/io_read_regotrack.f90 @@ -28,14 +28,14 @@ subroutine io_read_regotrack(user,surf) ! Internals integer(I4B), parameter :: LUN=7 - integer(I4B), parameter :: LUM=8 - integer(I4B), parameter :: LUP=9 - integer(I4B), parameter :: LUC=10 - integer(I4B), parameter :: LUA=11 + integer(I4B), parameter :: FMELT = 10 + integer(I4B), parameter :: FREGO = 11 + integer(I4B), parameter :: FCOMP = 12 + integer(I4B), parameter :: FAGE = 13 real(DP),dimension(user%gridsize,user%gridsize) :: regotop,melt,comp - real(SP),dimension(user%gridsize,user%gridsize,60) :: age + real(SP),dimension(user%gridsize,user%gridsize,MAXAGEBINS) :: age integer(I4B),dimension(user%gridsize,user%gridsize) :: stacks_num - real(DP), dimension(:),allocatable :: regotopi,melti,compi,agei + real(DP), dimension(:), allocatable :: regotopi,melti,compi,agei type(regodatatype) :: newsurfi integer(I4B) :: ioerr,i,j,k,q,itmp integer(kind=8) :: recsize @@ -46,34 +46,34 @@ subroutine io_read_regotrack(user,surf) ! Open a file for obtaining the number of stacks that is stored in each linked list ioerr = 0 recsize = sizeof(itmp) * user%gridsize * user%gridsize - open(LUP,file=STACKNUMFILE,status='old',form='unformatted',recl=recsize,access='direct',iostat=ioerr) + open(LUN,file=STACKNUMFILE,status='old',form='unformatted',recl=recsize,access='direct',iostat=ioerr) if (ioerr/=0) then write(*,*) 'Error! Cannot read file ',trim(adjustl(STACKNUMFILE)) stop end if - read(LUP,rec=1) stacks_num - close(LUP) + read(LUN,rec=1) stacks_num + close(LUN) ! Open files for regolith thickness/melt fraction - open(LUN,file=REGOFILE,status='old',form='unformatted',iostat=ioerr) + open(FREGO,file=REGOFILE,status='old',form='unformatted',iostat=ioerr) if (ioerr/=0) then write(*,*) 'Error! Cannot read file ',trim(adjustl(REGOFILE)) stop end if - open(LUC,file=COMPFILE,status='old',form='unformatted',iostat=ioerr) + open(FCOMP,file=COMPFILE,status='old',form='unformatted',iostat=ioerr) if (ioerr/=0) then write(*,*) 'Error! Cannot read file ',trim(adjustl(COMPFILE)) stop end if - open(LUM,file=MELTFILE,status='old',form='unformatted',iostat=ioerr) + open(FMELT,file=MELTFILE,status='old',form='unformatted',iostat=ioerr) if (ioerr/=0) then write(*,*) 'Error! Cannot read file ',trim(adjustl(MELTFILE)) stop end if - open(LUA,file=AGEFILE,status='old',form='unformatted',iostat=ioerr) + open(FAGE,file=AGEFILE,status='old',form='unformatted',iostat=ioerr) if (ioerr/=0) then write(*,*) 'Error! Cannot read file ',trim(adjustl(MELTFILE)) stop @@ -89,18 +89,18 @@ subroutine io_read_regotrack(user,surf) allocate(regotopi(stacks_num(i,j))) allocate(compi(stacks_num(i,j))) allocate(melti(stacks_num(i,j))) - allocate(agei(60 * stacks_num(i,j))) + allocate(agei(MAXAGEBINS * stacks_num(i,j))) do k=1,stacks_num(i,j) - read(LUN) regotop(i,j) + read(FREGO) regotop(i,j) regotopi(k) = regotop(i,j) - read(LUC) comp(i,j) + read(FCOMP) comp(i,j) compi(k) = comp(i,j) - read(LUM) melt(i,j) + read(FMELT) melt(i,j) melti(k) = melt(i,j) - read(LUA) age(i,j,:) - do q=1,60 - agei(60*k - (60-q)) = age(i,j,q) + read(FAGE) age(i,j,:) + do q=1,MAXAGEBINS + agei(MAXAGEBINS*k - (MAXAGEBINS-q)) = age(i,j,q) end do end do @@ -108,8 +108,8 @@ subroutine io_read_regotrack(user,surf) newsurfi%thickness = regotopi(k) newsurfi%comp = compi(k) newsurfi%meltfrac = melti(k) - do q=1,60 - newsurfi%age(q) = agei(60*k-(60-q)) + do q=1,MAXAGEBINS + newsurfi%age(q) = agei(MAXAGEBINS*k-(MAXAGEBINS-q)) end do call util_push(surf(i,j)%regolayer,newsurfi) end do @@ -118,8 +118,9 @@ subroutine io_read_regotrack(user,surf) end do end do - close(LUN) - close(LUC) - close(LUM) + close(FMELT) + close(FREGO) + close(FCOMP) + close(FAGE) return end subroutine io_read_regotrack diff --git a/src/io/io_write_regotrack.f90 b/src/io/io_write_regotrack.f90 index a6f3c2e1..5f522321 100644 --- a/src/io/io_write_regotrack.f90 +++ b/src/io/io_write_regotrack.f90 @@ -27,14 +27,12 @@ subroutine io_write_regotrack(user,surf) ! Regotrack Internals integer(I4B) :: i,j,k - integer(I4B), parameter :: LUN=7 - integer(I4B), parameter :: LUM=8 - integer(I4B), parameter :: LUC=9 - integer(I4B), parameter :: LUA=10 + integer(I4B), parameter :: LUN = 7 + integer(I4B), parameter :: FMELT = 10 + integer(I4B), parameter :: FREGO = 11 + integer(I4B), parameter :: FCOMP = 12 + integer(I4B), parameter :: FAGE = 13 type(regolisttype),pointer :: current => null() - real(DP),dimension(user%gridsize,user%gridsize) :: regotop,comp,melt - real(SP),dimension(user%gridsize,user%gridsize) :: agetop - real(SP),dimension(user%gridsize,user%gridsize,60) :: age integer(I4B),dimension(user%gridsize,user%gridsize) :: stacks_num integer(kind=8) :: recsize real(DP) :: dtmp @@ -42,104 +40,35 @@ subroutine io_write_regotrack(user,surf) integer(I4B) :: itmp real(DP),dimension(user%gridsize,user%gridsize) :: comptop, rego real(DP),dimension(:),allocatable :: marehisto - real(DP) :: mare, z - ! Mixing - !real(DP),parameter :: zmix = 0.0_DP - !real(DP) :: z, zmare - - ! Output multiple "comphisto" files - character(len=255) :: fname - character(len=255), parameter :: clockfile = 'tic-toc.dat' - integer(I4B) :: tictoc - logical :: exist - ! Executable code - ! Output mulitple "comphisto" files - inquire(file=clockfile, exist=exist) - if (exist) then - open(LUN,file=clockfile,status='old') - read(LUN,*) tictoc - else - write(*,*) clockfile,' is missing!' - end if - tictoc = tictoc + 1 - close(LUN) - open(LUN,file=clockfile,status='replace') - write(LUN,*) tictoc - close(LUN) - - write(fname,'(a,i4.4)') 'surface_melt',tictoc - open(LUN,file=fname,status='replace',form='unformatted') - write(fname,'(a,i4.4)') 'surface_rego',tictoc - open(LUM,file=fname,status='replace',form='unformatted') - write(fname,'(a,i4.4)') 'surface_comp',tictoc - open(LUC,file=fname,status='replace',form='unformatted') - write(fname,'(a,i4.4)') 'surface_age',tictoc - open(LUA,file=fname,status='replace',form='unformatted') + open(FMELT,file=MELTFILE,status='replace',form='unformatted') + open(FREGO,file=REGOFILE,status='replace',form='unformatted') + open(FCOMP,file=COMPFILE,status='replace',form='unformatted') + open(FAGE,file=AGEFILE,status='replace',form='unformatted') + stacks_num(:,:) = 0 do j=1,user%gridsize do i=1,user%gridsize - stacks_num(i,j) = 0 current => surf(i,j)%regolayer - comptop(i,j) = current%regodata%comp - rego(i,j) = current%regodata%thickness - agetop(i,j) = current%regodata%age(1) do - if (.not. associated(current)) exit - stacks_num(i,j) = stacks_num(i,j) + 1 - regotop(i,j) = current%regodata%thickness - comp(i,j) = current%regodata%comp - melt(i,j) = current%regodata%meltfrac - age(i,j,:)= current%regodata%age(:) - write(LUM) regotop(i,j) - write(LUC) comp(i,j) - write(LUN) melt(i,j) - write(LUA) age(i,j,:) - current => current%next + if (.not. associated(current)) exit ! We've reached the bottom of the linked list + stacks_num(i,j) = stacks_num(i,j) + 1 + write(FMELT) current%regodata%meltfrac + write(FREGO) current%regodata%comp + write(FCOMP) current%regodata%meltfrac + write(FAGE) current%regodata%age(:) + current => current%next end do end do end do - close(LUN) - close(LUM) - close(LUC) - close(LUA) - - recsize = sizeof(dtmp) * user%gridsize * user%gridsize - open(LUN,file='agetop.dat',status='replace',form='unformatted',recl=recsize,access='direct') - write(LUN,rec=1) agetop - close(LUN) - - recsize = sizeof(dtmp) * user%gridsize * user%gridsize - open(LUN,file='comptop.dat',status='replace',form='unformatted',recl=recsize,access='direct') - write(LUN,rec=1) comptop - close(LUN) - - allocate(marehisto(user%gridsize)) - !open(LUN,file='comphisto',status='replace') ! Output comphisto one time - write(fname,'(a,i4.4)') 'comphisto',tictoc - open(LUN,file=fname,status='replace') - - do i=1,user%gridsize - marehisto(i) = 0.0_DP - mare = 0.0_DP - do j=1,user%gridsize - mare = mare + comptop(i,j) - end do - marehisto(i) = mare/real(user%gridsize) - write(LUN,*) real(i-user%gridsize/2)*user%pix/1000.0,marehisto(i)*100.0 - end do - close(LUN) - deallocate(marehisto) - - recsize = sizeof(dtmp) * user%gridsize * user%gridsize - open(LUN,file='regotop.dat',status='replace',form='unformatted',recl=recsize,access='direct') - write(LUN,rec=1) rego - close(LUN) + close(FMELT) + close(FREGO) + close(FCOMP) + close(FAGE) recsize = sizeof(itmp) * user%gridsize * user%gridsize - write(fname,'(a,i4.4)') 'surface_stacknum',tictoc - open(LUN,file=fname,status='replace',form='unformatted',recl=recsize,access='direct') + open(LUN,file=STACKNUMFILE,status='replace',form='unformatted',recl=recsize,access='direct') write(LUN,rec=1) stacks_num close(LUN)