Observation operators

This documentation summarises the observation operator in HARMONIE and the use of the HOP_DRIVER tool. The test harness, HOP_DRIVER, calls the observation operator and generates FG departures without calling any model code or initialising any model modules. Firstly, the IFS is used to dump a single-observation gom_plus to file from the 1st trajectory of an experiment. Dumping multiple observations would require a more complex and full-featured dump (good file format, multi-process parallel). For code refactoring HOP_DRIVER can be used to test changes to the observation operator of a particular observation type.

HARMONIE and HOP_DRIVER

The HOP_DRIVER program was first added to CY42R2 code. The tool was initially implemented to test refactoring of the IFS observation operator code src/arpifs/op_obs/hop.F90. Instructions on how to prepare the code and run HOP_DRIVER using HARMONIE are outlined below. Presentation made at [wiki:HirlamMeetings/ModelMeetings/ObOpWorkshop OOPS Observation Operator Workshop] may provide some useful background information.

Comments on the branch

  • Code changes were required in order to compile cy42r2bf.04 + mods (provided by MF/ECMWF) in the HARMONIE system: [14312], [14325], [14326], [14330], [14331], [14332], [14333], [14334].
  • Changes were made to makeup in order to compile HOP_DRIVER correctly: [14310], [14327], [14328], [14329], [14335], [14362], [14382], [14392].
  • Included in [14362] is a change to ODBSQLFLAGS which is set to ODBSQLFLAGS=-O3 -C -UCANARI -DECMWF $(ODBEXTRAFLAGS) in order to use ECMWF flavoured ODB used by HOP_DRIVER
  • On cca GNU compilers 4.9 are not fully supported, ie I had to build GRIB-API and NetCDF locally using gcc/gfortran 4.9 on cca
  • An environment variable, HOPDIR, is used to define the location of necessary input data for HOP_DRIVER
  • An environment variable, HOPCOMPILER, is used by the HOP_driver script to define the compiler used. This is used to compare results.

HOPOBS: amsua

Currently there is only one observation type, AMSU-A (HOPOBS=amsua), available for testing with HOP_DRIVER. Alan Geer (ECMWF) has already carried out the refactoring of the HOP code related to AMSU-A observations. A single observation is provided in the ECMA and is used to test the refactoring of the HOP code. To carry out the testing of the amsua refactoring HOPOBS should be set to amsua in ecf/config_exp.h.

reportype@hdrobstype@hdrsensor@hdrstatid@hdrstalt@hdrdate@hdrtime@hdrdegrees(lat)degrees(lon)report_status@hdrdatum_status@bodyobsvalue@bodyvarno@bodyvertco_type@body
100773' 4'832800!20140131215914-29.59060.3113112173.281193
100773' 4'832800!20140131215914-29.59060.3113112158.861193
100773' 4'832800!20140131215914-29.59060.311313227.401193
100773' 4'832800!20140131215914-29.59060.311313260.821193
100773' 4'832800!20140131215914-29.59060.311311256.901193
100773' 4'832800!20140131215914-29.59060.311311239.601193
100773' 4'832800!20140131215914-29.59060.3113112NULL1193
100773' 4'832800!20140131215914-29.59060.311313217.691193
100773' 4'832800!20140131215914-29.59060.311311209.391193
100773' 4'832800!20140131215914-29.59060.311311214.051193
100773' 4'832800!20140131215914-29.59060.311311223.021193
100773' 4'832800!20140131215914-29.59060.311311234.421193
100773' 4'832800!20140131215914-29.59060.311311245.141193
100773' 4'832800!20140131215914-29.59060.311311257.181193
100773' 4'832800!20140131215914-29.59060.3113112227.911193

HOP_DRIVER

Using HOP_DRIVER

With LHOP_RESULTS=.TRUE. HOP_DRIVER will write results to a file called hop_results${MYPROC} for comparison between online and offline results. (The results file is opened by src/arpifs/var/taskob.F90. HOP_DRIVER results are written to hop_results${MYPROC} in src/arpifs/op_obs/hop.F90:

 :
 :
IF(LHOP_RESULTS) THEN
!$OMP CRITICAL
  ! Output for comparison between online and offline results:
  WRITE(CFILENAME,'("hop_results",I4.4)') MYPROC
  OPEN(NEWUNIT=IU,FILE=CFILENAME,POSITION='APPEND',ACTION='WRITE',FORM='FORMATTED')
  DO JOBS = 1,KDLEN
    DO JBODY=1,IMXBDY
      IF (JBODY>ICMBDY(JOBS)) CYCLE
      IBODY = ROBODY%MLNKH2B(JOBS)+(JBODY-1)
      WRITE(IU,'(6I8,2F30.14)') MYPROC, KSET, JOBS, NINT(ROBHDR%DATA(JOBS,ROBHDR%SEQNO_AT_HDR)),&
        & NINT(ROBODY%DATA(IBODY,ROBODY%VERTCO_REFERENCE_1_AT_BODY)), &
        & NINT(ROBODY%DATA(IBODY,ROBODY%VARNO_AT_BODY)), ZHOFX(JOBS,JBODY), ZXPPB(JOBS,JBODY)

    ENDDO
  ENDDO
  CLOSE(IU)
!$OMP END CRITICAL
ENDIF
 :
 :

The HOPdriver script (based a script provided by MF) sorts the contents of the `hopresults0001` file for comparison with some results made available by ECMWF/MF:

 :
 :
#
# Check HOP_DRIVER results (available for gfotran and intel)
#
ln -s $HOPDIR/${HOPOBS}/results.$HOPCOMPILER .
cat hop_results* | sort -k1,1n -k2,2n -k3,3n -k5,5n -k6,6n > results.driver
echo
cmp -s results.$HOPCOMPILER results.driver
if [ $? -eq 0] ; then
  echo "RESULTS ARE STRICTLY IDENTICAL TO THE REFERENCE FOR HOPCOMPILER=$HOPCOMPILER :-)"
else
  echo Compare exactly against the results dumped from hop:
  echo "xxdiff results.$HOPCOMPILER results.driver &"
  diff results.$HOPCOMPILER results.driver
  exit 1
fi
 :
 :

On cca you will find useful output from HOP_DRIVER in cca:$TEMP/hm_home/rfexp/archive/HOPDRIVEROUT:

fort.4
NODE.001_01
hop_results0001
results.gfortran
results.driver

The code

HOP_DRIVER is a short program written by Deborah Salmond (ECMWF) to test code changes made to the observation operator. The program src/arpifs/programs/hop_driver.F90 is summarised here.

  • The program sets up the model geometry and observations:
 :
 :
CALL GEOMETRY_SET(YRGEOMETRY)
CALL MODEL_SET(YRMODEL)

CALL IFS_INIT('gc7a')

CALL SUINTDYN

CALL SUGEOMETRY(YRGEOMETRY)        !From GEOMETRY_SETUP

CALL SURIP(YRGEOMETRY%YRDIM)             !From MODEL_CREATE

! Set up Observations, Sets
CALL SUDIMO(YRGEOMETRY,NULOUT)     !From SU0YOMB
CALL SUOAF              !From SU0YOMB
CALL SUALOBS            !From SU0YOMB
CALL SURINC             !From SU0YOMB
CALL SETUP_TESTVAR      !From SU0YOMB
CALL SUOBS(YRGEOMETRY)              !From CNT1
CALL ECSET(-1,NOBTOT,0) !From OBSV
CALL SUPHEC(YRGEOMETRY,NULOUT)

! Setup varbc (from cnt1.F90) and read VARBC.cycle
CALL YVARBC%SETUP_TRAJ
 :
 :
  • HOP_DRIVER then loops over the number of observation sets (NSETOT) and reads a GOM PLUS for each observation set. HRETR and HOP are then called:
 :
 :
DO ISET=1,NSETOT
  IDLEN   = MLNSET(ISET)
  IMXBDY = MAX(MMXBDY(ISET),1)

  ALLOCATE(ZHOFX(IDLEN,IMXBDY))
  ZHOFX=RMDI

  ! READ GOM_PLUS FROM DUMP
  CALL GOM_PLUS_READ_DUMP(YGP5,ISET)

  IF(IDLEN /= YGP5%NDLEN) THEN
    CALL ABOR1('Sets are incompatible')
  ENDIF

  :
  :
  :

  CALL HRETR(YRGEOMETRY%YRDIMV,IDLEN,IMXBDY,ISET,1,YGP5,YVARBC)

  CALL HOP(YRGEOMETRY%YRDIMV,YGP5,YVARBC,IDLEN,IMXBDY,ISET,1,LDOOPS=.TRUE.,PHOFX=ZHOFX)

  !write(0,*)'ZHOFX',ZHOFX
  DEALLOCATE(ZHOFX)

  CALL GOM_PLUS_DESTROY(YGP5)

ENDDO

 :
 :