Changeset 5682 for branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS

Timestamp:

2015-08-12T17:46:45+02:00 (9 years ago)

Author:

mattmartin

Message:

OBS simplification changes committed to branch after running SETTE tests to make sure we get the same results as the trunk for ORCA2_LIM_OBS.

Location:

branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS

Files:

• diaobs.F90 (modified) (29 diffs)
• obs_grid.F90 (modified) (4 diffs)
• obs_mpp.F90 (modified) (1 diff)
• obs_oper.F90 (modified) (32 diffs)
• obs_prep.F90 (modified) (22 diffs)
• obs_profiles_def.F90 (modified) (5 diffs)
• obs_read_altbias.F90 (modified) (4 diffs)
• obs_read_prof.F90 (modified) (33 diffs)
• obs_read_surf.F90 (modified) (22 diffs)
• obs_readmdt.F90 (modified) (8 diffs)
• obs_surf_def.F90 (modified) (5 diffs)
• obs_types.F90 (modified) (1 diff)
• obs_write.F90 (modified) (31 diffs)

branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

@@ -6 +6 @@
    !!======================================================================
 
-   !!----------------------------------------------------------------------
-   !!   'key_diaobs' : Switch on the observation diagnostic computation
    !!----------------------------------------------------------------------
    !!   dia_obs_init : Reading and prepare observations
@@ -15 +13 @@
    !!   fin_date     : Compute the final date YYYYMMDD.HHMMSS
    !!----------------------------------------------------------------------
-   !! * Modules used   
+   !! * Modules used
    USE wrk_nemo                 ! Memory Allocation
    USE par_kind                 ! Precision variables
@@ -21 +19 @@
    USE par_oce
    USE dom_oce                  ! Ocean space and time domain variables
-   USE obs_fbm, ONLY: ln_cl4    ! Class 4 diagnostic switch
-   USE obs_read_prof            ! Reading and allocation of observations (Coriolis)
-   USE obs_read_surf            ! Reading and allocation of SLA observations  
+   USE obs_read_prof            ! Reading and allocation of profile obs
+   USE obs_read_surf            ! Reading and allocation of surface obs
    USE obs_readmdt              ! Reading and allocation of MDT for SLA.
    USE obs_prep                 ! Preparation of obs. (grid search etc).
@@ -45 +42 @@
       &   dia_obs_dealloc  ! Deallocate dia_obs data
 
-   !! * Shared Module variables
-   LOGICAL, PUBLIC, PARAMETER :: &
-#if defined key_diaobs
-      & lk_diaobs = .TRUE.   !: Logical switch for observation diangostics
-#else
-      & lk_diaobs = .FALSE.  !: Logical switch for observation diangostics
-#endif
-
    !! * Module variables
-   LOGICAL, PUBLIC :: ln_t3d         !: Logical switch for temperature profiles
-   LOGICAL, PUBLIC :: ln_s3d         !: Logical switch for salinity profiles
-   LOGICAL, PUBLIC :: ln_sla         !: Logical switch for sea level anomalies 
-   LOGICAL, PUBLIC :: ln_sst         !: Logical switch for sea surface temperature
-   LOGICAL, PUBLIC :: ln_seaice      !: Logical switch for sea ice concentration
-   LOGICAL, PUBLIC :: ln_vel3d       !: Logical switch for velocity component (u,v) observations
-   LOGICAL, PUBLIC :: ln_ssh         !: Logical switch for sea surface height
-   LOGICAL, PUBLIC :: ln_sss         !: Logical switch for sea surface salinity
-   LOGICAL, PUBLIC :: ln_sstnight    !: Logical switch for night mean SST observations
-   LOGICAL, PUBLIC :: ln_nea         !: Remove observations near land
-   LOGICAL, PUBLIC :: ln_altbias     !: Logical switch for altimeter bias  
-   LOGICAL, PUBLIC :: ln_ignmis      !: Logical switch for ignoring missing files
-   LOGICAL, PUBLIC :: ln_s_at_t      !: Logical switch to compute model S at T observations
-
-   REAL(KIND=dp), PUBLIC :: dobsini   !: Observation window start date YYYYMMDD.HHMMSS
-   REAL(KIND=dp), PUBLIC :: dobsend   !: Observation window end date YYYYMMDD.HHMMSS
-  
-   INTEGER, PUBLIC :: numobtypes   !: Number of observation types to read in.
-   INTEGER, PUBLIC :: n1dint       !: Vertical interpolation method
-   INTEGER, PUBLIC :: n2dint       !: Horizontal interpolation method 
-   INTEGER, DIMENSION(:), ALLOCATABLE :: nvarsprof !Number of profile variables
-   INTEGER, DIMENSION(:), ALLOCATABLE :: nextrprof !Number of profile extra variables
-   INTEGER, DIMENSION(:), ALLOCATABLE :: nvarssurf !Number of surface variables
-   INTEGER, DIMENSION(:), ALLOCATABLE :: nextrsurf !Number of surface extra variables
+   LOGICAL, PUBLIC :: ln_diaobs   !: Logical switch for the obs operator
+   LOGICAL :: ln_sstnight         !: Logical switch for night mean SST obs
+
+   INTEGER :: nn_1dint       !: Vertical interpolation method
+   INTEGER :: nn_2dint       !: Horizontal interpolation method
    INTEGER, DIMENSION(imaxavtypes) :: &
-      & dailyavtypes !: Data types which are daily average
-
-   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: surfdata   ! Initial surface data
-   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: surfdataqc ! Surface data after quality control
-   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: profdata   ! Initial profile data
-   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: profdataqc ! Profile data after quality control
-
-   CHARACTER(len=6),   PUBLIC, DIMENSION(:),   ALLOCATABLE :: obstypesprof
-   CHARACTER(len=6),   PUBLIC, DIMENSION(:),   ALLOCATABLE :: obstypessurf
-
-
-      
-   INTEGER, PARAMETER :: MaxNumFiles = 1000
-   
-   LOGICAL, DIMENSION(MaxNumFiles) :: &
-      & ln_profb_ena, & !: Is the feedback files from ENACT data ?
-   !                    !: If so use dailyavtypes
-      & ln_profb_enatim !: Change tim for 820 enact data set.
-   
-   LOGICAL, DIMENSION(MaxNumFiles) :: &
-      & ln_velfb_av   !: Is the velocity feedback files daily average?
-   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
-      & ld_enact     !: Profile data is ENACT so use dailyavtypes
-   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
-      & ld_velav     !: Velocity data is daily averaged
-   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
-      & ld_sstnight  !: SST observation corresponds to night mean
+      & nn_profdavtypes      !: Profile data types representing a daily average
+   INTEGER :: nproftypes     !: Number of profile obs types
+   INTEGER :: nsurftypes     !: Number of surface obs types
+   INTEGER, DIMENSION(:), ALLOCATABLE :: &
+      & nvarsprof, &         !: Number of profile variables
+      & nvarssurf            !: Number of surface variables
+   INTEGER, DIMENSION(:), ALLOCATABLE :: &
+      & nextrprof, &         !: Number of profile extra variables
+      & nextrsurf            !: Number of surface extra variables
+
+   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: &
+      & surfdata, &          !: Initial surface data
+      & surfdataqc           !: Surface data after quality control
+   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: &
+      & profdata, &          !: Initial profile data
+      & profdataqc           !: Profile data after quality control
+
+   CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: &
+      & cobstypesprof, &     !: Profile obs types
+      & cobstypessurf        !: Surface obs types
 
    !!----------------------------------------------------------------------
@@ -136 +99 @@
 
       !! * Local declarations
-      CHARACTER(len=128) :: profbfiles(MaxNumFiles)
-      CHARACTER(len=128) :: sstfbfiles(MaxNumFiles)
-      CHARACTER(len=128) :: slafbfiles(MaxNumFiles)
-      CHARACTER(len=128) :: seaicefbfiles(MaxNumFiles)
-      CHARACTER(len=128) :: velfbfiles(MaxNumFiles)
-      CHARACTER(LEN=128) :: bias_file
-      CHARACTER(LEN=20)  :: datestr=" ", timestr=" "
-
-      NAMELIST/namobs/ln_t3d, ln_s3d, ln_sla, ln_sss, ln_ssh,         &
-         &            ln_sst, ln_seaice, ln_vel3d,                    &
+      INTEGER, PARAMETER :: &
+         & jpmaxnfiles = 1000    ! Maximum number of files for each obs type
+      INTEGER, DIMENSION(:), ALLOCATABLE :: &
+         & ifilesprof, &         ! Number of profile files
+         & ifilessurf            ! Number of surface files
+      INTEGER :: ios             ! Local integer output status for namelist read
+      INTEGER :: jtype           ! Counter for obs types
+      INTEGER :: jvar            ! Counter for variables
+      INTEGER :: jfile           ! Counter for files
+
+      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: &
+         & cn_profbfiles, &      ! T/S profile input filenames
+         & cn_sstfbfiles, &      ! Sea surface temperature input filenames
+         & cn_slafbfiles, &      ! Sea level anomaly input filenames
+         & cn_sicfbfiles, &      ! Seaice concentration input filenames
+         & cn_velfbfiles         ! Velocity profile input filenames
+      CHARACTER(LEN=128) :: &
+         & cn_altbiasfile        ! Altimeter bias input filename
+      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &
+         & clproffiles, &        ! Profile filenames
+         & clsurffiles           ! Surface filenames
+
+      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles
+      LOGICAL :: ln_s3d          ! Logical switch for salinity profiles
+      LOGICAL :: ln_sla          ! Logical switch for sea level anomalies 
+      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature
+      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration
+      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs
+      LOGICAL :: ln_nea          ! Logical switch to remove obs near land
+      LOGICAL :: ln_altbias      ! Logical switch for altimeter bias
+      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files
+      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs
+      LOGICAL :: llvar1          ! Logical for profile variable 1
+      LOGICAL :: llvar2          ! Logical for profile variable 1
+      LOGICAL :: llnightav       ! Logical for calculating night-time averages
+
+      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS
+      REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS
+      REAL(wp), DIMENSION(jpi,jpj) :: &
+         & zglam1, &             ! Model longitudes for profile variable 1
+         & zglam2                ! Model longitudes for profile variable 2
+      REAL(wp), DIMENSION(jpi,jpj) :: &
+         & zgphi1, &             ! Model latitudes for profile variable 1
+         & zgphi2                ! Model latitudes for profile variable 2
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
+         & zmask1, &             ! Model land/sea mask associated with variable 1
+         & zmask2                ! Model land/sea mask associated with variable 2
+
+      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              &
+         &            ln_sst, ln_sic, ln_vel3d,                       &
          &            ln_altbias, ln_nea, ln_grid_global,             &
-         &            ln_grid_search_lookup, ln_cl4,                  &
+         &            ln_grid_search_lookup,                          &
          &            ln_ignmis, ln_s_at_t, ln_sstnight,              &
-         &            ln_profb_ena, ln_profb_enatim,                  &
-         &            profbfiles, slafbfiles, sssfbfiles,             &
-         &            sshfbfiles, sstfbfiles, seaicefbfiles,          &
-         &            velfbfiles, bias_file, grid_search_file,        &
-         &            dobsini, dobsend, n1dint, n2dint,               &
-         &            nmsshc, mdtcorr, mdtcutoff,                     &
-         &            grid_search_res, dailyavtypes
-
-      INTEGER :: jtype
-      INTEGER :: ios                 ! Local integer output status for namelist read
-      INTEGER, DIMENSION(:), ALLOCATABLE :: jnumfilesprof
-      INTEGER, DIMENSION(:), ALLOCATABLE :: jnumfilessurf
-      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: obsfilesprof
-      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: obsfilessurf
-      LOGICAL :: lmask(MaxNumFiles)
+         &            cn_profbfiles, cn_slafbfiles,                   &
+         &            cn_sstfbfiles, cn_sicfbfiles,                   &
+         &            cn_velfbfiles, cn_altbiasfile,                  &
+         &            cn_gridsearchfile, rn_gridsearchres,            &
+         &            rn_dobsini, rn_dobsend, nn_1dint, nn_2dint,     &
+         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             &
+         &            nn_profdavtypes
+
       !-----------------------------------------------------------------------
       ! Read namelist parameters
       !-----------------------------------------------------------------------
 
-      profbfiles(:) = ''
-      slafbfiles(:) = ''
-      sstfbfiles(:) = ''
-      seaicefbfiles(:) = ''
-      velfbfiles(:) = ''
-      dailyavtypes(:) = -1
-      dailyavtypes(1) = 820
-      ln_profb_ena(:) = .FALSE.
-      ln_profb_enatim(:) = .TRUE.
-      ln_velfb_av(:) = .FALSE.
-      ln_ignmis = .FALSE.
-
-      CALL ini_date( dobsini )
-      CALL fin_date( dobsend )
-
-      ! Read Namelist namobs : control observation diagnostics
-      REWIND( numnam_ref )              ! Namelist namobs in reference namelist : Diagnostic: control observation
+      ! Some namelist arrays need initialising
+      cn_profbfiles(:) = ''
+      cn_slafbfiles(:) = ''
+      cn_sstfbfiles(:) = ''
+      cn_sicfbfiles(:) = ''
+      cn_velfbfiles(:) = ''
+      nn_profdavtypes(:) = -1
+
+      CALL ini_date( rn_dobsini )
+      CALL fin_date( rn_dobsend )
+
+      ! Read namelist namobs : control observation diagnostics
+      REWIND( numnam_ref )   ! Namelist namobs in reference namelist
       READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
 
-      REWIND( numnam_cfg )              ! Namelist namobs in configuration namelist : Diagnostic: control observation
+      REWIND( numnam_cfg )   ! Namelist namobs in configuration namelist
       READ  ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
       IF(lwm) WRITE ( numond, namobs )
 
-      !Set up list of observation types to be used
-      numproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) )
-      numsurftypes = COUNT( (/ln_sla, ln_sss, ln_sst, ln_seaice /) )
-      IF ( numproftypes > 0 ) THEN
-      
-         ALLOCATE( obstypesprof(numproftypes) )
-         ALLOCATE( jnumfilesprof(numproftypes) )
-         ALLOCATE( obsfilesprof(numproftypes, MaxNumFiles) )
-      
-         DO jtype = 1, numproftypes
-            IF (ln_t3d .OR. ln_s3d) THEN
-               obsfilesprof(:,jtype) = profbfiles(:)
-               obstypesprof(jtype) = 'prof'
-            ENDIF
-            IF (ln_vel3d) THEN
-               obsfilesprof(:,jtype) = velfbfiles(:)
-               obstypesprof(jtype) = 'vel'
-            ENDIF
-         
-            lmask(:) = .FALSE.
-            WHERE (obsfilesprof(jtype,:) /= '') lmask(:) = .TRUE.
-            jnumfilesprof(jtype) = COUNT(lmask)
-         END DO
-         
-      ENDIF
-      
-      IF ( numsurftypes > 0 ) THEN
-      
-         ALLOCATE( obstypessurf(numsurftypes) )
-         ALLOCATE( jnumfilessurf(numproftypes) )
-         ALLOCATE( obsfilessurf(numsurftypes, MaxNumFiles) )
-         
-         DO jtype = 1, numsurftypes
-            IF (ln_sla) THEN
-               obsfilessurf(:,jtype) = slafbfiles(:)
-               obstypessurf(jtype) = 'sla'
-            ENDIF
-            IF (ln_sss) THEN
-               obsfilessurf(:,jtype) = sssfbfiles(:)
-               obstypessurf(jtype) = 'sss'
-            ENDIF
-            IF (ln_sst) THEN
-               obsfilessurf(:,jtype) = sstfbfiles(:)
-               obstypessurf(jtype) = 'sst'
-            ENDIF
+      IF ( .NOT. ln_diaobs ) THEN
+         IF(lwp) WRITE(numout,cform_war)
+         IF(lwp) WRITE(numout,*)' ln_diaobs is set to false so not calling dia_obs'
+         RETURN
+      ENDIF
+
+      !-----------------------------------------------------------------------
+      ! Set up list of observation types to be used
+      ! and the files associated with each type
+      !-----------------------------------------------------------------------
+
+      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) )
+      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic /) )
+
+      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN
+         IF(lwp) WRITE(numout,cform_war)
+         IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', &
+            &                    ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', &
+            &                    ' are set to .FALSE. so turning off calls to dia_obs'
+         nwarn = nwarn + 1
+         ln_diaobs = .FALSE.
+         RETURN
+      ENDIF
+
+      IF ( nproftypes > 0 ) THEN
+
+         ALLOCATE( cobstypesprof(nproftypes) )
+         ALLOCATE( ifilesprof(nproftypes) )
+         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) )
+
+         jtype = 0
+         IF (ln_t3d .OR. ln_s3d) THEN
+            jtype = jtype + 1
+            clproffiles(jtype,:) = cn_profbfiles(:)
+            cobstypesprof(jtype) = 'prof  '
+            ifilesprof(jtype) = 0
+            DO jfile = 1, jpmaxnfiles
+               IF ( trim(clproffiles(jtype,jfile)) /= '' ) &
+                  ifilesprof(jtype) = ifilesprof(jtype) + 1
+            END DO
+         ENDIF
+         IF (ln_vel3d) THEN
+            jtype = jtype + 1
+            clproffiles(jtype,:) = cn_velfbfiles(:)
+            cobstypesprof(jtype) = 'vel   '
+            ifilesprof(jtype) = 0
+            DO jfile = 1, jpmaxnfiles
+               IF ( trim(clproffiles(jtype,jfile)) /= '' ) &
+                  ifilesprof(jtype) = ifilesprof(jtype) + 1
+            END DO
+         ENDIF
+
+      ENDIF
+
+      IF ( nsurftypes > 0 ) THEN
+
+         ALLOCATE( cobstypessurf(nsurftypes) )
+         ALLOCATE( ifilessurf(nsurftypes) )
+         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) )
+
+         jtype = 0
+         IF (ln_sla) THEN
+            jtype = jtype + 1
+            clsurffiles(jtype,:) = cn_slafbfiles(:)
+            cobstypessurf(jtype) = 'sla   '
+            ifilessurf(jtype) = 0
+            DO jfile = 1, jpmaxnfiles
+               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) &
+                  ifilessurf(jtype) = ifilessurf(jtype) + 1
+            END DO
+         ENDIF
+         IF (ln_sst) THEN
+            jtype = jtype + 1
+            clsurffiles(jtype,:) = cn_sstfbfiles(:)
+            cobstypessurf(jtype) = 'sst   '
+            ifilessurf(jtype) = 0
+            DO jfile = 1, jpmaxnfiles
+               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) &
+                  ifilessurf(jtype) = ifilessurf(jtype) + 1
+            END DO
+         ENDIF
 #if defined key_lim2 || defined key_lim3
-            IF (ln_seaice) THEN
-               obsfilessurf(:,jtype) = seaicefbfiles(:)
-               obstypessurf(jtype) = 'seaice'
-            ENDIF
+         IF (ln_sic) THEN
+            jtype = jtype + 1
+            clsurffiles(jtype,:) = cn_sicfbfiles(:)
+            cobstypessurf(jtype) = 'sic   '
+            ifilessurf(jtype) = 0
+            DO jfile = 1, jpmaxnfiles
+               IF ( trim(clsurffiles(jtype,jfile)) /= '' ) &
+                  ifilessurf(jtype) = ifilessurf(jtype) + 1
+            END DO
+         ENDIF
 #endif
 
-            lmask(:) = .FALSE.
-            WHERE (obsfilessurf(jtype,:) /= '') lmask(:) = .TRUE.
-            jnumfilessurf(jtype) = COUNT(lmask)
-         
-         END DO
-         
       ENDIF
 
@@ -259 +287 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '          Namelist namobs : set observation diagnostic parameters' 
-         WRITE(numout,*) '             Logical switch for T profile observations          ln_t3d = ', ln_t3d
-         WRITE(numout,*) '             Logical switch for S profile observations          ln_s3d = ', ln_s3d
-         WRITE(numout,*) '             Logical switch for SLA observations                ln_sla = ', ln_sla
-         WRITE(numout,*) '             Logical switch for SSH observations                ln_ssh = ', ln_ssh
-         WRITE(numout,*) '             Logical switch for SST observations                ln_sst = ', ln_sst
-         WRITE(numout,*) '             Logical switch for night-time SST obs         ln_sstnight = ', ln_sstnight
-         WRITE(numout,*) '             Logical switch for SSS observations                ln_sss = ', ln_sss
-         WRITE(numout,*) '             Logical switch for Sea Ice observations         ln_seaice = ', ln_seaice
-         WRITE(numout,*) '             Logical switch for velocity observations         ln_vel3d = ', ln_vel3d
-         WRITE(numout,*) '             Global distribution of observations        ln_grid_global = ',ln_grid_global
-         WRITE(numout,*) &
-   '             Logical switch for obs grid search w/lookup table  ln_grid_search_lookup = ',ln_grid_search_lookup
+         WRITE(numout,*) '             Logical switch for T profile observations                ln_t3d = ', ln_t3d
+         WRITE(numout,*) '             Logical switch for S profile observations                ln_s3d = ', ln_s3d
+         WRITE(numout,*) '             Logical switch for SLA observations                      ln_sla = ', ln_sla
+         WRITE(numout,*) '             Logical switch for SST observations                      ln_sst = ', ln_sst
+         WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic
+         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d
+         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ',ln_grid_global
+         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ',ln_grid_search_lookup
          IF (ln_grid_search_lookup) &
-            WRITE(numout,*) '             Grid search lookup file header       grid_search_file = ', grid_search_file
-         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS        dobsini = ', dobsin
-         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS          dobsend = ', dobsend
-         WRITE(numout,*) '             Type of vertical interpolation method          n1dint = ', n1dint
-         WRITE(numout,*) '             Type of horizontal interpolation method        n2dint = ', n2dint
-         WRITE(numout,*) '             Rejection of observations near land swithch    ln_nea = ', ln_nea
-         WRITE(numout,*) '             MSSH correction scheme                         nmsshc = ', nmsshc
-         WRITE(numout,*) '             MDT  correction                               mdtcorr = ', mdtcorr
-         WRITE(numout,*) '             MDT cutoff for computed correction          mdtcutoff = ', mdtcutoff
-         WRITE(numout,*) '             Logical switch for alt bias                ln_altbias = ', ln_altbias
-         WRITE(numout,*) '             Logical switch for ignoring missing files   ln_ignmis = ', ln_ignmis
-         WRITE(numout,*) '             Daily average types                                   = ', dailyavtypes
-
-         IF ( numproftypes > 0 ) THEN
-            DO jtype = 1, numproftypes
-               DO ji = 1, jnumfilesprof(jtype)
-                  WRITE(numout,'(1X,2A)') '             '//obstypesprof(jtype)//' input observation file names  = ', &
-                     TRIM(obsfilesprof(jtype,ji))
-                  IF ( TRIM(obstypesprof(jtype)) == 'prof' ) &
-                     WRITE(numout,'(1X,2A)') '       Enact feedback input time setting switch    ln_profb_enatim = ', ln_profb_enatim(ji)
+            WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile
+         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS               rn_dobsini = ', rn_dobsini
+         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend
+         WRITE(numout,*) '             Type of vertical interpolation method                  nn_1dint = ', nn_1dint
+         WRITE(numout,*) '             Type of horizontal interpolation method                nn_2dint = ', nn_2dint
+         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea
+         WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc
+         WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr
+         WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff
+         WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias
+         WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis
+         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes
+         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight
+         WRITE(numout,*) '          Number of profile obs types: ',nproftypes
+
+         IF ( nproftypes > 0 ) THEN
+            DO jtype = 1, nproftypes
+               DO jfile = 1, ifilesprof(jtype)
+                  WRITE(numout,'(1X,2A)') '             '//cobstypesprof(jtype)//' input observation file names  = ', &
+                     TRIM(clproffiles(jtype,jfile))
                END DO
             END DO
          ENDIF
-         
-         IF ( numsurftypes > 0 ) THEN
-            DO jtype = 1, numsurftypes
-               DO ji = 1, jnumfilessurf(jtype)
-                  WRITE(numout,'(1X,2A)') '             '//obstypessurf(jtype)//' input observation file names  = ', &
-                     TRIM(obsfilessurf(jtype,ji))
+
+         WRITE(numout,*)'          Number of surface obs types: ',nsurftypes
+         IF ( nsurftypes > 0 ) THEN
+            DO jtype = 1, nsurftypes
+               DO jfile = 1, ifilessurf(jtype)
+                  WRITE(numout,'(1X,2A)') '             '//cobstypessurf(jtype)//' input observation file names  = ', &
+                     TRIM(clsurffiles(jtype,jfile))
                END DO
             END DO
          ENDIF
-
-      ENDIF
-      
+         WRITE(numout,*) '~~~~~~~~~~~~'
+
+      ENDIF
+
+      !-----------------------------------------------------------------------
+      ! Obs operator parameter checking and initialisations
+      !-----------------------------------------------------------------------
+
       IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
          CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
@@ -312 +342 @@
       ENDIF
 
-      CALL obs_typ_init
-      
-      CALL mppmap_init
-      
-      ! Parameter control
-#if defined key_diaobs
-      IF ( numobtypes == 0 ) THEN
-         IF(lwp) WRITE(numout,cform_war)
-         IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', &
-            &                    ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d are all set to .FALSE.'
-         nwarn = nwarn + 1
-      ENDIF
-#endif
-
-      CALL obs_grid_setup( )
-      IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN
+      IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN
          CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
             &                    ' is not available')
       ENDIF
-      IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN
+
+      IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 4 ) ) THEN
          CALL ctl_stop(' Choice of horizontal (2D) interpolation method', &
             &                    ' is not available')
       ENDIF
 
+      CALL obs_typ_init
+
+      CALL mppmap_init
+
+      CALL obs_grid_setup( )
+
       !-----------------------------------------------------------------------
       ! Depending on switches read the various observation types
       !-----------------------------------------------------------------------
-      
-      IF ( numproftypes > 0 ) THEN
-      
-         ALLOCATE(profdata(numproftypes))
-         ALLOCATE(profdataqc(numproftypes))
-         ALLOCATE(nvarsprof(numproftypes))
-         ALLOCATE(nextrprof(numproftypes))
-            
-         DO jtype = 1, numproftypes
-      
+
+      IF ( nproftypes > 0 ) THEN
+
+         ALLOCATE(profdata(nproftypes))
+         ALLOCATE(profdataqc(nproftypes))
+         ALLOCATE(nvarsprof(nproftypes))
+         ALLOCATE(nextrprof(nproftypes))
+
+         DO jtype = 1, nproftypes
+
             nvarsprof(jtype) = 2
-            IF ( TRIM(obstypesprof(jtype)) == 'prof' ) nextrprof(jtype) = 1
-            IF ( TRIM(obstypesprof(jtype)) == 'vel' )  nextrprof(jtype) = 2
-
-            !Read in profile or velocity obs types
-            CALL obs_rea_prof( profdata(jtype),          &
-               &               jnumfilesprof(jtype),       &
-               &               obsfilesprof(jtype,1:jnumfilesprof(jtype)), &
-               &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2,             &
-               &               dobsini, dobsend, ln_t3d, ln_s3d, &
-               &               ln_ignmis, ln_s_at_t, .TRUE., .FALSE., &
-               &               kdailyavtypes = dailyavtypes )
-            
-            DO jvar = 1, nvars
+            IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN
+               nextrprof(jtype) = 1
+               llvar1 = ln_t3d
+               llvar2 = ln_s3d
+               zglam1 = glamt
+               zgphi1 = gphit
+               zmask1 = tmask
+               zglam2 = glamt
+               zgphi2 = gphit
+               zmask2 = tmask
+            ENDIF
+            IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN
+               nextrprof(jtype) = 2
+               llvar1 = ln_vel3d
+               llvar2 = ln_vel3d
+               zglam1 = glamu
+               zgphi1 = gphiu
+               zmask1 = umask
+               zglam2 = glamv
+               zgphi2 = gphiv
+               zmask2 = vmask
+            ENDIF
+
+            !Read in profile or profile obs types
+            CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype),       &
+               &               clproffiles(jtype,1:ifilesprof(jtype)), &
+               &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, &
+               &               rn_dobsini, rn_dobsend, llvar1, llvar2, &
+               &               ln_ignmis, ln_s_at_t, .FALSE., &
+               &               kdailyavtypes = nn_profdavtypes )
+
+            DO jvar = 1, nvarsprof(jtype)
                CALL obs_prof_staend( profdata(jtype), jvar )
             END DO
-         
-            CALL obs_pre_prof( profdata(jtype), profdataqc(jtype),   &
-               &              ln_t3d, ln_s3d, ln_nea, &
-               &              kdailyavtypes = dailyavtypes )
+
+            CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), &
+               &               llvar1, llvar2, &
+               &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2,  &
+               &               ln_nea, kdailyavtypes = nn_profdavtypes )
 
          END DO
 
-         DEALLOCATE( jnumfilesprof, obsfilesprof )
-
-      ENDIF
-
-      IF ( numsurftypes > 0 ) THEN
-      
-         ALLOCATE(surfdata(numsurftypes))
-         ALLOCATE(surfdatatqc(numsurftypes))
-         ALLOCATE(nvarssurf(numsurftypes))
-         ALLOCATE(nextrsurf(numsurftypes))
-         
-         DO jtype = 1, numsurftypes
-      
+         DEALLOCATE( ifilesprof, clproffiles )
+
+      ENDIF
+
+      IF ( nsurftypes > 0 ) THEN
+
+         ALLOCATE(surfdata(nsurftypes))
+         ALLOCATE(surfdataqc(nsurftypes))
+         ALLOCATE(nvarssurf(nsurftypes))
+         ALLOCATE(nextrsurf(nsurftypes))
+
+         DO jtype = 1, nsurftypes
+
             nvarssurf(jtype) = 1
             nextrsurf(jtype) = 0
-            IF ( TRIM(obstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2
+            llnightav = .FALSE.
+            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2
+            IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav = ln_sstnight
 
             !Read in surface obs types
-            CALL obs_rea_surf( surfdata(jtype), jnumfilessurf(jtype), &
-               &               obsfilessurf(jtype,1:jnumfilessurf(jtype)), &
+            CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), &
+               &               clsurffiles(jtype,1:ifilessurf(jtype)), &
                &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, &
-               &               dobsini, dobsend, ln_ignmis, .FALSE. )
+               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav )
 
             CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea )
-         
-            IF ( TRIM(obstypessurf(jtype)) == 'sla' ) THEN
-               CALL obs_rea_mdt( surfdataqc(jtype), n2dint )
-               IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), n2dint, bias_file )
+
+            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
+               CALL obs_rea_mdt( surfdataqc(jtype), nn_2dint )
+               IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), nn_2dint, cn_altbiasfile )
             ENDIF
 
-            DEALLOCATE( jnumfilessurf, obsfilessurf )
-
          END DO
+
+         DEALLOCATE( ifilessurf, clsurffiles )
+
+      ENDIF
 
    END SUBROUTINE dia_obs_init
@@ -415 +461 @@
       !!
       !! ** Method  : Call the observation operators on each time step to
-      !!              compute the model equivalent of the following date:
-      !!               - T profiles
-      !!               - S profiles
-      !!               - Sea surface height (referenced to a mean)
-      !!               - Sea surface temperature
-      !!               - Sea surface salinity
-      !!               - Velocity component (U,V) profiles
-      !!
-      !! ** Action  : 
+      !!              compute the model equivalent of the following data:
+      !!               - Profile data, currently T/S or U/V
+      !!               - Surface data, currently SST, SLA or sea-ice concentration.
+      !!
+      !! ** Action  :
       !!
       !! History :
@@ -432 +474 @@
       !!        !  07-04  (G. Smith) Generalized surface operators
       !!        !  08-10  (M. Valdivieso) obs operator for velocity profiles
+      !!        !  15-08  (M. Martin) Combined surface/profile routines.
       !!----------------------------------------------------------------------
       !! * Modules used
-      USE dom_oce, ONLY : &             ! Ocean space and time domain variables
-         & rdt,           &                       
-         & gdept_1d,       &             
-         & tmask, umask, vmask                            
-      USE phycst, ONLY : &              ! Physical constants
-         & rday                         
-      USE oce, ONLY : &                 ! Ocean dynamics and tracers variables
-         & tsn,  &             
-         & un, vn,  &
+      USE phycst, ONLY : &         ! Physical constants
+         & rday
+      USE oce, ONLY : &            ! Ocean dynamics and tracers variables
+         & tsn,       &
+         & un,        &
+         & vn,        &
          & sshn
 #if defined  key_lim3
-      USE ice, ONLY : &                     ! LIM Ice model variables
+      USE ice, ONLY : &            ! LIM3 Ice model variables
          & frld
 #endif
 #if defined key_lim2
-      USE ice_2, ONLY : &                     ! LIM Ice model variables
+      USE ice_2, ONLY : &          ! LIM2 Ice model variables
          & frld
 #endif
@@ -455 +495 @@
 
       !! * Arguments
-      INTEGER, INTENT(IN) :: kstp                         ! Current timestep
+      INTEGER, INTENT(IN) :: kstp  ! Current timestep
       !! * Local declarations
-      INTEGER :: idaystp                ! Number of timesteps per day
-      INTEGER :: jtype                  ! data loop variable
-      INTEGER :: jvar                   ! Variable number    
+      INTEGER :: idaystp           ! Number of timesteps per day
+      INTEGER :: jtype             ! Data loop variable
+      INTEGER :: jvar              ! Variable number
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
+         & zprofvar1, &            ! Model values for 1st variable in a prof ob
+         & zprofvar2               ! Model values for 2nd variable in a prof ob
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
+         & zprofmask1, &           ! Mask associated with zprofvar1
+         & zprofmask2              ! Mask associated with zprofvar2
+      REAL(wp), DIMENSION(jpi,jpj)     :: &
+         & zsurfvar                ! Model values equivalent to surface ob.
+      REAL(wp), DIMENSION(jpi,jpj) :: &
+         & zglam1,    &            ! Model longitudes for prof variable 1
+         & zglam2,    &            ! Model longitudes for prof variable 2
+         & zgphi1,    &            ! Model latitudes for prof variable 1
+         & zgphi2                  ! Model latitudes for prof variable 2
 #if ! defined key_lim2 && ! defined key_lim3
-      REAL(wp), POINTER, DIMENSION(:,:) :: frld   
+      REAL(wp), POINTER, DIMENSION(:,:) :: frld
 #endif
-      CHARACTER(LEN=20) :: datestr=" ",timestr=" "
- 
+      LOGICAL :: llnightav        ! Logical for calculating night-time average
+
 #if ! defined key_lim2 && ! defined key_lim3
       CALL wrk_alloc(jpi,jpj,frld) 
@@ -473 +526 @@
          WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
          WRITE(numout,*) '~~~~~~~'
+         CALL FLUSH(numout)
       ENDIF
 
@@ -484 +538 @@
 #endif
       !-----------------------------------------------------------------------
-      ! Depending on switches call various observation operators
-      !-----------------------------------------------------------------------
-
-      IF ( numproftypes > 0 ) THEN
-         DO jtype = 1, numproftypes
-      
-            SELECT CASE ( TRIM(obstypesprof(jtype)) )
+      ! Call the profile and surface observation operators
+      !-----------------------------------------------------------------------
+
+      IF ( nproftypes > 0 ) THEN
+
+         DO jtype = 1, nproftypes
+
+            SELECT CASE ( TRIM(cobstypesprof(jtype)) )
             CASE('prof')
-               profvar1(:,:,:) = tsn(:,:,:,jp_tem)
-               profvar2(:,:,:) = tsn(:,:,:,jp_sal)
-               profmask1(:,:,:) = tmask(:,:,:)
-               profmask2(:,:,:) = tmask(:,:,:)
+               zprofvar1(:,:,:) = tsn(:,:,:,jp_tem)
+               zprofvar2(:,:,:) = tsn(:,:,:,jp_sal)
+               zprofmask1(:,:,:) = tmask(:,:,:)
+               zprofmask2(:,:,:) = tmask(:,:,:)
+               zglam1(:,:) = glamt(:,:)
+               zglam2(:,:) = glamt(:,:)
+               zgphi1(:,:) = gphit(:,:)
+               zgphi2(:,:) = gphit(:,:)
             CASE('vel')
-               profvar1(:,:,:) = un(:,:,:)
-               profvar2(:,:,:) = vn(:,:,:)
-               profmask1(:,:,:) = umask(:,:,:)
-               profmask2(:,:,:) = vmask(:,:,:)
+               zprofvar1(:,:,:) = un(:,:,:)
+               zprofvar2(:,:,:) = vn(:,:,:)
+               zprofmask1(:,:,:) = umask(:,:,:)
+               zprofmask2(:,:,:) = vmask(:,:,:)
+               zglam1(:,:) = glamu(:,:)
+               zglam2(:,:) = glamv(:,:)
+               zgphi1(:,:) = gphiu(:,:)
+               zgphi2(:,:) = gphiv(:,:)
             END SELECT
-            
-            CALL obs_prof_opt( profdataqc(jtype),                       &
-               &               kstp, jpi, jpj, jpk, nit000, idaystp,   &
-               &               profvar1, profvar2,   &
-               &               gdept_1d, profmask1, profmask2, n1dint, n2dint,        &
-               &               kdailyavtypes = dailyavtypes )
-            
+
+            CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
+               &               nit000, idaystp,                         &
+               &               zprofvar1, zprofvar2,                    &
+               &               gdept_1d, zprofmask1, zprofmask2,        &
+               &               zglam1, zglam2, zgphi1, zgphi2,          &
+               &               nn_1dint, nn_2dint,                      &
+               &               kdailyavtypes = nn_profdavtypes )
+
          END DO
 
       ENDIF
 
-      IF ( numsurftypes > 0 ) THEN
-         DO jtype = 1, numsurftypes
-      
-            SELECT CASE ( TRIM(obstypessurf(jtype)) )
+      IF ( nsurftypes > 0 ) THEN
+
+         DO jtype = 1, nsurftypes
+
+            SELECT CASE ( TRIM(cobstypessurf(jtype)) )
             CASE('sst')
-               surfvar(:,:) = tsn(:,:,1,jp_tem)
+               zsurfvar(:,:) = tsn(:,:,1,jp_tem)
+               llnightav = ln_sstnight
             CASE('sla')
-               surfvar(:,:) = sshn(:,:)
-            CASE('sss')
-               surfvar(:,:) = tsn(:,:,1,jp_sal)
+               zsurfvar(:,:) = sshn(:,:)
+               llnightav = .FALSE.
 #if defined key_lim2 || defined key_lim3
-            CASE('seaice')
-               surfvar(:,:) = 1._wp - frld(:,:)
+            CASE('sic')
+               zsurfvar(:,:) = 1._wp - frld(:,:)
+               llnightav = .FALSE.
 #endif
             END SELECT
-         
-            CALL obs_surf_opt( surfdatqc(jtype),             &
-               &               kstp, jpi, jpj, nit000, surfvar, &
-               &               tmask(:,:,1), n2dint, ld_sstnight )
-            
+
+            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &
+               &               nit000, idaystp, zsurfvar, tmask(:,:,1), &
+               &               nn_2dint, llnightav )
+
          END DO
-         
-      ENDIF
-      
+
+      ENDIF
+
 #if ! defined key_lim2 && ! defined key_lim3
       CALL wrk_dealloc(jpi,jpj,frld) 
@@ -542 +609 @@
 
    END SUBROUTINE dia_obs
-  
+
    SUBROUTINE dia_obs_wri 
       !!----------------------------------------------------------------------
@@ -559 +626 @@
       !!        !  07-03  (K. Mogensen) General handling of profiles
       !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles
+      !!        !  15-08  (M. Martin) Combined writing for prof and surf types
       !!----------------------------------------------------------------------
       IMPLICIT NONE
 
       !! * Local declarations
-
       INTEGER :: jtype                    ! Data set loop variable
+
       !-----------------------------------------------------------------------
       ! Depending on switches call various observation output routines
       !-----------------------------------------------------------------------
 
-      IF ( numproftypes > 0 ) THEN
-         DO jtype = 1, numproftypes
-      
+      IF ( nproftypes > 0 ) THEN
+
+         DO jtype = 1, nproftypes
+
             CALL obs_prof_decompress( profdataqc(jtype), &
                &                      profdata(jtype), .TRUE., numout )
 
-            CALL obs_wri_prof( obstypesprof(jtype), profdata(jtype), n2dint )
-         
+            CALL obs_wri_prof( profdata(jtype), nn_2dint )
+
          END DO
-         
-      ENDIF
-
-      IF ( numsurftypes > 0 ) THEN
-         DO jtype = 1, numsurftypes
-      
-            CALL obs_surf_decompress( surfdatqc(jtype), &
+
+      ENDIF
+
+      IF ( nsurftypes > 0 ) THEN
+
+         DO jtype = 1, nsurftypes
+
+            CALL obs_surf_decompress( surfdataqc(jtype), &
                &                      surfdata(jtype), .TRUE., numout )
 
-            CALL obs_wri_surf( obstypessurf(jtype), surfdata(jtype), n2dint )
+            CALL obs_wri_surf( surfdata(jtype) )
 
          END DO
-         
-      ENDIF
-
+
+      ENDIF
 
    END SUBROUTINE dia_obs_wri
@@ -609 +678 @@
       !!
       !!----------------------------------------------------------------------
-      !! obs_grid deallocation
+      ! obs_grid deallocation
       CALL obs_grid_deallocate
 
-      !! diaobs deallocation
-      IF ( numproftypes > 0 ) DEALLOCATE(profdata, profdataqc, nvarsprof, nextrprof)
-      IF ( numsurftypes > 0 ) DEALLOCATE(surfdata, surfdataqc, nvarssurf, nextrsurf)
-      
+      ! diaobs deallocation
+      IF ( nproftypes > 0 ) &
+         &   DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof )
+
+      IF ( nsurftypes > 0 ) &
+         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf )
+
    END SUBROUTINE dia_obs_dealloc
 
@@ -621 +693 @@
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ini_date  ***
-      !!          
-      !! ** Purpose : Get initial data in double precision YYYYMMDD.HHMMSS format
-      !!
-      !! ** Method  : Get initial data in double precision YYYYMMDD.HHMMSS format
-      !!
-      !! ** Action  : Get initial data in double precision YYYYMMDD.HHMMSS format
+      !!
+      !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format
+      !!
+      !! ** Method  : Get initial date in double precision YYYYMMDD.HHMMSS format
+      !!
+      !! ** Action  : Get initial date in double precision YYYYMMDD.HHMMSS format
       !!
       !! History :
@@ -637 +709 @@
       USE phycst, ONLY : &            ! Physical constants
          & rday
-!      USE daymod, ONLY : &            ! Time variables
-!         & nmonth_len           
       USE dom_oce, ONLY : &           ! Ocean space and time domain variables
          & rdt
@@ -645 +715 @@
 
       !! * Arguments
-      REAL(KIND=dp), INTENT(OUT) :: ddobsini                         ! Initial date in YYYYMMDD.HHMMSS
+      REAL(dp), INTENT(OUT) :: ddobsini  ! Initial date in YYYYMMDD.HHMMSS
 
       !! * Local declarations
@@ -653 +723 @@
       INTEGER :: ihou
       INTEGER :: imin
-      INTEGER :: imday         ! Number of days in month.
-      REAL(KIND=wp) :: zdayfrc ! Fraction of day
-
-      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
-
-      !!----------------------------------------------------------------------
-      !! Initial date initialization (year, month, day, hour, minute)
-      !! (This assumes that the initial date is for 00z))
-      !!----------------------------------------------------------------------
+      INTEGER :: imday       ! Number of days in month.
+      INTEGER, DIMENSION(12) :: &
+         &       imonth_len  ! Length in days of the months of the current year
+      REAL(wp) :: zdayfrc    ! Fraction of day
+
+      !----------------------------------------------------------------------
+      ! Initial date initialization (year, month, day, hour, minute)
+      ! (This assumes that the initial date is for 00z))
+      !----------------------------------------------------------------------
       iyea =   ndate0 / 10000
       imon = ( ndate0 - iyea * 10000 ) / 100
@@ -668 +738 @@
       imin = 0
 
-      !!----------------------------------------------------------------------
-      !! Compute number of days + number of hours + min since initial time
-      !!----------------------------------------------------------------------
+      !----------------------------------------------------------------------
+      ! Compute number of days + number of hours + min since initial time
+      !----------------------------------------------------------------------
       iday = iday + ( nit000 -1 ) * rdt / rday
       zdayfrc = ( nit000 -1 ) * rdt / rday
@@ -677 +747 @@
       imin = int( (zdayfrc * 24 - ihou) * 60 )
 
-      !!-----------------------------------------------------------------------
-      !! Convert number of days (iday) into a real date
-      !!----------------------------------------------------------------------
+      !-----------------------------------------------------------------------
+      ! Convert number of days (iday) into a real date
+      !----------------------------------------------------------------------
 
       CALL calc_month_len( iyea, imonth_len )
-      
+
       DO WHILE ( iday > imonth_len(imon) )
          iday = iday - imonth_len(imon)
@@ -693 +763 @@
       END DO
 
-      !!----------------------------------------------------------------------
-      !! Convert it into YYYYMMDD.HHMMSS format.
-      !!----------------------------------------------------------------------
+      !----------------------------------------------------------------------
+      ! Convert it into YYYYMMDD.HHMMSS format.
+      !----------------------------------------------------------------------
       ddobsini = iyea * 10000_dp + imon * 100_dp + &
          &       iday + ihou * 0.01_dp + imin * 0.0001_dp
@@ -705 +775 @@
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE fin_date  ***
-      !!          
-      !! ** Purpose : Get final data in double precision YYYYMMDD.HHMMSS format
-      !!
-      !! ** Method  : Get final data in double precision YYYYMMDD.HHMMSS format
-      !!
-      !! ** Action  : Get final data in double precision YYYYMMDD.HHMMSS format
+      !!
+      !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format
+      !!
+      !! ** Method  : Get final date in double precision YYYYMMDD.HHMMSS format
+      !!
+      !! ** Action  : Get final date in double precision YYYYMMDD.HHMMSS format
       !!
       !! History :
@@ -720 +790 @@
       USE phycst, ONLY : &            ! Physical constants
          & rday
-!      USE daymod, ONLY : &            ! Time variables
-!         & nmonth_len                
       USE dom_oce, ONLY : &           ! Ocean space and time domain variables
          & rdt
@@ -728 +796 @@
 
       !! * Arguments
-      REAL(KIND=dp), INTENT(OUT) :: ddobsfin                   ! Final date in YYYYMMDD.HHMMSS
+      REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
 
       !! * Local declarations
@@ -736 +804 @@
       INTEGER :: ihou
       INTEGER :: imin
-      INTEGER :: imday         ! Number of days in month.
-      REAL(KIND=wp) :: zdayfrc       ! Fraction of day
-         
-      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
-            
+      INTEGER :: imday       ! Number of days in month.
+      INTEGER, DIMENSION(12) :: &
+         &       imonth_len  ! Length in days of the months of the current year
+      REAL(wp) :: zdayfrc    ! Fraction of day
+
       !-----------------------------------------------------------------------
       ! Initial date initialization (year, month, day, hour, minute)
@@ -750 +818 @@
       ihou = 0
       imin = 0
-      
+
       !-----------------------------------------------------------------------
       ! Compute number of days + number of hours + min since initial time
@@ -765 +833 @@
 
       CALL calc_month_len( iyea, imonth_len )
-      
+
       DO WHILE ( iday > imonth_len(imon) )
          iday = iday - imonth_len(imon)
@@ -783 +851 @@
 
     END SUBROUTINE fin_date
-    
+
 END MODULE diaobs

branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_grid.F90

@@ -52 +52 @@
 
    !! Default values
-   REAL, PUBLIC :: grid_search_res = 0.5    ! Resolution of grid
+   REAL, PUBLIC :: rn_gridsearchres = 0.5   ! Resolution of grid
    INTEGER, PRIVATE :: gsearch_nlons_def    ! Num of longitudes
    INTEGER, PRIVATE :: gsearch_nlats_def    ! Num of latitudes
@@ -83 +83 @@
    LOGICAL, PUBLIC :: ln_grid_global         ! Use global distribution of observations
    CHARACTER(LEN=44), PUBLIC :: &
-      & grid_search_file    ! file name head for grid search lookup 
+      & cn_gridsearchfile    ! file name head for grid search lookup 
 
    !!----------------------------------------------------------------------
@@ -690 +690 @@
          
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*)'Grid search resolution : ', grid_search_res
-         
-         gsearch_nlons_def  = NINT( 360.0_wp / grid_search_res ) 
-         gsearch_nlats_def  = NINT( 180.0_wp / grid_search_res )
-         gsearch_lonmin_def = -180.0_wp + 0.5_wp * grid_search_res
-         gsearch_latmin_def =  -90.0_wp + 0.5_wp * grid_search_res
-         gsearch_dlon_def   = grid_search_res
-         gsearch_dlat_def   = grid_search_res
+         IF(lwp) WRITE(numout,*)'Grid search resolution : ', rn_gridsearchres
+         
+         gsearch_nlons_def  = NINT( 360.0_wp / rn_gridsearchres ) 
+         gsearch_nlats_def  = NINT( 180.0_wp / rn_gridsearchres )
+         gsearch_lonmin_def = -180.0_wp + 0.5_wp * rn_gridsearchres
+         gsearch_latmin_def =  -90.0_wp + 0.5_wp * rn_gridsearchres
+         gsearch_dlon_def   = rn_gridsearchres
+         gsearch_dlat_def   = rn_gridsearchres
          
          IF (lwp) THEN
@@ -710 +710 @@
          IF ( ln_grid_global ) THEN
             WRITE(cfname, FMT="(A,'_',A)") &
-               &          TRIM(grid_search_file), 'global.nc'
+               &          TRIM(cn_gridsearchfile), 'global.nc'
          ELSE
             WRITE(cfname, FMT="(A,'_',I4.4,'of',I4.4,'by',I4.4,'.nc')") &
-               &          TRIM(grid_search_file), nproc, jpni, jpnj
+               &          TRIM(cn_gridsearchfile), nproc, jpni, jpnj
          ENDIF
 

branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_mpp.F90

@@ -115 +115 @@
       !!               ***  ROUTINE obs_mpp_find_obs_proc ***
       !!          
-      !! ** Purpose : From the array kobsp containing the results of the grid
+      !! ** Purpose : From the array kobsp containing the results of the
       !!              grid search on each processor the processor return a
       !!              decision of which processors should hold the observation.

branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

@@ -11 +11 @@
    !!----------------------------------------------------------------------
 
-   !! * Modules used   
+   !! * Modules used
    USE par_kind, ONLY : &         ! Precision variables
       & wp
    USE in_out_manager             ! I/O manager
    USE obs_inter_sup              ! Interpolation support
-   USE obs_inter_h2d, ONLY : &    ! Horizontal interpolation to the observation pt
+   USE obs_inter_h2d, ONLY : &    ! Horizontal interpolation to the obs pt
       & obs_int_h2d, &
       & obs_int_h2d_init
-   USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the observation pt
+   USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the obs pt
       & obs_int_z1d,    &
       & obs_int_z1d_spl
-   USE obs_const,  ONLY :     &
-      & obfillflt      ! Fillvalue   
-   USE dom_oce,       ONLY : &
-      & glamt, glamu, glamv, &
-      & gphit, gphiu, gphiv
-   USE lib_mpp,       ONLY : &
+   USE obs_const,  ONLY :    &    ! Obs fill value
+      & obfillflt
+   USE dom_oce,       ONLY : &    ! Model lats/lons
+      & glamt, &
+      & gphit
+   USE lib_mpp,       ONLY : &    ! Warning and stopping routines
       & ctl_warn, ctl_stop
+   USE sbcdcy,        ONLY : &    ! For calculation of where it is night-time
+      & sbc_dcy, nday_qsr
 
    IMPLICIT NONE
@@ -35 +37 @@
    PRIVATE
 
-   PUBLIC obs_prof_opt, &  ! Compute the model counterpart of profile observations
-      &   obs_surf_opt     ! Compute the model counterpart of surface observations
-
-   INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types
+   PUBLIC obs_prof_opt, &  ! Compute the model counterpart of profile obs
+      &   obs_surf_opt     ! Compute the model counterpart of surface obs
+
+   INTEGER, PARAMETER, PUBLIC :: &
+      & imaxavtypes = 20   ! Max number of daily avgd obs types
 
    !!----------------------------------------------------------------------
@@ -48 +51 @@
 CONTAINS
 
-   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
-      &                     pvar1, pvar2, pgdept, pmask1, k1dint, k2dint, &
-      &                     kdailyavtypes )
+   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk,          &
+      &                     kit000, kdaystp,                      &
+      &                     pvar1, pvar2, pgdept, pmask1, pmask2, &
+      &                     plam1, plam2, pphi1, pphi2,           &
+      &                     k1dint, k2dint, kdailyavtypes )
+
       !!-----------------------------------------------------------------------
       !!
@@ -99 +105 @@
       !!      ! 15-02 (M. Martin) Combined routine for all profile types
       !!-----------------------------------------------------------------------
-  
+
       !! * Modules used
       USE obs_profiles_def ! Definition of storage space for profile obs.
@@ -106 +112 @@
 
       !! * Arguments
-      TYPE(obs_prof), INTENT(INOUT) :: prodatqc  ! Subset of profile data not failing screening
-      INTEGER, INTENT(IN) :: kt        ! Time step
-      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
+      TYPE(obs_prof), INTENT(INOUT) :: &
+         & prodatqc                  ! Subset of profile data passing QC
+      INTEGER, INTENT(IN) :: kt      ! Time step
+      INTEGER, INTENT(IN) :: kpi     ! Model grid parameters
       INTEGER, INTENT(IN) :: kpj
       INTEGER, INTENT(IN) :: kpk
-      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step 
-                                       !   (kit000-1 = restart time)
-      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)
-      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
-      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                    
+      INTEGER, INTENT(IN) :: kit000  ! Number of the first time step
+                                     !   (kit000-1 = restart time)
+      INTEGER, INTENT(IN) :: k1dint  ! Vertical interpolation type (see header)
+      INTEGER, INTENT(IN) :: k2dint  ! Horizontal interpolation type (see header)
+      INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
-         & pvar1,    &    ! Model field 1
-         & pvar2,    &    ! Model field 2
-         & pmask1,   &    ! Land-sea mask 1
-         & pmask2         ! Land-sea mask 2
+         & pvar1,    &               ! Model field 1
+         & pvar2,    &               ! Model field 2
+         & pmask1,   &               ! Land-sea mask 1
+         & pmask2                    ! Land-sea mask 2
+      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
+         & plam1,    &               ! Model longitudes for variable 1
+         & plam2,    &               ! Model longitudes for variable 2
+         & pphi1,    &               ! Model latitudes for variable 1
+         & pphi2                     ! Model latitudes for variable 2
       REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
-         & pgdept       ! Model array of depth levels
+         & pgdept                    ! Model array of depth levels
       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
-         & kdailyavtypes! Types for daily averages
+         & kdailyavtypes             ! Types for daily averages
+
       !! * Local declarations
       INTEGER ::   ji
@@ -138 +151 @@
       INTEGER, DIMENSION(imaxavtypes) :: &
          & idailyavtypes
+      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
+         & igrdi1, &
+         & igrdi2, &
+         & igrdj1, &
+         & igrdj2
       REAL(KIND=wp) :: zlam
       REAL(KIND=wp) :: zphi
@@ -144 +162 @@
          & zobsmask1, &
          & zobsmask2, &
-         & zobsmask, &
          & zobsk,    &
          & zobs2k
@@ -162 +179 @@
          & zgphi1, &
          & zgphi2
-      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
-         & igrdi1, &
-         & igrdi2, &
-         & igrdj1, &
-         & igrdj2
+      LOGICAL :: ld_dailyav
 
       !------------------------------------------------------------------------
       ! Local initialization 
       !------------------------------------------------------------------------
-      ! ... Record and data counters
+      ! Record and data counters
       inrc = kt - kit000 + 2
       ipro = prodatqc%npstp(inrc)
- 
+
       ! Daily average types
+      ld_dailyav = .FALSE.
       IF ( PRESENT(kdailyavtypes) ) THEN
          idailyavtypes(:) = kdailyavtypes(:)
+         IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE.
       ELSE
          idailyavtypes(:) = -1
       ENDIF
 
-      ! Initialize daily mean for first timestep
+      ! Daily means are calculated for values over timesteps:
+      !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ...
       idayend = MOD( kt - kit000 + 1, kdaystp )
 
-      ! Added kt == 0 test to catch restart case 
-      IF ( idayend == 1 .OR. kt == 0) THEN
-         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
+      IF ( ld_dailyav ) THEN
+
+         ! Initialize daily mean for first timestep of the day
+         IF ( idayend == 1 .OR. kt == 0 ) THEN
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     prodatqc%vdmean(ji,jj,jk,1) = 0.0
+                     prodatqc%vdmean(ji,jj,jk,2) = 0.0
+                  END DO
+               END DO
+            END DO
+         ENDIF
+
          DO jk = 1, jpk
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  prodatqc%vdmean(ji,jj,jk,1) = 0.0
-                  prodatqc%vdmean(ji,jj,jk,2) = 0.0
+                  ! Increment field 1 for computing daily mean
+                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
+                     &                        + pvar1(ji,jj,jk)
+                  ! Increment field 2 for computing daily mean
+                  prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
+                     &                        + pvar2(ji,jj,jk)
                END DO
             END DO
          END DO
-      ENDIF
-
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ! Increment field 1 for computing daily mean
-               prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
-                  &                        + pvar1(ji,jj,jk)
-               ! Increment field 2 for computing daily mean
-               prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
-                  &                        + pvar2(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-   
-      ! Compute the daily mean at the end of day
-      zdaystp = 1.0 / REAL( kdaystp )
-      IF ( idayend == 0 ) THEN
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
-                     &                        * zdaystp
-                  prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
-                  &                           * zdaystp
+
+         ! Compute the daily mean at the end of day
+         zdaystp = 1.0 / REAL( kdaystp )
+         IF ( idayend == 0 ) THEN
+            IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt
+            CALL FLUSH(numout)
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
+                        &                        * zdaystp
+                     prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
+                        &                        * zdaystp
+                  END DO
                END DO
             END DO
-         END DO
+         ENDIF
+
       ENDIF
 
@@ -262 +283 @@
       END DO
 
-      CALL obs_int_comm_2d( 2, 2, ipro, igrdi1, igrdj1, glamt1, zglam1 )
-      CALL obs_int_comm_2d( 2, 2, ipro, igrdi1, igrdj1, gphit1, zgphi1 )
+      CALL obs_int_comm_2d( 2, 2, ipro, igrdi1, igrdj1, plam1, zglam1 )
+      CALL obs_int_comm_2d( 2, 2, ipro, igrdi1, igrdj1, pphi1, zgphi1 )
       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi1, igrdj1, pmask1, zmask1 )
       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi1, igrdj1, pvar1,   zint1 )
       
-      CALL obs_int_comm_2d( 2, 2, ipro, igrdi2, igrdj2, glamt2, zglam2 )
-      CALL obs_int_comm_2d( 2, 2, ipro, igrdi2, igrdj2, gphit2, zgphi2 )
-      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, pmask1, zmask2 )
+      CALL obs_int_comm_2d( 2, 2, ipro, igrdi2, igrdj2, plam2, zglam2 )
+      CALL obs_int_comm_2d( 2, 2, ipro, igrdi2, igrdj2, pphi2, zgphi2 )
+      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, pmask2, zmask2 )
       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, pvar2,   zint2 )
 
       ! At the end of the day also get interpolated means
-      IF ( idayend == 0 ) THEN
+      IF ( ld_dailyav .AND. idayend == 0 ) THEN
 
          ALLOCATE( &
@@ -292 +313 @@
 
          IF ( kt /= prodatqc%mstp(jobs) ) THEN
-            
+
             IF(lwp) THEN
                WRITE(numout,*)
@@ -307 +328 @@
             CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
          ENDIF
-         
+
          zlam = prodatqc%rlam(jobs)
          zphi = prodatqc%rphi(jobs)
-         
+
          ! Horizontal weights and vertical mask
 
@@ -333 +354 @@
             zobsk(:) = obfillflt
 
-       IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
+            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
 
                IF ( idayend == 0 )  THEN
-                  
                   ! Daily averaged data
-                  
                   CALL obs_int_h2d( kpk, kpk,      &
                      &              zweig1, zinm1(:,:,:,iobs), zobsk )
-                  
-                  
-               ELSE
-               
-                  CALL ctl_stop( ' A nonzero' //     &
-                     &           ' number of average profile data should' // &
-                     &           ' only occur at the end of a given day' )
 
                ENDIF
-          
+
             ELSE 
-               
+
                ! Point data
-
                CALL obs_int_h2d( kpk, kpk,      &
                   &              zweig1, zint1(:,:,:,iobs), zobsk )
@@ -364 +375 @@
             ! polynomial at obs points
             !-------------------------------------------------------------
-            
+
             IF ( k1dint == 1 ) THEN
                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   &
-                  &                  pgdept, zobsmask )
+                  &                  pgdept, zobsmask1 )
             ENDIF
-            
+
             !-----------------------------------------------------------------
             !  Vertical interpolation to the observation point
@@ -381 +392 @@
                & zobsk, zobs2k,                   &
                & prodatqc%var(1)%vmod(ista:iend), &
-               & pgdept, zobsmask )
+               & pgdept, zobsmask1 )
 
          ENDIF
@@ -394 +405 @@
 
                   ! Daily averaged data
-                  
                   CALL obs_int_h2d( kpk, kpk,      &
                      &              zweig2, zinm2(:,:,:,iobs), zobsk )
-                  
-               ELSE
-
-                  CALL ctl_stop( ' A nonzero' //     &
-                     &           ' number of average profile data should' // &
-                     &           ' only occur at the end of a given day' )
 
                ENDIF
 
             ELSE
-               
+
                ! Point data
-
                CALL obs_int_h2d( kpk, kpk,      &
                   &              zweig2, zint2(:,:,:,iobs), zobsk )
@@ -420 +423 @@
             ! polynomial at obs points
             !-------------------------------------------------------------
-            
+
             IF ( k1dint == 1 ) THEN
                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
-                  &                  pgdept, zobsmask )
+                  &                  pgdept, zobsmask2 )
             ENDIF
-            
+
             !----------------------------------------------------------------
             !  Vertical interpolation to the observation point
@@ -437 +440 @@
                & zobsk, zobs2k, &
                & prodatqc%var(2)%vmod(ista:iend),&
-               & pgdept, zobsmask )
+               & pgdept, zobsmask2 )
 
          ENDIF
 
       END DO
- 
+
       ! Deallocate the data for interpolation
       DEALLOCATE( &
@@ -458 +461 @@
          & zint2   &
          & )
+
       ! At the end of the day also get interpolated means
-      IF ( idayend == 0 ) THEN
+      IF ( ld_dailyav .AND. idayend == 0 ) THEN
          DEALLOCATE( &
             & zinm1,  &
@@ -467 +471 @@
 
       prodatqc%nprofup = prodatqc%nprofup + ipro 
-      
+
    END SUBROUTINE obs_prof_opt
 
-   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, kit000, &
-      &                    psurf, psurfmask, k2dint, ld_nightav )
+   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,         &
+      &                    kit000, kdaystp, psurf, psurfmask, &
+      &                    k2dint, ldnightav )
+
       !!-----------------------------------------------------------------------
       !!
@@ -491 +497 @@
       !!        - bilinear (quadrilateral grid)    (k2dint = 3)
       !!        - polynomial (quadrilateral grid)  (k2dint = 4)
-      !!  
+      !!
       !!
       !! ** Action  :
@@ -499 +505 @@
       !!      ! 15-02 (M. Martin) Combined routine for surface types
       !!-----------------------------------------------------------------------
-  
+
       !! * Modules used
       USE obs_surf_def  ! Definition of storage space for surface observations
@@ -506 +512 @@
 
       !! * Arguments
-      TYPE(obs_surf), INTENT(INOUT) :: surfdataqc     ! Subset of surface data not failing screening
-      INTEGER, INTENT(IN) :: kt      ! Time step
-      INTEGER, INTENT(IN) :: kpi     ! Model grid parameters
+      TYPE(obs_surf), INTENT(INOUT) :: &
+         & surfdataqc                  ! Subset of surface data passing QC
+      INTEGER, INTENT(IN) :: kt        ! Time step
+      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
       INTEGER, INTENT(IN) :: kpj
-      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step 
-                                      !   (kit000-1 = restart time)
-      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
-      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
-         & psurf,  &    ! Model surface field
-         & psurfmask    ! Land-sea mask
-         
+      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step 
+                                       !   (kit000-1 = restart time)
+      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day
+      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
+      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
+         & psurf,  &                   ! Model surface field
+         & psurfmask                   ! Land-sea mask
+      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
+
       !! * Local declarations
       INTEGER :: ji
@@ -524 +533 @@
       INTEGER :: isurf
       INTEGER :: iobs
-      REAL(KIND=wp) :: zlam
-      REAL(KIND=wp) :: zphi
-      REAL(KIND=wp) :: zext(1), zobsmask(1)
-      REAL(kind=wp), DIMENSION(2,2,1) :: &
-         & zweig
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
-         & zmask, &
-         & zsurf, &
-         & zglam, &
-         & zgphi
+      INTEGER :: idayend
       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
          & igrdi, &
          & igrdj
+      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
+         & icount_night,      &
+         & imask_night
+      REAL(wp) :: zlam
+      REAL(wp) :: zphi
+      REAL(wp), DIMENSION(1) :: zext, zobsmask
+      REAL(wp) :: zdaystp
+      REAL(wp), DIMENSION(2,2,1) :: &
+         & zweig
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
+         & zmask,  &
+         & zsurf,  &
+         & zsurfm, &
+         & zglam,  &
+         & zgphi
+      REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
+         & zintmp,  &
+         & zouttmp, &
+         & zmeanday    ! to compute model sst in region of 24h daylight (pole)
 
       !------------------------------------------------------------------------
       ! Local initialization 
       !------------------------------------------------------------------------
-      ! ... Record and data counters
+      ! Record and data counters
       inrc = kt - kit000 + 2
       isurf = surfdataqc%nsstp(inrc)
-      
-      IF ( ld_nightav ) THEN
+
+      IF ( ldnightav ) THEN
 
       ! Initialize array for night mean
-         IF ( kt .EQ. 0 ) THEN
+         IF ( kt == 0 ) THEN
             ALLOCATE ( icount_night(kpi,kpj) )
             ALLOCATE ( imask_night(kpi,kpj) )
@@ -557 +576 @@
          ENDIF
 
-         ! Initialize daily mean for first timestep
+         ! Night-time means are calculated for night-time values over timesteps:
+         !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
          idayend = MOD( kt - kit000 + 1, kdaystp )
 
-         ! Added kt == 0 test to catch restart case 
-         IF ( idayend == 1 .OR. kt == 0) THEN
-            IF (lwp) WRITE(numout,*) 'Reset surfdataqc%vdmean on time-step: ',kt
+         ! Initialize night-time mean for first timestep of the day
+         IF ( idayend == 1 .OR. kt == 0 ) THEN
             DO jj = 1, jpj
                DO ji = 1, jpi
@@ -580 +599 @@
                ! Increment the temperature field for computing night mean and counter
                surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj)  &
-                      &                        + psurf(ji,jj)*imask_night(ji,jj)
-               zmeanday(ji,jj)        = zmeanday(ji,jj) + psurf(ji,jj)
-               icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj)
+                      &                    + psurf(ji,jj) * REAL( imask_night(ji,jj) )
+               zmeanday(ji,jj)          = zmeanday(ji,jj) + psurf(ji,jj)
+               icount_night(ji,jj)      = icount_night(ji,jj) + imask_night(ji,jj)
             END DO
          END DO
-   
-         ! Compute the daily mean at the end of day
-
+
+         ! Compute the night-time mean at the end of the day
          zdaystp = 1.0 / REAL( kdaystp )
-
-         IF ( idayend == 0 ) THEN 
+         IF ( idayend == 0 ) THEN
+            IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt
             DO jj = 1, jpj
                DO ji = 1, jpi
                   ! Test if "no night" point
-                  IF ( icount_night(ji,jj) .NE. 0 ) THEN
+                  IF ( icount_night(ji,jj) > 0 ) THEN
                      surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
-                       &                        / icount_night(ji,jj) 
+                       &                        / REAL( icount_night(ji,jj) )
                   ELSE
+                     !At locations where there is no night (e.g. poles),
+                     ! calculate daily mean instead of night-time mean.
                      surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
                   ENDIF
@@ -616 +636 @@
          & zsurf(2,2,isurf)  &
          & )
-      
+
       DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
          iobs = jobs - surfdataqc%nsurfup
@@ -639 +659 @@
 
       ! At the end of the day get interpolated means
-      IF ( idayend == 0 .AND. ld_nightav ) THEN
+      IF ( idayend == 0 .AND. ldnightav ) THEN
 
          ALLOCATE( &
@@ -656 +676 @@
 
          IF ( kt /= surfdataqc%mstp(jobs) ) THEN
-            
+
             IF(lwp) THEN
                WRITE(numout,*)
@@ -669 +689 @@
                   &            ' ntyp    = ', surfdataqc%ntyp(jobs)
             ENDIF
-            CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' )
-            
-         ENDIF
-         
+            CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
+
+         ENDIF
+
          zlam = surfdataqc%rlam(jobs)
          zphi = surfdataqc%rphi(jobs)
@@ -680 +700 @@
             &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
             &                   zmask(:,:,iobs), zweig, zobsmask )
-         
 
          ! Interpolate the model field to the observation point
-         IF ( ld_nightav ) THEN
-
-           IF ( idayend == 0 )  THEN
-               ! Daily averaged data
-               CALL obs_int_h2d( 1, 1,      & 
-                     &              zweig, zsurfm(:,:,iobs), zext )
-            ELSE 
-               CALL ctl_stop( ' ld_nightav is set to true: a nonzero' //     &
-                     &           ' number of night data should' // &
-                     &           ' only occur at the end of a given day' )
-            ENDIF
-
+         IF ( ldnightav .AND. idayend == 0 ) THEN
+            ! Night-time averaged data
+            CALL obs_int_h2d( 1, 1, zweig, zsurfm(:,:,iobs), zext )
          ELSE
-
-            CALL obs_int_h2d( 1, 1,      &
-               &              zweig, zsurf(:,:,iobs),  zext )
-
-         ENDIF
-         
-         IF ( surfdataqc%nextr == 2 ) THEN
+            CALL obs_int_h2d( 1, 1, zweig, zsurf(:,:,iobs),  zext )
+         ENDIF
+
+         IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
             ! ... Remove the MDT from the SSH at the observation point to get the SLA
             surfdataqc%rext(jobs,1) = zext(1)
@@ -722 +729 @@
          & )
 
-      ! At the end of the day also get interpolated means
-      IF ( idayend == 0 .AND. ld_nightav ) THEN
+      ! At the end of the day also deallocate night-time mean array
+      IF ( idayend == 0 .AND. ldnightav ) THEN
          DEALLOCATE( &
             & zsurfm  &
@@ -734 +741 @@
 
 END MODULE obs_oper
-

branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

@@ -71 +71 @@
          & nproc
       !! * Arguments
-      TYPE(obs_surf), INTENT(INOUT) :: surfdata    ! Full set of SLA data
-      TYPE(obs_surf), INTENT(INOUT) :: surfdataqc   ! Subset of SLA data not failing screening
+      TYPE(obs_surf), INTENT(INOUT) :: surfdata    ! Full set of surface data
+      TYPE(obs_surf), INTENT(INOUT) :: surfdataqc   ! Subset of surface data not failing screening
       LOGICAL, INTENT(IN) :: ld_nea         ! Switch for rejecting observation near land
       !! * Local declarations
@@ -99 +99 @@
       INTEGER :: inrc         ! Time index variable
 
-      IF(lwp) WRITE(numout,*)'obs_pre_sla : Preparing the SLA observations...'
-
+      IF(lwp) WRITE(numout,*)'obs_pre_surf : Preparing the surface observations...'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      
       ! Initial date initialization (year, month, day, hour, minute)
       iyea0 =   ndate0 / 10000
@@ -185 +186 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'obs_pre_surf :'
-         WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*)
-         WRITE(numout,*) ' Surface data outside time domain               = ', &
+         WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data outside time domain                  = ', &
             &            iotdobsmpp
-         WRITE(numout,*) ' Remaining surface data that failed grid search    = ', &
+         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data that failed grid search    = ', &
             &            igrdobsmpp
-         WRITE(numout,*) ' Remaining surface data outside space domain       = ', &
+         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data outside space domain       = ', &
             &            iosdsobsmpp
-         WRITE(numout,*) ' Remaining surface data at land points             = ', &
+         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data at land points             = ', &
             &            ilansobsmpp
          IF (ld_nea) THEN
-            WRITE(numout,*) ' Remaining surface data near land points (removed) = ', &
+            WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (removed) = ', &
                &            inlasobsmpp
          ELSE
-            WRITE(numout,*) ' Remaining surface data near land points (kept)    = ', &
+            WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (kept)    = ', &
                &            inlasobsmpp
          ENDIF
-         WRITE(numout,*) ' surface data accepted                             = ', &
+         WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted                             = ', &
             &            surfdataqc%nsurfmpp
 
@@ -209 +207 @@
          WRITE(numout,*) ' Number of observations per time step :'
          WRITE(numout,*)
-         WRITE(numout,1997)
-         WRITE(numout,1998)
+         WRITE(numout,'(10X,A,10X,A)')'Time step',surfdataqc%cvars(1)
+         WRITE(numout,'(10X,A,5X,A)')'---------','-----------------'
+         CALL FLUSH(numout)
       ENDIF
       
@@ -225 +224 @@
             inrc = jstp - nit000 + 2
             WRITE(numout,1999) jstp, surfdataqc%nsstpmpp(inrc)
+            CALL FLUSH(numout)
          END DO
       ENDIF
 
-1997  FORMAT(10X,'Time step',5X,'Sea level anomaly')
-1998  FORMAT(10X,'---------',5X,'-----------------')
 1999  FORMAT(10X,I9,5X,I17)
 
@@ -235 +233 @@
 
 
-   SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_vel3d, ld_nea, ld_dailyav )
-      !!----------------------------------------------------------------------
+   SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_var1, ld_var2, &
+      &                     zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2,  &
+      &                     ld_nea, kdailyavtypes )
+
+!!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_pre_prof  ***
       !!
@@ -246 +247 @@
       !!        !  2007-06  (K. Mogensen) original : T and S profile data
       !!        !  2008-09  (M. Valdivieso) : TAO velocity data
-      !!        !  2009-01  (K. Mogensen) : New feedback strictuer
+      !!        !  2009-01  (K. Mogensen) : New feedback stricture
       !!        !  2015-02  (M. Martin) : Combined profile routine.
       !!
@@ -254 +255 @@
       USE par_oce             ! Ocean parameters
       USE dom_oce, ONLY : &   ! Geographical information
-         & glamt, glamu, glamv,    &
-         & gphit, gphiu, gphiv,    &
          & gdept_1d,             &
-         & tmask, umask, vmask,  &
          & nproc
+
       !! * Arguments
       TYPE(obs_prof), INTENT(INOUT) :: profdata   ! Full set of profile data
       TYPE(obs_prof), INTENT(INOUT) :: prodatqc   ! Subset of profile data not failing screening
-      LOGICAL, INTENT(IN) :: ld_vel3d      ! Switch for zonal and meridional velocity components
-      LOGICAL, INTENT(IN) :: ld_nea        ! Switch for rejecting observation near land
-      LOGICAL, INTENT(IN) :: ld_dailyav    ! Switch for daily average data
+      LOGICAL, INTENT(IN) :: ld_var1              ! Observed variables switches
+      LOGICAL, INTENT(IN) :: ld_var2
+      LOGICAL, INTENT(IN) :: ld_nea               ! Switch for rejecting observation near land
+      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
+         & kdailyavtypes                          ! Types for daily averages
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj,jpk) :: &
+         & zmask1, &
+         & zmask2
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: &
+         & pglam1, &
+         & pglam2, &
+         & pgphi1, &
+         & pgphi2
+
       !! * Local declarations
       INTEGER :: iyea0        ! Initial date
@@ -272 +282 @@
       INTEGER :: imin0
       INTEGER :: icycle       ! Current assimilation cycle
-                              ! Counters for observations that
+                              ! Counters for observations that are
       INTEGER :: iotdobs      !  - outside time domain
-      INTEGER :: iosduobs     !  - outside space domain (zonal velocity component)
-      INTEGER :: iosdvobs     !  - outside space domain (meridional velocity component)
-      INTEGER :: ilanuobs     !  - within a model land cell (zonal velocity component)
-      INTEGER :: ilanvobs     !  - within a model land cell (meridional velocity component)
-      INTEGER :: inlauobs     !  - close to land (zonal velocity component)
-      INTEGER :: inlavobs     !  - close to land (meridional velocity component)
+      INTEGER :: iosdv1obs    !  - outside space domain (variable 1)
+      INTEGER :: iosdv2obs    !  - outside space domain (variable 2)
+      INTEGER :: ilanv1obs    !  - within a model land cell (variable 1)
+      INTEGER :: ilanv2obs    !  - within a model land cell (variable 2)
+      INTEGER :: inlav1obs    !  - close to land (variable 1)
+      INTEGER :: inlav2obs    !  - close to land (variable 2)
       INTEGER :: igrdobs      !  - fail the grid search
       INTEGER :: iuvchku      !  - reject u if v rejected and vice versa
       INTEGER :: iuvchkv      !
-                              ! Global counters for observations that
+                              ! Global counters for observations that are
       INTEGER :: iotdobsmpp   !  - outside time domain
-      INTEGER :: iosduobsmpp  !  - outside space domain (zonal velocity component)
-      INTEGER :: iosdvobsmpp  !  - outside space domain (meridional velocity component)
-      INTEGER :: ilanuobsmpp  !  - within a model land cell (zonal velocity component)
-      INTEGER :: ilanvobsmpp  !  - within a model land cell (meridional velocity component)
-      INTEGER :: inlauobsmpp  !  - close to land (zonal velocity component)
-      INTEGER :: inlavobsmpp  !  - close to land (meridional velocity component)
+      INTEGER :: iosdv1obsmpp !  - outside space domain (variable 1)
+      INTEGER :: iosdv2obsmpp !  - outside space domain (variable 2)
+      INTEGER :: ilanv1obsmpp !  - within a model land cell (variable 1)
+      INTEGER :: ilanv2obsmpp !  - within a model land cell (variable 2)
+      INTEGER :: inlav1obsmpp !  - close to land (variable 1)
+      INTEGER :: inlav2obsmpp !  - close to land (variable 2)
       INTEGER :: igrdobsmpp   !  - fail the grid search
-      INTEGER :: iuvchkumpp   !  - reject u if v rejected and vice versa
+      INTEGER :: iuvchkumpp   !  - reject var1 if var2 rejected and vice versa
       INTEGER :: iuvchkvmpp   !
       TYPE(obs_prof_valid) ::  llvalid      ! Profile selection 
       TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
-         & llvvalid           ! U,V selection 
+         & llvvalid           ! var1,var2 selection 
       INTEGER :: jvar         ! Variable loop variable
       INTEGER :: jobs         ! Obs. loop variable
@@ -302 +312 @@
       INTEGER :: inrc         ! Time index variable
 
-      IF(lwp) WRITE(numout,*)'obs_pre_prof: Preparing the profile data'
+      IF(lwp) WRITE(numout,*)'obs_pre_prof: Preparing the profile data...'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
 
       ! Initial date initialization (year, month, day, hour, minute)
@@ -317 +328 @@
       iotdobs  = 0
       igrdobs  = 0
-      iosduobs = 0
-      iosdvobs = 0
-      ilanuobs = 0
-      ilanvobs = 0
-      inlauobs = 0
-      inlavobs = 0
+      iosdv1obs = 0
+      iosdv2obs = 0
+      ilanv1obs = 0
+      ilanv2obs = 0
+      inlav1obs = 0
+      inlav2obs = 0
       iuvchku  = 0
       iuvchkv = 0
@@ -330 +341 @@
       ! -----------------------------------------------------------------------
 
-      CALL obs_coo_tim_prof( icycle, &
-         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
-         &              profdata%nprof,   profdata%nyea, profdata%nmon, &
-         &              profdata%nday,    profdata%nhou, profdata%nmin, &
-         &              profdata%ntyp,    profdata%nqc,  profdata%mstp, &
-         &              iotdobs, ld_dailyav = ld_dailyav        )
-    
+      IF ( PRESENT(kdailyavtypes) ) THEN
+         CALL obs_coo_tim_prof( icycle, &
+            &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
+            &              profdata%nprof,   profdata%nyea, profdata%nmon, &
+            &              profdata%nday,    profdata%nhou, profdata%nmin, &
+            &              profdata%ntyp,    profdata%nqc,  profdata%mstp, &
+            &              iotdobs, kdailyavtypes = kdailyavtypes )
+      ELSE
+         CALL obs_coo_tim_prof( icycle, &
+            &              iyea0,   imon0,   iday0,   ihou0,   imin0,      &
+            &              profdata%nprof,   profdata%nyea, profdata%nmon, &
+            &              profdata%nday,    profdata%nhou, profdata%nmin, &
+            &              profdata%ntyp,    profdata%nqc,  profdata%mstp, &
+            &              iotdobs )
+      ENDIF
+
       CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
       
@@ -343 +363 @@
       ! -----------------------------------------------------------------------
 
-      CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,1), profdata%mj(:,1), &
-         &              profdata%nqc,     igrdobs                         )
-      CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,2), profdata%mj(:,2), &
-         &              profdata%nqc,     igrdobs                         )
+      CALL obs_coo_grd( profdata%nprof,   profdata%mi, profdata%mj, &
+         &              profdata%nqc,     igrdobs                )
 
       CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
@@ -361 +379 @@
       ! -----------------------------------------------------------------------
 
-      ! Zonal Velocity Component
-
+      ! Variable 1
       CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(1),   &
          &                 profdata%npvsta(:,1),  profdata%npvend(:,1), &
          &                 jpi,                   jpj,                  &
          &                 jpk,                                         &
-         &                 profdata%mi,           profdata%mj,          & 
+         &                 profdata%mi,           profdata%mj,          &
          &                 profdata%var(1)%mvk,                         &
          &                 profdata%rlam,         profdata%rphi,        &
          &                 profdata%var(1)%vdep,                        &
-         &                 glamu,                 gphiu,                &
-         &                 gdept_1d,              umask,                &
+         &                 pglam1,                pgphi1,               &
+         &                 gdept_1d,              zmask1,               &
          &                 profdata%nqc,          profdata%var(1)%nvqc, &
-         &                 iosduobs,              ilanuobs,             &
-         &                 inlauobs,              ld_nea                )
-
-      CALL obs_mpp_sum_integer( iosduobs, iosduobsmpp )
-      CALL obs_mpp_sum_integer( ilanuobs, ilanuobsmpp )
-      CALL obs_mpp_sum_integer( inlauobs, inlauobsmpp )
-
-      ! Meridional Velocity Component
-
+         &                 iosdv1obs,              ilanv1obs,           &
+         &                 inlav1obs,              ld_nea                )
+
+      CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp )
+      CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp )
+      CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp )
+
+      ! Variable 2
       CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(2),   &
          &                 profdata%npvsta(:,2),  profdata%npvend(:,2), &
@@ -391 +407 @@
          &                 profdata%rlam,         profdata%rphi,        &
          &                 profdata%var(2)%vdep,                        &
-         &                 glamv,                 gphiv,                &
-         &                 gdept_1d,              vmask,                &
+         &                 pglam2,                pgphi2,               &
+         &                 gdept_1d,              zmask2,               &
          &                 profdata%nqc,          profdata%var(2)%nvqc, &
-         &                 iosdvobs,              ilanvobs,             &
-         &                 inlavobs,              ld_nea                )
-
-      CALL obs_mpp_sum_integer( iosdvobs, iosdvobsmpp )
-      CALL obs_mpp_sum_integer( ilanvobs, ilanvobsmpp )
-      CALL obs_mpp_sum_integer( inlavobs, inlavobsmpp )
+         &                 iosdv2obs,              ilanv2obs,           &
+         &                 inlav2obs,              ld_nea                )
+
+      CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp )
+      CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp )
+      CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp )
 
       ! -----------------------------------------------------------------------
@@ -405 +421 @@
       ! -----------------------------------------------------------------------
 
-      CALL obs_uv_rej( profdata, iuvchku, iuvchkv )
-      CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp )
-      CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp )
+      IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN
+         CALL obs_uv_rej( profdata, iuvchku, iuvchkv )
+         CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp )
+         CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp )
+      ENDIF
 
       ! -----------------------------------------------------------------------
@@ -446 +464 @@
       
       IF(lwp) THEN
+      
          WRITE(numout,*)
-         WRITE(numout,*) 'obs_pre_vel :'
-         WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*)
-         WRITE(numout,*) ' Profiles outside time domain                = ', &
+         WRITE(numout,*) ' Profiles outside time domain                     = ', &
             &            iotdobsmpp
-         WRITE(numout,*) ' Remaining profiles that failed grid search  = ', &
+         WRITE(numout,*) ' Remaining profiles that failed grid search       = ', &
             &            igrdobsmpp
-         WRITE(numout,*) ' Remaining U data outside space domain       = ', &
-            &            iosduobsmpp
-         WRITE(numout,*) ' Remaining U data at land points             = ', &
-            &            ilanuobsmpp
+         WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data outside space domain       = ', &
+            &            iosdv1obsmpp
+         WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data at land points             = ', &
+            &            ilanv1obsmpp
          IF (ld_nea) THEN
-            WRITE(numout,*) ' Remaining U data near land points (removed) = ',&
-               &            inlauobsmpp
+            WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (removed) = ',&
+               &            inlav1obsmpp
          ELSE
-            WRITE(numout,*) ' Remaining U data near land points (kept)    = ',&
-               &            inlauobsmpp
-         ENDIF
-         WRITE(numout,*) ' U observation rejected since V rejected     = ', &
-            &            iuvchku     
-         WRITE(numout,*) ' U data accepted                             = ', &
+            WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (kept)    = ',&
+               &            inlav1obsmpp
+         ENDIF
+         IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN
+            WRITE(numout,*) ' U observation rejected since V rejected     = ', &
+               &            iuvchku
+         ENDIF
+         WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted                             = ', &
             &            prodatqc%nvprotmpp(1)
-         WRITE(numout,*) ' Remaining V data outside space domain       = ', &
-            &            iosdvobsmpp
-         WRITE(numout,*) ' Remaining V data at land points             = ', &
-            &            ilanvobsmpp
+         WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data outside space domain       = ', &
+            &            iosdv2obsmpp
+         WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data at land points             = ', &
+            &            ilanv2obsmpp
          IF (ld_nea) THEN
-            WRITE(numout,*) ' Remaining V data near land points (removed) = ',&
-               &            inlavobsmpp
+            WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (removed) = ',&
+               &            inlav2obsmpp
          ELSE
-            WRITE(numout,*) ' Remaining V data near land points (kept)    = ',&
-               &            inlavobsmpp
-         ENDIF
-         WRITE(numout,*) ' V observation rejected since U rejected     = ', &
-            &            iuvchkv     
-         WRITE(numout,*) ' V data accepted                             = ', &
+            WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (kept)    = ',&
+               &            inlav2obsmpp
+         ENDIF
+         IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN
+            WRITE(numout,*) ' V observation rejected since U rejected     = ', &
+               &            iuvchkv
+         ENDIF
+         WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted                             = ', &
             &            prodatqc%nvprotmpp(2)
 
@@ -488 +508 @@
          WRITE(numout,*) ' Number of observations per time step :'
          WRITE(numout,*)
-         WRITE(numout,997)
+         WRITE(numout,'(10X,A,5X,A,5X,A,A)')'Time step','Profiles', &
+            &                               '     '//prodatqc%cvars(1)//'     ', &
+            &                               '     '//prodatqc%cvars(2)//'     '
          WRITE(numout,998)
       ENDIF
@@ -522 +544 @@
       ENDIF
 
-997   FORMAT(10X,'Time step',5X,'Profiles',5X,'Zonal Comp.',5X,'Meridional Comp.')
 998   FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'----------------')
 999   FORMAT(10X,I9,5X,I8,5X,I11,5X,I8)
@@ -728 +749 @@
       &                    kobsno,                                        &
       &                    kobsyea, kobsmon, kobsday, kobshou, kobsmin,   &
-      &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes, &
-      &                    ld_dailyav )
+      &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE obs_coo_tim ***
@@ -773 +793 @@
       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
          & kdailyavtypes    ! Types for daily averages
-      LOGICAL, OPTIONAL :: ld_dailyav    ! All types are daily averages
       !! * Local declarations
       INTEGER :: jobs
@@ -807 +826 @@
       ENDIF
 
-      !------------------------------------------------------------------------
-      ! If ld_dailyav is set then all data assumed to be daily averaged
-      !------------------------------------------------------------------------
-      
-      IF ( PRESENT( ld_dailyav) ) THEN
-         IF (ld_dailyav) THEN
-            DO jobs = 1, kobsno
-               
-               IF ( kobsqc(jobs) <= 10 ) THEN
-                  
-                  IF ( kobsstp(jobs) == (nit000 - 1) ) THEN
-                     kobsqc(jobs) = kobsqc(jobs) + 14
-                     kotdobs      = kotdobs + 1
-                     CYCLE
-                  ENDIF
-                  
-               ENDIF
-            END DO
-         ENDIF
-      ENDIF
 
    END SUBROUTINE obs_coo_tim_prof

branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90

@@ -104 +104 @@
       ! Bookkeeping arrays with sizes equal to number of variables
 
+      CHARACTER(len=6), POINTER, DIMENSION(:) :: &
+         & cvars          !: Variable names
+
       INTEGER, POINTER, DIMENSION(:) :: &
          & nvprot,   &    !: Local total number of profile T data
@@ -237 +240 @@
 
       ALLOCATE( &
+         & prof%cvars(kvar),    &
          & prof%nvprot(kvar),   &
          & prof%nvprotmpp(kvar) &
@@ -242 +246 @@
          
       DO jvar = 1, kvar
+         prof%cvars    (jvar) = "NotSet"
          prof%nvprot   (jvar) = ko3dt(jvar)
          prof%nvprotmpp(jvar) = 0
@@ -452 +457 @@
 
       DEALLOCATE( &
-         & prof%nvprot,  &
+         & prof%cvars,    &
+         & prof%nvprot,   &
          & prof%nvprotmpp &
          )
@@ -770 +776 @@
       newprof%npj      = prof%npj
       newprof%npk      = prof%npk
+      newprof%cvars(:) = prof%cvars(:)
  
       ! Deallocate temporary data

branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_altbias.F90

@@ -50 +50 @@
 CONTAINS
 
-   SUBROUTINE obs_rea_altbias( kslano, sladata, k2dint, bias_file )
+   SUBROUTINE obs_rea_altbias( sladata, k2dint, bias_file )
       !!---------------------------------------------------------------------
       !!
@@ -70 +70 @@
       !
       !! * Arguments
-      INTEGER, INTENT(IN) :: kslano      ! Number of SLA Products
-      TYPE(obs_surf), DIMENSION(kslano), INTENT(INOUT) :: &
+      TYPE(obs_surf), INTENT(INOUT) :: &
          & sladata       ! SLA data
       INTEGER, INTENT(IN) :: k2dint
@@ -80 +79 @@
       CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_rea_altbias'
 
-      INTEGER :: jslano       ! Data set loop variable
       INTEGER :: jobs         ! Obs loop variable
       INTEGER :: jpialtbias   ! Number of grid point in latitude for the bias
@@ -144 +142 @@
       ! Intepolate the bias already on the model grid at the observation point
   
-      DO jslano = 1, kslano
-
-         ALLOCATE( &
-            & igrdi(2,2,sladata(jslano)%nsurf), &
-            & igrdj(2,2,sladata(jslano)%nsurf), &
-            & zglam(2,2,sladata(jslano)%nsurf), &
-            & zgphi(2,2,sladata(jslano)%nsurf), &
-            & zmask(2,2,sladata(jslano)%nsurf), &
-            & zbias(2,2,sladata(jslano)%nsurf)  &
-            & )
-         
-         DO jobs = 1, sladata(jslano)%nsurf
-
-            igrdi(1,1,jobs) = sladata(jslano)%mi(jobs)-1
-            igrdj(1,1,jobs) = sladata(jslano)%mj(jobs)-1
-            igrdi(1,2,jobs) = sladata(jslano)%mi(jobs)-1
-            igrdj(1,2,jobs) = sladata(jslano)%mj(jobs)
-            igrdi(2,1,jobs) = sladata(jslano)%mi(jobs)
-            igrdj(2,1,jobs) = sladata(jslano)%mj(jobs)-1
-            igrdi(2,2,jobs) = sladata(jslano)%mi(jobs)
-            igrdj(2,2,jobs) = sladata(jslano)%mj(jobs)
-
-         END DO
-
-         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, &
-            &                  igrdi, igrdj, glamt, zglam )
-         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, &
-            &                  igrdi, igrdj, gphit, zgphi )
-         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, &
-            &                  igrdi, igrdj, tmask(:,:,1), zmask )
-         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, &
-            &                  igrdi, igrdj, z_altbias, zbias )
-
-         DO jobs = 1, sladata(jslano)%nsurf
-
-            zlam = sladata(jslano)%rlam(jobs)
-            zphi = sladata(jslano)%rphi(jobs)
-            iico = sladata(jslano)%mi(jobs)
-            ijco = sladata(jslano)%mj(jobs)
+      ALLOCATE( &
+         & igrdi(2,2,sladata%nsurf), &
+         & igrdj(2,2,sladata%nsurf), &
+         & zglam(2,2,sladata%nsurf), &
+         & zgphi(2,2,sladata%nsurf), &
+         & zmask(2,2,sladata%nsurf), &
+         & zbias(2,2,sladata%nsurf)  &
+         & )
+         
+      DO jobs = 1, sladata%nsurf
+
+         igrdi(1,1,jobs) = sladata%mi(jobs)-1
+         igrdj(1,1,jobs) = sladata%mj(jobs)-1
+         igrdi(1,2,jobs) = sladata%mi(jobs)-1
+         igrdj(1,2,jobs) = sladata%mj(jobs)
+         igrdi(2,1,jobs) = sladata%mi(jobs)
+         igrdj(2,1,jobs) = sladata%mj(jobs)-1
+         igrdi(2,2,jobs) = sladata%mi(jobs)
+         igrdj(2,2,jobs) = sladata%mj(jobs)
+
+      END DO
+
+      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, &
+         &                  igrdi, igrdj, glamt, zglam )
+      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, &
+         &                  igrdi, igrdj, gphit, zgphi )
+      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, &
+         &                  igrdi, igrdj, tmask(:,:,1), zmask )
+      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, &
+         &                  igrdi, igrdj, z_altbias, zbias )
+
+      DO jobs = 1, sladata%nsurf
+
+         zlam = sladata%rlam(jobs)
+         zphi = sladata%rphi(jobs)
+         iico = sladata%mi(jobs)
+         ijco = sladata%mj(jobs)
             
-            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
-               &                   zglam(:,:,jobs), zgphi(:,:,jobs), &
-               &                   zmask(:,:,jobs), zweig, zobsmask )
+         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
+            &                   zglam(:,:,jobs), zgphi(:,:,jobs), &
+            &                   zmask(:,:,jobs), zweig, zobsmask )
             
-            CALL obs_int_h2d( 1, 1,      &
-               &              zweig, zbias(:,:,jobs),  zext )
-
-            ! adjust mdt with bias field
-            sladata(jslano)%rext(jobs,2) = &
-               sladata(jslano)%rext(jobs,2) - zext(1)
+         CALL obs_int_h2d( 1, 1,      &
+            &              zweig, zbias(:,:,jobs),  zext )
+
+         ! adjust mdt with bias field
+         sladata%rext(jobs,2) = sladata%rext(jobs,2) - zext(1)
             
-         END DO
-
-         DEALLOCATE( &
-            & igrdi, &
-            & igrdj, &
-            & zglam, &
-            & zgphi, &
-            & zmask, &
-            & zbias  &
-            & )
-         
       END DO
 
+      DEALLOCATE( &
+         & igrdi, &
+         & igrdj, &
+         & zglam, &
+         & zgphi, &
+         & zmask, &
+         & zbias  &
+         & )
+         
       CALL wrk_dealloc(jpi,jpj,z_altbias) 
 

branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90

@@ -25 +25 @@
    USE netcdf                   ! NetCDF library
    USE obs_oper                 ! Observation operators
-   USE obs_prof_io              ! Profile files I/O (non-FB files)
    USE lib_mpp                  ! For ctl_warn/stop
+   USE obs_fbm                  ! Feedback routines
 
    IMPLICIT NONE
@@ -42 +42 @@
 
 CONTAINS
- 
-   SUBROUTINE obs_rea_prof( profdata, knumfiles, cfilenames, &
+
+   SUBROUTINE obs_rea_prof( profdata, knumfiles, cdfilenames, &
       &                     kvars, kextr, kstp, ddobsini, ddobsend, &
-      &                     ldt3d, lds3d, ldignmis, ldsatt, ldavtimset, &
+      &                     ldvar1, ldvar2, ldignmis, ldsatt, &
       &                     ldmod, kdailyavtypes )
       !!---------------------------------------------------------------------
@@ -62 +62 @@
       !! History :  
       !!      ! :  2009-09 (K. Mogensen) : New merged version of old routines
+      !!      ! :  2015-08 (M. Martin) : Merged profile and velocity routines
       !!----------------------------------------------------------------------
-      !! * Modules used
-   
+
       !! * Arguments
-      TYPE(obs_prof), INTENT(OUT) ::  profdata     ! Profile data to be read
-      INTEGER, INTENT(IN) :: knumfiles      ! Number of files to read in
+      TYPE(obs_prof), INTENT(OUT) :: &
+         & profdata                     ! Profile data to be read
+      INTEGER, INTENT(IN) :: knumfiles  ! Number of files to read
       CHARACTER(LEN=128), INTENT(IN) ::  &
-         & cfilenames(knumfiles)  ! File names to read in
+         & cdfilenames(knumfiles)        ! File names to read in
       INTEGER, INTENT(IN) :: kvars      ! Number of variables in profdata
-      INTEGER, INTENT(IN) :: kextr      ! Number of extra fields for each var in profdata
-      INTEGER, INTENT(IN) :: kstp        ! Ocean time-step index
-      LOGICAL, INTENT(IN) :: ldt3d       ! Observed variables switches
-      LOGICAL, INTENT(IN) :: lds3d
-      LOGICAL, INTENT(IN) :: ldignmis    ! Ignore missing files
-      LOGICAL, INTENT(IN) :: ldsatt      ! Compute salinity at all temperature points
-      LOGICAL, INTENT(IN) :: ldavtimset  ! Correct time for daily averaged data
-      LOGICAL, INTENT(IN) :: ldmod       ! Initialize model from input data
-      REAL(KIND=dp), INTENT(IN) :: ddobsini    ! Obs. ini time in YYYYMMDD.HHMMSS
-      REAL(KIND=dp), INTENT(IN) :: ddobsend    ! Obs. end time in YYYYMMDD.HHMMSS
+      INTEGER, INTENT(IN) :: kextr      ! Number of extra fields for each var
+      INTEGER, INTENT(IN) :: kstp       ! Ocean time-step index
+      LOGICAL, INTENT(IN) :: ldvar1     ! Observed variables switches
+      LOGICAL, INTENT(IN) :: ldvar2
+      LOGICAL, INTENT(IN) :: ldignmis   ! Ignore missing files
+      LOGICAL, INTENT(IN) :: ldsatt     ! Compute salinity at all temperature points
+      LOGICAL, INTENT(IN) :: ldmod      ! Initialize model from input data
+      REAL(dp), INTENT(IN) :: ddobsini  ! Obs. ini time in YYYYMMDD.HHMMSS
+      REAL(dp), INTENT(IN) :: ddobsend  ! Obs. end time in YYYYMMDD.HHMMSS
       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
-         & kdailyavtypes
+         & kdailyavtypes                ! Types of daily average observations
 
       !! * Local declarations
       CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof'
+      CHARACTER(len=8) :: clrefdate
+      CHARACTER(len=6), DIMENSION(:), ALLOCATABLE :: clvars
       INTEGER :: jvar
       INTEGER :: ji
@@ -101 +103 @@
       INTEGER :: imin
       INTEGER :: isec
+      INTEGER :: iprof
+      INTEGER :: iproftot
+      INTEGER :: ivar1t0
+      INTEGER :: ivar2t0
+      INTEGER :: ivar1t
+      INTEGER :: ivar2t
+      INTEGER :: ip3dt
+      INTEGER :: ios
+      INTEGER :: ioserrcount
+      INTEGER :: ivar1tmpp
+      INTEGER :: ivar2tmpp
+      INTEGER :: ip3dtmpp
+      INTEGER :: itype
       INTEGER, DIMENSION(knumfiles) :: &
          & irefdate
       INTEGER, DIMENSION(ntyp1770+1) :: &
-         & itypt,    &
-         & ityptmpp, &
-         & ityps,    &
-         & itypsmpp 
-      INTEGER :: it3dtmpp
-      INTEGER :: is3dtmpp
-      INTEGER :: ip3dtmpp
+         & itypvar1,    &
+         & itypvar1mpp, &
+         & itypvar2,    &
+         & itypvar2mpp 
       INTEGER, DIMENSION(:), ALLOCATABLE :: &
          & iobsi,    &
@@ -118 +130 @@
          & ifileidx, &
          & iprofidx
-      INTEGER :: itype
       INTEGER, DIMENSION(imaxavtypes) :: &
          & idailyavtypes
+      INTEGER, DIMENSION(kvars) :: &
+         & iv3dt
       REAL(wp), DIMENSION(:), ALLOCATABLE :: &
          & zphi, &
          & zlam
-      real(wp), DIMENSION(:), ALLOCATABLE :: &
+      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
          & zdat
+      REAL(wp), DIMENSION(knumfiles) :: &
+         & djulini, &
+         & djulend
       LOGICAL :: llvalprof
+      LOGICAL :: lldavtimset
       TYPE(obfbdata), POINTER, DIMENSION(:) :: &
          & inpfiles
-      real(wp), DIMENSION(knumfiles) :: &
-         & djulini, &
-         & djulend
-      INTEGER :: iprof
-      INTEGER :: iproftot
-      INTEGER :: it3dt0
-      INTEGER :: is3dt0
-      INTEGER :: it3dt
-      INTEGER :: is3dt
-      INTEGER :: ip3dt
-      INTEGER :: ios
-      INTEGER :: ioserrcount
-      INTEGER, DIMENSION(kvars) :: &
-         & iv3dt
-      CHARACTER(len=8) :: cl_refdate
-   
+
       ! Local initialization
       iprof = 0
-      it3dt0 = 0
-      is3dt0 = 0
+      ivar1t0 = 0
+      ivar2t0 = 0
       ip3dt = 0
 
       ! Daily average types
+      lldavtimset = .FALSE.
       IF ( PRESENT(kdailyavtypes) ) THEN
          idailyavtypes(:) = kdailyavtypes(:)
+         IF ( ANY (idailyavtypes(:) /= -1) ) lldavtimset = .TRUE.
       ELSE
          idailyavtypes(:) = -1
@@ -159 +163 @@
 
       !-----------------------------------------------------------------------
-      ! Check data the model part is just with feedback data files
-      !-----------------------------------------------------------------------
-      IF ( ldmod .AND. ( kformat /= 0 ) ) THEN
-         CALL ctl_stop( 'Model can only be read from feedback data' )
-         RETURN
-      ENDIF
-
-      !-----------------------------------------------------------------------
       ! Count the number of files needed and allocate the obfbdata type
       !-----------------------------------------------------------------------
-      
+
       inobf = knumfiles
-      
+
       ALLOCATE( inpfiles(inobf) )
 
       prof_files : DO jj = 1, inobf
-          
+
          !---------------------------------------------------------------------
          ! Prints
@@ -182 +178 @@
             WRITE(numout,*)
             WRITE(numout,*) ' obs_rea_pro_dri : Reading from file = ', &
-               & TRIM( TRIM( cfilenames(jj) ) )
+               & TRIM( TRIM( cdfilenames(jj) ) )
             WRITE(numout,*) ' ~~~~~~~~~~~~~~~'
             WRITE(numout,*)
@@ -190 +186 @@
          !  Initialization: Open file and get dimensions only
          !---------------------------------------------------------------------
-         
-         iflag = nf90_open( TRIM( cfilenames(jj) ), nf90_nowrite, &
+
+         iflag = nf90_open( TRIM( cdfilenames(jj) ), nf90_nowrite, &
             &                      i_file_id )
-         
+
          IF ( iflag /= nf90_noerr ) THEN
 
             IF ( ldignmis ) THEN
                inpfiles(jj)%nobs = 0
-               CALL ctl_warn( 'File ' // TRIM( cfilenames(jj) ) // &
+               CALL ctl_warn( 'File ' // TRIM( cdfilenames(jj) ) // &
                   &           ' not found' )
             ELSE 
-               CALL ctl_stop( 'File ' // TRIM( cfilenames(jj) ) // &
+               CALL ctl_stop( 'File ' // TRIM( cdfilenames(jj) ) // &
                   &           ' not found' )
             ENDIF
 
          ELSE 
-            
+
             !------------------------------------------------------------------
-            !  Close the file since it is opened in read_proffile
+            !  Close the file since it is opened in read_obfbdata
             !------------------------------------------------------------------
-            
+
             iflag = nf90_close( i_file_id )
 
@@ -217 +213 @@
             !------------------------------------------------------------------
             CALL init_obfbdata( inpfiles(jj) )
-            IF(lwp) THEN
-               WRITE(numout,*)
-               WRITE(numout,*)'Reading from feedback file :', &
-                  &           TRIM( cfilenames(jj) )
-            ENDIF
-            CALL read_obfbdata( TRIM( cfilenames(jj) ), inpfiles(jj), &
+            CALL read_obfbdata( TRIM( cdfilenames(jj) ), inpfiles(jj), &
                &                ldgrid = .TRUE. )
-               
+
             IF ( inpfiles(jj)%nvar < 2 ) THEN
-               CALL ctl_stop( 'Feedback format error' )
-               RETURN
-            ENDIF
+               CALL ctl_stop( 'Feedback format error: ', &
+                  &           ' less than 2 vars in profile file' )
+            ENDIF
+
             IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN
                CALL ctl_stop( 'Model not in input data' )
-               RETURN
-            ENDIF
-            
+            ENDIF
+
+            IF ( jj == 1 ) THEN
+               ALLOCATE( clvars( inpfiles(jj)%nvar ) )
+               DO ji = 1, inpfiles(jj)%nvar
+                 clvars(ji) = inpfiles(jj)%cname(ji)
+               END DO
+            ELSE
+               DO ji = 1, inpfiles(jj)%nvar
+                  IF ( inpfiles(jj)%cname(ji) /= clvars(ji) ) THEN
+                     CALL ctl_stop( 'Feedback file variables not consistent', &
+                        &           ' with previous files for this type' )
+                  ENDIF
+               END DO
+            ENDIF
+
             !------------------------------------------------------------------
             !  Change longitude (-180,180)
@@ -251 +256 @@
             !  Calculate the date  (change eventually)
             !------------------------------------------------------------------
-            cl_refdate=inpfiles(jj)%cdjuldref(1:8)
-            READ(cl_refdate,'(I8)') irefdate(jj)
-            
+            clrefdate=inpfiles(jj)%cdjuldref(1:8)
+            READ(clrefdate,'(I8)') irefdate(jj)
+
             CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec )
             CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), &
@@ -262 +267 @@
 
             ioserrcount=0
-            IF ( ldavtimset ) THEN
+            IF ( lldavtimset ) THEN
+
+               IF ( ANY ( idailyavtypes(:) /= -1 ) .AND. lwp) THEN
+                  WRITE(numout,*)' Resetting time of daily averaged', &
+                     &           ' observations to the end of the day'
+               ENDIF
+
                DO ji = 1, inpfiles(jj)%nobs
-                  ! 
-                  !  for daily averaged data for example
-                  !  MRB data (itype==820) force the time
-                  !  to be the  end of the day
-                  !
                   READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 900 ) itype
 900               IF ( ios /= 0 ) THEN
-                     itype = 0         ! Set type to zero if there is a problem in the string conversion
+                     ! Set type to zero if there is a problem in the string conversion
+                     itype = 0
                   ENDIF
-                  IF ( ANY (idailyavtypes == itype ) ) THEN
-                     inpfiles(jj)%ptim(ji) = &
-                     & INT(inpfiles(jj)%ptim(ji)) + 1
+
+                  IF ( ANY ( idailyavtypes(:) == itype ) ) THEN
+                  !  for daily averaged data force the time
+                  !  to be the last time-step of the day, but still within the day.
+                     IF ( inpfiles(jj)%ptim(ji) >= 0. ) THEN
+                        inpfiles(jj)%ptim(ji) = &
+                           & INT(inpfiles(jj)%ptim(ji)) + 0.9999
+                     ELSE
+                        inpfiles(jj)%ptim(ji) = &
+                           & INT(inpfiles(jj)%ptim(ji)) - 0.0001
+                     ENDIF
                   ENDIF
+
                END DO
-            ENDIF
-            
+
+            ENDIF
+
             IF ( inpfiles(jj)%nobs > 0 ) THEN
                inpfiles(jj)%iproc = -1
@@ -342 +359 @@
                   ENDIF
                   llvalprof = .FALSE.
-                  IF ( ldt3d ) THEN
+                  IF ( ldvar1 ) THEN
                      loop_t_count : DO ij = 1,inpfiles(jj)%nlev
                         IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
@@ -348 +365 @@
                         IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
                            & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
-                           it3dt0 = it3dt0 + 1
+                           ivar1t0 = ivar1t0 + 1
                         ENDIF
                      END DO loop_t_count
                   ENDIF
-                  IF ( lds3d ) THEN
+                  IF ( ldvar2 ) THEN
                      loop_s_count : DO ij = 1,inpfiles(jj)%nlev
                         IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
@@ -358 +375 @@
                         IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
                            & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
-                           is3dt0 = is3dt0 + 1
+                           ivar2t0 = ivar2t0 + 1
                         ENDIF
                      END DO loop_s_count
@@ -367 +384 @@
                      IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
                         &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
-                        &     ldt3d ) .OR. &
+                        &     ldvar1 ) .OR. &
                         & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
                         &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
-                        &     lds3d ) ) THEN
+                        &     ldvar2 ) ) THEN
                         ip3dt = ip3dt + 1
                         llvalprof = .TRUE.
@@ -384 +401 @@
 
       END DO prof_files
-      
+
       !-----------------------------------------------------------------------
       ! Get the time ordered indices of the input data
@@ -425 +442 @@
          &               zdat,     &
          &               iindx   )
-      
+
       iv3dt(:) = -1
       IF (ldsatt) THEN
@@ -431 +448 @@
          iv3dt(2) = ip3dt
       ELSE
-         iv3dt(1) = it3dt0
-         iv3dt(2) = is3dt0
+         iv3dt(1) = ivar1t0
+         iv3dt(2) = ivar2t0
       ENDIF
       CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, &
          &                 kstp, jpi, jpj, jpk )
-      
+
       ! * Read obs/positions, QC, all variable and assign to profdata
 
       profdata%nprof     = 0
       profdata%nvprot(:) = 0
-
+      profdata%cvars(:)  = clvars(:)
       iprof = 0
 
       ip3dt = 0
-      it3dt = 0
-      is3dt = 0
-      itypt   (:) = 0
-      ityptmpp(:) = 0
-      
-      ityps   (:) = 0
-      itypsmpp(:) = 0
-      
-      ioserrcount = 0      
+      ivar1t = 0
+      ivar2t = 0
+      itypvar1   (:) = 0
+      itypvar1mpp(:) = 0
+
+      itypvar2   (:) = 0
+      itypvar2mpp(:) = 0
+
+      ioserrcount = 0
       DO jk = 1, iproftot
-         
+
          jj = ifileidx(iindx(jk))
          ji = iprofidx(iindx(jk))
@@ -465 +482 @@
          IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
             & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
-            
+
             IF ( nproc == 0 ) THEN
                IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
@@ -471 +488 @@
                IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
             ENDIF
-            
+
             llvalprof = .FALSE.
 
@@ -480 +497 @@
 
             loop_prof : DO ij = 1, inpfiles(jj)%nlev
-               
+
                IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
                   & CYCLE
-               
+
                IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
                   & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
-                  
+
                   llvalprof = .TRUE. 
                   EXIT loop_prof
-                  
+
                ENDIF
-               
+
                IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
                   & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
-                  
+
                   llvalprof = .TRUE. 
                   EXIT loop_prof
-                  
+
                ENDIF
-               
+
             END DO loop_prof
-            
+
             ! Set profile information
-            
+
             IF ( llvalprof ) THEN
-               
+
                iprof = iprof + 1
 
@@ -524 +541 @@
                profdata%nhou(iprof) = ihou
                profdata%nmin(iprof) = imin
-               
+
                ! Profile space coordinates
                profdata%rlam(iprof) = inpfiles(jj)%plam(ji)
@@ -532 +549 @@
                profdata%mi  (iprof,:) = inpfiles(jj)%iobsi(ji,1)
                profdata%mj  (iprof,:) = inpfiles(jj)%iobsj(ji,1)
-               
+
                ! Profile WMO number
                profdata%cwmo(iprof) = inpfiles(jj)%cdwmo(ji)
-               
+
                ! Instrument type
                READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype
@@ -543 +560 @@
                   itype = 0
                ENDIF
-               
+
                profdata%ntyp(iprof) = itype
-               
+
                ! QC stuff
 
@@ -564 +581 @@
                profdata%nqc(iprof)  = 0 !TODO
 
-               loop_p : DO ij = 1, inpfiles(jj)%nlev            
-                  
+               loop_p : DO ij = 1, inpfiles(jj)%nlev
+
                   IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
                      & CYCLE
@@ -573 +590 @@
                      IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
                         &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
-                        &     ldt3d ) .OR. &
+                        &     ldvar1 ) .OR. &
                         & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
                         &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
-                        &     lds3d ) ) THEN
+                        &     ldvar2 ) ) THEN
                         ip3dt = ip3dt + 1
                      ELSE
                         CYCLE
                      ENDIF
-                     
+
                   ENDIF
 
                   IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
                      &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
-                     &       ldt3d ) .OR. ldsatt ) THEN
-                     
+                     &       ldvar1 ) .OR. ldsatt ) THEN
+
                      IF (ldsatt) THEN
 
-                        it3dt = ip3dt
+                        ivar1t = ip3dt
 
                      ELSE
 
-                        it3dt = it3dt + 1
-                        
+                        ivar1t = ivar1t + 1
+
                      ENDIF
 
-                     ! Depth of T observation
-                     profdata%var(1)%vdep(it3dt) = &
+                     ! Depth of var1 observation
+                     profdata%var(1)%vdep(ivar1t) = &
                         &                inpfiles(jj)%pdep(ij,ji)
-                     
-                     ! Depth of T observation QC
-                     profdata%var(1)%idqc(it3dt) = &
+
+                     ! Depth of var1 observation QC
+                     profdata%var(1)%idqc(ivar1t) = &
                         &                inpfiles(jj)%idqc(ij,ji)
-                     
-                     ! Depth of T observation QC flags
-                     profdata%var(1)%idqcf(:,it3dt) = &
+
+                     ! Depth of var1 observation QC flags
+                     profdata%var(1)%idqcf(:,ivar1t) = &
                         &                inpfiles(jj)%idqcf(:,ij,ji)
-                     
+
                      ! Profile index
-                     profdata%var(1)%nvpidx(it3dt) = iprof
-                     
+                     profdata%var(1)%nvpidx(ivar1t) = iprof
+
                      ! Vertical index in original profile
-                     profdata%var(1)%nvlidx(it3dt) = ij
-
-                     ! Profile potential T value
+                     profdata%var(1)%nvlidx(ivar1t) = ij
+
+                     ! Profile potential var1 value
                      IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &
                         & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
-                        profdata%var(1)%vobs(it3dt) = &
+                        profdata%var(1)%vobs(ivar1t) = &
                            &                inpfiles(jj)%pob(ij,ji,1)
                         IF ( ldmod ) THEN
-                           profdata%var(1)%vmod(it3dt) = &
+                           profdata%var(1)%vmod(ivar1t) = &
                               &                inpfiles(jj)%padd(ij,ji,1,1)
                         ENDIF
-                        ! Count number of profile T data as function of type
-                        itypt( profdata%ntyp(iprof) + 1 ) = &
-                           & itypt( profdata%ntyp(iprof) + 1 ) + 1
+                        ! Count number of profile var1 data as function of type
+                        itypvar1( profdata%ntyp(iprof) + 1 ) = &
+                           & itypvar1( profdata%ntyp(iprof) + 1 ) + 1
                      ELSE
-                        profdata%var(1)%vobs(it3dt) = fbrmdi
+                        profdata%var(1)%vobs(ivar1t) = fbrmdi
                      ENDIF
 
-                     ! Profile T qc
-                     profdata%var(1)%nvqc(it3dt) = &
+                     ! Profile var1 qc
+                     profdata%var(1)%nvqc(ivar1t) = &
                         & inpfiles(jj)%ivlqc(ij,ji,1)
 
-                     ! Profile T qc flags
-                     profdata%var(1)%nvqcf(:,it3dt) = &
+                     ! Profile var1 qc flags
+                     profdata%var(1)%nvqcf(:,ivar1t) = &
                         & inpfiles(jj)%ivlqcf(:,ij,ji,1)
 
                      ! Profile insitu T value
-                     profdata%var(1)%vext(it3dt,1) = &
+                     profdata%var(1)%vext(ivar1t,1) = &
                         &                inpfiles(jj)%pext(ij,ji,1)
-                     
+
                   ENDIF
-                  
+
                   IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
                      &   ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &
-                     &   lds3d ) .OR. ldsatt ) THEN
-                     
+                     &   ldvar2 ) .OR. ldsatt ) THEN
+
                      IF (ldsatt) THEN
 
-                        is3dt = ip3dt
+                        ivar2t = ip3dt
 
                      ELSE
 
-                        is3dt = is3dt + 1
-                        
+                        ivar2t = ivar2t + 1
+
                      ENDIF
 
-                     ! Depth of S observation
-                     profdata%var(2)%vdep(is3dt) = &
+                     ! Depth of var2 observation
+                     profdata%var(2)%vdep(ivar2t) = &
                         &                inpfiles(jj)%pdep(ij,ji)
-                     
-                     ! Depth of S observation QC
-                     profdata%var(2)%idqc(is3dt) = &
+
+                     ! Depth of var2 observation QC
+                     profdata%var(2)%idqc(ivar2t) = &
                         &                inpfiles(jj)%idqc(ij,ji)
-                     
-                     ! Depth of S observation QC flags
-                     profdata%var(2)%idqcf(:,is3dt) = &
+
+                     ! Depth of var2 observation QC flags
+                     profdata%var(2)%idqcf(:,ivar2t) = &
                         &                inpfiles(jj)%idqcf(:,ij,ji)
-                     
+
                      ! Profile index
-                     profdata%var(2)%nvpidx(is3dt) = iprof
-                     
+                     profdata%var(2)%nvpidx(ivar2t) = iprof
+
                      ! Vertical index in original profile
-                     profdata%var(2)%nvlidx(is3dt) = ij
-
-                     ! Profile S value
+                     profdata%var(2)%nvlidx(ivar2t) = ij
+
+                     ! Profile var2 value
                      IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &
                         & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN
-                        profdata%var(2)%vobs(is3dt) = &
+                        profdata%var(2)%vobs(ivar2t) = &
                            &                inpfiles(jj)%pob(ij,ji,2)
                         IF ( ldmod ) THEN
-                           profdata%var(2)%vmod(is3dt) = &
+                           profdata%var(2)%vmod(ivar2t) = &
                               &                inpfiles(jj)%padd(ij,ji,1,2)
                         ENDIF
-                        ! Count number of profile S data as function of type
-                        ityps( profdata%ntyp(iprof) + 1 ) = &
-                           & ityps( profdata%ntyp(iprof) + 1 ) + 1
+                        ! Count number of profile var2 data as function of type
+                        itypvar2( profdata%ntyp(iprof) + 1 ) = &
+                           & itypvar2( profdata%ntyp(iprof) + 1 ) + 1
                      ELSE
-                        profdata%var(2)%vobs(is3dt) = fbrmdi
+                        profdata%var(2)%vobs(ivar2t) = fbrmdi
                      ENDIF
-                     
-                     ! Profile S qc
-                     profdata%var(2)%nvqc(is3dt) = &
+
+                     ! Profile var2 qc
+                     profdata%var(2)%nvqc(ivar2t) = &
                         & inpfiles(jj)%ivlqc(ij,ji,2)
 
-                     ! Profile S qc flags
-                     profdata%var(2)%nvqcf(:,is3dt) = &
+                     ! Profile var2 qc flags
+                     profdata%var(2)%nvqcf(:,ivar2t) = &
                         & inpfiles(jj)%ivlqcf(:,ij,ji,2)
 
                   ENDIF
-            
+
                END DO loop_p
 
@@ -715 +732 @@
       ! Sum up over processors
       !-----------------------------------------------------------------------
-      
-      CALL obs_mpp_sum_integer ( it3dt0, it3dtmpp )
-      CALL obs_mpp_sum_integer ( is3dt0, is3dtmpp )
-      CALL obs_mpp_sum_integer ( ip3dt, ip3dtmpp )
-      
-      CALL obs_mpp_sum_integers( itypt, ityptmpp, ntyp1770 + 1 )
-      CALL obs_mpp_sum_integers( ityps, itypsmpp, ntyp1770 + 1 )
-      
+
+      CALL obs_mpp_sum_integer ( ivar1t0, ivar1tmpp )
+      CALL obs_mpp_sum_integer ( ivar2t0, ivar2tmpp )
+      CALL obs_mpp_sum_integer ( ip3dt,   ip3dtmpp  )
+
+      CALL obs_mpp_sum_integers( itypvar1, itypvar1mpp, ntyp1770 + 1 )
+      CALL obs_mpp_sum_integers( itypvar2, itypvar2mpp, ntyp1770 + 1 )
+
       !-----------------------------------------------------------------------
       ! Output number of observations.
@@ -728 +745 @@
       IF(lwp) THEN
          WRITE(numout,*) 
-         WRITE(numout,'(1X,A)') 'Profile data'
+         WRITE(numout,'(A)') ' Profile data'
          WRITE(numout,'(1X,A)') '------------'
          WRITE(numout,*) 
-         WRITE(numout,'(1X,A)') 'Profile T data'
-         WRITE(numout,'(1X,A)') '--------------'
+         WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(1) )
+         WRITE(numout,'(1X,A)') '------------------------'
          DO ji = 0, ntyp1770
-            IF ( ityptmpp(ji+1) > 0 ) THEN
+            IF ( itypvar1mpp(ji+1) > 0 ) THEN
                WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
                   & cwmonam1770(ji)(1:52),' = ', &
-                  & ityptmpp(ji+1)
+                  & itypvar1mpp(ji+1)
             ENDIF
          END DO
@@ -743 +760 @@
             & '---------------------------------------------------------------'
          WRITE(numout,'(1X,A55,I8)') &
-            & 'Total profile T data                                 = ',&
-            & it3dtmpp
+            & 'Total profile data for variable '//TRIM( profdata%cvars(1) )// &
+            & '             = ', ivar1tmpp
          WRITE(numout,'(1X,A)') &
             & '---------------------------------------------------------------'
          WRITE(numout,*) 
-         WRITE(numout,'(1X,A)') 'Profile S data'
-         WRITE(numout,'(1X,A)') '--------------'
+         WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(2) )
+         WRITE(numout,'(1X,A)') '------------------------'
          DO ji = 0, ntyp1770
-            IF ( itypsmpp(ji+1) > 0 ) THEN
+            IF ( itypvar2mpp(ji+1) > 0 ) THEN
                WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
                   & cwmonam1770(ji)(1:52),' = ', &
-                  & itypsmpp(ji+1)
+                  & itypvar2mpp(ji+1)
             ENDIF
          END DO
@@ -760 +777 @@
             & '---------------------------------------------------------------'
          WRITE(numout,'(1X,A55,I8)') &
-            & 'Total profile S data                                 = ',&
-            & is3dtmpp
+            & 'Total profile data for variable '//TRIM( profdata%cvars(2) )// &
+            & '             = ', ivar2tmpp
          WRITE(numout,'(1X,A)') &
             & '---------------------------------------------------------------'
          WRITE(numout,*) 
       ENDIF
-      
+
       IF (ldsatt) THEN
          profdata%nvprot(1)    = ip3dt
@@ -773 +790 @@
          profdata%nvprotmpp(2) = ip3dtmpp
       ELSE
-         profdata%nvprot(1)    = it3dt
-         profdata%nvprot(2)    = is3dt
-         profdata%nvprotmpp(1) = it3dtmpp
-         profdata%nvprotmpp(2) = is3dtmpp
+         profdata%nvprot(1)    = ivar1t
+         profdata%nvprot(2)    = ivar2t
+         profdata%nvprotmpp(1) = ivar1tmpp
+         profdata%nvprotmpp(2) = ivar2tmpp
       ENDIF
       profdata%nprof        = iprof
@@ -783 +800 @@
       ! Model level search
       !-----------------------------------------------------------------------
-      IF ( ldt3d ) THEN
+      IF ( ldvar1 ) THEN
          CALL obs_level_search( jpk, gdept_1d, &
             & profdata%nvprot(1), profdata%var(1)%vdep, &
             & profdata%var(1)%mvk )
       ENDIF
-      IF ( lds3d ) THEN
+      IF ( ldvar2 ) THEN
          CALL obs_level_search( jpk, gdept_1d, &
             & profdata%nvprot(2), profdata%var(2)%vdep, &
             & profdata%var(2)%mvk )
       ENDIF
-      
+
       !-----------------------------------------------------------------------
       ! Set model equivalent to 99999
@@ -805 +822 @@
       ! Deallocate temporary data
       !-----------------------------------------------------------------------
-      DEALLOCATE( ifileidx, iprofidx, zdat )
+      DEALLOCATE( ifileidx, iprofidx, zdat, clvars )
 
       !-----------------------------------------------------------------------

branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90

@@ -9 +9 @@
    !!----------------------------------------------------------------------
 
-   !! * Modules used   
+   !! * Modules used
    USE par_kind                 ! Precision variables
    USE in_out_manager           ! I/O manager
@@ -20 +20 @@
    USE obs_surf_def             ! Surface observation definitions
    USE obs_types                ! Observation type definitions
+   USE obs_fbm                  ! Feedback routines
    USE netcdf                   ! NetCDF library
 
@@ -28 +29 @@
 
    PUBLIC obs_rea_surf      ! Read the surface observations from the point data
-   
+
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -37 +38 @@
 CONTAINS
 
-   SUBROUTINE obs_rea_surf( surfdata, knumfiles, cfilenames, &
+   SUBROUTINE obs_rea_surf( surfdata, knumfiles, cdfilenames, &
       &                     kvars, kextr, kstp, ddobsini, ddobsend, &
-      &                     ldignmis, ldmod )
+      &                     ldignmis, ldmod, ldnightav )
       !!---------------------------------------------------------------------
       !!
@@ -59 +60 @@
 
       !! * Arguments
-      TYPE(obs_surf), INTENT(INOUT) :: surfdata   ! Surface data to be read
-      INTEGER, INTENT(IN) :: knumfiles    ! Number of corio format files to read in
-      CHARACTER(LEN=128), INTENT(IN) :: cfilenames(knumfiles) ! File names to read in
+      TYPE(obs_surf), INTENT(INOUT) :: &
+         & surfdata                     ! Surface data to be read
+      INTEGER, INTENT(IN) :: knumfiles  ! Number of corio format files to read
+      CHARACTER(LEN=128), INTENT(IN) :: &
+         & cdfilenames(knumfiles)       ! File names to read in
       INTEGER, INTENT(IN) :: kvars      ! Number of variables in surfdata
-      INTEGER, INTENT(IN) :: kextr      ! Number of extra fields for each var in surfdata
+      INTEGER, INTENT(IN) :: kextr      ! Number of extra fields for each var
       INTEGER, INTENT(IN) :: kstp       ! Ocean time-step index
       LOGICAL, INTENT(IN) :: ldignmis   ! Ignore missing files
       LOGICAL, INTENT(IN) :: ldmod      ! Initialize model from input data
-      REAL(KIND=dp), INTENT(IN) :: ddobsini   ! Obs. ini time in YYYYMMDD.HHMMSS
-      REAL(KIND=dp), INTENT(IN) :: ddobsend   ! Obs. end time in YYYYMMDD.HHMMSS
-         
+      LOGICAL, INTENT(IN) :: ldnightav  ! Observations represent a night-time average
+      REAL(dp), INTENT(IN) :: ddobsini   ! Obs. ini time in YYYYMMDD.HHMMSS
+      REAL(dp), INTENT(IN) :: ddobsend   ! Obs. end time in YYYYMMDD.HHMMSS
+
       !! * Local declarations
       CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_surf'
+      CHARACTER(len=8) :: clrefdate
+      CHARACTER(len=6), DIMENSION(:), ALLOCATABLE :: clvars
       INTEGER :: ji
       INTEGER :: jj
@@ -85 +91 @@
       INTEGER :: imin
       INTEGER :: isec
+      INTEGER :: itype
+      INTEGER :: iobsmpp
+      INTEGER :: iobs
+      INTEGER :: iobstot
+      INTEGER :: ios
+      INTEGER :: ioserrcount
+      INTEGER, PARAMETER :: jpsurfmaxtype = 1024
       INTEGER, DIMENSION(knumfiles) :: irefdate
-      INTEGER :: iobsmpp
-      INTEGER, PARAMETER :: isurfmaxtype = 1024
-      INTEGER, DIMENSION(0:isurfmaxtype) :: &
+      INTEGER, DIMENSION(jpsurfmaxtype+1) :: &
          & ityp, &
          & itypmpp
@@ -98 +109 @@
          & ifileidx, &
          & isurfidx
-      INTEGER :: itype
       REAL(wp), DIMENSION(:), ALLOCATABLE :: &
          & zphi, &
          & zlam
-      real(wp), DIMENSION(:), ALLOCATABLE :: &
+      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
          & zdat
+      REAL(wp), DIMENSION(knumfiles) :: &
+         & djulini, &
+         & djulend
       LOGICAL :: llvalprof
       TYPE(obfbdata), POINTER, DIMENSION(:) :: &
          & inpfiles
-      real(wp), DIMENSION(knumfiles) :: &
-         & djulini, &
-         & djulend
-      INTEGER :: iobs
-      INTEGER :: iobstot
-      INTEGER :: ios
-      INTEGER :: ioserrcount
-      CHARACTER(len=8) :: cl_refdate
-   
+
       ! Local initialization
       iobs = 0
- 
-      !-----------------------------------------------------------------------
-      ! Check data the model part is just with feedback data files
-      !-----------------------------------------------------------------------
-      IF ( ldmod .AND. ( kformat /= 0 ) ) THEN
-         CALL ctl_stop( 'Model can only be read from feedback data' )
-         RETURN
-      ENDIF
 
       !-----------------------------------------------------------------------
       ! Count the number of files needed and allocate the obfbdata type
       !-----------------------------------------------------------------------
-      
+
       inobf = knumfiles
-      
+
       ALLOCATE( inpfiles(inobf) )
 
       surf_files : DO jj = 1, inobf
 
-         CALL init_obfbdata( inpfiles(jj) )
-          
          !---------------------------------------------------------------------
          ! Prints
@@ -145 +140 @@
             WRITE(numout,*)
             WRITE(numout,*) ' obs_rea_surf : Reading from file = ', &
-               & TRIM( TRIM( cfilenames(jj) ) )
+               & TRIM( TRIM( cdfilenames(jj) ) )
             WRITE(numout,*) ' ~~~~~~~~~~~'
             WRITE(numout,*)
@@ -153 +148 @@
          !  Initialization: Open file and get dimensions only
          !---------------------------------------------------------------------
-         
-         iflag = nf90_open( TRIM( TRIM( cfilenames(jj) ) ), nf90_nowrite, &
+
+         iflag = nf90_open( TRIM( TRIM( cdfilenames(jj) ) ), nf90_nowrite, &
             &                      i_file_id )
-         
+
          IF ( iflag /= nf90_noerr ) THEN
 
             IF ( ldignmis ) THEN
                inpfiles(jj)%nobs = 0
-               CALL ctl_warn( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
+               CALL ctl_warn( 'File ' // TRIM( TRIM( cdfilenames(jj) ) ) // &
                   &           ' not found' )
             ELSE 
-               CALL ctl_stop( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
+               CALL ctl_stop( 'File ' // TRIM( TRIM( cdfilenames(jj) ) ) // &
                   &           ' not found' )
             ENDIF
 
          ELSE 
-            
-            !------------------------------------------------------------------
-            !  Close the file since it is opened in read_proffile
-            !------------------------------------------------------------------
-            
+
+            !------------------------------------------------------------------
+            !  Close the file since it is opened in read_obfbdata
+            !------------------------------------------------------------------
+
             iflag = nf90_close( i_file_id )
 
@@ -179 +174 @@
             !  Read the profile file into inpfiles
             !------------------------------------------------------------------
-            IF(lwp) THEN
-               WRITE(numout,*)
-               WRITE(numout,*)'Reading from feedback file :', &
-                  &           TRIM( cfilenames(jj) )
-            ENDIF
-            CALL read_obfbdata( TRIM( cfilenames(jj) ), inpfiles(jj), &
+            CALL init_obfbdata( inpfiles(jj) )
+            CALL read_obfbdata( TRIM( cdfilenames(jj) ), inpfiles(jj), &
                &                ldgrid = .TRUE. )
+
             IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN
                CALL ctl_stop( 'Model not in input data' )
@@ -191 +183 @@
             ENDIF
 
+            IF ( jj == 1 ) THEN
+               ALLOCATE( clvars( inpfiles(jj)%nvar ) )
+               DO ji = 1, inpfiles(jj)%nvar
+                 clvars(ji) = inpfiles(jj)%cname(ji)
+               END DO
+            ELSE
+               DO ji = 1, inpfiles(jj)%nvar
+                  IF ( inpfiles(jj)%cname(ji) /= clvars(ji) ) THEN
+                     CALL ctl_stop( 'Feedback file variables not consistent', &
+                        &           ' with previous files for this type' )
+                  ENDIF
+               END DO
+            ENDIF
+
             !------------------------------------------------------------------
             !  Change longitude (-180,180)
             !------------------------------------------------------------------
 
-            DO ji = 1, inpfiles(jj)%nobs 
+            DO ji = 1, inpfiles(jj)%nobs
 
                IF ( inpfiles(jj)%plam(ji) < -180. ) &
@@ -208 +214 @@
             !  Calculate the date  (change eventually)
             !------------------------------------------------------------------
-            cl_refdate=inpfiles(jj)%cdjuldref(1:8)
-            READ(cl_refdate,'(I8)') irefdate(jj)
-            
+            clrefdate=inpfiles(jj)%cdjuldref(1:8)
+            READ(clrefdate,'(I8)') irefdate(jj)
+
             CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec )
             CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), &
@@ -217 +223 @@
             CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulend(jj), &
                &           krefdate = irefdate(jj) )
+
+            IF ( ldnightav ) THEN
+
+               IF ( lwp ) THEN
+                  WRITE(numout,*)'Resetting time of night-time averaged observations', &
+                     &             ' to the end of the day'
+               ENDIF
+
+               DO ji = 1, inpfiles(jj)%nobs
+                  !  for night-time averaged data force the time
+                  !  to be the last time-step of the day, but still within the day.
+                  IF ( inpfiles(jj)%ptim(ji) >= 0. ) THEN
+                     inpfiles(jj)%ptim(ji) = &
+                        & INT(inpfiles(jj)%ptim(ji)) + 0.9999
+                  ELSE
+                     inpfiles(jj)%ptim(ji) = &
+                        & INT(inpfiles(jj)%ptim(ji)) - 0.0001
+                  ENDIF
+               END DO
+            ENDIF
+
             IF ( inpfiles(jj)%nobs > 0 ) THEN
                inpfiles(jj)%iproc = -1
@@ -312 +339 @@
          &               zdat,     &
          &               iindx   )
-      
+
       CALL obs_surf_alloc( surfdata, iobs, kvars, kextr, kstp, jpi, jpj )
-      
-      ! * Read obs/positions, QC, all variable and assign to surfdata
- 
+
+      ! Read obs/positions, QC, all variable and assign to surfdata
+
       iobs = 0
+
+      surfdata%cvars(:)  = clvars(:)
 
       ityp   (:) = 0
       itypmpp(:) = 0
-      
+
       ioserrcount = 0
-      
+
       DO jk = 1, iobstot
-         
+
          jj = ifileidx(iindx(jk))
          ji = isurfidx(iindx(jk))
          IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
             & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
-            
+
             IF ( nproc == 0 ) THEN
                IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
@@ -336 +365 @@
                IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
             ENDIF
-            
+
             ! Set observation information
-            
+
             IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. &
                & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN
@@ -360 +389 @@
                surfdata%nhou(iobs) = ihou
                surfdata%nmin(iobs) = imin
-               
+
                ! Surface space coordinates
                surfdata%rlam(iobs) = inpfiles(jj)%plam(ji)
@@ -368 +397 @@
                surfdata%mi  (iobs) = inpfiles(jj)%iobsi(ji,1)
                surfdata%mj  (iobs) = inpfiles(jj)%iobsj(ji,1)
-               
+
                ! Instrument type
                READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype
 901            IF ( ios /= 0 ) THEN
-                  IF (ioserrcount == 0) CALL ctl_warn ( 'Problem converting an instrument type to integer. Setting type to zero' ) 
+                  IF (ioserrcount == 0) THEN
+                     CALL ctl_warn ( 'Problem converting an instrument type ', &
+                        &            'to integer. Setting type to zero' )
+                  ENDIF
                   ioserrcount = ioserrcount + 1
                   itype = 0
                ENDIF
                surfdata%ntyp(iobs) = itype
-               IF ( itype < isurfmaxtype + 1 ) THEN
+               IF ( itype < jpsurfmaxtype + 1 ) THEN
                   ityp(itype+1) = ityp(itype+1) + 1
                ELSE
-                  IF(lwp)WRITE(numout,*)'WARNING:Increase isurfmaxtype in ',&
+                  IF(lwp)WRITE(numout,*)'WARNING:Increase jpsurfmaxtype in ',&
                      &                  cpname
                ENDIF
@@ -398 +430 @@
                IF ( ldmod ) THEN
                   surfdata%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1)
-               ELSE
+                  IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN
+                     surfdata%rext(iobs,1) = inpfiles(jj)%padd(1,ji,2,1)
+                     surfdata%rext(iobs,2) = inpfiles(jj)%pext(1,ji,1)
+                  ENDIF
+                ELSE
                   surfdata%rmod(iobs,1) = fbrmdi
+                  IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) surfdata%rext(iobs,:) = fbrmdi
                ENDIF
             ENDIF
@@ -409 +446 @@
       ! Sum up over processors
       !-----------------------------------------------------------------------
-      
+
       CALL obs_mpp_sum_integer( iobs, iobsmpp )
-      
+      CALL obs_mpp_sum_integers( ityp, itypmpp, jpsurfmaxtype + 1 )
+
       !-----------------------------------------------------------------------
       ! Output number of observations.
@@ -418 +456 @@
 
          WRITE(numout,*)
-         WRITE(numout,'(1X,A)')'Surface data types'
+         WRITE(numout,'(1X,A)')TRIM( surfdata%cvars(1) )//' data'
          WRITE(numout,'(1X,A)')'--------------'
          DO jj = 1,8
@@ -425 +463 @@
             ENDIF
          END DO
-         WRITE(numout,'(1X,A50)')'--------------------------------------------------'
-         WRITE(numout,'(1X,A40,I10)')'Total                                 = ',iobsmpp
+         WRITE(numout,'(1X,A)') &
+            & '---------------------------------------------------------------'
+         WRITE(numout,'(1X,A,I8)') &
+            & 'Total data for variable '//TRIM( surfdata%cvars(1) )// &
+            & '           = ', iobsmpp
+         WRITE(numout,'(1X,A)') &
+            & '---------------------------------------------------------------'
          WRITE(numout,*)
 
@@ -434 +477 @@
       ! Deallocate temporary data
       !-----------------------------------------------------------------------
-      DEALLOCATE( ifileidx, isurfidx, zdat )
+      DEALLOCATE( ifileidx, isurfidx, zdat, clvars )
 
       !-----------------------------------------------------------------------

branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90

@@ -31 +31 @@
    PRIVATE
    
-   PUBLIC   obs_rea_mdt     ! called by ?
-   PUBLIC   obs_offset_mdt  ! called by ?
-
-   INTEGER , PUBLIC ::   nmsshc    = 1         ! MDT correction scheme
-   REAL(wp), PUBLIC ::   mdtcorr   = 1.61_wp   ! User specified MDT correction
-   REAL(wp), PUBLIC ::   mdtcutoff = 65.0_wp   ! MDT cutoff for computed correction
+   PUBLIC   obs_rea_mdt     ! called by dia_obs_init
+   PUBLIC   obs_offset_mdt  ! called by obs_rea_mdt
+
+   INTEGER , PUBLIC :: nn_msshc    = 1         ! MDT correction scheme
+   REAL(wp), PUBLIC :: rn_mdtcorr   = 1.61_wp  ! User specified MDT correction
+   REAL(wp), PUBLIC :: rn_mdtcutoff = 65.0_wp  ! MDT cutoff for computed correction
 
    !!----------------------------------------------------------------------
@@ -45 +45 @@
 CONTAINS
 
-   SUBROUTINE obs_rea_mdt( kslano, sladata, k2dint )
+   SUBROUTINE obs_rea_mdt( sladata, k2dint )
       !!---------------------------------------------------------------------
       !!
@@ -58 +58 @@
       USE iom
       !
-      INTEGER                          , INTENT(IN)    ::   kslano    ! Number of SLA Products
-      TYPE(obs_surf), DIMENSION(kslano), INTENT(inout) ::   sladata   ! SLA data
-      INTEGER                          , INTENT(in)    ::   k2dint    ! ?
+      TYPE(obs_surf), INTENT(inout) ::   sladata   ! SLA data
+      INTEGER       , INTENT(in)    ::   k2dint    ! ?
       !
       CHARACTER(LEN=12), PARAMETER ::   cpname  = 'obs_rea_mdt'
       CHARACTER(LEN=20), PARAMETER ::   mdtname = 'slaReferenceLevel.nc'
 
-      INTEGER ::   jslano              ! Data set loop variable
       INTEGER ::   jobs                ! Obs loop variable
       INTEGER ::   jpimdt, jpjmdt      ! Number of grid point in lat/lon for the MDT
@@ -88 +86 @@
       IF(lwp)WRITE(numout,*) ' obs_rea_mdt : Read MDT for referencing altimeter anomalies'
       IF(lwp)WRITE(numout,*) ' ------------- '
+      CALL FLUSH(numout)
 
       CALL iom_open( mdtname, nummdt )       ! Open the file
@@ -109 +108 @@
 
       ! Remove the offset between the MDT used with the sla and the model MDT
-      IF( nmsshc == 1 .OR. nmsshc == 2 )   CALL obs_offset_mdt( z_mdt, zfill )
+      IF( nn_msshc == 1 .OR. nn_msshc == 2 )   CALL obs_offset_mdt( z_mdt, zfill )
 
       ! Intepolate the MDT already on the model grid at the observation point
   
-      DO jslano = 1, kslano
-         ALLOCATE( &
-            & igrdi(2,2,sladata(jslano)%nsurf), &
-            & igrdj(2,2,sladata(jslano)%nsurf), &
-            & zglam(2,2,sladata(jslano)%nsurf), &
-            & zgphi(2,2,sladata(jslano)%nsurf), &
-            & zmask(2,2,sladata(jslano)%nsurf), &
-            & zmdtl(2,2,sladata(jslano)%nsurf)  &
-            & )
+      ALLOCATE( &
+         & igrdi(2,2,sladata%nsurf), &
+         & igrdj(2,2,sladata%nsurf), &
+         & zglam(2,2,sladata%nsurf), &
+         & zgphi(2,2,sladata%nsurf), &
+         & zmask(2,2,sladata%nsurf), &
+         & zmdtl(2,2,sladata%nsurf)  &
+         & )
          
-         DO jobs = 1, sladata(jslano)%nsurf
-
-            igrdi(1,1,jobs) = sladata(jslano)%mi(jobs)-1
-            igrdj(1,1,jobs) = sladata(jslano)%mj(jobs)-1
-            igrdi(1,2,jobs) = sladata(jslano)%mi(jobs)-1
-            igrdj(1,2,jobs) = sladata(jslano)%mj(jobs)
-            igrdi(2,1,jobs) = sladata(jslano)%mi(jobs)
-            igrdj(2,1,jobs) = sladata(jslano)%mj(jobs)-1
-            igrdi(2,2,jobs) = sladata(jslano)%mi(jobs)
-            igrdj(2,2,jobs) = sladata(jslano)%mj(jobs)
-
-         END DO
-
-         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, igrdi, igrdj, glamt  , zglam )
-         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, igrdi, igrdj, gphit  , zgphi )
-         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, igrdi, igrdj, mdtmask, zmask )
-         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, igrdi, igrdj, z_mdt  , zmdtl )
-
-         DO jobs = 1, sladata(jslano)%nsurf
+      DO jobs = 1, sladata%nsurf
+
+         igrdi(1,1,jobs) = sladata%mi(jobs)-1
+         igrdj(1,1,jobs) = sladata%mj(jobs)-1
+         igrdi(1,2,jobs) = sladata%mi(jobs)-1
+         igrdj(1,2,jobs) = sladata%mj(jobs)
+         igrdi(2,1,jobs) = sladata%mi(jobs)
+         igrdj(2,1,jobs) = sladata%mj(jobs)-1
+         igrdi(2,2,jobs) = sladata%mi(jobs)
+         igrdj(2,2,jobs) = sladata%mj(jobs)
+
+      END DO
+
+      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, igrdi, igrdj, glamt  , zglam )
+      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, igrdi, igrdj, gphit  , zgphi )
+      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, igrdi, igrdj, mdtmask, zmask )
+      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, igrdi, igrdj, z_mdt  , zmdtl )
+
+      DO jobs = 1, sladata%nsurf
             
-            zlam = sladata(jslano)%rlam(jobs)
-            zphi = sladata(jslano)%rphi(jobs)
-
-            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
-               &                   zglam(:,:,jobs), zgphi(:,:,jobs), &
-               &                   zmask(:,:,jobs), zweig, zobsmask )
+         zlam = sladata%rlam(jobs)
+         zphi = sladata%rphi(jobs)
+
+         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
+            &                   zglam(:,:,jobs), zgphi(:,:,jobs), &
+            &                   zmask(:,:,jobs), zweig, zobsmask )
             
-            CALL obs_int_h2d( 1, 1, zweig, zmdtl(:,:,jobs),  zext )
+         CALL obs_int_h2d( 1, 1, zweig, zmdtl(:,:,jobs),  zext )
  
-            sladata(jslano)%rext(jobs,2) = zext(1)
+         sladata%rext(jobs,2) = zext(1)
 
 ! mark any masked data with a QC flag
-            IF( zobsmask(1) == 0 )   sladata(jslano)%nqc(jobs) = 11
+         IF( zobsmask(1) == 0 )   sladata%nqc(jobs) = 11
 
          END DO
          
-         DEALLOCATE( &
-            & igrdi, &
-            & igrdj, &
-            & zglam, &
-            & zgphi, &
-            & zmask, &
-            & zmdtl  &
-            & )
-
-      END DO
+      DEALLOCATE( &
+         & igrdi, &
+         & igrdj, &
+         & zglam, &
+         & zgphi, &
+         & zmask, &
+         & zmdtl  &
+         & )
 
       CALL wrk_dealloc(jpi,jpj,z_mdt,mdtmask) 
+      IF(lwp)WRITE(numout,*) ' ------------- '
+      CALL FLUSH(numout)
       !
    END SUBROUTINE obs_rea_mdt
@@ -205 +203 @@
         DO jj = 1, jpj
            zpromsk(ji,jj) = tmask_i(ji,jj)
-           IF (    ( gphit(ji,jj) .GT.  mdtcutoff ) &
-              &.OR.( gphit(ji,jj) .LT. -mdtcutoff ) &
+           IF (    ( gphit(ji,jj) .GT.  rn_mdtcutoff ) &
+              &.OR.( gphit(ji,jj) .LT. -rn_mdtcutoff ) &
               &.OR.( mdt(ji,jj) .EQ. zfill ) ) &
               &        zpromsk(ji,jj) = 0.0
@@ -212 +210 @@
       END DO 
 
-      ! Compute MSSH mean over [0,360] x [-mdtcutoff,mdtcutoff]
+      ! Compute MSSH mean over [0,360] x [-rn_mdtcutoff,rn_mdtcutoff]
 
       zarea = 0.0
@@ -240 +238 @@
       !  Correct spatial mean of the MSSH
 
-      IF( nmsshc == 1 )   mdt(:,:) = mdt(:,:) - zcorr  
+      IF( nn_msshc == 1 )   mdt(:,:) = mdt(:,:) - zcorr  
 
       ! User defined value : 1.6 m for the Rio MDT compared to ORCA2 MDT
 
-      IF( nmsshc == 2 )   mdt(:,:) = mdt(:,:) - mdtcorr
+      IF( nn_msshc == 2 )   mdt(:,:) = mdt(:,:) - rn_mdtcorr
 
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) ' obs_readmdt : mdtcutoff     = ', mdtcutoff
+         WRITE(numout,*) ' obs_readmdt : rn_mdtcutoff     = ', rn_mdtcutoff
          WRITE(numout,*) ' -----------   zcorr_mdt     = ', zcorr_mdt
          WRITE(numout,*) '               zcorr_bcketa  = ', zcorr_bcketa
          WRITE(numout,*) '               zcorr         = ', zcorr
-         WRITE(numout,*) '               nmsshc        = ', nmsshc
+         WRITE(numout,*) '               nn_msshc        = ', nn_msshc
       ENDIF
 
-      IF ( nmsshc == 0 ) WRITE(numout,*) '           MSSH correction is not applied'
-      IF ( nmsshc == 1 ) WRITE(numout,*) '           MSSH correction is applied'
-      IF ( nmsshc == 2 ) WRITE(numout,*) '           User defined MSSH correction' 
+      IF ( nn_msshc == 0 ) WRITE(numout,*) '           MSSH correction is not applied'
+      IF ( nn_msshc == 1 ) WRITE(numout,*) '           MSSH correction is applied'
+      IF ( nn_msshc == 2 ) WRITE(numout,*) '           User defined MSSH correction' 
 
       CALL wrk_dealloc( jpi,jpj, zpromsk )

branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90

@@ -67 +67 @@
          & ntyp           !: Type of surface observation product
 
+      CHARACTER(len=6), POINTER, DIMENSION(:) :: &
+         & cvars          !: Variable names
+
       CHARACTER(LEN=8), POINTER, DIMENSION(:) :: &
          & cwmo           !: WMO indentifier
@@ -130 +133 @@
       !!* Local variables
       INTEGER :: ji
+      INTEGER :: jvar
 
       ! Set bookkeeping variables
@@ -140 +144 @@
       surf%npi      = kpi
       surf%npj      = kpj
+
+      ! Allocate arrays of size number of variables
+
+      ALLOCATE( &
+         & surf%cvars(kvar)    &
+         & )
+
+      DO jvar = 1, kvar
+         surf%cvars(jvar) = "NotSet"
+      END DO
       
       ! Allocate arrays of number of surface data size
@@ -271 +285 @@
          & )
 
+      ! Dellocate arrays of size number of variables
+
+      DEALLOCATE( &
+         & surf%cvars     &
+         & )
+
    END SUBROUTINE obs_surf_dealloc
 
@@ -392 +412 @@
       ! Set book keeping variables which do not depend on number of obs.
 
-      newsurf%nstp  = surf%nstp
+      newsurf%nstp     = surf%nstp
+      newsurf%cvars(:) = surf%cvars(:)
  
       ! Deallocate temporary data

branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_types.F90

@@ -117 +117 @@
 
          cwmonam1770(ji) = 'Not defined'
-         ctypshort(ji) = 'XBT'
+         ctypshort(ji) = '---'
 
 !         IF ( ji < 1000 ) THEN

branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90

@@ -55 +55 @@
 CONTAINS
 
-   SUBROUTINE obs_wri_prof( cobstype, profdata, padd, pext )
+   SUBROUTINE obs_wri_prof( profdata, k2dint, padd, pext )
       !!-----------------------------------------------------------------------
       !!
@@ -76 +76 @@
       !!-----------------------------------------------------------------------
 
-      !! * Modules used
-
       !! * Arguments
-      CHARACTER(LEN=*), INTENT(IN) :: cobstype        ! Prefix for output files
       TYPE(obs_prof), INTENT(INOUT) :: profdata      ! Full set of profile data
+      INTEGER, INTENT(IN)        :: k2dint           ! Horizontal interpolation method
       TYPE(obswriinfo), OPTIONAL :: padd             ! Additional info for each variable
       TYPE(obswriinfo), OPTIONAL :: pext             ! Extra info
-      
+
       !! * Local declarations
       TYPE(obfbdata) :: fbdata
-      CHARACTER(LEN=40) :: cfname
+      CHARACTER(LEN=40) :: clfname
+      CHARACTER(LEN=6) :: clfiletype
       INTEGER :: ilevel
       INTEGER :: jvar
@@ -94 +93 @@
       INTEGER :: ja
       INTEGER :: je
+      INTEGER :: iadd
+      INTEGER :: iext
       REAL(wp) :: zpres
-      INTEGER :: nadd
-      INTEGER :: next
+      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
+         & zu, &
+         & zv
 
       IF ( PRESENT( padd ) ) THEN
-         nadd = padd%inum
+         iadd = padd%inum
       ELSE
-         nadd = 0
+         iadd = 0
       ENDIF
 
       IF ( PRESENT( pext ) ) THEN
-         next = pext%inum
+         iext = pext%inum
       ELSE
-         next = 0
-      ENDIF
-      
+         iext = 0
+      ENDIF
+
       CALL init_obfbdata( fbdata )
 
@@ -118 +120 @@
       END DO
 
-      SELECT CASE ( TRIM(cobstype) )
-      CASE('prof')
-
+      SELECT CASE ( TRIM(profdata%cvars(1)) )
+      CASE('POTM')
+
+         clfiletype='profb'
          CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, &
-            &                 1 + nadd, 1 + next, .TRUE. )
-         fbdata%cname(1)      = 'POTM'
-         fbdata%cname(2)      = 'PSAL'
+            &                 1 + iadd, 1 + iext, .TRUE. )
+         fbdata%cname(1)      = profdata%cvars(1)
+         fbdata%cname(2)      = profdata%cvars(2)
          fbdata%coblong(1)    = 'Potential temperature'
          fbdata%coblong(2)    = 'Practical salinity'
@@ -137 +140 @@
          fbdata%caddunit(1,2) = 'PSU'
          fbdata%cgrid(:)      = 'T'
-         DO je = 1, next
+         DO je = 1, iext
             fbdata%cextname(1+je) = pext%cdname(je)
             fbdata%cextlong(1+je) = pext%cdlong(je,1)
             fbdata%cextunit(1+je) = pext%cdunit(je,1)
          END DO
-         DO ja = 1, nadd
+         DO ja = 1, iadd
             fbdata%caddname(1+ja) = padd%cdname(ja)
             DO jvar = 1, 2
@@ -149 +152 @@
             END DO
          END DO
-      
-      CASE('vel')
-      
+
+      CASE('UVEL')
+
+         clfiletype='velfb'
          CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, 2, 0, .TRUE. )
-         fbdata%cname(1)      = 'UVEL'
-         fbdata%cname(2)      = 'VVEL'
+         fbdata%cname(1)      = profdata%cvars(1)
+         fbdata%cname(2)      = profdata%cvars(2)
          fbdata%coblong(1)    = 'Zonal velocity'
          fbdata%coblong(2)    = 'Meridional velocity'
          fbdata%cobunit(1)    = 'm/s'
          fbdata%cobunit(2)    = 'm/s'
-         DO je = 1, next
+         DO je = 1, iext
             fbdata%cextname(je) = pext%cdname(je)
             fbdata%cextlong(je) = pext%cdlong(je,1)
@@ -175 +179 @@
          fbdata%cgrid(1)      = 'U' 
          fbdata%cgrid(2)      = 'V'
-         DO ja = 1, nadd
+         DO ja = 1, iadd
             fbdata%caddname(2+ja) = padd%cdname(ja)
             fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1)
@@ -187 +191 @@
 
       END SELECT
-         
+
       fbdata%caddname(1)   = 'Hx'
-         
-      WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cobstype), nproc
+
+      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
 
       IF(lwp) THEN
@@ -196 +200 @@
          WRITE(numout,*)'obs_wri_prof :'
          WRITE(numout,*)'~~~~~~~~~~~~~'
-         WRITE(numout,*)'Writing profile feedback file : ',TRIM(cfname)
+         WRITE(numout,*)'Writing '//TRIM(clfiletype)//' feedback file : ',TRIM(clfname)
       ENDIF
 
@@ -242 +246 @@
             DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar)
                ik = profdata%var(jvar)%nvlidx(jk)
-               IF ( TRIM(cobstype) == 'prof' ) THEN
+               IF ( TRIM(profdata%cvars(1)) == 'POTM' ) THEN
                   fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk)
-               ELSE IF ( TRIM(cobstype) == 'vel' ) THEN
+               ELSE IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN
                   IF ( jvar == 1 ) THEN
                      fbdata%padd(ik,jo,1,jvar) = zu(jk)
@@ -265 +269 @@
                ENDIF
                fbdata%iobsk(ik,jo,jvar)  = profdata%var(jvar)%mvk(jk)
-               DO ja = 1, nadd
+               DO ja = 1, iadd
                   fbdata%padd(ik,jo,1+ja,jvar) = &
                      & profdata%var(jvar)%vext(jk,padd%ipoint(ja))
                END DO
-               DO je = 1, next
+               DO je = 1, iext
                   fbdata%pext(ik,jo,1+je) = &
                      & profdata%var(jvar)%vext(jk,pext%ipoint(je))
@@ -280 +284 @@
       END DO
 
-      IF ( TRIM(cobstype) == 'prof' ) THEN
+      IF ( TRIM(profdata%cvars(1)) == 'POTM' ) THEN
          ! Convert insitu temperature to potential temperature using the model
          ! salinity if no potential temperature
@@ -301 +305 @@
          END DO
       ENDIF
-      
+
       ! Write the obfbdata structure
-      CALL write_obfbdata( cfname, fbdata )
+      CALL write_obfbdata( clfname, fbdata )
 
       ! Output some basic statistics
@@ -309 +313 @@
 
       CALL dealloc_obfbdata( fbdata )
-     
+
    END SUBROUTINE obs_wri_prof
 
-   SUBROUTINE obs_wri_surf( cobstype, surfdata, padd, pext )
+   SUBROUTINE obs_wri_surf( surfdata, padd, pext )
       !!-----------------------------------------------------------------------
       !!
@@ -332 +336 @@
 
       !! * Arguments
-      CHARACTER(LEN=*), INTENT(IN) :: cobstype          ! Prefix for output files
       TYPE(obs_surf), INTENT(INOUT) :: surfdata         ! Full set of surface data
       TYPE(obswriinfo), OPTIONAL :: padd               ! Additional info for each variable
@@ -339 +342 @@
       !! * Local declarations
       TYPE(obfbdata) :: fbdata
-      CHARACTER(LEN=40) :: cfname         ! netCDF filename
+      CHARACTER(LEN=40) :: clfname         ! netCDF filename
+      CHARACTER(LEN=6)  :: clfiletype
       CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf'
       INTEGER :: jo
       INTEGER :: ja
       INTEGER :: je
-      INTEGER :: nadd
-      INTEGER :: next
+      INTEGER :: iadd
+      INTEGER :: iext
 
       IF ( PRESENT( padd ) ) THEN
-         nadd = padd%inum
+         iadd = padd%inum
       ELSE
-         nadd = 0
+         iadd = 0
       ENDIF
 
       IF ( PRESENT( pext ) ) THEN
-         next = pext%inum
+         iext = pext%inum
       ELSE
-         next = 0
+         iext = 0
       ENDIF
 
@@ -362 +366 @@
 
       CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
-         &                 2 + nadd, 1 + next, .TRUE. )
-
-      SELECT CASE ( TRIM(cobstype) )
-      CASE('sla')
-
-         fbdata%cname(1)      = 'SLA'
+         &                 2 + iadd, 1 + iext, .TRUE. )
+
+      SELECT CASE ( TRIM(surfdata%cvars(1)) )
+      CASE('SLA')
+
+         clfiletype = 'slafb'
+         fbdata%cname(1)      = surfdata%cvars(1)
          fbdata%coblong(1)    = 'Sea level anomaly'
          fbdata%cobunit(1)    = 'Metres'
@@ -373 +378 @@
          fbdata%cextlong(1)   = 'Mean dynamic topography'
          fbdata%cextunit(1)   = 'Metres'
-         DO je = 1, next
+         DO je = 1, iext
             fbdata%cextname(je) = pext%cdname(je)
             fbdata%cextlong(je) = pext%cdlong(je,1)
@@ -384 +389 @@
          fbdata%caddunit(2,1) = 'Metres'
          fbdata%cgrid(1)      = 'T'
-         DO ja = 1, nadd
+         DO ja = 1, iadd
             fbdata%caddname(2+ja) = padd%cdname(ja)
             fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1)
@@ -390 +395 @@
          END DO
 
-      CASE('sst')
-
-         fbdata%cname(1)      = 'SST'
+      CASE('SST')
+
+         clfiletype = 'sstfb'
+         fbdata%cname(1)      = surfdata%cvars(1)
          fbdata%coblong(1)    = 'Sea surface temperature'
          fbdata%cobunit(1)    = 'Degree centigrade'
-         DO je = 1, next
+         DO je = 1, iext
             fbdata%cextname(je) = pext%cdname(je)
             fbdata%cextlong(je) = pext%cdlong(je,1)
@@ -403 +409 @@
          fbdata%caddunit(1,1) = 'Degree centigrade'
          fbdata%cgrid(1)      = 'T'
-         DO ja = 1, nadd
+         DO ja = 1, iadd
             fbdata%caddname(1+ja) = padd%cdname(ja)
             fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
@@ -409 +415 @@
          END DO
 
-      CASE('sss')
-
-         fbdata%cname(1)      = 'SSS'
-         fbdata%coblong(1)    = 'Sea surface salinity'
-         fbdata%cobunit(1)    = 'PSU'
-         DO je = 1, next
+      CASE('SEAICE')
+
+         clfiletype = 'sicfb'
+         fbdata%cname(1)      = surfdata%cvars(1)
+         fbdata%coblong(1)    = 'Sea ice'
+         fbdata%cobunit(1)    = 'Fraction'
+         DO je = 1, iext
             fbdata%cextname(je) = pext%cdname(je)
             fbdata%cextlong(je) = pext%cdlong(je,1)
             fbdata%cextunit(je) = pext%cdunit(je,1)
          END DO
-         fbdata%caddlong(1,1) = 'Model interpolated SSS'
-         fbdata%caddunit(1,1) = 'PSU'
+         fbdata%caddlong(1,1) = 'Model interpolated ICE'
+         fbdata%caddunit(1,1) = 'Fraction'
          fbdata%cgrid(1)      = 'T'
-         DO ja = 1, nadd
+         DO ja = 1, iadd
             fbdata%caddname(1+ja) = padd%cdname(ja)
             fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
@@ -428 +435 @@
          END DO
 
-      CASE('seaice')
-
-         fbdata%cname(1)      = 'SEAICE'
-         fbdata%coblong(1)    = 'Sea ice'
-         fbdata%cobunit(1)    = 'Fraction'
-         DO je = 1, next
-            fbdata%cextname(je) = pext%cdname(je)
-            fbdata%cextlong(je) = pext%cdlong(je,1)
-            fbdata%cextunit(je) = pext%cdunit(je,1)
-         END DO
-         fbdata%caddlong(1,1) = 'Model interpolated ICE'
-         fbdata%caddunit(1,1) = 'Fraction'
-         fbdata%cgrid(1)      = 'T'
-         DO ja = 1, nadd
-            fbdata%caddname(1+ja) = padd%cdname(ja)
-            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
-            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
-         END DO
-
       END SELECT
 
       fbdata%caddname(1)   = 'Hx'
 
-      WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cobstype), nproc
+      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
 
       IF(lwp) THEN
@@ -457 +445 @@
          WRITE(numout,*)'obs_wri_surf :'
          WRITE(numout,*)'~~~~~~~~~~~~~'
-         WRITE(numout,*)'Writing surface feedback file : ',TRIM(cfname)
+         WRITE(numout,*)'Writing '//TRIM(surfdata%cvars(1))//' feedback file : ',TRIM(clfname)
       ENDIF
 
@@ -498 +486 @@
             &           krefdate = 19500101 )
          fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1)
-         IF ( TRIM(cobstype) == 'sla' ) fbdata%padd(1,jo,2,1) = surfdata%rext(jo,1)
+         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%padd(1,jo,2,1) = surfdata%rext(jo,1)
          fbdata%pob(1,jo,1)    = surfdata%robs(jo,1) 
          fbdata%pdep(1,jo)     = 0.0
@@ -514 +502 @@
          ENDIF
          fbdata%iobsk(1,jo,1)  = 0
-         IF ( TRIM(cobstype) == 'sla' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)
-         DO ja = 1, nadd
+         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)
+         DO ja = 1, iadd
             fbdata%padd(1,jo,2+ja,1) = &
                & surfdata%rext(jo,padd%ipoint(ja))
          END DO
-         DO je = 1, next
+         DO je = 1, iext
             fbdata%pext(1,jo,1+je) = &
                & surfdata%rext(jo,pext%ipoint(je))
@@ -526 +514 @@
 
       ! Write the obfbdata structure
-      CALL write_obfbdata( cfname, fbdata )
+      CALL write_obfbdata( clfname, fbdata )
 
       ! Output some basic statistics
@@ -556 +544 @@
       INTEGER :: jo
       INTEGER :: jk
-
-      INTEGER :: numgoodobs
-      INTEGER :: numgoodobsmpp
+      INTEGER :: inumgoodobs
+      INTEGER :: inumgoodobsmpp
       REAL(wp) :: zsumx
       REAL(wp) :: zsumx2
@@ -566 +553 @@
          WRITE(numout,*) ''
          WRITE(numout,*) 'obs_wri_stats :'
-         WRITE(numout,*) '~~~~~~~~~~~~~~~' 
+         WRITE(numout,*) '~~~~~~~~~~~~~~~'
       ENDIF
 
@@ -572 +559 @@
          zsumx=0.0_wp
          zsumx2=0.0_wp
-         numgoodobs=0
+         inumgoodobs=0
          DO jo = 1, fbdata%nobs
             DO jk = 1, fbdata%nlev
@@ -578 +565 @@
                   & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
                   & ( fbdata%padd(jk,jo,1,jvar) < 9999.0 ) ) THEN
-       
-             zomb=fbdata%pob(jk, jo, jvar)-fbdata%padd(jk, jo, 1, jvar)
+
+                  zomb=fbdata%pob(jk, jo, jvar)-fbdata%padd(jk, jo, 1, jvar)
                   zsumx=zsumx+zomb
                   zsumx2=zsumx2+zomb**2
-                  numgoodobs=numgoodobs+1
-          ENDIF
+                  inumgoodobs=inumgoodobs+1
+               ENDIF
             ENDDO
          ENDDO
 
-         CALL obs_mpp_sum_integer( numgoodobs, numgoodobsmpp )
+         CALL obs_mpp_sum_integer( inumgoodobs, inumgoodobsmpp )
          CALL mpp_sum(zsumx)
          CALL mpp_sum(zsumx2)
 
          IF (lwp) THEN
-       WRITE(numout,*) 'Type: ',fbdata%cname(jvar),'  Total number of good observations: ',numgoodobsmpp 
-       WRITE(numout,*) 'Overall mean obs minus model of the good observations: ',zsumx/numgoodobsmpp
-            WRITE(numout,*) 'Overall RMS obs minus model of the good observations: ',sqrt( zsumx2/numgoodobsmpp )
-       WRITE(numout,*) ''
+            WRITE(numout,*) 'Type: ',fbdata%cname(jvar),'  Total number of good observations: ',inumgoodobsmpp 
+            WRITE(numout,*) 'Overall mean obs minus model of the good observations: ',zsumx/inumgoodobsmpp
+            WRITE(numout,*) 'Overall RMS obs minus model of the good observations: ',sqrt( zsumx2/inumgoodobsmpp )
+            WRITE(numout,*) ''
          ENDIF
- 
+
       ENDDO