EnKF/conversion_uf/TMP/p_prep_obs_ORCA1.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', 'siconc' if model forced by itself)
!!!!!TAG-DEV-AD : change mask-ORCA1.nc into mask-eORCA1.nc
!  character(len=80), parameter :: maskfile = 'mask-ORCA1.nc' !hc!
  character(len=80), parameter :: maskfile = 'mask-eORCA1.nc' !hc!
! Data vars
!!!!!TAG-DEV-AD : adapted coordinate for eORCA1, mlon mlat mx my for model dim
  integer, parameter           :: mx=360, my=331 ! Unfortunately, still hard coded.
  real, dimension(mx,my)       :: mlons, mlats
!!!!!TAG-DEV-AD : add OSI-SAF dim = obs dim = ox,oy
  integer, parameter           :: ox=432, oy=432 ! Unfortunately, still hard coded.
  real, dimension(ox,oy)       :: olons, olats, errorfld, obsfld, obsfldu, obsfldv, errorfldu, errorfldv

  REAL                         :: obsmin, obsmax, errmin, errmax
  REAL                         :: latmin_NH = 40.0 !76.0 Assim only circle more or less with Svlbard latitudes
  REAL                         :: latmax_NH = 90.0 !82.0
  REAL                         :: latmin_SH = 40.0 ! Same for SH (sign will be added)
  REAL                         :: latmax_SH = 90.0
  integer, dimension(mx,my)    :: mask
  integer                      :: obscnt, obscnt2, cpt, cptemin, cptemax, cptomin, cptomax       ! 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, mlons)       ; 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, mlats)       ; 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) .eq. 'rfb' ) THEN
    WRITE(*,*) "Handling ", TRIM(varname)
    ! Min and max values for error used to screen the data (any data with
    ! standard deviation in between those values will be selected
    obsmin = 0.01
    obsmax = 10.0
    errmin = 0.01
    errmax = 1.0
  ELSEIF ( TRIM(varname) .eq. 'vt_i') THEN
    WRITE(*,*) "Handling ", TRIM(varname)
    obsmin = 0.01
    obsmax = 50.0
    errmin = 0.01
    errmax = 1.0
  ELSEIF ( TRIM(varname) .eq. 'at_i') THEN
    WRITE(*,*) "Handling ", TRIM(varname)
    obsmin = 0.001
    obsmax = 1.0
    errmin = 0.001
    errmax = 0.5
  ELSEIF ( TRIM(varname) .eq. 'siconc') THEN
    WRITE(*,*) "Handling ", TRIM(varname)
    obsmin = 0.001
    obsmax = 1.0
    errmin = 0.001
    errmax = 0.5
  ELSEIF ( TRIM(varname) .eq. 'ice_conc') THEN ! OSI-SAF-450
    WRITE(*,*) "Handling ", TRIM(varname)
    obsmin = 0.0
    obsmax = 101.0
    errmin = 0.0001
    errmax = 60.0
  ELSE
    WRITE(*,*) " (prep_obs) Currently only the variables rfb  (sea ice freeboard),"
    WRITE(*,*) "                                         vt_i (total sea ice volume)"
    WRITE(*,*) "                                         at_i (total sea ice concentration)"
    WRITE(*,*) "            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 longitude & get vals from obs
  ierr = nf90_inq_varid(ncid, 'lon', varID)           ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "inquiring varID lon obs")
  ierr = nf90_get_var(ncid, varID, olons)             ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "getting variable lon obs")
  ! Find VarID of latitude & get vals
  ierr = nf90_inq_varid(ncid, 'lat', varID)           ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "inquiring varID lat obs")
  ierr = nf90_get_var(ncid, varID, olats)             ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "getting variable lat obs")
  ! Find VarID of varname (variable to assim)
  WRITE(*,*) "varname="//trim(varname)
  ierr = nf90_inq_varid(ncid, trim(varname), varID)   ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "inquiring varID obs")
  ! get values
  ierr = nf90_get_var(ncid, varID, obsfld)            ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "getting var in obs")
  WRITE(*,*) "OK var"//trim(varname)
  ! Find VarID of varname (error)
  cvarerror=TRIM(varname)//'_sd'
  
  WRITE(*,*) "cvarerror="//trim(varname)//"_sd"
  ierr = nf90_inq_varid(ncid, cvarerror, varID)       ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "inquiring varID obs")
  !! get values
  ierr = nf90_get_var(ncid, varID, errorfld)          ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "getting error var obs")
  ! 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"

  obscnt2=0
  cpt=0
  cptemin=0
  cptemax=0
  cptomin=0
  cptomax=0

  WRITE(*,*) 'obs min:',obsmin,'obs max:',obsmax,'err min:',errmin,'err max:',errmax
  ! Loop over space
  DO i = 1, SIZE(olats, 1)
    DO j = 1, SIZE(olats, 2)
     IF (obsfld(i,j) .GT. 0.0 ) THEN
      cpt = cpt + 1
      !WRITE(*,*) 'error(i,j): ', i, j, errorfld(i,j)
      !WRITE(*,*) 'obs(i,j): ',i,j, obsfld(i,j)
      !WRITE(*,*) 'obs: ', (obsfld(i,j))
     END IF


      ! Pick out ocean points where data is existent
      !IF (       (errorfld(i,j)) .GT. errmin  &
      !     .AND. (errorfld(i,j)) .LT. errmax  &
      !     .AND. obsfld(i,j) .GT. obsmin                &
       !    .AND. obsfld(i,j) .LT. obsmax &
       !    ) THEN
       
      IF ( (errorfld(i,j)**2) .GT. (errmin**2) ) THEN
              cptemin = cptemin + 1
              WRITE(*,*) 'error: ', errmin**2, (errorfld(i,j)**2), errmax**2
              PAUSE

              IF ( (errorfld(i,j)**2) .LT. errmax**2 ) THEN
                cptemax = cptemax + 1
                !IF ( (obsfld(i,j)) .GT. obsmin ) THEN
                !        cptomin = cptomin + 1
                !        IF ( (obsfld(i,j)) .LT. obsmax ) THEN
                !                cptomax = cptomax + 1

      !     .AND. (errorfld(i,j)) .LT. errmax  &
      !     .AND. obsfld(i,j) .GT. obsmin                &
       !    .AND. obsfld(i,j) .LT. obsmax &
       !    ) THEN 
       
       obscnt2 = obscnt2 + 1
        ! Limit 'obs' spatially
        IF (     (       olats(i,j) .GE. latmin_NH        & 
                  .AND.  olats(i,j) .LE. latmax_NH      ) &
             .OR.(       olats(i,j) .LE. (-latmin_SH)     &
                  .AND.  olats(i,j) .GE. (-latmax_SH)   ) &
           )  THEN 


          obscnt = obscnt + 1
          obs(obscnt)%d    = obsfld(i,j)
          obs(obscnt)%lon  = olons(i,j)
          obs(obscnt)%lat  = olats(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))**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 ! err max
        END IF ! Spatial selection
      END IF ! if valid point
    END DO ! j
  END DO ! i
  WRITE(*,*) 'Nb obscnt: ', obscnt
  WRITE(*,*) 'Nb obscnt2: ', obscnt2
  WRITE(*,*) 'Compteur obs > 0: ', cpt
  WRITE(*,*) 'Compteur t err min: ', cptemin
  WRITE(*,*) 'Compteur t err max: ', cptemax
  WRITE(*,*) 'Compteur t obs min: ', cptomin
  WRITE(*,*) 'Compteur t obs max: ', cptomax

  !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)

  ! 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)

  ! 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