EnKF/EnKF-MPI-TOPAZ/m_io_mod_fld.F90

module m_io_mod_fld
! Get or put one of the fields of a restart file, specified by
! ensemble number, field name and type, and vertical level. The 
! time level is currently not used (restart files have only one)
! but who knows. Grid dimension is also needed, as well as if you
! want to 'get' or 'put'.
!
! This replaces the two routines 'm_get_mod_fld.F90' and m_put_mod_fld.F90'.
! There was so much overlap that it became easier to merge the two. I think.
!
! (c) July 2009, Christof.KonigBeatty@uclouvain.be

  use netcdf
#if defined (QMPI)
  use qmpi
#else
  use qmpi_fake
#endif

  private handle_err


contains 
  subroutine io_mod_fld(fld,k,enslist,cfld,type,vlevel,tlevel,nx,ny,gorp,rdate_obs)

  implicit none

  ! In/out
  real,dimension(nx,ny),intent(inout):: fld    ! output fl
  integer,                intent(in) :: k      ! Index to enslist
  integer,dimension(:),   intent(in) :: enslist! List of existing ensemble members
  character(len=*),       intent(in) :: cfld   ! name of fld
  integer,                intent(in) :: type   ! which file to use
  integer,                intent(in) :: vlevel ! vertical level (ignored)
  integer,                intent(in) :: tlevel ! time level (ignored)
  integer,                intent(in) :: nx,ny  ! Grid dimension
  character(len=3),       intent(in) :: gorp   ! 'get' or 'put' (sorry, couldn't come up with anything better)
  real(kind=8),           intent(in) :: rdate_obs
  ! NetCDF vars
  integer           :: iens             ! Ensemble member to read
  character(len=99) :: fcfile, anafile, cfile
  integer           :: ncid, varID, error
  logical           :: exfc, exan
  ! Other
  character(len=3)  :: cmem
  integer           :: zvlevel   ! for i/o ocean variable
  real(kind=8)      :: rdate_mod
  ! Find iens withing enslist
  iens = enslist(k)

  ! Create filename dep. on type of variable/parameter requested
  write(cmem,'(i3.3)') 100+iens  ! iens=1 gives cmem = 101
  select case( type )
  case(1) ! ice variable
     fcfile ='forecast_ice_'//cmem//'.nc'
     anafile='analysis_ice_'//cmem//'.nc'
  case(2) ! ocean variable
     fcfile ='forecast_oce_'//cmem//'.nc'
     anafile='analysis_oce_'//cmem//'.nc'
     zvlevel = max(vlevel,1)
  case(3) ! ice namelist parameter
     if (master) print *, '(io_mod_fld): ice parameter writing not implemented yet!'
     call stop_mpi()
  case(4) ! ocean namelist parameter
     if (master) print *, '(io_mod_fld): ocean parameter writing not implemented yet!'
     call stop_mpi()
  case default
     if (master) print *, '(io_mod_fld): variable type not understood!'
     call stop_mpi()
  end select


  ! If the fc file exists we turn it into the analysis file (unless that's already there).
  inquire(file=fcfile,  exist=exfc)
  inquire(file=anafile, exist=exan)
  if ((.not.exfc).and.(.not.exan)) then       ! Neither file is there
    if (master) print *, '(io_mod_fld): Restart file '//cmem//' missing!'
    call stop_mpi()
  elseif (exfc.and.(.not.exan)) then          ! fcfile here but no anafile
!     call system('mv '//trim(fcfile)//' '//trim(anafile) )  ! "operational" to save space
     call system('cp '//trim(fcfile)//' '//trim(anafile) )  ! for debugging
  end if
  
  ! Decide on which file to use
  if (gorp=='get') cfile=fcfile
  if (gorp=='put') cfile=anafile

  ! ckb prefers only one file at the time, so take care of this special case
  inquire(file=fcfile,  exist=exfc)
  if (.not.exfc) cfile=anafile

!!$  !XXX:
!!$  write(*,*) "XXX: "
!!$  write(*,*) "XXX: iens           : ", iens
!!$  write(*,*) "XXX: cfld           : ", cfld
!!$  write(*,*) "XXX: type           : ", type
!!$  write(*,*) "XXX: nx, ny, zvlevel: ", nx, ny, zvlevel
!!$  write(*,*) "XXX: fcfile         : ", trim(fcfile)
!!$  write(*,*) "XXX: anafile        : ", trim(anafile)
!!$  write(*,*) "XXX: shape(fldIO)   : ", shape(fldIO)
!!$  write(*,*) "XXX: "
!!$  !:XXX
  
  ! open the netCDF file
  error = nf90_open(trim(cfile),nf90_Write,ncid); if (error.ne.nf90_noerr) call handle_err(error, "opening")

  ! Find VarID of cfld
  error = nf90_inq_varid(ncid, trim(cfld), varID); if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID")

  ! Put/Get
  select case( type )
  case(1, 3, 4) ! 2D
     if     (gorp=='get') then
        error = nf90_get_var(ncid, varID, fld); if (error.ne.nf90_noerr) call handle_err(error, "getting 2D variable")
     elseif (gorp=='put') then
        error = nf90_put_var(ncid, varID, fld); if (error.ne.nf90_noerr) call handle_err(error, "putting 2D variable")
     else
        if (master) print *, "(io_mod_fld): Either 'put' or 'get'!"
        call stop_mpi()
     endif
  case(2) ! 3D ocean variable
     if     (gorp=='get') then
        error = nf90_get_var(ncid, varID, fld, start=(/1,1,zvlevel/), count=(/nx,ny,1/))
        if (error.ne.nf90_noerr) call handle_err(error, "getting ocean variable")
     elseif (gorp=='put') then
        error = nf90_put_var(ncid, varID, fld, start=(/1,1,zvlevel/), count=(/nx,ny,1/))
        if (error.ne.nf90_noerr) call handle_err(error, "putting ocean variable")
     else
        if (master) print *, "(io_mod_fld): Either 'put' or 'get'!"
        call stop_mpi()
     endif
  end select

  !if (master) PRINT *, " Find VarID of cfld "
  error = nf90_inq_varid(ncid, 'time_counter', varID); if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID")
  error = nf90_get_var(ncid, varID, rdate_mod); if (error.ne.nf90_noerr) call handle_err(error, "getting ocean variable")

  ! Close file
  error = nf90_close(ncid); if (error.ne.nf90_noerr) call handle_err(error, "closing")

  ! Check date mode and date obs
  IF (INT(rdate_mod) .NE. INT(rdate_obs)) THEN
     !PRINT *, 'date mod not egal to date obs, stop, (',INT(rdate_mod),' ',INT(rdate_obs),')'
     !STOP 1
  END IF
end subroutine  io_mod_fld


subroutine handle_err(status, infomsg) 
  integer,            intent ( in) :: status 
  character(len = *), intent ( in), optional :: infomsg
  if(status /= nf90_noerr) then
     if (master) then
        if (present(infomsg)) then
           print *, 'Error while '//infomsg//' - '//trim(nf90_strerror(status)) 
        else
           print *, trim(nf90_strerror(status)) 
        endif ! opt arg
        print *,'(io_mod_fld)'
     endif ! only master outputs
     call stop_mpi()
  end if ! check error status
end subroutine handle_err

end module  m_io_mod_fld