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-ORCA1.nc' !hc! ! Data vars integer, parameter :: nx=362, ny=292 ! Unfortunately, still hard coded. real, dimension(nx,ny) :: lons, lats, obsfld, errorfld, obsfldu, obsfldv, errorfldu, errorfldv REAL :: obsmin, obsmax, errmin, errmax REAL :: latmin_NH = 40.0 REAL :: latmax_NH = 90.0 REAL :: latmin_SH = 40.0 ! Same for SH (sign will be added) REAL :: latmax_SH = 90.0 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) .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 ELSE WRITE(*,*) " (prep_obs) Currently only the variables rfb (sea ice freeboard)," WRITE(*,*) " at_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 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=TRIM(varname)//'_sd' 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)) ** 2 .GT. errmin ** 2 & .AND. (errorfld(i,j)) ** 2 .LT. errmax ** 2 & .AND. obsfld(i,j) .GT. obsmin & .AND. obsfld(i,j) .LT. obsmax ) THEN ! Limit 'obs' spatially IF ( ( lats(i,j) .GE. latmin_NH & .AND. lats(i,j) .LE. latmax_NH ) & .OR.( lats(i,j) .LE. (-latmin_SH) & .AND. lats(i,j) .GE. (-latmax_SH) ) & ) THEN obscnt = obscnt + 1 obs(obscnt)%d = obsfld(i,j) 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))**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 ! Spatial selection END IF ! if valid 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) ! 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