history.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

! Arguments
!
      integer, intent(in) :: nrg            ! unit number
      integer, intent(in) :: luhrest        ! unit number
!
! Local workspace
!
      integer t, f                     ! tape, field indices
      integer c                        ! chunk or lat index
      integer lenc                     ! length of useful character data
      integer numlev                   ! number of vertical levels (dimension and loop)
      integer begver                   ! on-node vert start index
      integer endver                   ! on-node vert end index
      integer ioerr                    ! error code from read()
      integer coldimin                 ! column dimension of model array
      integer begdim3                  ! on-node chunk or lat start index
      integer enddim3                  ! on-node chunk or lat end index
      integer ncol                     ! number of active columns per chunk
      integer lenarr                   ! global size of array to be read

      character(len=80)  :: locfn       ! Local filename

      integer, pointer :: nacs(:,:)    ! accumulation counter

      type (hbuffer_3d) :: hbuf             ! history buffer
      integer,  pointer :: fullnacs(:,:)    ! accumulation buffer
      type (dim_index_3d) :: dimind    ! 3-D dimension index

      if (masterproc) then

         read (nrg,iostat=ioerr) rgnht, mtapes, varid, fincl, fexcl
         if (ioerr /= 0 ) then
            write (6,*) 'READ ioerror ',ioerr,' on i/o unit = ',nrg
            call endrun
         end if

         do t=1,mtapes
            read (nrg,iostat=ioerr) nhtfrq(t), nflds(t), nfils(t), mfilt(t), &
                                    nfpath(t), cpath(t), nhfil(t), &
                                    nhstpr(t), ndens(t), ncprec(t), beg_time(t)
            if (ioerr /= 0) then
               write (6,*) 'READ ioerror ',ioerr,' on i/o unit = ',nrg
               call endrun
            end if

            do f=1,nflds(t)
               read (nrg,iostat=ioerr) tape(t)%hlist(f)%field%name,        &
                                       tape(t)%hlist(f)%field%long_name,   &
                                       tape(t)%hlist(f)%field%units,       &
                                       tape(t)%hlist(f)%field%numlev,      &
                                       tape(t)%hlist(f)%field%decomp_type, &
                                       tape(t)%hlist(f)%avgflag,           &
                                       tape(t)%hlist(f)%time_op,           &
                                       tape(t)%hlist(f)%hbuf_prec,         &
                                       tape(t)%hlist(f)%hwrt_prec
               if (ioerr /= 0) then
                  write(6,*)'READ_RESTART_HISTORY: ', &
                     'End or error condition reading history restart field ', &
                     f,' from tape ',t
                  write(6,*)'ioerr=',ioerr
                  call endrun
               end if
            end do
            if (rgnht(t)) then
               read (nrg,iostat=ioerr) hrestpath(t)
               if (ioerr /= 0) then
                  write (6,*) 'READ ioerror on read of filename ',ioerr,' on i/o unit = ',nrg
                  call endrun
               end if
            end if
         end do

      end if

#if ( defined SPMD )
      call mpibcast (rgnht  ,ptapes   ,mpilog ,0,mpicom)

      call mpibcast (mtapes ,1        ,mpiint ,0,mpicom)
      call mpibcast (nhtfrq ,ptapes   ,mpiint ,0,mpicom)
      call mpibcast (nflds  ,ptapes   ,mpiint ,0,mpicom)
      call mpibcast (nfils  ,ptapes   ,mpiint ,0,mpicom)
      call mpibcast (mfilt  ,ptapes   ,mpiint ,0,mpicom)

      do t=1,mtapes
         do f=1,nflds(t)
            call mpibcast (tape(t)%hlist(f)%field%numlev,     1,        mpiint, 0,mpicom)
            call mpibcast (tape(t)%hlist(f)%field%decomp_type,1,        mpiint, 0,mpicom)
            call mpibcast (tape(t)%hlist(f)%field%name,       8,        mpichar,0,mpicom)
            call mpibcast (tape(t)%hlist(f)%field%units,      max_chars,mpichar,0,mpicom)
            call mpibcast (tape(t)%hlist(f)%avgflag,          1,        mpichar,0,mpicom)
            call mpibcast (tape(t)%hlist(f)%hbuf_prec,        1,        mpiint, 0,mpicom)
            call mpibcast (tape(t)%hlist(f)%hwrt_prec,        1,        mpiint, 0,mpicom)
         end do
      end do
#endif
!
! Allocate space for history buffers and initialize
!
      do t=1,mtapes
         do f=1,nflds(t)
            numlev   = tape(t)%hlist(f)%field%numlev
            select case (tape(t)%hlist(f)%field%decomp_type)
            case (phys_decomp)
               tape(t)%hlist(f)%field%begdim3 = begchunk
               tape(t)%hlist(f)%field%enddim3 = endchunk
               tape(t)%hlist(f)%field%begver = 1
               tape(t)%hlist(f)%field%endver = numlev
               allocate (tape(t)%hlist(f)%field%colperdim3(begchunk:endchunk))
               do c=begchunk,endchunk
                  ncol = get_ncols_p(c)
                  tape(t)%hlist(f)%field%colperdim3(c) = ncol
               end do
               tape(t)%hlist(f)%field%coldimin = pcols
            case (dyn_decomp)
               tape(t)%hlist(f)%field%begdim3 = beglat
               tape(t)%hlist(f)%field%enddim3 = endlat
               if ( dycore_is('LR') )then
# if ( defined STAGGERED )
                  select case (numlev)
                  case (1)
                     tape(t)%hlist(f)%field%begver = 1
                     tape(t)%hlist(f)%field%endver = 1
                  case (plev)
                     tape(t)%hlist(f)%field%begver = beglev
                     tape(t)%hlist(f)%field%endver = endlev
                  case (plevp)
                     tape(t)%hlist(f)%field%begver = beglev
                     tape(t)%hlist(f)%field%endver = endlevp
                  case default
                     write(6,*)'READ_RESTART_HISTORY: invalid number of levels=', numlev
                     call endrun ()
                  end select
# endif
               else
                  tape(t)%hlist(f)%field%begver = 1
                  tape(t)%hlist(f)%field%endver = numlev
               endif
               allocate (tape(t)%hlist(f)%field%colperdim3(beglat:endlat))
               do c=beglat,endlat
                  tape(t)%hlist(f)%field%colperdim3(c) = nlon(c)
               end do
               tape(t)%hlist(f)%field%coldimin = plon
            case default
               write(6,*)'READ_RESTART_HISTORY: bad decomp_type=',tape(t)%hlist(f)%field%decomp_type
               call endrun ()
            end select

            coldimin = tape(t)%hlist(f)%field%coldimin
            begdim3  = tape(t)%hlist(f)%field%begdim3
            enddim3  = tape(t)%hlist(f)%field%enddim3
            begver   = tape(t)%hlist(f)%field%begver
            endver   = tape(t)%hlist(f)%field%endver

            dimind = dim_index_3d (1,coldimin,begver,endver,begdim3,enddim3)
            call allocate_hbuf (tape(t)%hlist(f)%hbuf,dimind,tape(t)%hlist(f)%hbuf_prec)
            tape(t)%hlist(f)%hbuf = 0._r8
            allocate (tape(t)%hlist(f)%nacs(coldimin,begdim3:enddim3))
            tape(t)%hlist(f)%nacs(:coldimin,begdim3:enddim3) = 0
         end do
      end do
!
!-----------------------------------------------------------------------
! Read history restart files
!-----------------------------------------------------------------------
!
! Loop over the total number of history files declared and
! read the pathname for any history restart files
! that are present (if any). Test to see if the run is a restart run
! AND if any history buffer regen files exist (rgnht=.T.). Note, rgnht 
! is preset to false, reset to true in routine WSDS if hbuf restart files
! are written and saved in the master restart file. Each history buffer
! restart file is then obtained.
! Note: some f90 compilers (e.g. SGI) complain about I/O of 
! derived types which have pointer components, so explicitly read each one.
! 
      do t=1,mtapes
         if (rgnht(t)) then
!
! Open history restart file
!
            if (masterproc) then
               call getfil (hrestpath(t), locfn)
               call opnfil (locfn, luhrest, 'u')
            end if
!
! Read history restart file
!
            do f=1,nflds(t)
               coldimin   =  tape(t)%hlist(f)%field%coldimin
               begdim3    =  tape(t)%hlist(f)%field%begdim3
               enddim3    =  tape(t)%hlist(f)%field%enddim3
               numlev     =  tape(t)%hlist(f)%field%numlev
               nacs       => tape(t)%hlist(f)%nacs(:coldimin,begdim3:enddim3)
#ifdef SPMD
               if (masterproc) then
                  dimind = dim_index_3d (1,plon,1,numlev,1,plat)
                  call allocate_hbuf (hbuf,dimind,tape(t)%hlist(f)%hbuf_prec)
                  allocate (fullnacs(plon,plat))
               else
                  call assoc_hbuf_with_nothing (hbuf,tape(t)%hlist(f)%hbuf_prec)
                  fullnacs => nothing_int
               end if

               select case (tape(t)%hlist(f)%field%decomp_type)
               case (phys_decomp)
                  if (masterproc) then
                     call read_hbuf (hbuf,luhrest,ioerr)
                     read (luhrest) fullnacs
                  end if
                  call scatter_field_to_chunk_hbuf (1, numlev, 1, plon, hbuf, tape(t)%hlist(f)%hbuf)
                  call scatter_field_to_chunk_int (1, 1, 1, plon, fullnacs, tape(t)%hlist(f)%nacs)
               case (dyn_decomp)
                  if ( dycore_is('LR') )then
# if ( defined STAGGERED )
! NEW LR CODING
                     if (tape(t)%hlist(f)%hbuf_prec == 8) then
                        lenarr = plon*numlev*plat
                        select case (numlev)
                        case (1)
                           call lrreadin(luhrest, strip2d, tape(t)%hlist(f)%hbuf%buf8,  &
                                         lenarr, 2)
                        case (plev)
                           call lrreadin(luhrest, strip3dxzy, tape(t)%hlist(f)%hbuf%buf8,  &
                                         lenarr, 3)
                        case (plevp)
                           call lrreadin(luhrest, strip3dxzyp, tape(t)%hlist(f)%hbuf%buf8,  &
                                         lenarr, 3)
                        case default
                           write(6,*)'READ_RESTART_HISTORY: invalid number of levels=', numlev
                           call endrun ()
                        end select
                        lenarr = plon*plat
                        call lrreadini(luhrest, strip2d, nacs,  &
                                      lenarr, 2)
                     else
                        lenarr = plon*numlev*plat
                        select case (numlev)
                        case (1)
                           call lrreadin4(luhrest, strip2d, tape(t)%hlist(f)%hbuf%buf4,  &
                                         lenarr, 2)
                        case (plev)
                           call lrreadin4(luhrest, strip3dxzy, tape(t)%hlist(f)%hbuf%buf4,  &
                                         lenarr, 3)
                        case (plevp)
                           call lrreadin4(luhrest, strip3dxzyp, tape(t)%hlist(f)%hbuf%buf4,  &
                                         lenarr, 3)
                        case default
                           write(6,*)'READ_RESTART_HISTORY: invalid number of levels=', numlev
                           call endrun ()
                        end select
                        lenarr = plon*plat
                        call lrreadini(luhrest, strip2d, nacs,  &
                                      lenarr, 2)
                     endif
# endif
                  else
                     call readin_hbuf(luhrest, tape(t)%hlist(f)%hbuf, coldimin*numlev)
                     call readin_int (luhrest, nacs, coldimin)
                  endif
               case default
                  write(6,*)'READ_RESTART_HISTORY: bad decomp_type=',tape(t)%hlist(f)%field%decomp_type
                  call endrun ()
               end select

               if (masterproc) then
                  call deallocate_hbuf (hbuf)
                  deallocate (fullnacs)
               else
                  call nullify_hbuf (hbuf)
                  nullify (fullnacs)
               end if
#else
               select case (tape(t)%hlist(f)%field%decomp_type)
               case (phys_decomp)
                  dimind = dim_index_3d (1,plon,1,numlev,1,plat)
                  call allocate_hbuf (hbuf,dimind,tape(t)%hlist(f)%hbuf_prec)
                  call read_hbuf (hbuf,luhrest,ioerr)
                  call scatter_field_to_chunk_hbuf (1, numlev, 1, plon, hbuf, tape(t)%hlist(f)%hbuf)
                  call deallocate_hbuf (hbuf)
                  allocate (fullnacs(plon,plat))
                  read (luhrest) fullnacs
                  call scatter_field_to_chunk_int (1, 1, 1, plon, fullnacs, tape(t)%hlist(f)%nacs)
                  deallocate (fullnacs)
               case (dyn_decomp)
                  call readin_hbuf(luhrest, tape(t)%hlist(f)%hbuf, coldimin*numlev)
                  call readin_int (luhrest, nacs, coldimin)
               case default
                  write(6,*)'READ_RESTART_HISTORY: bad decomp_type=',tape(t)%hlist(f)%field%decomp_type
                  call endrun ()
               end select
#endif
            end do
!          
! Done reading this history restart file
!
            if (masterproc) close (luhrest)

         end if  ! rgnht(t)
      end do     ! end of do mtapes loop
!
! If the history files are partially complete (contain less than
! mfilt(t) time samples, then get the files and open them.
!
      do t=1,mtapes
         if (masterproc .and. nfils(t) > 0) then
            call getfil (cpath(t), locfn)
            call wrap_open (locfn, NF_WRITE, nfid(t))
            call h_inquire (t)
         end if
!
! If the history file is full, close the current unit
!
         if (nfils(t) >= mfilt(t)) then
            if (masterproc) then
               write(6,*)'READ_RESTART_HISTORY: nf_close(',nfid(t),')=',nhfil(t)
               call wrap_close (nfid(t))
            end if
            nfils(t) = 0
         end if
      end do
!
! set flag indicating h-tape contents are now defined (needed by addfld)
!
      htapes_defined = .true.      

      return
   end subroutine read_restart_history

!#######################################################################

   character(len=nlen) function get_hfilepath( tape )
!----------------------------------------------------------------------- 
! 
! Purpose: Return full filepath of history file for given tape number
! This allows public read access to the filenames without making
! the filenames public data.
!
!----------------------------------------------------------------------- 
  integer, intent(in) :: tape  ! Tape number

  get_hfilepath = cpath( tape )
  end function get_hfilepath