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,1) :: errorfld, obsfld, obsfldu, obsfldv, errorfldu, errorfldv 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,1) :: errorfld, obsfld, obsfldu, obsfldv, errorfldu, errorfldv ! real, dimension(ox,oy) :: olons, olats 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, varID_lon, varID_lat, varID_err, varID_var, 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") WRITE(*,*)'tmask varID:',varID 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") WRITE(*,*)'nav_lon varID:',varID 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 = 100.0 errmin = 0.001 errmax = 40.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. errorfld=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_lon) ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "inq varID lon obs") WRITE(*,*)'lon obs varID:',varID_lon ierr = nf90_get_var(ncid, varID_lon, olons) ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "getting var lon obs") ! Find VarID of latitude & get vals ierr = nf90_inq_varid(ncid, 'lat', varID_lat) ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "inq varID lat obs") ierr = nf90_get_var(ncid, varID_lat, olats) ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "getting var lat obs") ! Find VarID of varname (variable to assim) WRITE(*,*) "varname="//trim(varname) ierr = nf90_inq_varid(ncid, trim(varname), varID_var) ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "inq varID obs") ! get values ierr = nf90_get_var(ncid, varID_var, obsfld) ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "getting var obs") obsfld=obsfld/100. 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_err) ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "inq varID obs") !! get values ierr = nf90_get_var(ncid, varID_err, errorfld) ; if (ierr.ne.nf90_noerr) call handle_err(ierr, "getting error var obs") errorfld=errorfld/100. ! 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 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) ! Pick out ocean points where data is existent IF ( obsfld(i,j,1) .LT. obsmax & .AND. obsfld(i,j,1) .GT. obsmin & .AND. ( (errorfld(i,j,1)) ** 2 ) .GT. ( errmin ** 2 ) & .AND. ( (errorfld(i,j,1)) ** 2 ) .LT. ( errmax ** 2 ) & ) 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,1) 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,1))**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(*,*) 'Nb obscnt: ', obscnt WRITE(*,*) 'Nb obscnt2: ', obscnt2 !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,f9.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