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