Cleaned up file io for rego tracking. File names now conform to the s…

…tandard CTEM style, removed redundant set of top files, and added rego files (age, melt, comp, and stack) to the redirected save files when doregotrack is turned on.

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')
 
 
 

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
 

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'

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,27 +89,27 @@ 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
 
          do k=max(stacks_num(i,j)-1,1),1,-1
             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

src/io/io_write_regotrack.f90

@@ -27,119 +27,48 @@ 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
    real(SP) :: stmp
    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)