c Copyright (C) Stichting Deltares, 2005-2024. c c This file is part of iMOD. c c This program is free software: you can redistribute it and/or modify c it under the terms of the GNU General Public License as published by c the Free Software Foundation, either version 3 of the License, or c (at your option) any later version. c c This program is distributed in the hope that it will be useful, c but WITHOUT ANY WARRANTY; without even the implied warranty of c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the c GNU General Public License for more details. c c You should have received a copy of the GNU General Public License c along with this program. If not, see . c c Contact: imod.support@deltares.nl c Stichting Deltares c P.O. Box 177 c 2600 MH Delft, The Netherlands. c c state save for Modflow2005 c c Method: c during forward run c ------------------ c 1: when initialising packages c - subscribe all input and output files for each package c 2: before reading stress/timestep information c - save the file pointers of all input files c - save the file pointers of all output files c - save the timestep information kper,kstp,delt c - save HNEW, IBOUND, VBVL to file c 3: when reading stress information c - save if new information is read or not for every input file c This has to be implemented in every RP-routine of every package c restoring a timestep c -------------------- c - find the last timestep at or before the the timestep to rollback to c where data was available c - set the file pointers of the input files to these corresponding records c - call the stress data read routines c - restore the timestep information kper,kstp,delt c - restore HNEW, IBOUND, VBVL from file c - set the file pointers of the input files to the wanted timestep c - set the file pointers of the output files to the wanted timestep module m_sts2 ! description: ! State Save module implicit none ! parameters to define filedat() array with ! REMARK: the order of the numbering of these parameters is important, the value itself isn't. ! In routine sts_setdata() the maximum value of the already defined ! value and the new value will be used. So, the highest value of these parameters ! has always precedence over the lower value. integer, parameter :: p_dataunk=0 ! unknown whether data is avilable for this timestep integer, parameter :: p_datano=1 ! NO data avilable for this timestep integer, parameter :: p_data=2 ! data avilable for this timestep ! parameter to set restore2 or not integer, parameter :: p_res2_yes=1 ! restore for unit has to be done in sts_restore2 too integer, parameter :: p_res2_no =0 ! restore for unit has NOT to be done in sts_restore2 logical, save :: usests=.false. ! .true.: sts-package is active integer, save :: nluns=0 ! number of units defined integer, save :: nts =0 ! number of timesteps stored double precision, save:: currenttime ! current time value logical, save :: currentTimeStepCalculated=.false. integer, save :: lastts=-1 ! last saved timestep logical, save :: restorets ! restore timestep integer, save :: stslun=0 ! unit number to save state data to type sts_tsinfo ! integer :: kper ! not used ! integer :: kstp ! not used ! real :: delt ! not used double precision :: timevalue ! start time of current time step (MJD) end type type(sts_tsinfo),dimension(:),pointer,save :: tsinfo type(sts_tsinfo),dimension(:),pointer,save :: tsinfo_dum type sts_luninfo ! information of unit integer :: lun ! unit number for which next information holds integer :: res2 ! Value is set in sts_restore1 routine ! 0: don't restore this unit in routine sts_restore2 ! 1: restore in restore2 integer(kind=8), dimension(:), pointer :: filepsb=>null() ! file pointer for corresponding lun at ! the beginning of this timestep ! -1: not set yet ! >=0: set integer(kind=8), dimension(:), pointer :: filepse=>null() ! file pointer for corresponding lun at ! the end of this timestep ! -1: not set yet ! >=0: set integer, dimension(:), pointer :: filedat=>null() ! Type of data filepointer points to ! p_dataunk: unknown whether data is available ! p_datano : data of former timestep reused ! p_data : new data available at this timestep ! This value is set by sts2data or sts2nodata end type type(sts_luninfo),dimension(:),pointer,save :: luninfo type(sts_luninfo),dimension(:),pointer,save :: luninfo_dum end module m_sts2 c ****************************************************************************** subroutine sts2init(lsts,lun) c ****************************************************************************** c Init State Save use m_sts2 implicit none logical, intent(in) :: lsts ! .true. switch sts on integer, intent(in) :: lun ! unit number to store state save integer :: ios,n character(len=52) :: fname integer :: osd_open2,cfn_getlun c init if (.not.lsts) then usests=.false.; write(*,*) '>>> State Save will not be used. <<<' return endif usests=.true.; write(*,*) '>>> State Save in effect. <<<' c allocate ! init timestep currenttime=0.0d0 restorets=.false. stslun =lun tsinfo=>null(); tsinfo_dum=>null() luninfo=>null(); luninfo_dum=>null() ! allocate timestep info pointer call ssts2increaseStructures(1,1) ! call ssts2increaseStructures(2000,200) ! unit number if (stslun.le.0)stslun=cfn_getlun(10,99) ! open state save file fname='StateSaveBin.dat' n=0; call pks7mpifname(fname,n) ios=osd_open2(stslun,0,fname,'unformatted') !'StateSaveBin.dat','unformatted') call sts2subscribe(stslun) end subroutine sts2init c ****************************************************************************** subroutine sts2start(timevalue) c ****************************************************************************** use m_sts2 implicit none double precision, intent(in) :: timevalue ! start date/time of current time step c check or sts is in effect if (.not. usests) return c init timestep currenttime=timevalue end subroutine sts2start c ****************************************************************************** subroutine ssts2nextts(timevalue) c ****************************************************************************** c Init next timestep c If timestep> currentts then continue run c timestep<=currentts data of former timestep needs to be restored c datainformation of units will be cleared use m_sts2 implicit none double precision, intent(in) :: timevalue ! start date/time of current time step integer ilun,its integer timestep,currentts integer(kind=8), pointer :: filepsb(:),filepse(:) integer, pointer :: filedat(:) c check or sts is in effect if (.not. usests) return c find idices for time values call ssts2timeposition(timevalue ,timestep, .false.) call ssts2timeposition(currenttime,currentts, .false.) c check c timestep has to be older then or equal to currentts c and c timestep has to be older then or equal to last saved timestep if (timestep.le.0 .or. currentts.le.0) then restorets=.false. else if (timestep.lt.currentts) then restorets=.true. else if (timestep.eq.currentts.and. 1 currentTimeStepCalculated) then restorets=.true. else restorets=.false. endif end if write(*,*)' sts2nextts: ',timevalue,currenttime,timestep, 1 currentts,restorets if (restorets) then ! restore timestep is wanted ! reset unit information from timestep+2...currentts do ilun=1,nluns filedat=>luninfo(ilun)%filedat filepsb=>luninfo(ilun)%filepsb filepse=>luninfo(ilun)%filepse do its=timestep+1,currentts filedat(its)=p_dataunk filepsb(its)=-1 filepse(its)=-1 enddo enddo endif write(*,*) ' ***** State Save nextts: ',restorets,currentts, 1 timestep c set currentts to the new value currentTimeStepCalculated=.false. end subroutine ssts2nextts c ****************************************************************************** subroutine sts2subscribe(lun) c ****************************************************************************** c description: c subscribe a logical unit number to the State Save module c This routine can be called by everyone who wants to subscribe a unit number use m_sts2 implicit none integer, intent(in) :: lun ! logical unit number to subscribe integer :: i,ilun character(len=256) :: fname c check or sts is in effect if (.not. usests) return c check or lun is already defined call ssts2getilun(lun,ilun) c add when not already defined if (ilun.le.0) then nluns=nluns+1 ! increase number of units call ssts2increaseStructures(0,nluns) ! assign luninfo(nluns)%lun=lun ! write(*,*) '>>> State Save Subscribe <<< ' inquire(unit=lun,name=fname) write(*,'(1x,a,i5,a)') 'State Save Subscribe unit=',lun, 1'; file='//trim(fname(index(fname,'\',.true.)+1:)) endif end subroutine sts2subscribe c ****************************************************************************** subroutine ssts2save(timevalue,phase,inoc) c ****************************************************************************** c save state for timestep c - save file-pointers c Check array sizes and re-allocate when necessary c Only save when no restore has to be done use m_sts2 implicit none double precision, intent(in) :: timevalue ! start date/time of current time step integer , intent(in) :: phase ! 1: at begin of timestep loop ! 2: at end of timestep loop integer,intent(in) :: inoc integer ilun integer ts integer(kind=8), pointer :: filepsb(:),filepse(:),filepos(:) integer, pointer :: filedat(:) character(len=256) :: fname c check or sts is in effect if (.not. usests) return c find idices for time values call ssts2timeposition(timevalue,ts,.false.) c if last saved timestep equal to ts, not necessary to save if (ts.eq.lastts .and. phase.eq.1) return call ssts2timeposition(timevalue,ts,.true.) c check array sizes call ssts2increaseStructures(ts,0) c save file pointers if (.not.restorets) then write(*,*) ' ***** State Save save: ',ts do ilun=1,nluns ! use pointer filedat => luninfo(ilun)%filedat filepsb => luninfo(ilun)%filepsb filepse => luninfo(ilun)%filepse ! set position depending on phase if (phase.eq.1) then filepos=>luninfo(ilun)%filepsb else filepos=>luninfo(ilun)%filepse endif ! save file position call osd_ftell(luninfo(ilun)%lun,filepos(ts)) !## do not increase byte-pos with 2 for oc-package if(luninfo(ilun)%lun.eq.inoc)then filepos(ts)=filepos(ts)-2 endif ! set filedat to unknown !filedat(ts)=p_dataunk inquire(unit=luninfo(ilun)%lun,name=fname) write(*,'(i5,i16,1x,a)') luninfo(ilun)%lun,filepos(ts), 1trim(fname(index(fname,'\',.true.)+1:)) ! write(*,*) ' lun: ',luninfo(ilun)%lun,filepos(ts) enddo ! save last save timestep lastts=ts ! save time value tsinfo(ts)%timevalue=timevalue ! set new number of saved timesteps nts=max(nts,ts) endif end subroutine ssts2save c ****************************************************************************** subroutine ssts2restore1(timevalue) c ****************************************************************************** c restore timestep c set file-pointers to the first timestep at or before ts where data was available use m_sts2 implicit none double precision, intent(in) :: timevalue ! start date/time of current time step integer ilun,its integer ts integer(kind=8), pointer :: filepsb(:),filepse(:) integer, pointer :: filedat(:) character(len=256) :: fname c check or sts is in effect if (.not. usests) return c find idices for time values call ssts2timeposition(timevalue,ts,.true.) c for each lun if (restorets) then write(*,*) ' ***** State Save Restore 1: ',ts do ilun=1,nluns ! use pointer filedat => luninfo(ilun)%filedat filepsb => luninfo(ilun)%filepsb filepse => luninfo(ilun)%filepse ! find position where data was found its=ts do while (its.gt.1 .and. filedat(its).ne.p_data) its=its-1 enddo ! check or pinfo%filedat(its).eq.p_data if (filedat(its).ne.p_data) then ! see or timestep ts had some data to read if (filepsb(ts).ne.filepse(ts)) then its=ts else ! try again, look for a timestep where the fileposition changes its=ts do while (its.gt.1 .and. filepsb(its).eq.filepsb(ts)) its=its-1 enddo endif endif ! check for res2 value if (its.eq.ts) then ! no restore needed in routine sts2restore2() luninfo(ilun)%res2=p_res2_no else luninfo(ilun)%res2=p_res2_yes endif ! restore file position if (filepsb(its).ge.0) then call osd_fseek(luninfo(ilun)%lun,filepsb(its),0) inquire(unit=luninfo(ilun)%lun,name=fname) write(*,'(i5,i16,1x,a)') luninfo(ilun)%lun,filepsb(its), 1trim(fname(index(fname,'\',.true.)+1:)) ! write(*,'(a)') trim(fname) ! write(*,'(6x,a,6i8)') 'restore lun,psb,dat,res2,ts,its:', ! 1 luninfo(ilun)%lun,filepsb(its),filedat(its), ! 1 luninfo(ilun)%res2,ts,its else call osd_fseek(luninfo(ilun)%lun,0,0) endif if (luninfo(ilun)%res2.ne.p_res2_yes) then ! reset data info to unknown for this timestep filedat(ts)=p_dataunk endif enddo endif end subroutine ssts2restore1 c ****************************************************************************** subroutine ssts2restore2(timevalue) c ****************************************************************************** use m_sts2 implicit none double precision, intent(in) :: timevalue ! start date/time of current time step integer ilun integer ts integer(kind=8), pointer :: filepsb(:),filepse(:) integer, pointer :: filedat(:) character(len=256) :: fname c check or sts is in effect if (.not. usests) return c find idices for time values call ssts2timeposition(timevalue,ts,.true.) c for each lun if (restorets) then write(*,*) ' ***** State Save Restore 2: ',ts do ilun=1,nluns ! use pointer filedat => luninfo(ilun)%filedat filepsb => luninfo(ilun)%filepsb filepse => luninfo(ilun)%filepse if (luninfo(ilun)%res2.eq.p_res2_yes) then inquire(unit=luninfo(ilun)%lun,name=fname) ! restore file position if (filepse(ts).ge.0) then call osd_fseek(luninfo(ilun)%lun,filepse(ts),0) write(*,'(i5,i16,1x,a)') luninfo(ilun)%lun,filepse(ts), 1trim(fname(index(fname,'\',.true.)+1:)) ! write(*,'(6x,a,5i8)') 'restore lun,pse,dat,res2,ts:', ! 1 luninfo(ilun)%lun,filepse(ts),filedat(ts), ! 1 luninfo(ilun)%res2,ts else write(*,'(3a,i4,/,9x,a,i10,/,9x,a)') 1 ' WARNING sts_restore2: need to restore file ', 1 trim(fname),', unit ',luninfo(ilun)%lun, 1 'at an unknown fileposition ',filepse(ts), 1 'File pointer will be restored to position 0' call osd_fseek(luninfo(ilun)%lun,0,0) call exit(13) endif ! reset res2 luninfo(ilun)%res2=p_res2_no endif enddo endif ! reset restorets restorets=.false. end subroutine ssts2restore2 c ****************************************************************************** subroutine sts2saverestore(timevalue,stsave,strestore,phase,inoc) c ****************************************************************************** c state save/restore c - argument strestore has no function yet, may be removed later on c - Depending on 'timevalue' this routines determines itself or a StateRestore c has to be done (When timevalue<=lastSavedTime) c - when stsave==TRUE state save will be don, unless a stateRestore is executed use m_sts2 implicit none double precision, intent(in) :: timevalue ! start date/time of current time step logical, intent(in) :: stsave ! .true. save state logical, intent(out) :: strestore ! .true. restore state integer, intent(in) :: phase ! 1: at begin of timestep loop integer,intent(in) :: inoc logical restore c check or sts is in effect if (.not. usests) return c check or restore is necessary if (phase.eq.1) call ssts2nextts(timevalue) c save if wanted write(*,*) ' ===== stsave : ',stsave write(*,*) ' ===== restore : ',restorets write(*,*) ' ===== timevalue: ',timevalue write(*,*) ' ===== phase : ',phase if (stsave .and. .not.restorets) then if (phase.eq.1) currenttime = timevalue call ssts2save(timevalue,phase,inoc) end if c restore if (restorets) then ! select phase if (phase.eq.1) then ! phase 1 - before reading call ssts2restore1(timevalue) ! restore timestep (1) else ! phase 2 - after reading call ssts2restore2(timevalue) ! restore timestep (2) endif endif c set currentTimeStepCalculated to true c phase 2 is at the end of the timestep loop if (phase.eq.2) then currentTimeStepCalculated=.true. endif c set argument value strestore=restorets end subroutine sts2saverestore c ****************************************************************************** subroutine sts2data(lun) c ****************************************************************************** c save for unit lun for current timestep data has been read use m_sts2 implicit none integer, intent(in) :: lun ! logical unit number to subscribe c check or sts is in effect if (.not. usests) return c set value call ssts2setdata(lun,p_data) end subroutine sts2data c ****************************************************************************** subroutine sts2nodata(lun) c ****************************************************************************** c save for unit lun for current timestep data has not been read use m_sts2 implicit none integer, intent(in) :: lun ! logical unit number to subscribe c check or sts is in effect if (.not. usests) return c set value call ssts2setdata(lun,p_datano) end subroutine sts2nodata c ****************************************************************************** subroutine ssts2setdata(lun,value) c ****************************************************************************** c save for unit lun for current timestep data has been read or not use m_sts2 implicit none integer, intent(in) :: lun ! logical unit number to subscribe integer, intent(in) :: value ! value to set to filedat() integer :: ilun,currentts integer, pointer :: filedat(:) c only save data when currentts>0 call ssts2timeposition(currenttime,currentts,.true.) if (currentts.le.0) return c check or lun is defined call ssts2getilun(lun,ilun) c set data if (ilun.gt.0) then ! increase strucure call ssts2increaseStructures(currentts,0) ! check res2, ! if res2 is set to p_res2_yes then do not change the ! value of filedat() ! reason: this unit number is in a restore state and gets data from a different ! timestep then the current timestep filedat => luninfo(ilun)%filedat if (luninfo(ilun)%res2.ne.p_res2_yes) then ! set data ! see m_sts2 for the defined order of 'value' filedat(currentts)=max(value,filedat(currentts)) endif endif end subroutine ssts2setdata c ****************************************************************************** subroutine ssts2increaseStructures(ts,nlun) c ****************************************************************************** use m_sts2 implicit none integer, intent(in) :: ts ! minimum number of time steps in structures integer, intent(in) :: nlun ! minimum number of luns in structures integer :: i,j,mxts,mxln mxts=0 mxln=0 !## increase time array if needed if(ts.gt.0)then if(associated(tsinfo))mxts=size(tsinfo) if(ts.gt.mxts)then allocate(tsinfo_dum(ts)) do i=1,mxts; tsinfo_dum(i)=tsinfo(i); enddo if(associated(tsinfo))deallocate(tsinfo); tsinfo=>tsinfo_dum endif endif !## increase luninfo if needed if(nlun.gt.0.or.ts.gt.0)then mxln=0; if(associated(luninfo))mxln=size(luninfo) !## resize array if(nlun.gt.mxln.or.ts.gt.mxts)then allocate(luninfo_dum(max(nlun,mxln))) mxts=size(tsinfo) do i=1,max(nlun,mxln) allocate(luninfo_dum(i)%filepsb(mxts)) luninfo_dum(i)%filepsb=-1 allocate(luninfo_dum(i)%filepse(mxts)) luninfo_dum(i)%filepse=-1 allocate(luninfo_dum(i)%filedat(mxts)) luninfo_dum(i)%filedat=p_dataunk enddo !## process all unit number stored already do i=1,mxln luninfo_dum(i)%lun =luninfo(i)%lun luninfo_dum(i)%res2=luninfo(i)%res2 !## copy if already part existed if(associated(luninfo(i)%filepsb))then mxts=size(luninfo(i)%filepsb) do j=1,mxts luninfo_dum(i)%filepsb(j)=luninfo(i)%filepsb(j) luninfo_dum(i)%filepse(j)=luninfo(i)%filepse(j) luninfo_dum(i)%filedat(j)=luninfo(i)%filedat(j) enddo endif enddo do i=1,mxln if(associated(luninfo(i)%filepsb))deallocate(luninfo(i)%filepsb) if(associated(luninfo(i)%filepse))deallocate(luninfo(i)%filepse) if(associated(luninfo(i)%filedat))deallocate(luninfo(i)%filedat) enddo if(associated(luninfo))deallocate(luninfo); luninfo=>luninfo_dum endif endif end subroutine ssts2increaseStructures c ****************************************************************************** subroutine ssts2getilun(lun,ilun) c ****************************************************************************** use m_sts2 implicit none integer, intent(in) :: lun ! logical unit number to find integer, intent(out) :: ilun ! >0: position in luninfo structure integer :: i ilun=-1 ! search unitnumber do i=1,nluns if (luninfo(i)%lun.eq.lun) then ilun=i; exit endif enddo end subroutine ssts2getilun !c ****************************************************************************** ! subroutine sts2dump(lun) !c ****************************************************************************** ! use m_sts2 ! implicit none ! integer, intent(in) :: lun ! <=0: all units are displayed ! ! >0: only unit lun will be displayed ! integer i,j,k,n ! integer currentts ! !c check or sts is in effect ! if (.not. usests) return ! !c find idices for time values ! call ssts2timeposition(currenttime,currentts,.true.) ! !c dump ! do i=1,nluns ! ! if (lun.le.0 .or. lun.eq.luninfo(i)%lun) then ! write(*,'(/,i4,'':'')') luninfo(i)%lun ! do j=1,currentts,10 ! n=min(j+9,currentts) ! write(*,'(i5,a,i5)') j,'...',n ! write(*,'('' dat'',10I8)') (luninfo(i)%filedat(k),k=j,n) ! write(*,'('' psb'',10I8)') (luninfo(i)%filepsb(k),k=j,n) ! write(*,'('' pse'',10I8)') (luninfo(i)%filepse(k),k=j,n) ! enddo ! endif ! ! enddo ! ! end subroutine sts2dump !c ****************************************************************************** ! subroutine sts2test(cla,stsave,strestore,timevalue, ! 1 ts1,ts2,tsstep,deltats,option) !c ****************************************************************************** !c routine to test save/restore !c option 1: init test save/restore for timestep 1 !c 2: determine action for NEXT timestep (timestep+1) ! implicit none ! character (len=*), intent(in) :: cla ! command line argument of sts ! logical, intent(inout) :: stsave ! state save ! logical, intent(inout) :: strestore ! state restore ! double precision, intent(in) :: timevalue ! start date/time of current time step ! ! integer, intent(in) :: timestep ! current timestep ! integer, intent(inout) :: ts1,ts2,tsstep ! save/restore/step ! integer, intent(inout) :: deltats ! delta timestep for next timestep ! integer, intent(in) :: option ! option ! integer :: tsplus1,ios,timestep ! logical, save :: usedummy=.true. ! !c check usage of routine ! if (.not.usedummy) return ! !c find idices for time values ! call ssts2timeposition(timevalue ,timestep,.true.) ! !c execute action ! tsplus1=timestep+1 ! select case( option ) ! ! case( 1 ) ! ! init test save/restore ! if (cla(1:5).eq.'test') then ! ! set parameters for test save/restore ! write(*,*) ' MESSAGE: STS-test switched on.' ! ! try to read data from cla ! read(cla(6:),*,iostat=ios) ts1,tsstep ! if (ios.ne.0) then ! tsstep=5 ! ts1=1 ! endif ! ts2=ts1 ! if (timestep.eq.ts1) then ! stsave =.true. ! ts2 =ts2+tsstep ! ! strestore=.false. ! endif ! else ! ! don't use this routine anymore ! usedummy=.false. ! endif ! ! case( 2 ) ! ! check what to do ! if (timestep.eq.ts2) then ! ! when timestep==ts2: restore ! !stsave =.false. ! strestore=.true. ! deltats =min(deltats,ts1-timestep) ! restore data at ts1 ! ts1 =ts1+tsstep ! else ! if (tsplus1.eq.ts1) then ! ! when tsplus1==ts1: save ! stsave =.true. ! !strestore=.false. ! ts2 =ts2+tsstep ! endif ! endif ! ! case default ! ! ! write(*,*) ' ERROR: unknown option in sts_test: ',option ! call exit(10) ! ! end select ! write(*,'(a,2l4,5i8)') ! 1 ' STS TEST: stsave,strestore,timestep,ts1,ts2,tsstep,deltats:', ! 1 stsave,strestore,timestep,ts1,ts2,tsstep,deltats ! ! end subroutine sts2test c ****************************************************************************** subroutine ssts2timeposition(timevalue,timestep,new) c ****************************************************************************** use m_sts2 implicit none double precision, intent(in) :: timevalue !> time value integer, intent(out) :: timestep !> timestep index !! >0: index number !! <=0: not found logical, intent(in) :: new integer :: i double precision d,cfn_mjd_delta c find index number for timevalue timestep=0 do i=1,nts d=cfn_mjd_delta(timevalue,tsinfo(i)%timevalue) ! write(*,'(a,4(1x,f15.7))') 'd,timevalue,timevalue(),i',d, ! 1timevalue,tsinfo(i)%ptsinfo%timevalue,i if (d.eq.0.d0) then timestep=i exit endif enddo ! may be a new number needed if (timestep.eq.0 .and. new) then if (nts.le.0) then ! first usage of these data structures timestep=1 else if (timevalue.gt.tsinfo(nts)%timevalue) then ! new index timestep=nts+1 endif endif endif end subroutine ssts2timeposition c ****************************************************************************** subroutine sts2getlun(lun) c ****************************************************************************** use m_sts2 implicit none integer, intent(out) :: lun !> unit number c get unit number if (usests) then lun=stslun else lun=-1 endif end subroutine sts2getlun c ****************************************************************************** subroutine utl_write_unf_int(lun,x,ndata) c ****************************************************************************** use m_sts2 implicit none integer, intent(in) :: lun ! logical unit number integer, intent(in) :: ndata ! number of data elements integer, intent(in) :: x(ndata) ! data integer :: i write(lun) (x(i),i=1,ndata) end subroutine utl_write_unf_int c ****************************************************************************** subroutine utl_write_unf_real(lun,x,ndata) c ****************************************************************************** use m_sts2 implicit none integer, intent(in) :: lun ! logical unit number integer, intent(in) :: ndata ! number of data elements real , intent(in) :: x(ndata) ! data integer :: i write(lun) (x(i),i=1,ndata) end subroutine utl_write_unf_real c ****************************************************************************** subroutine utl_write_unf_dble(lun,x,ndata) c ****************************************************************************** use m_sts2 implicit none integer, intent(in) :: lun ! logical unit number integer, intent(in) :: ndata ! number of data elements double precision, intent(in) :: x(ndata) ! data integer :: i write(lun) (x(i),i=1,ndata) end subroutine utl_write_unf_dble c ****************************************************************************** subroutine utl_read_unf_int(lun,x,ndata) c ****************************************************************************** use m_sts2 implicit none integer, intent(in) :: lun ! logical unit number integer, intent(in) :: ndata ! number of data elements integer, intent(out) :: x(ndata) ! data integer :: i read(lun) (x(i),i=1,ndata) end subroutine utl_read_unf_int c ****************************************************************************** subroutine utl_read_unf_real(lun,x,ndata) c ****************************************************************************** use m_sts2 implicit none integer, intent(in) :: lun ! logical unit number integer, intent(in) :: ndata ! number of data elements real , intent(out) :: x(ndata) ! data integer :: i read(lun) (x(i),i=1,ndata) end subroutine utl_read_unf_real c ****************************************************************************** subroutine utl_read_unf_dble(lun,x,ndata) c ****************************************************************************** use m_sts2 implicit none integer, intent(in) :: lun ! logical unit number integer, intent(in) :: ndata ! number of data elements double precision, intent(out) :: x(ndata) ! data integer :: i read(lun) (x(i),i=1,ndata) end subroutine utl_read_unf_dble