! subroutine read_oldsrc ! ! Search for and read old cray-blocked style source history: ! (this routine executes only on unicos machines, i.e., cannot ! read old source histories on non-unicos platforms) ! use input_module,only: tempdir,source,source_start,output,start use hist_module,only: nsource,h,nhist implicit none include "params.h" include "cons.h" include "buff.h" include "fcom.h" include "fgcom.h" include "index.h" include "dynphi.h" ! ! Local: integer(kind=8) :: longint integer :: isearch,istrt,nwds,j,ixk,iyk, | iend_iyd,iend_day,iend_hr,iend_min,idelmin,nxk,i,k, | iendonly,ibeg_iyd,ibeg_day,ibeg_hr,ibeg_min,nphik,nb integer :: isorc,mtime(3),ios integer(kind=8) :: iheader(512) ! kind=8 for the equiv's real :: summary(100),dum,fmin,fmax,fminj,fmaxj character(len=80) :: dskfile,mssfile real :: hdr_mag(4),hdr_sdtide(10),hdr_dtide(2),hdr_hp, | hdr_cp,hdr_byimf,hdr_colfac,hdr_f107d,hdr_f107a ! ! Equivalences for old history header: equivalence (hdr_mag ,iheader(126)) equivalence (hdr_sdtide ,iheader(132)) equivalence (hdr_dtide ,iheader(154)) equivalence (hdr_hp ,iheader(166)) equivalence (hdr_cp ,iheader(167)) equivalence (hdr_byimf ,iheader(168)) equivalence (hdr_colfac ,iheader(230)) ! assumes naur=3 equivalence (hdr_f107d ,iheader(231)) ! assumes naur=3 equivalence (hdr_f107a ,iheader(232)) ! assumes naur=3 ! ! External: integer,external :: readirec,readfrec real,external :: vsum ! if (nsource==1) then mssfile = source mtime = source_start else mssfile = output(1) mtime(:) = start(:,1) endif isorc = 2 call mkdiskflnm(mssfile,dskfile) call getms(mssfile,'SOURCE',tempdir,dskfile) open(unit=isorc,file='SOURCE',status='OLD',form='UNFORMATTED', | iostat=ios) if (ios == 0) then write(6,"('read_oldsrc: opened file SOURCE with isorc=',i3)") | isorc else write(6,"(/,'>>> read_oldsrc: error opening file SOURCE')") stop 'read_oldsrc' endif c c In this case, search 1st given OUTPUT volume. c (acquires and assigns were done in sub input) c longint = -1 istrt = isorc write(6,"('Searching for start time ',i3,',',i2,',',i2, + ' on SOURCE vol ',a)") mtime,mssfile ! ! Read header: isearch = 1 75 continue nwds = readirec(istrt,iheader,512,3,0) if (nwds /= 512) then write(6,"('>>> START: ERROR reading history header from', + ' unit ',i2,' nwds=',i5)") istrt,nwds stop 'i/o' endif ! ! Read summary: nwds = readfrec(istrt,summary,100,2,0) if (nwds /= 100) then write(6,"('>>> Error reading summary: istrt=',i3,' nwds=',i6)") + istrt,nwds stop 'i/o' endif ! ! Report time read: write(6,"('Reading time ',i3,':',i2,':',i2,' Seeking time ', + i3,':',i2,':',i2)") iheader(2:4),mtime c if (iheader(2)==mtime(1).and.iheader(3)==mtime(2).and. | iheader(4)==mtime(3)) then ! history found ! ! Update gpi vars: h%f107d = hdr_f107d h%f107a = hdr_f107a h%hpower = hdr_hp h%ctpoten = hdr_cp ! write(6,"('read_oldsrc: h%f107d,a=',2e12.4, ! | ' h%hpower,ctpoten=',2e12.4)") h%f107d,h%f107a,h%hpower, ! | h%ctpoten goto 140 endif ! ! Read past unwanted history: do j=1,jmax nwds = readfrec(istrt,dum,1,2,0) enddo isearch = isearch+1 go to 75 140 continue ! ! Set nhist (hist_mod.f): nhist = isearch ! ! History found -- read into global shared fg-array: ! (retain moulo, init, and flip temporarilly, for testing, ! see also putf call below) ! do 240 j=1,jmax call flip ! keep flipping f-array temporarilly for testing nwds = readfrec(istrt,fg(1,1,j,ixtimep),lbscm,2,0) if (nwds <= 0) then write(6,"('>>> start: error reading latitude slice j=', | i3)") j stop 'i/o' endif ! ! Change from "external" to "internal" format: nxk = ncols+1 ixk = lbtap do nb=1,ncols iyk = imaxp4 nxk = nxk-1 do i=1,imaxp2 fg(iyk,nxk,j,ixtimep) = fg(ixk,1,j,ixtimep) ixk = ixk-1 iyk = iyk-1 enddo enddo C **** ADD PERIODIC POINTS do nxk = 1,ncols-1 fg( 1,nxk,j,ixtimep) = fg(imax+1,nxk,j,ixtimep) fg( 2,nxk,j,ixtimep) = fg(imax+2,nxk,j,ixtimep) fg(imax+3,nxk,j,ixtimep) = fg( 3,nxk,j,ixtimep) fg(imax+4,nxk,j,ixtimep) = fg( 4,nxk,j,ixtimep) enddo ! ! Model tgcm13 accepted histories that did not have dynamo potential. ! In this case, the following conditional was true, and the data was ! necessarilly shifted to make room for the potential. Because now ! "all" histories currently used contain potential, this ability is ! removed from tgcm13mt. ! if (nwds < lbtap) then write(6,"(/'>>> start: nwds < lbtap: lenh=',i5,' lbtap=',i5)") | nwds,lbtap write(6,"(4x,'tgcm13mt does not accept histories without', | ' dynamo potential.',/,4x,'Please use tgcm13.')") stop 'start' endif ! ! Define dynpot: nphik = nphi do k = 1,kmaxp1 nphik = nphik+1 do i = 1,imax+1 dynpot(i,j,k) = fg(i+2,nphik,j,ixtimep) enddo enddo ! ! End j read loop: 240 CONTINUE ! ! Insert polar points in dynpot: ! do k = 1,kmaxp1 dynpot(1,0,k) = (9.*vsum(imaxg,dynpot(1,1,k),1)- | vsum(imaxg,dynpot(1,2,k),1))/(8.*float(imaxg)) dynpot(1,jmaxgp,k) = (9.*vsum(imaxg,dynpot(1,jmaxg,k),1) | -vsum(imaxg,dynpot(1,jmaxg-1,k),1))/(8.*float(imaxg)) do i = 2,imaxgp dynpot(i,0,k) = dynpot(1,0,k) dynpot(i,jmaxgp,k) = dynpot(1,jmaxgp,k) enddo enddo ! ! Top level of NOPNM is not defined in old cray-blocked histories. ! Define it here with spval so it won't crash in min,max routines: do j=1,jmax fg(:,nopnm+zkmxp,j,ixtimep) = spval enddo ! ! Debug lat min,max: ! write(6,"('read_oldsrc: call fgmnmx after reading old', ! | ' cray-blocked source: ixtimep=',i2)") ixtimep ! do j=1,jmax ! call fgmnmx(0,'DUM',1,ndisk,j,ixtimep) ! enddo ! ! Debug 3d min,max: ! do i=1,ndisk ! fmin = 1.e36 ! fmax = -1.e36 ! do j=1,jmax ! call fminmaxspv(fg(1,ndexa(i+1)+1,j,ixtimep),zimxp*zkmxp, ! | fminj,fmaxj,spval) ! if (fminj < fmin) fmin = fminj ! if (fmaxj > fmax) fmax = fmaxj ! enddo ! write(6,"('read_oldsrc: Read field ',a,' 3d min,max = ', ! | 2e12.4)") nflds_lab(i),fmin,fmax ! enddo end subroutine read_oldsrc !----------------------------------------------------------------------- integer function readirec(lu,iaddr,nwds,iconv,irew) implicit none ! ! Read an integer record from lu into iaddr for nwds words, ! with conversion flag iconv (iconv used only if SGI). ! Return number of words actually read. ! If irew==1, rewind unit prior to reading. ! ! Args: integer,intent(in) :: lu,nwds,iconv,irew integer,intent(in) :: iaddr(nwds) ! ! Locals: real :: stat integer :: ier ! ! Externals: ! if (irew==1) rewind(lu) ! ! This is not executed unless we are on unicos. This pre-proc ! conditional allows this file to be compiled on non-unicos machines. ! #ifdef UNICOS buffer in (lu,1) (iaddr,iaddr(nwds)) stat = unit(lu) #else stat = -1. #endif if (stat.eq.0.) then write(6,"('>>> readirec: EOF on unit ',i2)") lu readirec = int(stat) elseif (stat.eq.1.) then write(6,"('>>> readirec: ERROR reading from unit ',i2)") lu readirec = int(stat) else readirec = nwds endif end function readirec !------------------------------------------------------------------- integer function readfrec(lu,faddr,nwds,iconv,irew) implicit none ! ! Read a floating record from lu into faddr for nwds words, ! with conversion flag iconv (iconv used only if SGI). ! Return number of words actually read. ! If irew==1, rewind unit prior to reading. ! ! Args: integer,intent(in) :: lu,nwds,iconv,irew real,intent(in) :: faddr(nwds) ! ! Locals: real :: stat integer :: ier ! ! Externals: ! if (irew==1) rewind(lu) ! ! This is not executed unless we are on unicos. This pre-proc ! conditional allows this file to be compiled on non-unicos machines. ! #ifdef UNICOS buffer in (lu,1) (faddr,faddr(nwds)) ! ! return codes for unit(): -1.0 -> success ! 0.0 -> EOF ! 1.0 -> error stat = unit(lu) #endif if (stat.eq.0.) then write(6,"('>>> readfrec: EOF on unit ',i2,' nwds=',i6)") | lu,nwds readfrec = int(stat) elseif (stat.eq.1.) then write(6,"('>>> readfrec: ERROR reading from unit ',i2)") lu readfrec = int(stat) else readfrec = nwds endif end function readfrec