EnKF/src_conversion/p_prep_obs_ORCA25.F90

program prep_obs
! This takes a model restart file, extracts the desired variable
! and brings into a format that the EnKF V2 can read & treat ('observations.uf').
!
! !!! AGAIN: THIS USES EnKF VERSION 2 !!!
!
! Two command line arguments are expected:
! 1. Path to the nc file of which the data is to be extracted.
! 2. Variable name that can be found in there.
!
! Output is written to 'observations.uf'. If lots of files are to be treated,
! use a shell script to call me and rename output.
!
! Warning: Currently only the data from the first time step are being read.
!          (No problem with restart files, normally)
!
!
! TO DO: Add possibility of treating several obs types.
!
!
! (c) September 2008, Christof König Beatty, Christof.KonigBeatty@uclouvain.be
! (c) May 2009, generalized by same.
! (c) Jun 2011, simplified by F. Massonnet
! (c) April 2016, cleaned by F. Massonnet


  use mod_measurement
  use netcdf

  implicit none


  ! NetCDF vars
  character(len=99)            :: filename, cvarerror, cvarerroru, cvarerrorv
  integer                      :: ncid, ierr, nvar, narg
  logical                      :: ex
  character(len=16)            :: varname, varnameu, varnamev      ! Type of measurement ('a_i_htcX', 'u_ice', 'v_ice', maybe 'v_i_htcX')
  character(len=80), parameter :: maskfile = 'mask-ORCA25.nc' !hc!
  ! Data vars
  integer, parameter           :: nx=1442, ny=1050 ! Unfortunately, still hard coded.
  real, dimension(nx,ny)       :: lons, lats, obsfld, errorfld, obsfldu, obsfldv, errorfldu, errorfldv
  integer, dimension(nx,ny)    :: mask
  integer                      :: obscnt       ! Counts nr of obs's.
  ! Other vars
  character(len=99)               dummy        ! To read cmd line args
  ! for loops (haha)
  integer                      :: i, j, varID, icomp
  ! Ice thickness category stuff
   integer, parameter :: nhtc=5     !hc! nr of ice thickn. cat's
   real(KIND=8)      :: rdate

  ! Obs stuff
  type (measurement), allocatable :: obs(:)
  ! Need to fill:
  ! d (val), var (error var), id (iiceconc etc.), lon, lat, depths,
  ! ipic, jpic (i/j-pivot point in grid), ns (support, 0: point meas.), 
  ! a1-4 (bilin. coeff), status (not used)
  narg= iargc()

  PRINT *, narg
  
  if (narg<=1) then
     ! Write info
     write(*,*)
     write(*,*) " (prep_obs) takes a real obs, extracts the desired variable and outputs"
     write(*,*) " it in a format that the EnKF can read & treat ('observations.uf')."
     write(*,*)
     write(*,*) " A file named mask.nc containing the variables tmaskutil, nav_lon and nav_lat"
     write(*,*) " is expected to be in the current directory (ORCA-file)"
     write(*,*) 
     write(*,*) " Three command line arguments are expected:"
     write(*,*) " 1. Path to the nc file of which the data is to be extracted."
     write(*,*) " 2. Variable name that can be found in there, 'h_i_htc1' or"
     write(*,*) "    'sic'. or dxdy_ice"
     write(*,*) " 3. A tag with the date, e.g. 19790520"
     write(*,*)
     write(*,*) " Hope to see you again soon."
     write(*,*)
     stop "(prep_obs): Stopped."
  end if
  ! Command line arguments
  nvar=narg-1  
  ! Get them
  call getarg(1, dummy)
  read(dummy,'(A)') filename
  ! Some parameter checking
  inquire(file=trim(filename), exist=ex)
  if (.not.ex) then
     print *, '(p_prep_obs): file does not exist: '// trim(filename)
     stop
  end if
  
  ! Get mask, lons & lats
  ! open the maskfile
  ierr = nf90_open(trim(maskfile),nf90_NoWrite,ncid);   if (ierr.ne.nf90_noerr) call handle_err(ierr, "opening mask file")

  ! Find VarID of tmask & get values
  ierr = nf90_inq_varid(ncid, 'tmask', varID); if (ierr.ne.nf90_noerr) call handle_err(ierr, "inquiring varID tmask")
  ierr = nf90_get_var(ncid, varID, mask) ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "getting variable tmaks")
  ! Find VarID of longitude & get vals
  ierr = nf90_inq_varid(ncid, 'nav_lon', varID); if (ierr.ne.nf90_noerr) call handle_err(ierr, "inquiring varID nav_lon")
  ierr = nf90_get_var(ncid, varID, lons)       ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "getting variable nav_lon")
  ! Find VarID of latitude & get vals
  ierr = nf90_inq_varid(ncid, 'nav_lat', varID); if (ierr.ne.nf90_noerr) call handle_err(ierr, "inquiring varID nav_lat")
  ierr = nf90_get_var(ncid, varID, lats)       ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "getting variable nav_lat")

  ! Close maskfile
  ierr = nf90_close(ncid)
  if (ierr.ne.nf90_noerr) call handle_err(ierr, "closing")
  
  allocate( obs(nvar*sum(mask)), STAT=ierr )
  if (ierr.ne.0) call handle_err(ierr, "allocating obs") !no netcdf-error! ohwell.
  
  obscnt = 0 ! Keeps track of nr of added obs's.
  
  call getarg(2, dummy)
  read(dummy,'(A)') varname
  call getarg(3, dummy)
  read(dummy,*) rdate
  if (      trim(varname).ne.'sic' ) then
    write(*,*) ' (prep_obs) Currently only the total sea ice concentration'
    write(*,*) '            sic    can be processed                      '
    stop " (prep_obs) Aborted"
  endif

  ! User info
  write(*,*) ' (prep_obs) Extracting '//trim(varname)//' from '//trim(filename)

  ! Some preparations
  obsfld=0.
  ! open the netCDF file
  ierr = nf90_open(trim(filename),nf90_NoWrite,ncid) ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "opening")

  ! Read observation data
  ! Find VarID of varname
  ierr = nf90_inq_varid(ncid, trim(varname), varID) ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "inquiring varID")
  ! get values
  ierr = nf90_get_var(ncid, varID, obsfld) ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "getting variable")

  ! Find VarID of varname

  cvarerror='error_'//varname

   ierr = nf90_inq_varid(ncid, cvarerror, varID) ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "inquiring varID")
   ! get values
   ierr = nf90_get_var(ncid, varID, errorfld) ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "getting variable")

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



  ! User info - ADAPT ACCORDINGLY
  write (*,*) " (prep_obs) Using data >40N and <45S"

  ! Loop over space
  do i=1,size(mask,1)
    do j=1,size(mask,2)

      ! Pick out ocean points where data is existent
      if ( errorfld(i,j) .GE. 1.0 ) then

        ! Limit 'obs' to >40°N or <45°S
        if ( (lats(i,j).ge.40) .or. (lats(i,j).le.-45) ) then
          obscnt = obscnt + 1
          obs(obscnt)%d    = obsfld(i,j)/100.0
          obs(obscnt)%lon  = lons(i,j)
          obs(obscnt)%lat  = lats(i,j)
          obs(obscnt)%ipiv = i ! the i-point of the grid of the model
          obs(obscnt)%jpiv = j ! the j-point of the grid of the model

          ! Put other data into obs type array
          obs(obscnt)%var    = (errorfld(i,j)/100.0)**2 ! The variance
          obs(obscnt)%id     = TRIM(varname)
          obs(obscnt)%depths = 0
          obs(obscnt)%ns     = 0
          obs(obscnt)%a1     = 1
          obs(obscnt)%a2     = 0
          obs(obscnt)%a3     = 0
          obs(obscnt)%a4     = 0
          obs(obscnt)%status      = .TRUE.
          obs(obscnt)%i_orig_grid = -1
          obs(obscnt)%j_orig_grid = -1
          obs(obscnt)%h           = -1.0
          obs(obscnt)%date        = rdate
          obs(obscnt)%orig_id     =  0

        end if ! >40°N or <50°S
      end if ! ocean point
    end do ! j
  end do ! i

  !Write data:
  inquire(iolength=i)obs(1)
  open (11, file='observations.uf', status='replace',form='unformatted', access='direct', recl=i)
  do j=1,obscnt
     write(11,rec=j)obs(j)
  enddo
  close(11)

PRINT *, rdate

!  !XXX: 
!  ! Write data to textfile, for visual checking
  open(12, file='observations.txt')
  do j=1,obscnt
     write(12,FMT=53) obs(j)
53   FORMAT(f8.4,X,f8.4,X,a8,X,2(f10.5,X),f4.2,X,2(I3,X),I1,X,4(f5.2,X),L,X,2(I3,X),f5.2,X,I8,X,I1)
  enddo
  close(12)
!  !:XXX

  ! Tell user how many obs there were
  write(*,*) " (prep_obs) Wrote out this many obs: ", obscnt
  write(*,*) " (prep_obs) Number of ocean points : ", sum(mask)

  ! Cleanup
  if (allocated(obs)) deallocate(obs,STAT=ierr)

  write(*,*) ' (prep_obs) End successfully reached'; write(*,*)


contains

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

end program prep_obs