New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 13758 – NEMO

Changeset 13758


Ignore:
Timestamp:
2020-11-09T17:36:09+01:00 (3 years ago)
Author:
dford
Message:

Initial implementation of changes following generic aspects of changesets 9306 and 11546.

Location:
NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/cfgs/SHARED/namelist_ref

    r13731 r13758  
    13451345   ln_sstnight = .false.             ! Logical switch for calculating night-time average for SST obs 
    13461346   ln_bound_reject  = .false.        ! Logical to remove obs near boundaries in LAMs. 
     1347   ln_default_fp_indegs = .true.     ! Logical: T=> averaging footprint is in degrees, F=> in metres 
    13471348   ln_sla_fp_indegs = .true.         ! Logical for SLA: T=> averaging footprint is in degrees, F=> in metres 
    13481349   ln_sst_fp_indegs = .true.         ! Logical for SST: T=> averaging footprint is in degrees, F=> in metres 
     
    13601361   cn_gridsearchfile ='gridsearch.nc' ! Grid search file name 
    13611362   rn_gridsearchres = 0.5            ! Grid search resolution 
     1363   rn_default_avglamscl = 0.         ! Default E/W diameter of observation footprint (metres/degrees) 
     1364   rn_default_avgphiscl = 0.         ! Default N/S diameter of observation footprint (metres/degrees) 
    13621365   rn_mdtcorr  = 1.61                ! MDT  correction 
    13631366   rn_mdtcutoff = 65.0               ! MDT cutoff for computed correction 
     
    13731376   rn_sic_avgphiscl = 0.             ! N/S diameter of SIC observation footprint (metres/degrees) 
    13741377   nn_1dint = 0                      ! Type of vertical interpolation method 
    1375    nn_2dint = 0                      ! Default horizontal interpolation method 
     1378   nn_2dint_default = 0              ! Default horizontal interpolation method 
    13761379   nn_2dint_sla = 0                  ! Horizontal interpolation method for SLA 
    13771380   nn_2dint_sst = 0                  ! Horizontal interpolation method for SST 
  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/OBS/diaobs.F90

    r13216 r13758  
    6464   LOGICAL         :: ln_sic_fp_indegs   !  T=> sea-ice obs footprint size specified in degrees, F=> in metres 
    6565 
    66    REAL(wp) ::   rn_sla_avglamscl   ! E/W diameter of SLA observation footprint (metres) 
    67    REAL(wp) ::   rn_sla_avgphiscl   ! N/S diameter of SLA observation footprint (metres) 
    68    REAL(wp) ::   rn_sst_avglamscl   ! E/W diameter of SST observation footprint (metres) 
    69    REAL(wp) ::   rn_sst_avgphiscl   ! N/S diameter of SST observation footprint (metres) 
    70    REAL(wp) ::   rn_sss_avglamscl   ! E/W diameter of SSS observation footprint (metres) 
    71    REAL(wp) ::   rn_sss_avgphiscl   ! N/S diameter of SSS observation footprint (metres) 
    72    REAL(wp) ::   rn_sic_avglamscl   ! E/W diameter of sea-ice observation footprint (metres) 
    73    REAL(wp) ::   rn_sic_avgphiscl   ! N/S diameter of sea-ice observation footprint (metres) 
    74  
    75    INTEGER :: nn_1dint       ! Vertical interpolation method 
    76    INTEGER :: nn_2dint       ! Default horizontal interpolation method 
    77    INTEGER :: nn_2dint_sla   ! SLA horizontal interpolation method  
    78    INTEGER :: nn_2dint_sst   ! SST horizontal interpolation method  
    79    INTEGER :: nn_2dint_sss   ! SSS horizontal interpolation method  
    80    INTEGER :: nn_2dint_sic   ! Seaice horizontal interpolation method  
     66   REAL(wp) ::   rn_default_avglamscl    ! E/W diameter of SLA observation footprint (metres) 
     67   REAL(wp) ::   rn_default_avgphiscl    ! N/S diameter of SLA observation footprint (metre 
     68   REAL(wp) ::   rn_sla_avglamscl        ! E/W diameter of SLA observation footprint (metres) 
     69   REAL(wp) ::   rn_sla_avgphiscl        ! N/S diameter of SLA observation footprint (metres) 
     70   REAL(wp) ::   rn_sst_avglamscl        ! E/W diameter of SST observation footprint (metres) 
     71   REAL(wp) ::   rn_sst_avgphiscl        ! N/S diameter of SST observation footprint (metres) 
     72   REAL(wp) ::   rn_sss_avglamscl        ! E/W diameter of SSS observation footprint (metres) 
     73   REAL(wp) ::   rn_sss_avgphiscl        ! N/S diameter of SSS observation footprint (metres) 
     74   REAL(wp) ::   rn_sic_avglamscl        ! E/W diameter of sea-ice observation footprint (metres) 
     75   REAL(wp) ::   rn_sic_avgphiscl        ! N/S diameter of sea-ice observation footprint (metres) 
     76 
     77   INTEGER :: nn_1dint                   ! Vertical interpolation method 
     78   INTEGER :: nn_2dint_default           ! Default horizontal interpolation method 
     79   INTEGER :: nn_2dint_sla               ! SLA horizontal interpolation method  
     80   INTEGER :: nn_2dint_sst               ! SST horizontal interpolation method  
     81   INTEGER :: nn_2dint_sss               ! SSS horizontal interpolation method  
     82   INTEGER :: nn_2dint_sic               ! Seaice horizontal interpolation method  
    8183   INTEGER, DIMENSION(imaxavtypes) ::   nn_profdavtypes   ! Profile data types representing a daily average 
    8284   INTEGER :: nproftypes     ! Number of profile obs types 
     
    9496   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) ::   profdataqc   !: Profile data after quality control 
    9597 
    96    CHARACTER(len=lca), PUBLIC, DIMENSION(:), ALLOCATABLE ::   cobstypesprof, cobstypessurf   !: Profile & surface obs types 
     98   CHARACTER(len=8), PUBLIC, DIMENSION(:), ALLOCATABLE ::   cobstypesprof, cobstypessurf   !: Profile & surface obs types 
    9799 
    98100   !!---------------------------------------------------------------------- 
     
    121123      INTEGER :: jvar            ! Counter for variables 
    122124      INTEGER :: jfile           ! Counter for files 
    123       INTEGER :: jnumsstbias 
     125      INTEGER :: jnumsstbias     ! Number of SST bias files to read and apply 
     126      INTEGER :: n2dint_type     ! Local version of nn_2dint* 
    124127      ! 
    125128      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 
     
    150153      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs 
    151154      LOGICAL :: ln_bound_reject ! Logical to remove obs near boundaries in LAMs. 
    152       LOGICAL :: llvar1          ! Logical for profile variable 1 
    153       LOGICAL :: llvar2          ! Logical for profile variable 1 
     155      LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs 
     156      LOGICAL :: ltype_night     ! Local version of ln_sstnight (false for other variables) 
     157      LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar   ! Logical for profile variable read 
    154158      LOGICAL, DIMENSION(jpmaxnfiles) :: lmask ! Used for finding number of sstbias files 
    155159      ! 
    156       REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS 
    157       REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS 
    158       REAL(wp), DIMENSION(jpi,jpj)     ::   zglam1, zglam2   ! Model longitudes for profile variable 1 & 2 
    159       REAL(wp), DIMENSION(jpi,jpj)     ::   zgphi1, zgphi2   ! Model latitudes  for profile variable 1 & 2 
    160       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmask1, zmask2   ! Model land/sea mask associated with variable 1 & 2 
     160      REAL(dp) :: rn_dobsini      ! Obs window start date YYYYMMDD.HHMMSS 
     161      REAL(dp) :: rn_dobsend      ! Obs window end date   YYYYMMDD.HHMMSS 
     162      REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl 
     163      REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl 
     164      REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: zglam   ! Model longitudes for profile variables 
     165      REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: zgphi   ! Model latitudes  for profile variables 
     166      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zmask   ! Model land/sea mask associated with variables 
    161167      !! 
    162168      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              & 
     
    165171         &            ln_grid_global, ln_grid_search_lookup,          & 
    166172         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          & 
    167          &            ln_sstnight,                                    & 
     173         &            ln_sstnight, ln_default_fp_indegs,              & 
    168174         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
    169175         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             & 
     
    174180         &            cn_gridsearchfile, rn_gridsearchres,            & 
    175181         &            rn_dobsini, rn_dobsend,                         & 
     182         &            rn_default_avglamscl, rn_default_avgphiscl,     & 
    176183         &            rn_sla_avglamscl, rn_sla_avgphiscl,             & 
    177184         &            rn_sst_avglamscl, rn_sst_avgphiscl,             & 
    178185         &            rn_sss_avglamscl, rn_sss_avgphiscl,             & 
    179186         &            rn_sic_avglamscl, rn_sic_avgphiscl,             & 
    180          &            nn_1dint, nn_2dint,                             & 
     187         &            nn_1dint, nn_2dint_default,                     & 
    181188         &            nn_2dint_sla, nn_2dint_sst,                     & 
    182189         &            nn_2dint_sss, nn_2dint_sic,                     & 
     
    234241         WRITE(numout,*) '      Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend 
    235242         WRITE(numout,*) '      Type of vertical interpolation method                  nn_1dint = ', nn_1dint 
    236          WRITE(numout,*) '      Type of horizontal interpolation method                nn_2dint = ', nn_2dint 
     243         WRITE(numout,*) '      Default horizontal interpolation method        nn_2dint_default = ', nn_2dint_default 
     244         WRITE(numout,*) '      Type of horizontal interpolation method for SLA    nn_2dint_sla = ', nn_2dint_sla 
     245         WRITE(numout,*) '      Type of horizontal interpolation method for SST    nn_2dint_sst = ', nn_2dint_sst 
     246         WRITE(numout,*) '      Type of horizontal interpolation method for SSS    nn_2dint_sss = ', nn_2dint_sss         
     247         WRITE(numout,*) '      Type of horizontal interpolation method for SIC    nn_2dint_sic = ', nn_2dint_sic 
     248         WRITE(numout,*) '      Default E/W diameter of obs footprint      rn_default_avglamscl = ', rn_default_avglamscl 
     249         WRITE(numout,*) '      Default N/S diameter of obs footprint      rn_default_avgphiscl = ', rn_default_avgphiscl 
     250         WRITE(numout,*) '      Default obs footprint in deg [T] or m [F]  ln_default_fp_indegs = ', ln_default_fp_indegs 
     251         WRITE(numout,*) '      SLA E/W diameter of obs footprint              rn_sla_avglamscl = ', rn_sla_avglamscl 
     252         WRITE(numout,*) '      SLA N/S diameter of obs footprint              rn_sla_avgphiscl = ', rn_sla_avgphiscl 
     253         WRITE(numout,*) '      SLA obs footprint in deg [T] or m [F]          ln_sla_fp_indegs = ', ln_sla_fp_indegs 
     254         WRITE(numout,*) '      SST E/W diameter of obs footprint              rn_sst_avglamscl = ', rn_sst_avglamscl 
     255         WRITE(numout,*) '      SST N/S diameter of obs footprint              rn_sst_avgphiscl = ', rn_sst_avgphiscl 
     256         WRITE(numout,*) '      SST obs footprint in deg [T] or m [F]          ln_sst_fp_indegs = ', ln_sst_fp_indegs 
     257         WRITE(numout,*) '      SIC E/W diameter of obs footprint              rn_sic_avglamscl = ', rn_sic_avglamscl 
     258         WRITE(numout,*) '      SIC N/S diameter of obs footprint              rn_sic_avgphiscl = ', rn_sic_avgphiscl 
     259         WRITE(numout,*) '      SIC obs footprint in deg [T] or m [F]          ln_sic_fp_indegs = ', ln_sic_fp_indegs 
    237260         WRITE(numout,*) '      Rejection of observations near land switch               ln_nea = ', ln_nea 
    238261         WRITE(numout,*) '      Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject 
     
    278301         IF( ln_t3d .OR. ln_s3d ) THEN 
    279302            jtype = jtype + 1 
    280             CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof  ', & 
    281                &                   cn_profbfiles, ifilesprof, cobstypesprof, clproffiles ) 
     303            cobstypesprof(jtype) = 'prof' 
     304            clproffiles(jtype,:) = cn_profbfiles 
    282305         ENDIF 
    283306         IF( ln_vel3d ) THEN 
    284307            jtype = jtype + 1 
    285             CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel   ', & 
    286                &                   cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles ) 
     308            cobstypesprof(jtype) = 'vel' 
     309            clproffiles(jtype,:) = cn_velfbfiles 
    287310         ENDIF 
     311         ! 
     312         CALL obs_settypefiles( nproftypes, jpmaxnfiles, ifilesprof, cobstypesprof, clproffiles ) 
    288313         ! 
    289314      ENDIF 
     
    303328         IF( ln_sla ) THEN 
    304329            jtype = jtype + 1 
    305             CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla   ', & 
    306                &                   cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
    307             CALL obs_setinterpopts( nsurftypes, jtype, 'sla   ',      & 
    308                &                  nn_2dint, nn_2dint_sla,             & 
    309                &                  rn_sla_avglamscl, rn_sla_avgphiscl, & 
    310                &                  ln_sla_fp_indegs, .FALSE.,          & 
    311                &                  n2dintsurf, zavglamscl, zavgphiscl, & 
    312                &                  lfpindegs, llnightav ) 
     330            cobstypessurf(jtype) = 'sla' 
     331            clsurffiles(jtype,:) = cn_slafbfiles 
    313332         ENDIF 
    314333         IF( ln_sst ) THEN 
    315334            jtype = jtype + 1 
    316             CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst   ', & 
    317                &                   cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
    318             CALL obs_setinterpopts( nsurftypes, jtype, 'sst   ',      & 
    319                &                  nn_2dint, nn_2dint_sst,             & 
    320                &                  rn_sst_avglamscl, rn_sst_avgphiscl, & 
    321                &                  ln_sst_fp_indegs, ln_sstnight,      & 
    322                &                  n2dintsurf, zavglamscl, zavgphiscl, & 
    323                &                  lfpindegs, llnightav ) 
     335            cobstypessurf(jtype) = 'sst' 
     336            clsurffiles(jtype,:) = cn_sstfbfiles 
    324337         ENDIF 
    325338#if defined key_si3 || defined key_cice 
    326339         IF( ln_sic ) THEN 
    327340            jtype = jtype + 1 
    328             CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic   ', & 
    329                &                   cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
    330             CALL obs_setinterpopts( nsurftypes, jtype, 'sic   ',      & 
    331                &                  nn_2dint, nn_2dint_sic,             & 
    332                &                  rn_sic_avglamscl, rn_sic_avgphiscl, & 
    333                &                  ln_sic_fp_indegs, .FALSE.,          & 
    334                &                  n2dintsurf, zavglamscl, zavgphiscl, & 
    335                &                  lfpindegs, llnightav ) 
     341            cobstypessurf(jtype) = 'sic' 
     342            clsurffiles(jtype,:) = cn_sicfbfiles 
    336343         ENDIF 
    337344#endif 
    338345         IF( ln_sss ) THEN 
    339346            jtype = jtype + 1 
    340             CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss   ', & 
    341                &                   cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
    342             CALL obs_setinterpopts( nsurftypes, jtype, 'sss   ',      & 
    343                &                  nn_2dint, nn_2dint_sss,             & 
    344                &                  rn_sss_avglamscl, rn_sss_avgphiscl, & 
    345                &                  ln_sss_fp_indegs, .FALSE.,          & 
    346                &                  n2dintsurf, zavglamscl, zavgphiscl, & 
    347                &                  lfpindegs, llnightav ) 
     347            cobstypessurf(jtype) = 'sss' 
     348            clsurffiles(jtype,:) = cn_sssfbfiles 
    348349         ENDIF 
     350         ! 
     351         CALL obs_settypefiles( nsurftypes, jpmaxnfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     352 
     353         DO jtype = 1, nsurftypes 
     354 
     355            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
     356               IF ( nn_2dint_sla == -1 ) THEN 
     357                  n2dint_type  = nn_2dint_default 
     358               ELSE 
     359                  n2dint_type  = nn_2dint_sla 
     360               ENDIF 
     361               ztype_avglamscl = rn_sla_avglamscl 
     362               ztype_avgphiscl = rn_sla_avgphiscl 
     363               ltype_fp_indegs = ln_sla_fp_indegs 
     364               ltype_night     = .FALSE. 
     365            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN 
     366               IF ( nn_2dint_sst == -1 ) THEN 
     367                  n2dint_type  = nn_2dint_default 
     368               ELSE 
     369                  n2dint_type  = nn_2dint_sst 
     370               ENDIF 
     371               ztype_avglamscl = rn_sst_avglamscl 
     372               ztype_avgphiscl = rn_sst_avgphiscl 
     373               ltype_fp_indegs = ln_sst_fp_indegs 
     374               ltype_night     = ln_sstnight 
     375            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN 
     376               IF ( nn_2dint_sic == -1 ) THEN 
     377                  n2dint_type  = nn_2dint_default 
     378               ELSE 
     379                  n2dint_type  = nn_2dint_sic 
     380               ENDIF 
     381               ztype_avglamscl = rn_sic_avglamscl 
     382               ztype_avgphiscl = rn_sic_avgphiscl 
     383               ltype_fp_indegs = ln_sic_fp_indegs 
     384               ltype_night     = .FALSE. 
     385            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 
     386               IF ( nn_2dint_sss == -1 ) THEN 
     387                  n2dint_type  = nn_2dint_default 
     388               ELSE 
     389                  n2dint_type  = nn_2dint_sss 
     390               ENDIF 
     391               ztype_avglamscl = rn_sss_avglamscl 
     392               ztype_avgphiscl = rn_sss_avgphiscl 
     393               ltype_fp_indegs = ln_sss_fp_indegs 
     394               ltype_night     = .FALSE. 
     395            ELSE 
     396               n2dint_type     = nn_2dint_default 
     397               ztype_avglamscl = rn_default_avglamscl 
     398               ztype_avgphiscl = rn_default_avgphiscl 
     399               ltype_fp_indegs = ln_default_fp_indegs 
     400               ltype_night     = .FALSE. 
     401            ENDIF 
     402             
     403            CALL obs_setinterpopts( nsurftypes, jtype, TRIM(cobstypessurf(jtype)), & 
     404               &                    nn_2dint_default, n2dint_type,                 & 
     405               &                    ztype_avglamscl, ztype_avgphiscl,              & 
     406               &                    ltype_fp_indegs, ltype_night,                  & 
     407               &                    n2dintsurf, zavglamscl, zavgphiscl,            & 
     408               &                    lfpindegs, llnightav ) 
     409 
     410         END DO 
    349411         ! 
    350412      ENDIF 
     
    368430      ENDIF 
    369431      ! 
    370       IF( nn_2dint < 0  .OR.  nn_2dint > 6  ) THEN 
    371          CALL ctl_stop('dia_obs_init: Choice of horizontal (2D) interpolation method is not available') 
     432      IF( nn_2dint_default < 0  .OR.  nn_2dint_default > 6  ) THEN 
     433         CALL ctl_stop('dia_obs_init: Choice of default horizontal (2D) interpolation method is not available') 
    372434      ENDIF 
    373435      ! 
     
    388450         DO jtype = 1, nproftypes 
    389451            ! 
    390             nvarsprof(jtype) = 2 
    391452            IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN 
     453               nvarsprof(jtype) = 2 
    392454               nextrprof(jtype) = 1 
    393                llvar1 = ln_t3d 
    394                llvar2 = ln_s3d 
    395                zglam1 = glamt 
    396                zgphi1 = gphit 
    397                zmask1 = tmask 
    398                zglam2 = glamt 
    399                zgphi2 = gphit 
    400                zmask2 = tmask 
    401             ENDIF 
    402             IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN 
     455               IF ( ln_output_clim ) ltype_clim = .TRUE.               
     456               ALLOCATE( llvar (nvarsprof(jtype)) ) 
     457               ALLOCATE( clvars(nvarsprof(jtype)) ) 
     458               ALLOCATE( zglam(jpi, jpj,      nvarsprof(jtype)) ) 
     459               ALLOCATE( zgphi(jpi, jpj,      nvarsprof(jtype)) ) 
     460               ALLOCATE( zmask(jpi, jpj, jpk, nvarsprof(jtype)) ) 
     461               llvar(1)       = ln_t3d 
     462               llvar(2)       = ln_s3d 
     463               clvars(1)      = 'POTM' 
     464               clvars(2)      = 'PSAL' 
     465               zglam(:,:,1)   = glamt(:,:) 
     466               zglam(:,:,2)   = glamt(:,:) 
     467               zgphi(:,:,1)   = gphit(:,:) 
     468               zgphi(:,:,2)   = gphit(:,:) 
     469               zmask(:,:,:,1) = tmask(:,:,:) 
     470               zmask(:,:,:,2) = tmask(:,:,:) 
     471            ELSE IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN 
     472               nvarsprof(jtype) = 2 
    403473               nextrprof(jtype) = 2 
    404                llvar1 = ln_vel3d 
    405                llvar2 = ln_vel3d 
    406                zglam1 = glamu 
    407                zgphi1 = gphiu 
    408                zmask1 = umask 
    409                zglam2 = glamv 
    410                zgphi2 = gphiv 
    411                zmask2 = vmask 
     474               ALLOCATE( llvar (nvarsprof(jtype)) ) 
     475               ALLOCATE( clvars(nvarsprof(jtype)) ) 
     476               ALLOCATE( zglam(jpi, jpj,      nvarsprof(jtype)) ) 
     477               ALLOCATE( zgphi(jpi, jpj,      nvarsprof(jtype)) ) 
     478               ALLOCATE( zmask(jpi, jpj, jpk, nvarsprof(jtype)) ) 
     479               llvar(1)       = ln_vel3d 
     480               llvar(2)       = ln_vel3d 
     481               clvars(1)      = 'UVEL' 
     482               clvars(2)      = 'VVEL' 
     483               zglam(:,:,1)   = glamu(:,:) 
     484               zglam(:,:,2)   = glamv(:,:) 
     485               zgphi(:,:,1)   = gphiu(:,:) 
     486               zgphi(:,:,2)   = gphiv(:,:) 
     487               zmask(:,:,:,1) = umask(:,:,:) 
     488               zmask(:,:,:,2) = vmask(:,:,:) 
     489            ELSE 
     490               nvarsprof(jtype) = 1 
     491               nextrprof(jtype) = 0 
     492               ALLOCATE( llvar (nvarsprof(jtype)) ) 
     493               ALLOCATE( clvars(nvarsprof(jtype)) ) 
     494               ALLOCATE( zglam(jpi, jpj,      nvarsprof(jtype)) ) 
     495               ALLOCATE( zgphi(jpi, jpj,      nvarsprof(jtype)) ) 
     496               ALLOCATE( zmask(jpi, jpj, jpk, nvarsprof(jtype)) ) 
     497               llvar(1)       = .TRUE. 
     498               zglam(:,:,1)   = glamt(:,:) 
     499               zgphi(:,:,1)   = gphit(:,:) 
     500               zmask(:,:,:,1) = tmask(:,:,:) 
    412501            ENDIF 
    413502            ! 
     
    416505               &               clproffiles(jtype,1:ifilesprof(jtype)), & 
    417506               &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & 
    418                &               rn_dobsini, rn_dobsend, llvar1, llvar2, & 
    419                &               ln_ignmis, ln_s_at_t, .FALSE., & 
     507               &               rn_dobsini, rn_dobsend, llvar, & 
     508               &               ln_ignmis, ln_s_at_t, .FALSE., clvars, & 
    420509               &               kdailyavtypes = nn_profdavtypes ) 
    421510               ! 
     
    425514            ! 
    426515            CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), & 
    427                &               llvar1, llvar2, & 
     516               &               llvar, & 
    428517               &               jpi, jpj, jpk, & 
    429                &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2, & 
     518               &               zmask, zglam, zgphi, & 
    430519               &               ln_nea, ln_bound_reject, Kmm, & 
    431520               &               kdailyavtypes = nn_profdavtypes ) 
     521            ! 
     522            DEALLOCATE( llvar, clvars, zglam, zgphi, zmask ) 
     523            ! 
    432524         END DO 
    433525         ! 
     
    449541            IF( TRIM(cobstypessurf(jtype)) == 'sst' )   llnightav(jtype) = ln_sstnight 
    450542            ! 
     543            ALLOCATE( clvars( nvarssurf(jtype) ) ) 
     544            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
     545               clvars(1) = 'SLA' 
     546            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN 
     547               clvars(1) = 'SST' 
     548            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN 
     549               clvars(1) = 'ICECONC' 
     550            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 
     551               clvars(1) = 'SSS' 
     552            ENDIF 
     553            ! 
    451554            ! Read in surface obs types 
    452555            CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & 
    453556               &               clsurffiles(jtype,1:ifilessurf(jtype)), & 
    454557               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 
    455                &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) 
     558               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype), * 
     559               &               clvars ) 
    456560               ! 
    457561            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 
     
    473577                  &                  jnumsstbias      , cn_sstbiasfiles(1:jnumsstbias) )  
    474578            ENDIF 
     579            ! 
     580            DEALLOCATE( clvars ) 
    475581         END DO 
    476582         ! 
     
    516622      INTEGER :: jvar              ! Variable number 
    517623      INTEGER :: ji, jj            ! Loop counters 
    518       REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 
    519          & zprofvar1, &            ! Model values for 1st variable in a prof ob 
    520          & zprofvar2               ! Model values for 2nd variable in a prof ob 
    521       REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 
    522          & zprofmask1, &           ! Mask associated with zprofvar1 
    523          & zprofmask2              ! Mask associated with zprofvar2 
     624      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
     625         & zprofvar,               ! Model values for variables in a prof ob 
     626      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
     627         & zprofmask               ! Mask associated with zprofvar 
    524628      REAL(wp), DIMENSION(jpi,jpj) :: & 
    525629         & zsurfvar, &             ! Model values equivalent to surface ob. 
    526630         & zsurfmask               ! Mask associated with surface variable 
    527       REAL(wp), DIMENSION(jpi,jpj) :: & 
    528          & zglam1,    &            ! Model longitudes for prof variable 1 
    529          & zglam2,    &            ! Model longitudes for prof variable 2 
    530          & zgphi1,    &            ! Model latitudes for prof variable 1 
    531          & zgphi2                  ! Model latitudes for prof variable 2 
     631      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
     632         & zglam,    &             ! Model longitudes for prof variables 
     633         & zgphi,                  ! Model latitudes for prof variables 
    532634 
    533635      !----------------------------------------------------------------------- 
     
    549651         DO jtype = 1, nproftypes 
    550652 
     653            ! Allocate local work arrays 
     654            ALLOCATE( zprofvar (jpi, jpj, jpk, profdataqc(jtype)%nvar) ) 
     655            ALLOCATE( zprofmask(jpi, jpj, jpk, profdataqc(jtype)%nvar) ) 
     656            ALLOCATE( zglam    (jpi, jpj,      profdataqc(jtype)%nvar) ) 
     657            ALLOCATE( zgphi    (jpi, jpj,      profdataqc(jtype)%nvar) )   
     658                               
     659            ! Defaults which might change 
     660            DO jvar = 1, profdataqc(jtype)%nvar 
     661               zprofmask(:,:,:,jvar) = tmask(:,:,:) 
     662               zglam(:,:,jvar)       = glamt(:,:) 
     663               zgphi(:,:,jvar)       = gphit(:,:) 
     664            END DO 
     665 
    551666            SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 
    552667            CASE('prof') 
    553668               zprofvar1(:,:,:) = ts(:,:,:,jp_tem,Kmm) 
    554669               zprofvar2(:,:,:) = ts(:,:,:,jp_sal,Kmm) 
    555                zprofmask1(:,:,:) = tmask(:,:,:) 
    556                zprofmask2(:,:,:) = tmask(:,:,:) 
    557                zglam1(:,:) = glamt(:,:) 
    558                zglam2(:,:) = glamt(:,:) 
    559                zgphi1(:,:) = gphit(:,:) 
    560                zgphi2(:,:) = gphit(:,:) 
    561670            CASE('vel') 
    562                zprofvar1(:,:,:) = uu(:,:,:,Kmm) 
    563                zprofvar2(:,:,:) = vv(:,:,:,Kmm) 
    564                zprofmask1(:,:,:) = umask(:,:,:) 
    565                zprofmask2(:,:,:) = vmask(:,:,:) 
    566                zglam1(:,:) = glamu(:,:) 
    567                zglam2(:,:) = glamv(:,:) 
    568                zgphi1(:,:) = gphiu(:,:) 
    569                zgphi2(:,:) = gphiv(:,:) 
     671               zprofvar(:,:,:,1) = uu(:,:,:,Kmm) 
     672               zprofvar(:,:,:,2) = vv(:,:,:,Kmm) 
     673               zprofmask(:,:,:,1) = umask(:,:,:) 
     674               zprofmask(:,:,:,2) = vmask(:,:,:) 
     675               zglam(:,:,1) = glamu(:,:) 
     676               zglam(:,:,2) = glamv(:,:) 
     677               zgphi(:,:,1) = gphiu(:,:) 
     678               zgphi(:,:,2) = gphiv(:,:) 
    570679            CASE DEFAULT 
    571680               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 
    572681            END SELECT 
    573682 
    574             CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  & 
    575                &               nit000, idaystp,                         & 
    576                &               zprofvar1, zprofvar2,                    & 
    577                &               gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm),      &  
    578                &               zprofmask1, zprofmask2,                  & 
    579                &               zglam1, zglam2, zgphi1, zgphi2,          & 
    580                &               nn_1dint, nn_2dint,                      & 
    581                &               kdailyavtypes = nn_profdavtypes ) 
     683            DO jvar = 1, profdataqc(jtype)%nvar 
     684               CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  & 
     685                  &               nit000, idaystp, jvar,                   & 
     686                  &               zprofvar(:,:,:,jvar),                    & 
     687                  &               gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm),      &  
     688                  &               zprofmask(:,:,:,jvar),                   & 
     689                  &               zglam(:,:,jvar), zgphi(:,:,jvar),        & 
     690                  &               nn_1dint, nn_2dint_default,              & 
     691                  &               kdailyavtypes = nn_profdavtypes ) 
     692            END DO 
     693             
     694            DEALLOCATE( zprofvar, zprofmask, zglam, zgphi ) 
    582695 
    583696         END DO 
     
    680793                  & ) 
    681794 
    682                CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv ) 
     795               CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv ) 
    683796 
    684797               DO jo = 1, profdataqc(jtype)%nprof 
     
    8961009   END SUBROUTINE fin_date 
    8971010    
    898     SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, jtype, ctypein, & 
    899        &                         cfilestype, ifiles, cobstypes, cfiles ) 
    900  
    901     INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types 
    902     INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type 
    903     INTEGER, INTENT(IN) :: jtype       ! Index of the current type of obs 
    904     INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 
    905        &                   ifiles      ! Out appended number of files for this type 
    906  
    907     CHARACTER(len=6), INTENT(IN) :: ctypein  
    908     CHARACTER(len=128), DIMENSION(jpmaxnfiles), INTENT(IN) :: & 
    909        &                   cfilestype  ! In list of files for this obs type 
    910     CHARACTER(len=6), DIMENSION(ntypes), INTENT(INOUT) :: & 
    911        &                   cobstypes   ! Out appended list of obs types 
    912     CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(INOUT) :: & 
    913        &                   cfiles      ! Out appended list of files for all types 
    914  
    915     !Local variables 
    916     INTEGER :: jfile 
    917  
    918     cfiles(jtype,:) = cfilestype(:) 
    919     cobstypes(jtype) = ctypein 
    920     ifiles(jtype) = 0 
    921     DO jfile = 1, jpmaxnfiles 
    922        IF ( trim(cfiles(jtype,jfile)) /= '' ) & 
    923                  ifiles(jtype) = ifiles(jtype) + 1 
    924     END DO 
    925  
    926     IF ( ifiles(jtype) == 0 ) THEN 
    927          CALL ctl_stop( 'Logical for observation type '//TRIM(ctypein)//   & 
    928             &           ' set to true but no files available to read' ) 
    929     ENDIF 
    930  
    931     IF(lwp) THEN     
    932        WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:' 
    933        DO jfile = 1, ifiles(jtype) 
    934           WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile)) 
    935        END DO 
    936     ENDIF 
    937  
    938     END SUBROUTINE obs_settypefiles 
    939  
    940     SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             & 
    941                &                  n2dint_default, n2dint_type,        & 
    942                &                  zavglamscl_type, zavgphiscl_type,   & 
    943                &                  lfp_indegs_type, lavnight_type,     & 
    944                &                  n2dint, zavglamscl, zavgphiscl,     & 
    945                &                  lfpindegs, lavnight ) 
    946  
    947     INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types 
    948     INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs 
    949     INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type 
    950     INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type 
    951     REAL(wp), INTENT(IN) :: & 
    952        &                    zavglamscl_type, & !E/W diameter of obs footprint for this type 
    953        &                    zavgphiscl_type    !N/S diameter of obs footprint for this type 
    954     LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres 
    955     LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average 
    956     CHARACTER(len=6), INTENT(IN) :: ctypein  
    957  
    958     INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 
    959        &                    n2dint  
    960     REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & 
    961        &                    zavglamscl, zavgphiscl 
    962     LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & 
    963        &                    lfpindegs, lavnight 
    964  
    965     lavnight(jtype) = lavnight_type 
    966  
    967     IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN 
    968        n2dint(jtype) = n2dint_type 
    969     ELSE 
    970        n2dint(jtype) = n2dint_default 
    971     ENDIF 
    972  
    973     ! For averaging observation footprints set options for size of footprint  
    974     IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN 
    975        IF ( zavglamscl_type > 0._wp ) THEN 
    976           zavglamscl(jtype) = zavglamscl_type 
    977        ELSE 
    978           CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
    979                          'scale (zavglamscl) for observation type '//TRIM(ctypein) )       
    980        ENDIF 
    981  
    982        IF ( zavgphiscl_type > 0._wp ) THEN 
    983           zavgphiscl(jtype) = zavgphiscl_type 
    984        ELSE 
    985           CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
    986                          'scale (zavgphiscl) for observation type '//TRIM(ctypein) )       
    987        ENDIF 
    988  
    989        lfpindegs(jtype) = lfp_indegs_type  
    990  
    991     ENDIF 
    992  
    993     ! Write out info  
    994     IF(lwp) THEN 
    995        IF ( n2dint(jtype) <= 4 ) THEN 
    996           WRITE(numout,*) '             '//TRIM(ctypein)// & 
    997              &            ' model counterparts will be interpolated horizontally' 
    998        ELSE IF ( n2dint(jtype) <= 6 ) THEN 
    999           WRITE(numout,*) '             '//TRIM(ctypein)// & 
    1000              &            ' model counterparts will be averaged horizontally' 
    1001           WRITE(numout,*) '             '//'    with E/W scale: ',zavglamscl(jtype) 
    1002           WRITE(numout,*) '             '//'    with N/S scale: ',zavgphiscl(jtype) 
    1003           IF ( lfpindegs(jtype) ) THEN 
    1004               WRITE(numout,*) '             '//'    (in degrees)' 
    1005           ELSE 
    1006               WRITE(numout,*) '             '//'    (in metres)' 
    1007           ENDIF 
    1008        ENDIF 
    1009     ENDIF 
    1010  
    1011     END SUBROUTINE obs_setinterpopts 
     1011   SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, ifiles, cobstypes, cfiles ) 
     1012 
     1013      INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types 
     1014      INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type 
     1015      INTEGER, DIMENSION(ntypes), INTENT(OUT) :: & 
     1016         &                   ifiles      ! Out number of files for each type 
     1017      CHARACTER(len=8), DIMENSION(ntypes), INTENT(IN) :: & 
     1018         &                   cobstypes   ! List of obs types 
     1019      CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(IN) :: & 
     1020         &                   cfiles      ! List of files for all types 
     1021 
     1022      !Local variables 
     1023      INTEGER :: jfile 
     1024      INTEGER :: jtype 
     1025 
     1026      DO jtype = 1, ntypes 
     1027 
     1028         ifiles(jtype) = 0 
     1029         DO jfile = 1, jpmaxnfiles 
     1030            IF ( trim(cfiles(jtype,jfile)) /= '' ) & 
     1031                      ifiles(jtype) = ifiles(jtype) + 1 
     1032         END DO 
     1033 
     1034         IF ( ifiles(jtype) == 0 ) THEN 
     1035              CALL ctl_stop( 'Logical for observation type '//TRIM(cobstypes(jtype))//  & 
     1036                 &           ' set to true but no files available to read' ) 
     1037         ENDIF 
     1038 
     1039         IF(lwp) THEN     
     1040            WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:' 
     1041            DO jfile = 1, ifiles(jtype) 
     1042               WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile)) 
     1043            END DO 
     1044         ENDIF 
     1045 
     1046      END DO 
     1047 
     1048   END SUBROUTINE obs_settypefiles 
     1049 
     1050   SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             & 
     1051              &                  n2dint_default, n2dint_type,        & 
     1052              &                  ravglamscl_type, ravgphiscl_type,   & 
     1053              &                  lfp_indegs_type, lavnight_type,     & 
     1054              &                  n2dint, ravglamscl, ravgphiscl,     & 
     1055              &                  lfpindegs, lavnight ) 
     1056 
     1057      INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types 
     1058      INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs 
     1059      INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type 
     1060      INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type 
     1061      REAL(wp), INTENT(IN) :: & 
     1062         &                    ravglamscl_type, & !E/W diameter of obs footprint for this type 
     1063         &                    ravgphiscl_type    !N/S diameter of obs footprint for this type 
     1064      LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres 
     1065      LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average 
     1066      CHARACTER(len=8), INTENT(IN) :: ctypein  
     1067 
     1068      INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1069         &                    n2dint  
     1070      REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & 
     1071         &                    ravglamscl, ravgphiscl 
     1072      LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1073         &                    lfpindegs, lavnight 
     1074 
     1075      lavnight(jtype) = lavnight_type 
     1076 
     1077      IF ( (n2dint_type >= 0) .AND. (n2dint_type <= 6) ) THEN 
     1078         n2dint(jtype) = n2dint_type 
     1079      ELSE IF ( n2dint_type == -1 ) THEN 
     1080         n2dint(jtype) = n2dint_default 
     1081      ELSE 
     1082         CALL ctl_stop(' Choice of '//TRIM(ctypein)//' horizontal (2D) interpolation method', & 
     1083           &                    ' is not available') 
     1084      ENDIF 
     1085 
     1086      ! For averaging observation footprints set options for size of footprint  
     1087      IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN 
     1088         IF ( ravglamscl_type > 0._wp ) THEN 
     1089            ravglamscl(jtype) = ravglamscl_type 
     1090         ELSE 
     1091            CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
     1092                           'scale (ravglamscl) for observation type '//TRIM(ctypein) )       
     1093         ENDIF 
     1094 
     1095         IF ( ravgphiscl_type > 0._wp ) THEN 
     1096            ravgphiscl(jtype) = ravgphiscl_type 
     1097         ELSE 
     1098            CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
     1099                           'scale (ravgphiscl) for observation type '//TRIM(ctypein) )       
     1100         ENDIF 
     1101 
     1102         lfpindegs(jtype) = lfp_indegs_type  
     1103 
     1104      ENDIF 
     1105 
     1106      ! Write out info  
     1107      IF(lwp) THEN 
     1108         IF ( n2dint(jtype) <= 4 ) THEN 
     1109            WRITE(numout,*) '             '//TRIM(ctypein)// & 
     1110               &            ' model counterparts will be interpolated horizontally' 
     1111         ELSE IF ( n2dint(jtype) <= 6 ) THEN 
     1112            WRITE(numout,*) '             '//TRIM(ctypein)// & 
     1113               &            ' model counterparts will be averaged horizontally' 
     1114            WRITE(numout,*) '             '//'    with E/W scale: ',ravglamscl(jtype) 
     1115            WRITE(numout,*) '             '//'    with N/S scale: ',ravgphiscl(jtype) 
     1116            IF ( lfpindegs(jtype) ) THEN 
     1117                WRITE(numout,*) '             '//'    (in degrees)' 
     1118            ELSE 
     1119                WRITE(numout,*) '             '//'    (in metres)' 
     1120            ENDIF 
     1121         ENDIF 
     1122      ENDIF 
     1123 
     1124   END SUBROUTINE obs_setinterpopts 
    10121125 
    10131126END MODULE diaobs 
  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/OBS/obs_oper.F90

    r13295 r13758  
    4040CONTAINS 
    4141 
    42    SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk,          & 
    43       &                     kit000, kdaystp,                      & 
    44       &                     pvar1, pvar2, pgdept, pgdepw,         & 
    45       &                     pmask1, pmask2,                       &   
    46       &                     plam1, plam2, pphi1, pphi2,           & 
     42   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 
     43      &                     kit000, kdaystp, kvar,       & 
     44      &                     pvar, pgdept, pgdepw,        & 
     45      &                     pmask,                       &   
     46      &                     plam, pphi,                  & 
    4747      &                     k1dint, k2dint, kdailyavtypes ) 
    4848      !!----------------------------------------------------------------------- 
     
    105105      INTEGER       , INTENT(in   ) ::   k2dint          ! Horizontal interpolation type (see header) 
    106106      INTEGER       , INTENT(in   ) ::   kdaystp         ! Number of time steps per day 
    107       REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj,kpk) ::   pvar1 , pvar2    ! Model field     1 and 2 
    108       REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj,kpk) ::   pmask1, pmask2   ! Land-sea mask   1 and 2 
    109       REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj)     ::   plam1 , plam2    ! Model longitude 1 and 2 
    110       REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj)     ::   pphi1 , pphi2    ! Model latitudes 1 and 2 
     107      INTEGER       , INTENT(in   ) ::   kvar            ! Number of variables in prodatqc 
     108      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj,kpk) ::   pvar             ! Model field 
     109      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj,kpk) ::   pmask            ! Land-sea mask 
     110      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj)     ::   plam             ! Model longitude 
     111      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj)     ::   pphi             ! Model latitudes 
    111112      REAL(KIND=wp) , INTENT(in   ), DIMENSION(kpi,kpj,kpk) ::   pgdept, pgdepw   ! depth of T and W levels  
    112113      INTEGER, DIMENSION(imaxavtypes), OPTIONAL ::   kdailyavtypes             ! Types for daily averages 
     
    128129         & idailyavtypes 
    129130      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    130          & igrdi1, & 
    131          & igrdi2, & 
    132          & igrdj1, & 
    133          & igrdj2 
     131         & igrdi, & 
     132         & igrdj 
    134133      INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 
    135134 
     
    138137      REAL(KIND=wp) :: zdaystp 
    139138      REAL(KIND=wp), DIMENSION(kpk) :: & 
    140          & zobsmask1, & 
    141          & zobsmask2, & 
    142          & zobsk,    & 
     139         & zobsk,  & 
    143140         & zobs2k 
    144141      REAL(KIND=wp), DIMENSION(2,2,1) :: & 
    145142         & zweig1, & 
    146          & zweig2, & 
    147143         & zweig 
    148144      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
    149          & zmask1, & 
    150          & zmask2, & 
    151          & zint1,  & 
    152          & zint2,  & 
    153          & zinm1,  & 
    154          & zinm2,  & 
     145         & zmask,  & 
     146         & zint,   & 
     147         & zinm,   & 
    155148         & zgdept, &  
    156149         & zgdepw 
    157150      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    158          & zglam1, & 
    159          & zglam2, & 
    160          & zgphi1, & 
    161          & zgphi2 
    162       REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2    
     151         & zglam,  & 
     152         & zgphi 
     153      REAL(KIND=wp), DIMENSION(1) :: zmsk 
    163154      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 
    164155 
     
    190181         IF ( idayend == 1 .OR. kt == 0 ) THEN 
    191182            DO_3D( 1, 1, 1, 1, 1, jpk ) 
    192                prodatqc%vdmean(ji,jj,jk,1) = 0.0 
    193                prodatqc%vdmean(ji,jj,jk,2) = 0.0 
     183               prodatqc%vdmean(ji,jj,jk,kvar) = 0.0 
    194184            END_3D 
    195185         ENDIF 
     
    197187         DO_3D( 1, 1, 1, 1, 1, jpk ) 
    198188            ! Increment field 1 for computing daily mean 
    199             prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
    200                &                        + pvar1(ji,jj,jk) 
    201             ! Increment field 2 for computing daily mean 
    202             prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
    203                &                        + pvar2(ji,jj,jk) 
     189            prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) & 
     190               &                           + pvar(ji,jj,jk) 
    204191         END_3D 
    205192 
     
    210197            CALL FLUSH(numout) 
    211198            DO_3D( 1, 1, 1, 1, 1, jpk ) 
    212                prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
    213                   &                        * zdaystp 
    214                prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
    215                   &                        * zdaystp 
     199               prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) & 
     200                  &                           * zdaystp 
    216201            END_3D 
    217202         ENDIF 
     
    221206      ! Get the data for interpolation 
    222207      ALLOCATE( & 
    223          & igrdi1(2,2,ipro),      & 
    224          & igrdi2(2,2,ipro),      & 
    225          & igrdj1(2,2,ipro),      & 
    226          & igrdj2(2,2,ipro),      & 
    227          & zglam1(2,2,ipro),      & 
    228          & zglam2(2,2,ipro),      & 
    229          & zgphi1(2,2,ipro),      & 
    230          & zgphi2(2,2,ipro),      & 
    231          & zmask1(2,2,kpk,ipro),  & 
    232          & zmask2(2,2,kpk,ipro),  & 
    233          & zint1(2,2,kpk,ipro),   & 
    234          & zint2(2,2,kpk,ipro),   & 
    235          & zgdept(2,2,kpk,ipro),  &  
    236          & zgdepw(2,2,kpk,ipro)   &  
     208         & igrdi(2,2,ipro),      & 
     209         & igrdj(2,2,ipro),      & 
     210         & zglam(2,2,ipro),      & 
     211         & zgphi(2,2,ipro),      & 
     212         & zmask(2,2,kpk,ipro),  & 
     213         & zint(2,2,kpk,ipro),   & 
     214         & zgdept(2,2,kpk,ipro), &  
     215         & zgdepw(2,2,kpk,ipro)  &  
    237216         & ) 
    238217 
    239218      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
    240219         iobs = jobs - prodatqc%nprofup 
    241          igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1 
    242          igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1 
    243          igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1 
    244          igrdj1(1,2,iobs) = prodatqc%mj(jobs,1) 
    245          igrdi1(2,1,iobs) = prodatqc%mi(jobs,1) 
    246          igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1 
    247          igrdi1(2,2,iobs) = prodatqc%mi(jobs,1) 
    248          igrdj1(2,2,iobs) = prodatqc%mj(jobs,1) 
    249          igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1 
    250          igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1 
    251          igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1 
    252          igrdj2(1,2,iobs) = prodatqc%mj(jobs,2) 
    253          igrdi2(2,1,iobs) = prodatqc%mi(jobs,2) 
    254          igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1 
    255          igrdi2(2,2,iobs) = prodatqc%mi(jobs,2) 
    256          igrdj2(2,2,iobs) = prodatqc%mj(jobs,2) 
     220         igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1 
     221         igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1 
     222         igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1 
     223         igrdj(1,2,iobs) = prodatqc%mj(jobs,1) 
     224         igrdi(2,1,iobs) = prodatqc%mi(jobs,1) 
     225         igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1 
     226         igrdi(2,2,iobs) = prodatqc%mi(jobs,1) 
     227         igrdj(2,2,iobs) = prodatqc%mj(jobs,1) 
    257228      END DO 
    258229 
     
    261232      zgdepw(:,:,:,:) = 0.0 
    262233 
    263       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 ) 
    264       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 ) 
    265       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pmask1, zmask1 ) 
    266       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1,   zint1 ) 
    267        
    268       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 ) 
    269       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 ) 
    270       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 ) 
    271       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2,   zint2 ) 
    272  
    273       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept )  
    274       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw )  
     234      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, plam, zglam ) 
     235      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 
     236      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pmask, zmask ) 
     237      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pvar,   zint ) 
     238 
     239      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept, zgdept )  
     240      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw )  
    275241 
    276242      ! At the end of the day also get interpolated means 
    277243      IF ( ld_dailyav .AND. idayend == 0 ) THEN 
    278244 
    279          ALLOCATE( & 
    280             & zinm1(2,2,kpk,ipro),  & 
    281             & zinm2(2,2,kpk,ipro)   & 
    282             & ) 
    283  
    284          CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, & 
    285             &                  prodatqc%vdmean(:,:,:,1), zinm1 ) 
    286          CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, & 
    287             &                  prodatqc%vdmean(:,:,:,2), zinm2 ) 
     245         ALLOCATE( zinm(2,2,kpk,ipro) ) 
     246 
     247         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 
     248            &                  prodatqc%vdmean(:,:,:,kvar), zinm ) 
    288249 
    289250      ENDIF 
     
    320281         ! Horizontal weights  
    321282         ! Masked values are calculated later.   
    322          IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 
     283         IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN 
    323284 
    324285            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
    325                &                   zglam1(:,:,iobs), zgphi1(:,:,iobs), & 
    326                &                   zmask1(:,:,1,iobs), zweig1, zmsk_1 ) 
    327  
    328          ENDIF 
    329  
    330          IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    331  
    332             CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
    333                &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), & 
    334                &                   zmask2(:,:,1,iobs), zweig2, zmsk_2) 
    335   
    336          ENDIF 
    337  
    338          IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 
     286               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     287               &                   zmask(:,:,1,iobs), zweig1, zmsk_1 ) 
     288 
     289         ENDIF 
     290 
     291         IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN 
    339292 
    340293            zobsk(:) = obfillflt 
     
    346299 
    347300                  ! vertically interpolate all 4 corners  
    348                   ista = prodatqc%npvsta(jobs,1)  
    349                   iend = prodatqc%npvend(jobs,1)  
     301                  ista = prodatqc%npvsta(jobs,kvar)  
     302                  iend = prodatqc%npvend(jobs,kvar)  
    350303                  inum_obs = iend - ista + 1  
    351304                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
     
    356309                        IF ( k1dint == 1 ) THEN  
    357310                           CALL obs_int_z1d_spl( kpk, &  
    358                               &     zinm1(iin,ijn,:,iobs), &  
     311                              &     zinm(iin,ijn,:,iobs), &  
    359312                              &     zobs2k, zgdept(iin,ijn,:,iobs), &  
    360                               &     zmask1(iin,ijn,:,iobs))  
     313                              &     zmask(iin,ijn,:,iobs))  
    361314                        ENDIF  
    362315        
    363316                        CALL obs_level_search(kpk, &  
    364317                           &    zgdept(iin,ijn,:,iobs), &  
    365                            &    inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
     318                           &    inum_obs, prodatqc%var(kvar)%vdep(ista:iend), &  
    366319                           &    iv_indic)  
    367320 
    368321                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
    369                            &    prodatqc%var(1)%vdep(ista:iend), &  
    370                            &    zinm1(iin,ijn,:,iobs), &  
     322                           &    prodatqc%var(kvar)%vdep(ista:iend), &  
     323                           &    zinm(iin,ijn,:,iobs), &  
    371324                           &    zobs2k, interp_corner(iin,ijn,:), &  
    372325                           &    zgdept(iin,ijn,:,iobs), &  
    373                            &    zmask1(iin,ijn,:,iobs))  
     326                           &    zmask(iin,ijn,:,iobs))  
    374327        
    375328                     ENDDO  
     
    383336      
    384337               ! vertically interpolate all 4 corners  
    385                ista = prodatqc%npvsta(jobs,1)  
    386                iend = prodatqc%npvend(jobs,1)  
     338               ista = prodatqc%npvsta(jobs,kvar)  
     339               iend = prodatqc%npvend(jobs,kvar)  
    387340               inum_obs = iend - ista + 1  
    388341               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     
    392345                     IF ( k1dint == 1 ) THEN  
    393346                        CALL obs_int_z1d_spl( kpk, &  
    394                            &    zint1(iin,ijn,:,iobs),&  
     347                           &    zint(iin,ijn,:,iobs),&  
    395348                           &    zobs2k, zgdept(iin,ijn,:,iobs), &  
    396                            &    zmask1(iin,ijn,:,iobs))  
     349                           &    zmask(iin,ijn,:,iobs))  
    397350   
    398351                     ENDIF  
     
    400353                     CALL obs_level_search(kpk, &  
    401354                         &        zgdept(iin,ijn,:,iobs),&  
    402                          &        inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
     355                         &        inum_obs, prodatqc%var(kvar)%vdep(ista:iend), &  
    403356                         &        iv_indic)  
    404357 
    405358                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
    406                          &          prodatqc%var(1)%vdep(ista:iend),     &  
    407                          &          zint1(iin,ijn,:,iobs),            &  
     359                         &          prodatqc%var(kvar)%vdep(ista:iend),     &  
     360                         &          zint(iin,ijn,:,iobs),            &  
    408361                         &          zobs2k,interp_corner(iin,ijn,:), &  
    409362                         &          zgdept(iin,ijn,:,iobs),         &  
    410                          &          zmask1(iin,ijn,:,iobs) )       
     363                         &          zmask(iin,ijn,:,iobs) )       
    411364          
    412365                  ENDDO  
     
    432385                  DO ijn=1,2  
    433386      
    434                      depth_loop1: DO ik=kpk,2,-1  
    435                         IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
     387                     depth_loop: DO ik=kpk,2,-1  
     388                        IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
    436389                             
    437390                           zweig(iin,ijn,1) = &   
    438391                              & zweig1(iin,ijn,1) * &  
    439392                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
    440                               &  - prodatqc%var(1)%vdep(iend)),0._wp)  
     393                              &  - prodatqc%var(kvar)%vdep(iend)),0._wp)  
    441394                             
    442                            EXIT depth_loop1  
     395                           EXIT depth_loop  
    443396 
    444397                        ENDIF  
    445398 
    446                      ENDDO depth_loop1  
     399                     ENDDO depth_loop 
    447400      
    448401                  ENDDO  
     
    454407                  ! Set QC flag for any observations found below the bottom 
    455408                  ! needed as the check here is more strict than that in obs_prep 
    456                IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4 
     409               IF (sum(zweig) == 0.0_wp) prodatqc%var(kvar)%nvqc(iend:iend)=4 
    457410  
    458411            ENDDO  
     
    460413            DEALLOCATE(interp_corner,iv_indic)  
    461414           
    462          ENDIF  
    463  
    464          ! For the second variable 
    465          IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    466  
    467             zobsk(:) = obfillflt 
    468  
    469             IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 
    470  
    471                IF ( idayend == 0 )  THEN 
    472                   ! Daily averaged data 
    473  
    474                   ! vertically interpolate all 4 corners  
    475                   ista = prodatqc%npvsta(jobs,2)  
    476                   iend = prodatqc%npvend(jobs,2)  
    477                   inum_obs = iend - ista + 1  
    478                   ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
    479  
    480                   DO iin=1,2  
    481                      DO ijn=1,2  
    482  
    483                         IF ( k1dint == 1 ) THEN  
    484                            CALL obs_int_z1d_spl( kpk, &  
    485                               &     zinm2(iin,ijn,:,iobs), &  
    486                               &     zobs2k, zgdept(iin,ijn,:,iobs), &  
    487                               &     zmask2(iin,ijn,:,iobs))  
    488                         ENDIF  
    489         
    490                         CALL obs_level_search(kpk, &  
    491                            &    zgdept(iin,ijn,:,iobs), &  
    492                            &    inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
    493                            &    iv_indic)  
    494  
    495                         CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
    496                            &    prodatqc%var(2)%vdep(ista:iend), &  
    497                            &    zinm2(iin,ijn,:,iobs), &  
    498                            &    zobs2k, interp_corner(iin,ijn,:), &  
    499                            &    zgdept(iin,ijn,:,iobs), &  
    500                            &    zmask2(iin,ijn,:,iobs))  
    501         
    502                      ENDDO  
    503                   ENDDO  
    504  
    505                ENDIF !idayend 
    506  
    507             ELSE    
    508  
    509                ! Point data  
    510       
    511                ! vertically interpolate all 4 corners  
    512                ista = prodatqc%npvsta(jobs,2)  
    513                iend = prodatqc%npvend(jobs,2)  
    514                inum_obs = iend - ista + 1  
    515                ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
    516                DO iin=1,2   
    517                   DO ijn=1,2  
    518                      
    519                      IF ( k1dint == 1 ) THEN  
    520                         CALL obs_int_z1d_spl( kpk, &  
    521                            &    zint2(iin,ijn,:,iobs),&  
    522                            &    zobs2k, zgdept(iin,ijn,:,iobs), &  
    523                            &    zmask2(iin,ijn,:,iobs))  
    524    
    525                      ENDIF  
    526         
    527                      CALL obs_level_search(kpk, &  
    528                          &        zgdept(iin,ijn,:,iobs),&  
    529                          &        inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
    530                          &        iv_indic)  
    531  
    532                      CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
    533                          &          prodatqc%var(2)%vdep(ista:iend),     &  
    534                          &          zint2(iin,ijn,:,iobs),            &  
    535                          &          zobs2k,interp_corner(iin,ijn,:), &  
    536                          &          zgdept(iin,ijn,:,iobs),         &  
    537                          &          zmask2(iin,ijn,:,iobs) )       
    538           
    539                   ENDDO  
    540                ENDDO  
    541               
    542             ENDIF  
    543  
    544             !-------------------------------------------------------------  
    545             ! Compute the horizontal interpolation for every profile level  
    546             !-------------------------------------------------------------  
    547               
    548             DO ikn=1,inum_obs  
    549                iend=ista+ikn-1 
    550                    
    551                zweig(:,:,1) = 0._wp  
    552     
    553                ! This code forces the horizontal weights to be   
    554                ! zero IF the observation is below the bottom of the   
    555                ! corners of the interpolation nodes, Or if it is in   
    556                ! the mask. This is important for observations near   
    557                ! steep bathymetry  
    558                DO iin=1,2  
    559                   DO ijn=1,2  
    560       
    561                      depth_loop2: DO ik=kpk,2,-1  
    562                         IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
    563                              
    564                            zweig(iin,ijn,1) = &   
    565                               & zweig2(iin,ijn,1) * &  
    566                               & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
    567                               &  - prodatqc%var(2)%vdep(iend)),0._wp)  
    568                              
    569                            EXIT depth_loop2  
    570  
    571                         ENDIF  
    572  
    573                      ENDDO depth_loop2  
    574       
    575                   ENDDO  
    576                ENDDO  
    577     
    578                CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
    579                   &              prodatqc%var(2)%vmod(iend:iend) )  
    580  
    581                   ! Set QC flag for any observations found below the bottom 
    582                   ! needed as the check here is more strict than that in obs_prep 
    583                IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 
    584   
    585             ENDDO  
    586   
    587             DEALLOCATE(interp_corner,iv_indic)  
    588            
    589          ENDIF  
     415         ENDIF 
    590416 
    591417      ENDDO 
    592418 
    593419      ! Deallocate the data for interpolation 
    594       DEALLOCATE( & 
    595          & igrdi1, & 
    596          & igrdi2, & 
    597          & igrdj1, & 
    598          & igrdj2, & 
    599          & zglam1, & 
    600          & zglam2, & 
    601          & zgphi1, & 
    602          & zgphi2, & 
    603          & zmask1, & 
    604          & zmask2, & 
    605          & zint1,  & 
    606          & zint2,  & 
     420      DEALLOCATE(  & 
     421         & igrdi,  & 
     422         & igrdj,  & 
     423         & zglam,  & 
     424         & zgphi,  & 
     425         & zmask,  & 
     426         & zint,   & 
    607427         & zgdept, & 
    608428         & zgdepw  & 
     
    611431      ! At the end of the day also get interpolated means 
    612432      IF ( ld_dailyav .AND. idayend == 0 ) THEN 
    613          DEALLOCATE( & 
    614             & zinm1,  & 
    615             & zinm2   & 
    616             & ) 
     433         DEALLOCATE( zinm ) 
    617434      ENDIF 
    618435 
    619       prodatqc%nprofup = prodatqc%nprofup + ipro  
     436      IF ( kvar == prodatqc%nvar ) THEN 
     437         prodatqc%nprofup = prodatqc%nprofup + ipro  
     438      ENDIF 
    620439 
    621440   END SUBROUTINE obs_prof_opt 
  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/OBS/obs_prep.F90

    r12489 r13758  
    241241 
    242242 
    243    SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_var1, ld_var2, & 
     243   SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_var, & 
    244244      &                     kpi, kpj, kpk, & 
    245       &                     zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2,  & 
     245      &                     zmask, pglam, pgphi,  & 
    246246      &                     ld_nea, ld_bound_reject, Kmm, kdailyavtypes,  kqc_cutoff ) 
    247247 
     
    269269      TYPE(obs_prof), INTENT(INOUT) :: profdata   ! Full set of profile data 
    270270      TYPE(obs_prof), INTENT(INOUT) :: prodatqc   ! Subset of profile data not failing screening 
    271       LOGICAL, INTENT(IN) :: ld_var1              ! Observed variables switches 
    272       LOGICAL, INTENT(IN) :: ld_var2 
     271      LOGICAL, DIMENSION(profdata%nvar), INTENT(IN) :: & 
     272         & ld_var                                 ! Observed variables switches 
    273273      LOGICAL, INTENT(IN) :: ld_nea               ! Switch for rejecting observation near land 
    274274      LOGICAL, INTENT(IN) :: ld_bound_reject      ! Switch for rejecting observations near the boundary 
     
    277277      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    278278         & kdailyavtypes                          ! Types for daily averages 
    279       REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
    280          & zmask1, & 
    281          & zmask2 
    282       REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    283          & pglam1, & 
    284          & pglam2, & 
    285          & pgphi1, & 
    286          & pgphi2 
     279      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,kpk,profdata%nvar) :: & 
     280         & zmask 
     281      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,profdata%nvar) :: & 
     282         & pglam, & 
     283         & pgphi 
    287284      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value 
    288285 
     
    295292      INTEGER :: imin0 
    296293      INTEGER :: icycle       ! Current assimilation cycle 
    297                               ! Counters for observations that are 
    298       INTEGER :: iotdobs      !  - outside time domain 
    299       INTEGER :: iosdv1obs    !  - outside space domain (variable 1) 
    300       INTEGER :: iosdv2obs    !  - outside space domain (variable 2) 
    301       INTEGER :: ilanv1obs    !  - within a model land cell (variable 1) 
    302       INTEGER :: ilanv2obs    !  - within a model land cell (variable 2) 
    303       INTEGER :: inlav1obs    !  - close to land (variable 1) 
    304       INTEGER :: inlav2obs    !  - close to land (variable 2) 
    305       INTEGER :: ibdyv1obs    !  - boundary (variable 1)  
    306       INTEGER :: ibdyv2obs    !  - boundary (variable 2)       
    307       INTEGER :: igrdobs      !  - fail the grid search 
    308       INTEGER :: iuvchku      !  - reject u if v rejected and vice versa 
    309       INTEGER :: iuvchkv      ! 
    310                               ! Global counters for observations that are 
    311       INTEGER :: iotdobsmpp   !  - outside time domain 
    312       INTEGER :: iosdv1obsmpp !  - outside space domain (variable 1) 
    313       INTEGER :: iosdv2obsmpp !  - outside space domain (variable 2) 
    314       INTEGER :: ilanv1obsmpp !  - within a model land cell (variable 1) 
    315       INTEGER :: ilanv2obsmpp !  - within a model land cell (variable 2) 
    316       INTEGER :: inlav1obsmpp !  - close to land (variable 1) 
    317       INTEGER :: inlav2obsmpp !  - close to land (variable 2) 
    318       INTEGER :: ibdyv1obsmpp !  - boundary (variable 1)  
    319       INTEGER :: ibdyv2obsmpp !  - boundary (variable 2)       
    320       INTEGER :: igrdobsmpp   !  - fail the grid search 
    321       INTEGER :: iuvchkumpp   !  - reject var1 if var2 rejected and vice versa 
    322       INTEGER :: iuvchkvmpp   ! 
     294                                                       ! Counters for observations that are 
     295      INTEGER                           :: iotdobs     !  - outside time domain 
     296      INTEGER, DIMENSION(profdata%nvar) :: iosdvobs    !  - outside space domain 
     297      INTEGER, DIMENSION(profdata%nvar) :: ilanvobs    !  - within a model land cell 
     298      INTEGER, DIMENSION(profdata%nvar) :: inlavobs    !  - close to land 
     299      INTEGER, DIMENSION(profdata%nvar) :: ibdyvobs    !  - boundary    
     300      INTEGER                           :: igrdobs     !  - fail the grid search 
     301      INTEGER                           :: iuvchku     !  - reject UVEL if VVEL rejected 
     302      INTEGER                           :: iuvchkv     !  - reject VVEL if UVEL rejected 
     303                                                       ! Global counters for observations that are 
     304      INTEGER                           :: iotdobsmpp  !  - outside time domain 
     305      INTEGER, DIMENSION(profdata%nvar) :: iosdvobsmpp !  - outside space domain 
     306      INTEGER, DIMENSION(profdata%nvar) :: ilanvobsmpp !  - within a model land cell 
     307      INTEGER, DIMENSION(profdata%nvar) :: inlavobsmpp !  - close to land 
     308      INTEGER, DIMENSION(profdata%nvar) :: ibdyvobsmpp !  - boundary 
     309      INTEGER :: igrdobsmpp                            !  - fail the grid search 
     310      INTEGER :: iuvchkumpp                            !  - reject UVEL if VVEL rejected 
     311      INTEGER :: iuvchkvmpp                            !  - reject VVEL if UVEL rejected 
    323312      TYPE(obs_prof_valid) ::  llvalid      ! Profile selection  
    324313      TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: & 
    325          & llvvalid           ! var1,var2 selection  
     314         & llvvalid           ! var selection  
    326315      INTEGER :: jvar         ! Variable loop variable 
    327316      INTEGER :: jobs         ! Obs. loop variable 
    328317      INTEGER :: jstp         ! Time loop variable 
    329318      INTEGER :: inrc         ! Time index variable 
     319      CHARACTER(LEN=256) :: cout1  ! Diagnostic output line 
     320      CHARACTER(LEN=256) :: cout2  ! Diagnostic output line 
    330321      !!---------------------------------------------------------------------- 
    331322 
     
    342333      icycle = nn_no     ! Assimilation cycle 
    343334 
    344       ! Diagnotics counters for various failures. 
    345  
    346       iotdobs   = 0 
    347       igrdobs   = 0 
    348       iosdv1obs = 0 
    349       iosdv2obs = 0 
    350       ilanv1obs = 0 
    351       ilanv2obs = 0 
    352       inlav1obs = 0 
    353       inlav2obs = 0 
    354       ibdyv1obs = 0 
    355       ibdyv2obs = 0 
    356       iuvchku   = 0 
    357       iuvchkv   = 0 
     335      ! Diagnostic counters for various failures. 
     336 
     337      iotdobs     = 0 
     338      igrdobs     = 0 
     339      iosdvobs(:) = 0 
     340      ilanvobs(:) = 0 
     341      inlavobs(:) = 0 
     342      ibdyvobs(:) = 0 
     343      iuvchku     = 0 
     344      iuvchkv     = 0 
    358345 
    359346 
     
    388375      ! ----------------------------------------------------------------------- 
    389376 
    390       CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,1), profdata%mj(:,1), & 
    391          &              profdata%nqc,     igrdobs                         ) 
    392       CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,2), profdata%mj(:,2), & 
    393          &              profdata%nqc,     igrdobs                         ) 
     377      DO jvar = 1, profdata%nvar 
     378         CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,jvar), profdata%mj(:,jvar), & 
     379            &              profdata%nqc,     igrdobs                         ) 
     380      END DO 
    394381 
    395382      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 
     
    406393      ! ----------------------------------------------------------------------- 
    407394 
    408       ! Variable 1 
    409       CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(1),   & 
    410          &                 profdata%npvsta(:,1),  profdata%npvend(:,1), & 
    411          &                 jpi,                   jpj,                  & 
    412          &                 jpk,                                         & 
    413          &                 profdata%mi,           profdata%mj,          & 
    414          &                 profdata%var(1)%mvk,                         & 
    415          &                 profdata%rlam,         profdata%rphi,        & 
    416          &                 profdata%var(1)%vdep,                        & 
    417          &                 pglam1,                pgphi1,               & 
    418          &                 gdept_1d,              zmask1,               & 
    419          &                 profdata%nqc,          profdata%var(1)%nvqc, & 
    420          &                 iosdv1obs,             ilanv1obs,            & 
    421          &                 inlav1obs,             ld_nea,               & 
    422          &                 ibdyv1obs,             ld_bound_reject,      & 
    423          &                 iqc_cutoff,            Kmm                 ) 
    424  
    425       CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp ) 
    426       CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp ) 
    427       CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp ) 
    428       CALL obs_mpp_sum_integer( ibdyv1obs, ibdyv1obsmpp ) 
    429  
    430       ! Variable 2 
    431       CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(2),   & 
    432          &                 profdata%npvsta(:,2),  profdata%npvend(:,2), & 
    433          &                 jpi,                   jpj,                  & 
    434          &                 jpk,                                         & 
    435          &                 profdata%mi,           profdata%mj,          &  
    436          &                 profdata%var(2)%mvk,                         & 
    437          &                 profdata%rlam,         profdata%rphi,        & 
    438          &                 profdata%var(2)%vdep,                        & 
    439          &                 pglam2,                pgphi2,               & 
    440          &                 gdept_1d,              zmask2,               & 
    441          &                 profdata%nqc,          profdata%var(2)%nvqc, & 
    442          &                 iosdv2obs,             ilanv2obs,            & 
    443          &                 inlav2obs,             ld_nea,               & 
    444          &                 ibdyv2obs,             ld_bound_reject,      & 
    445          &                 iqc_cutoff,            Kmm                 ) 
    446  
    447       CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp ) 
    448       CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp ) 
    449       CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp ) 
    450       CALL obs_mpp_sum_integer( ibdyv2obs, ibdyv2obsmpp ) 
     395      DO jvar = 1, profdata%nvar 
     396         CALL obs_coo_spc_3d( profdata%nprof,          profdata%nvprot(jvar),   & 
     397            &                 profdata%npvsta(:,jvar), profdata%npvend(:,jvar), & 
     398            &                 jpi,                     jpj,                     & 
     399            &                 jpk,                                              & 
     400            &                 profdata%mi,             profdata%mj,             & 
     401            &                 profdata%var(jvar)%mvk,                           & 
     402            &                 profdata%rlam,           profdata%rphi,           & 
     403            &                 profdata%var(jvar)%vdep,                          & 
     404            &                 pglam(:,:,jvar),         pgphi(:,:,jvar),         & 
     405            &                 gdept_1d,                zmask(:,:,:,jvar),       & 
     406            &                 profdata%nqc,            profdata%var(jvar)%nvqc, & 
     407            &                 iosdvobs(jvar),          ilanvobs(jvar),          & 
     408            &                 inlavobs(jvar),          ld_nea,                  & 
     409            &                 ibdyvobs(jvar),          ld_bound_reject,         & 
     410            &                 iqc_cutoff       ) 
     411 
     412         CALL obs_mpp_sum_integer( iosdvobs(jvar), iosdvobsmpp(jvar) ) 
     413         CALL obs_mpp_sum_integer( ilanvobs(jvar), ilanvobsmpp(jvar) ) 
     414         CALL obs_mpp_sum_integer( inlavobs(jvar), inlavobsmpp(jvar) ) 
     415         CALL obs_mpp_sum_integer( ibdyvobs(jvar), ibdyvobsmpp(jvar) ) 
     416      END DO 
    451417 
    452418      ! ----------------------------------------------------------------------- 
     
    499465       
    500466         WRITE(numout,*) 
    501          WRITE(numout,*) ' Profiles outside time domain                     = ', & 
     467         WRITE(numout,*) ' Profiles outside time domain                       = ', & 
    502468            &            iotdobsmpp 
    503          WRITE(numout,*) ' Remaining profiles that failed grid search       = ', & 
     469         WRITE(numout,*) ' Remaining profiles that failed grid search         = ', & 
    504470            &            igrdobsmpp 
    505          WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data outside space domain       = ', & 
    506             &            iosdv1obsmpp 
    507          WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data at land points             = ', & 
    508             &            ilanv1obsmpp 
    509          IF (ld_nea) THEN 
    510             WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (removed) = ',& 
    511                &            inlav1obsmpp 
    512          ELSE 
    513             WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (kept)    = ',& 
    514                &            inlav1obsmpp 
    515          ENDIF 
    516          IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 
    517             WRITE(numout,*) ' U observation rejected since V rejected     = ', & 
    518                &            iuvchku 
    519          ENDIF 
    520          WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near open boundary (removed) = ',& 
    521                &            ibdyv1obsmpp 
    522          WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted                             = ', & 
    523             &            prodatqc%nvprotmpp(1) 
    524          WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data outside space domain       = ', & 
    525             &            iosdv2obsmpp 
    526          WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data at land points             = ', & 
    527             &            ilanv2obsmpp 
    528          IF (ld_nea) THEN 
    529             WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (removed) = ',& 
    530                &            inlav2obsmpp 
    531          ELSE 
    532             WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (kept)    = ',& 
    533                &            inlav2obsmpp 
    534          ENDIF 
    535          IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 
    536             WRITE(numout,*) ' V observation rejected since U rejected     = ', & 
    537                &            iuvchkv 
    538          ENDIF 
    539          WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near open boundary (removed) = ',& 
    540                &            ibdyv2obsmpp 
    541          WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted                             = ', & 
    542             &            prodatqc%nvprotmpp(2) 
     471         DO jvar = 1, profdata%nvar 
     472            WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data outside space domain       = ', & 
     473               &            iosdvobsmpp(jvar) 
     474            WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data at land points             = ', & 
     475               &            ilanvobsmpp(jvar) 
     476            IF (ld_nea) THEN 
     477               WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near land points (removed) = ',& 
     478                  &            inlavobsmpp(jvar) 
     479            ELSE 
     480               WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near land points (kept)    = ',& 
     481                  &            inlavobsmpp(jvar) 
     482            ENDIF 
     483            IF ( TRIM(profdata%cvars(jvar)) == 'UVEL' ) THEN 
     484               WRITE(numout,*) ' U observation rejected since V rejected     = ', & 
     485                  &            iuvchku 
     486            ELSE IF ( TRIM(profdata%cvars(jvar)) == 'VVEL' ) THEN 
     487               WRITE(numout,*) ' V observation rejected since U rejected     = ', & 
     488                  &            iuvchkv 
     489            ENDIF 
     490            WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near open boundary (removed) = ',& 
     491                  &            ibdyvobsmpp(jvar) 
     492            WRITE(numout,*) ' '//prodatqc%cvars(jvar)//' data accepted                             = ', & 
     493               &            prodatqc%nvprotmpp(jvar) 
     494         END DO 
    543495 
    544496         WRITE(numout,*) 
    545497         WRITE(numout,*) ' Number of observations per time step :' 
    546498         WRITE(numout,*) 
    547          WRITE(numout,'(10X,A,5X,A,5X,A,A)')'Time step','Profiles', & 
    548             &                               '     '//prodatqc%cvars(1)//'     ', & 
    549             &                               '     '//prodatqc%cvars(2)//'     ' 
    550          WRITE(numout,998) 
     499         WRITE(cout1,'(10X,A9,5X,A8)') 'Time step', 'Profiles' 
     500         WRITE(cout2,'(10X,A9,5X,A8)') '---------', '--------' 
     501         DO jvar = 1, prodatqc%nvar 
     502            WRITE(cout1,'(A,5X,A11)') TRIM(cout1), TRIM(prodatqc%cvars(jvar)) 
     503            WRITE(cout2,'(A,5X,A11)') TRIM(cout2), '-----------' 
     504         END DO 
     505         WRITE(numout,*) cout1 
     506         WRITE(numout,*) cout2 
    551507      ENDIF 
    552508       
     
    575531         DO jstp = nit000 - 1, nitend 
    576532            inrc = jstp - nit000 + 2 
    577             WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), & 
    578                &                    prodatqc%nvstpmpp(inrc,1), & 
    579                &                    prodatqc%nvstpmpp(inrc,2) 
     533            WRITE(cout1,'(10X,I9,5X,I8)') jstp, prodatqc%npstpmpp(inrc) 
     534            DO jvar = 1, prodatqc%nvar 
     535               WRITE(cout1,'(A,5X,I11)') TRIM(cout1), prodatqc%nvstpmpp(inrc,jvar) 
     536            END DO 
     537            WRITE(numout,*) cout1 
    580538         END DO 
    581539      ENDIF 
    582  
    583 998   FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'----------------') 
    584 999   FORMAT(10X,I9,5X,I8,5X,I11,5X,I8) 
    585540 
    586541   END SUBROUTINE obs_pre_prof 
  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/OBS/obs_read_prof.F90

    r13226 r13758  
    4545   SUBROUTINE obs_rea_prof( profdata, knumfiles, cdfilenames, & 
    4646      &                     kvars, kextr, kstp, ddobsini, ddobsend, & 
    47       &                     ldvar1, ldvar2, ldignmis, ldsatt, & 
    48       &                     ldmod, kdailyavtypes ) 
     47      &                     ldvar, ldignmis, ldsatt, & 
     48      &                     ldmod, cdvars, kdailyavtypes ) 
    4949      !!--------------------------------------------------------------------- 
    5050      !! 
     
    7474      INTEGER, INTENT(IN) :: kextr      ! Number of extra fields for each var 
    7575      INTEGER, INTENT(IN) :: kstp       ! Ocean time-step index 
    76       LOGICAL, INTENT(IN) :: ldvar1     ! Observed variables switches 
    77       LOGICAL, INTENT(IN) :: ldvar2 
     76      LOGICAL, DIMENSION(kvars), INTENT(IN) :: ldvar     ! Observed variables switches 
    7877      LOGICAL, INTENT(IN) :: ldignmis   ! Ignore missing files 
    7978      LOGICAL, INTENT(IN) :: ldsatt     ! Compute salinity at all temperature points 
     
    8180      REAL(dp), INTENT(IN) :: ddobsini  ! Obs. ini time in YYYYMMDD.HHMMSS 
    8281      REAL(dp), INTENT(IN) :: ddobsend  ! Obs. end time in YYYYMMDD.HHMMSS 
     82      CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars 
    8383      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    8484         & kdailyavtypes                ! Types of daily average observations 
     
    8787      CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof' 
    8888      CHARACTER(len=8) :: clrefdate 
    89       CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 
     89      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvarsin 
    9090      INTEGER :: jvar 
    9191      INTEGER :: ji 
     
    105105      INTEGER :: iprof 
    106106      INTEGER :: iproftot 
    107       INTEGER :: ivar1t0 
    108       INTEGER :: ivar2t0 
    109       INTEGER :: ivar1t 
    110       INTEGER :: ivar2t 
     107      INTEGER, DIMENSION(kvars) :: ivart0 
     108      INTEGER, DIMENSION(kvars) :: ivart 
    111109      INTEGER :: ip3dt 
    112110      INTEGER :: ios 
    113111      INTEGER :: ioserrcount 
    114       INTEGER :: ivar1tmpp 
    115       INTEGER :: ivar2tmpp 
     112      INTEGER, DIMENSION(kvars) :: ivartmpp 
    116113      INTEGER :: ip3dtmpp 
    117114      INTEGER :: itype 
    118115      INTEGER, DIMENSION(knumfiles) :: & 
    119116         & irefdate 
    120       INTEGER, DIMENSION(ntyp1770+1) :: & 
    121          & itypvar1,    & 
    122          & itypvar1mpp, & 
    123          & itypvar2,    & 
    124          & itypvar2mpp  
     117      INTEGER, DIMENSION(ntyp1770+1,kvars) :: & 
     118         & itypvar,    & 
     119         & itypvarmpp 
     120      INTEGER, DIMENSION(:,:), ALLOCATABLE :: & 
     121         & iobsi,    & 
     122         & iobsj,    & 
     123         & iproc 
    125124      INTEGER, DIMENSION(:), ALLOCATABLE :: & 
    126          & iobsi1,    & 
    127          & iobsj1,    & 
    128          & iproc1,    & 
    129          & iobsi2,    & 
    130          & iobsj2,    & 
    131          & iproc2,    & 
    132125         & iindx,    & 
    133126         & ifileidx, & 
     
    147140      LOGICAL :: llvalprof 
    148141      LOGICAL :: lldavtimset 
     142      LOGICAL :: llcycle 
    149143      TYPE(obfbdata), POINTER, DIMENSION(:) :: & 
    150144         & inpfiles 
     
    152146      ! Local initialization 
    153147      iprof = 0 
    154       ivar1t0 = 0 
    155       ivar2t0 = 0 
     148      ivart0(:) = 0 
    156149      ip3dt = 0 
    157150 
     
    219212               &                ldgrid = .TRUE. ) 
    220213 
    221             IF ( inpfiles(jj)%nvar < 2 ) THEN 
     214            IF ( inpfiles(jj)%nvar /= kvars ) THEN 
    222215               CALL ctl_stop( 'Feedback format error: ', & 
    223                   &           ' less than 2 vars in profile file' ) 
     216                  &           ' unexpected number of vars in profile file' ) 
    224217            ENDIF 
    225218 
     
    229222 
    230223            IF ( jj == 1 ) THEN 
    231                ALLOCATE( clvars( inpfiles(jj)%nvar ) ) 
     224               ALLOCATE( clvarsin( inpfiles(jj)%nvar ) ) 
    232225               DO ji = 1, inpfiles(jj)%nvar 
    233                  clvars(ji) = inpfiles(jj)%cname(ji) 
     226                 clvarsin(ji) = inpfiles(jj)%cname(ji) 
     227                 IF ( clvarsin(ji) /= cdvars(ji) ) THEN 
     228                    CALL ctl_stop( 'Feedback file variables do not match', & 
     229                        &           ' expected variable names for this type' ) 
     230                 ENDIF 
    234231               END DO 
    235232            ELSE 
    236233               DO ji = 1, inpfiles(jj)%nvar 
    237                   IF ( inpfiles(jj)%cname(ji) /= clvars(ji) ) THEN 
     234                  IF ( inpfiles(jj)%cname(ji) /= clvarsin(ji) ) THEN 
    238235                     CALL ctl_stop( 'Feedback file variables not consistent', & 
    239236                        &           ' with previous files for this type' ) 
     
    308305            DO ji = 1, inpfiles(jj)%nobs 
    309306               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
    310                IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
    311                   & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
     307               llcycle = .TRUE. 
     308               DO jvar = 1, kvars 
     309                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 
     310                     llcycle = .FALSE. 
     311                     EXIT 
     312                  ENDIF 
     313               END DO 
     314               IF ( llcycle ) CYCLE 
    312315               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    313316                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    317320            ALLOCATE( zlam(inowin)  ) 
    318321            ALLOCATE( zphi(inowin)  ) 
    319             ALLOCATE( iobsi1(inowin) ) 
    320             ALLOCATE( iobsj1(inowin) ) 
    321             ALLOCATE( iproc1(inowin) ) 
    322             ALLOCATE( iobsi2(inowin) ) 
    323             ALLOCATE( iobsj2(inowin) ) 
    324             ALLOCATE( iproc2(inowin) ) 
     322            ALLOCATE( iobsi(inowin,kvars) ) 
     323            ALLOCATE( iobsj(inowin,kvars) ) 
     324            ALLOCATE( iproc(inowin,kvars) ) 
    325325            inowin = 0 
    326326            DO ji = 1, inpfiles(jj)%nobs 
    327327               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
    328                IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
    329                   & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
     328               llcycle = .TRUE. 
     329               DO jvar = 1, kvars 
     330                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 
     331                     llcycle = .FALSE. 
     332                     EXIT 
     333                  ENDIF 
     334               END DO 
     335               IF ( llcycle ) CYCLE 
    330336               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    331337                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    336342            END DO 
    337343 
    338             IF ( TRIM( inpfiles(jj)%cname(1) ) == 'POTM' ) THEN 
    339                CALL obs_grid_search( inowin, zlam, zphi, iobsi1, iobsj1, & 
    340                   &                  iproc1, 'T' ) 
    341                iobsi2(:) = iobsi1(:) 
    342                iobsj2(:) = iobsj1(:) 
    343                iproc2(:) = iproc1(:) 
    344             ELSEIF ( TRIM( inpfiles(jj)%cname(1) ) == 'UVEL' ) THEN 
    345                CALL obs_grid_search( inowin, zlam, zphi, iobsi1, iobsj1, & 
    346                   &                  iproc1, 'U' ) 
    347                CALL obs_grid_search( inowin, zlam, zphi, iobsi2, iobsj2, & 
    348                   &                  iproc2, 'V' ) 
     344            ! Assume anything other than velocity is on T grid 
     345            IF ( TRIM( inpfiles(jj)%cname(1) ) == 'UVEL' ) THEN 
     346               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), & 
     347                  &                  iproc(:,1), 'U' ) 
     348               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,2), iobsj(:,2), & 
     349                  &                  iproc(:,2), 'V' ) 
     350            ELSE 
     351               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), & 
     352                  &                  iproc(:,1), 'T' ) 
     353               IF ( kvars > 1 ) THEN 
     354                  DO jvar = 2, kvars 
     355                     iobsi(:,jvar) = iobsi(:,1) 
     356                     iobsj(:,jvar) = iobsj(:,1) 
     357                     iproc(:,jvar) = iproc(:,1) 
     358                  END DO 
     359               ENDIF 
    349360            ENDIF 
    350361 
     
    352363            DO ji = 1, inpfiles(jj)%nobs 
    353364               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
    354                IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
    355                   & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
     365               llcycle = .TRUE. 
     366               DO jvar = 1, kvars 
     367                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 
     368                     llcycle = .FALSE. 
     369                     EXIT 
     370                  ENDIF 
     371               END DO 
     372               IF ( llcycle ) CYCLE 
    356373               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    357374                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
    358375                  inowin = inowin + 1 
    359                   inpfiles(jj)%iproc(ji,1) = iproc1(inowin) 
    360                   inpfiles(jj)%iobsi(ji,1) = iobsi1(inowin) 
    361                   inpfiles(jj)%iobsj(ji,1) = iobsj1(inowin) 
    362                   inpfiles(jj)%iproc(ji,2) = iproc2(inowin) 
    363                   inpfiles(jj)%iobsi(ji,2) = iobsi2(inowin) 
    364                   inpfiles(jj)%iobsj(ji,2) = iobsj2(inowin) 
    365                   IF ( inpfiles(jj)%iproc(ji,1) /= & 
    366                      & inpfiles(jj)%iproc(ji,2) ) THEN 
    367                      CALL ctl_stop( 'Error in obs_read_prof:', & 
    368                         & 'var1 and var2 observation on different processors') 
     376                  DO jvar = 1, kvars 
     377                     inpfiles(jj)%iproc(ji,jvar) = iproc(inowin,jvar) 
     378                     inpfiles(jj)%iobsi(ji,jvar) = iobsi(inowin,jvar) 
     379                     inpfiles(jj)%iobsj(ji,jvar) = iobsj(inowin,jvar) 
     380                  END DO 
     381                  IF ( kvars > 1 ) THEN 
     382                     DO jvar = 2, kvars 
     383                        IF ( inpfiles(jj)%iproc(ji,jvar) /= & 
     384                           & inpfiles(jj)%iproc(ji,1) ) THEN 
     385                           CALL ctl_stop( 'Error in obs_read_prof:', & 
     386                              & 'observation on different processors for different vars') 
     387                        ENDIF 
     388                     END DO 
    369389                  ENDIF 
    370390               ENDIF 
    371391            END DO 
    372             DEALLOCATE( zlam, zphi, iobsi1, iobsj1, iproc1, iobsi2, iobsj2, iproc2 ) 
     392            DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc ) 
    373393 
    374394            DO ji = 1, inpfiles(jj)%nobs 
    375395               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
    376                IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
    377                   & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
     396               llcycle = .TRUE. 
     397               DO jvar = 1, kvars 
     398                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 
     399                     llcycle = .FALSE. 
     400                     EXIT 
     401                  ENDIF 
     402               END DO 
     403               IF ( llcycle ) CYCLE 
    378404               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    379405                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    384410                  ENDIF 
    385411                  llvalprof = .FALSE. 
    386                   IF ( ldvar1 ) THEN 
    387                      loop_t_count : DO ij = 1,inpfiles(jj)%nlev 
    388                         IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    389                            & CYCLE 
    390                         IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
    391                            & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
    392                            ivar1t0 = ivar1t0 + 1 
    393                         ENDIF 
    394                      END DO loop_t_count 
    395                   ENDIF 
    396                   IF ( ldvar2 ) THEN 
    397                      loop_s_count : DO ij = 1,inpfiles(jj)%nlev 
    398                         IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    399                            & CYCLE 
    400                         IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
    401                            & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
    402                            ivar2t0 = ivar2t0 + 1 
    403                         ENDIF 
    404                      END DO loop_s_count 
    405                   ENDIF 
    406                   loop_p_count : DO ij = 1,inpfiles(jj)%nlev 
     412                  DO jvar = 1, kvars 
     413                     IF ( ldvar(jvar) ) THEN 
     414                        DO ij = 1,inpfiles(jj)%nlev 
     415                           IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
     416                              & CYCLE 
     417                           IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 
     418                              & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
     419                              ivart0(jvar) = ivart0(jvar) + 1 
     420                           ENDIF 
     421                        END DO 
     422                     ENDIF 
     423                  END DO 
     424                  DO ij = 1,inpfiles(jj)%nlev 
    407425                     IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    408426                        & CYCLE 
    409                      IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
    410                         &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
    411                         &    ldvar1 ) .OR. & 
    412                         & ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
    413                         &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
    414                         &     ldvar2 ) ) THEN 
    415                         ip3dt = ip3dt + 1 
    416                         llvalprof = .TRUE. 
    417                      ENDIF 
    418                   END DO loop_p_count 
     427                     DO jvar = 1, kvars 
     428                        IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 
     429                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     430                           &    ldvar(jvar) ) ) THEN 
     431                           ip3dt = ip3dt + 1 
     432                           llvalprof = .TRUE. 
     433                           EXIT 
     434                        ENDIF 
     435                     END DO 
     436                  END DO 
    419437 
    420438                  IF ( llvalprof ) iprof = iprof + 1 
     
    438456         DO ji = 1, inpfiles(jj)%nobs 
    439457            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
    440             IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
    441                & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
     458            llcycle = .TRUE. 
     459            DO jvar = 1, kvars 
     460               IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 
     461                  llcycle = .FALSE. 
     462                  EXIT 
     463               ENDIF 
     464            END DO 
     465            IF ( llcycle ) CYCLE 
    442466            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    443467               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    453477         DO ji = 1, inpfiles(jj)%nobs 
    454478            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
    455             IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
    456                & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
     479            llcycle = .TRUE. 
     480            DO jvar = 1, kvars 
     481               IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 
     482                  llcycle = .FALSE. 
     483                  EXIT 
     484               ENDIF 
     485            END DO 
     486            IF ( llcycle ) CYCLE 
    457487            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    458488               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    470500      iv3dt(:) = -1 
    471501      IF (ldsatt) THEN 
    472          iv3dt(1) = ip3dt 
    473          iv3dt(2) = ip3dt 
     502         iv3dt(:) = ip3dt 
    474503      ELSE 
    475          iv3dt(1) = ivar1t0 
    476          iv3dt(2) = ivar2t0 
     504         iv3dt(:) = ivart0(:) 
    477505      ENDIF 
    478506      CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, & 
     
    483511      profdata%nprof     = 0 
    484512      profdata%nvprot(:) = 0 
    485       profdata%cvars(:)  = clvars(:) 
     513      profdata%cvars(:)  = clvarsin(:) 
    486514      iprof = 0 
    487515 
    488516      ip3dt = 0 
    489       ivar1t = 0 
    490       ivar2t = 0 
    491       itypvar1   (:) = 0 
    492       itypvar1mpp(:) = 0 
    493  
    494       itypvar2   (:) = 0 
    495       itypvar2mpp(:) = 0 
     517      ivart(:) = 0 
     518      itypvar   (:,:) = 0 
     519      itypvarmpp(:,:) = 0 
    496520 
    497521      ioserrcount = 0 
     
    501525         ji = iprofidx(iindx(jk)) 
    502526 
    503             IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
    504             IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
    505                & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
     527         IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     528         llcycle = .TRUE. 
     529         DO jvar = 1, kvars 
     530            IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 
     531               llcycle = .FALSE. 
     532               EXIT 
     533            ENDIF 
     534         END DO 
     535         IF ( llcycle ) CYCLE 
    506536 
    507537         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  & 
     
    519549 
    520550            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
    521             IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
    522                & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
     551            llcycle = .TRUE. 
     552            DO jvar = 1, kvars 
     553               IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 
     554                  llcycle = .FALSE. 
     555                  EXIT 
     556               ENDIF 
     557            END DO 
     558            IF ( llcycle ) CYCLE 
    523559 
    524560            loop_prof : DO ij = 1, inpfiles(jj)%nlev 
     
    527563                  & CYCLE 
    528564 
    529                IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
    530                   & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
    531  
    532                   llvalprof = .TRUE.  
    533                   EXIT loop_prof 
    534  
    535                ENDIF 
    536  
    537                IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
    538                   & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
    539  
    540                   llvalprof = .TRUE.  
    541                   EXIT loop_prof 
    542  
    543                ENDIF 
     565               DO jvar = 1, kvars 
     566                  IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 
     567                     & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
     568 
     569                     llvalprof = .TRUE.  
     570                     EXIT loop_prof 
     571 
     572                  ENDIF 
     573               END DO 
    544574 
    545575            END DO loop_prof 
     
    573603 
    574604               ! Coordinate search parameters 
    575                profdata%mi  (iprof,1) = inpfiles(jj)%iobsi(ji,1) 
    576                profdata%mj  (iprof,1) = inpfiles(jj)%iobsj(ji,1) 
    577                profdata%mi  (iprof,2) = inpfiles(jj)%iobsi(ji,2) 
    578                profdata%mj  (iprof,2) = inpfiles(jj)%iobsj(ji,2) 
     605               DO jvar = 1, kvars 
     606                  profdata%mi  (iprof,jvar) = inpfiles(jj)%iobsi(ji,jvar) 
     607                  profdata%mj  (iprof,jvar) = inpfiles(jj)%iobsj(ji,jvar) 
     608               END DO 
    579609 
    580610               ! Profile WMO number 
     
    616646                  IF (ldsatt) THEN 
    617647 
    618                      IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
    619                         &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
    620                         &    ldvar1 ) .OR. & 
    621                         & ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
    622                         &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
    623                         &   ldvar2 ) ) THEN 
    624                         ip3dt = ip3dt + 1 
    625                      ELSE 
    626                         CYCLE 
     648                     DO jvar = 1, kvars 
     649                        IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 
     650                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     651                           &    ldvar(jvar) ) ) THEN 
     652                           ip3dt = ip3dt + 1 
     653                           EXIT 
     654                        ELSE IF ( jvar == kvars ) THEN 
     655                           CYCLE loop_p 
     656                        ENDIF 
     657                     END DO 
     658 
     659                  ENDIF 
     660 
     661                  DO jvar = 1, kvars 
     662                   
     663                     IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 
     664                       &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     665                       &    ldvar(jvar) ) .OR. ldsatt ) THEN 
     666 
     667                        IF (ldsatt) THEN 
     668 
     669                           ivart(jvar) = ip3dt 
     670 
     671                        ELSE 
     672 
     673                           ivart(jvar) = ivart(jvar) + 1 
     674 
     675                        ENDIF 
     676 
     677                        ! Depth of jvar observation 
     678                        profdata%var(jvar)%vdep(ivart(jvar)) = & 
     679                           &                inpfiles(jj)%pdep(ij,ji) 
     680 
     681                        ! Depth of jvar observation QC 
     682                        profdata%var(jvar)%idqc(ivart(jvar)) = & 
     683                           &                inpfiles(jj)%idqc(ij,ji) 
     684 
     685                        ! Depth of jvar observation QC flags 
     686                        profdata%var(jvar)%idqcf(:,ivart(jvar)) = & 
     687                           &                inpfiles(jj)%idqcf(:,ij,ji) 
     688 
     689                        ! Profile index 
     690                        profdata%var(jvar)%nvpidx(ivart(jvar)) = iprof 
     691 
     692                        ! Vertical index in original profile 
     693                        profdata%var(jvar)%nvlidx(ivart(jvar)) = ij 
     694 
     695                        ! Profile jvar value 
     696                        IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 
     697                           & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
     698                           profdata%var(jvar)%vobs(ivart(jvar)) = & 
     699                              &                inpfiles(jj)%pob(ij,ji,jvar) 
     700                           IF ( ldmod ) THEN 
     701                              profdata%var(jvar)%vmod(ivart(jvar)) = & 
     702                                 &                inpfiles(jj)%padd(ij,ji,1,jvar) 
     703                           ENDIF 
     704                           IF ( profdata%lclim ) THEN 
     705                               profdata%var(jvar)%vclm(ivart(jvar)) = fbrmdi 
     706                           ENDIF                           
     707                           ! Count number of profile var1 data as function of type 
     708                           itypvar( profdata%ntyp(iprof) + 1, jvar ) = & 
     709                              & itypvar( profdata%ntyp(iprof) + 1, jvar ) + 1 
     710                        ELSE 
     711                           profdata%var(jvar)%vobs(ivart(jvar)) = fbrmdi 
     712                        ENDIF 
     713 
     714                        ! Profile jvar qc 
     715                        profdata%var(jvar)%nvqc(ivart(jvar)) = & 
     716                           & inpfiles(jj)%ivlqc(ij,ji,jvar) 
     717 
     718                        ! Profile jvar qc flags 
     719                        profdata%var(jvar)%nvqcf(:,ivart(jvar)) = & 
     720                           & inpfiles(jj)%ivlqcf(:,ij,ji,jvar) 
     721 
     722                        ! Profile insitu T value 
     723                        IF ( TRIM( inpfiles(jj)%cname(jvar) ) == 'POTM' ) THEN 
     724                           profdata%var(jvar)%vext(ivart(jvar),1) = & 
     725                              &                inpfiles(jj)%pext(ij,ji,1) 
     726                        ENDIF 
     727 
    627728                     ENDIF 
    628  
    629                   ENDIF 
    630  
    631                   IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
    632                     &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
    633                     &    ldvar1 ) .OR. ldsatt ) THEN 
    634  
    635                      IF (ldsatt) THEN 
    636  
    637                         ivar1t = ip3dt 
    638  
    639                      ELSE 
    640  
    641                         ivar1t = ivar1t + 1 
    642  
    643                      ENDIF 
    644  
    645                      ! Depth of var1 observation 
    646                      profdata%var(1)%vdep(ivar1t) = & 
    647                         &                inpfiles(jj)%pdep(ij,ji) 
    648  
    649                      ! Depth of var1 observation QC 
    650                      profdata%var(1)%idqc(ivar1t) = & 
    651                         &                inpfiles(jj)%idqc(ij,ji) 
    652  
    653                      ! Depth of var1 observation QC flags 
    654                      profdata%var(1)%idqcf(:,ivar1t) = & 
    655                         &                inpfiles(jj)%idqcf(:,ij,ji) 
    656  
    657                      ! Profile index 
    658                      profdata%var(1)%nvpidx(ivar1t) = iprof 
    659  
    660                      ! Vertical index in original profile 
    661                      profdata%var(1)%nvlidx(ivar1t) = ij 
    662  
    663                      ! Profile var1 value 
    664                      IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
    665                         & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
    666                         profdata%var(1)%vobs(ivar1t) = & 
    667                            &                inpfiles(jj)%pob(ij,ji,1) 
    668                         IF ( ldmod ) THEN 
    669                            profdata%var(1)%vmod(ivar1t) = & 
    670                               &                inpfiles(jj)%padd(ij,ji,1,1) 
    671                         ENDIF 
    672                         ! Count number of profile var1 data as function of type 
    673                         itypvar1( profdata%ntyp(iprof) + 1 ) = & 
    674                            & itypvar1( profdata%ntyp(iprof) + 1 ) + 1 
    675                      ELSE 
    676                         profdata%var(1)%vobs(ivar1t) = fbrmdi 
    677                      ENDIF 
    678  
    679                      ! Profile var1 qc 
    680                      profdata%var(1)%nvqc(ivar1t) = & 
    681                         & inpfiles(jj)%ivlqc(ij,ji,1) 
    682  
    683                      ! Profile var1 qc flags 
    684                      profdata%var(1)%nvqcf(:,ivar1t) = & 
    685                         & inpfiles(jj)%ivlqcf(:,ij,ji,1) 
    686  
    687                      ! Profile insitu T value 
    688                      IF ( TRIM( inpfiles(jj)%cname(1) ) == 'POTM' ) THEN 
    689                         profdata%var(1)%vext(ivar1t,1) = & 
    690                            &                inpfiles(jj)%pext(ij,ji,1) 
    691                      ENDIF 
    692  
    693                   ENDIF 
    694  
    695                   IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
    696                      &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2)    .AND. & 
    697                      &   ldvar2 ) .OR. ldsatt ) THEN 
    698  
    699                      IF (ldsatt) THEN 
    700  
    701                         ivar2t = ip3dt 
    702  
    703                      ELSE 
    704  
    705                         ivar2t = ivar2t + 1 
    706  
    707                      ENDIF 
    708  
    709                      ! Depth of var2 observation 
    710                      profdata%var(2)%vdep(ivar2t) = & 
    711                         &                inpfiles(jj)%pdep(ij,ji) 
    712  
    713                      ! Depth of var2 observation QC 
    714                      profdata%var(2)%idqc(ivar2t) = & 
    715                         &                inpfiles(jj)%idqc(ij,ji) 
    716  
    717                      ! Depth of var2 observation QC flags 
    718                      profdata%var(2)%idqcf(:,ivar2t) = & 
    719                         &                inpfiles(jj)%idqcf(:,ij,ji) 
    720  
    721                      ! Profile index 
    722                      profdata%var(2)%nvpidx(ivar2t) = iprof 
    723  
    724                      ! Vertical index in original profile 
    725                      profdata%var(2)%nvlidx(ivar2t) = ij 
    726  
    727                      ! Profile var2 value 
    728                   IF (  ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) ) .AND. & 
    729                     &   ( .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2)    )  ) THEN 
    730                         profdata%var(2)%vobs(ivar2t) = & 
    731                            &                inpfiles(jj)%pob(ij,ji,2) 
    732                         IF ( ldmod ) THEN 
    733                            profdata%var(2)%vmod(ivar2t) = & 
    734                               &                inpfiles(jj)%padd(ij,ji,1,2) 
    735                         ENDIF 
    736                         ! Count number of profile var2 data as function of type 
    737                         itypvar2( profdata%ntyp(iprof) + 1 ) = & 
    738                            & itypvar2( profdata%ntyp(iprof) + 1 ) + 1 
    739                      ELSE 
    740                         profdata%var(2)%vobs(ivar2t) = fbrmdi 
    741                      ENDIF 
    742  
    743                      ! Profile var2 qc 
    744                      profdata%var(2)%nvqc(ivar2t) = & 
    745                         & inpfiles(jj)%ivlqc(ij,ji,2) 
    746  
    747                      ! Profile var2 qc flags 
    748                      profdata%var(2)%nvqcf(:,ivar2t) = & 
    749                         & inpfiles(jj)%ivlqcf(:,ij,ji,2) 
    750  
    751                   ENDIF 
     729                   
     730                  END DO 
    752731 
    753732               END DO loop_p 
     
    763742      !----------------------------------------------------------------------- 
    764743 
    765       CALL obs_mpp_sum_integer ( ivar1t0, ivar1tmpp ) 
    766       CALL obs_mpp_sum_integer ( ivar2t0, ivar2tmpp ) 
     744      DO jvar = 1, kvars 
     745         CALL obs_mpp_sum_integer ( ivart0(jvar), ivartmpp(jvar) ) 
     746      END DO 
    767747      CALL obs_mpp_sum_integer ( ip3dt,   ip3dtmpp  ) 
    768748 
    769       CALL obs_mpp_sum_integers( itypvar1, itypvar1mpp, ntyp1770 + 1 ) 
    770       CALL obs_mpp_sum_integers( itypvar2, itypvar2mpp, ntyp1770 + 1 ) 
     749      DO jvar = 1, kvars 
     750         CALL obs_mpp_sum_integers( itypvar(:,jvar), itypvarmpp(:,jvar), ntyp1770 + 1 ) 
     751      END DO 
    771752 
    772753      !----------------------------------------------------------------------- 
     
    778759         WRITE(numout,'(1X,A)') '------------' 
    779760         WRITE(numout,*)  
    780          WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(1) ) 
    781          WRITE(numout,'(1X,A)') '------------------------' 
    782          DO ji = 0, ntyp1770 
    783             IF ( itypvar1mpp(ji+1) > 0 ) THEN 
    784                WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), & 
    785                   & cwmonam1770(ji)(1:52),' = ', & 
    786                   & itypvar1mpp(ji+1) 
    787             ENDIF 
     761         DO jvar = 1, kvars 
     762            WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(jvar) ) 
     763            WRITE(numout,'(1X,A)') '------------------------' 
     764            DO ji = 0, ntyp1770 
     765               IF ( itypvarmpp(ji+1,jvar) > 0 ) THEN 
     766                  WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), & 
     767                     & cwmonam1770(ji)(1:52),' = ', & 
     768                     & itypvarmpp(ji+1,jvar) 
     769               ENDIF 
     770            END DO 
     771            WRITE(numout,'(1X,A)') & 
     772               & '---------------------------------------------------------------' 
     773            WRITE(numout,'(1X,A55,I8)') & 
     774               & 'Total profile data for variable '//TRIM( profdata%cvars(jvar) )// & 
     775               & '             = ', ivartmpp(jvar) 
     776            WRITE(numout,'(1X,A)') & 
     777               & '---------------------------------------------------------------' 
     778            WRITE(numout,*)  
    788779         END DO 
    789          WRITE(numout,'(1X,A)') & 
    790             & '---------------------------------------------------------------' 
    791          WRITE(numout,'(1X,A55,I8)') & 
    792             & 'Total profile data for variable '//TRIM( profdata%cvars(1) )// & 
    793             & '             = ', ivar1tmpp 
    794          WRITE(numout,'(1X,A)') & 
    795             & '---------------------------------------------------------------' 
    796          WRITE(numout,*)  
    797          WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(2) ) 
    798          WRITE(numout,'(1X,A)') '------------------------' 
    799          DO ji = 0, ntyp1770 
    800             IF ( itypvar2mpp(ji+1) > 0 ) THEN 
    801                WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), & 
    802                   & cwmonam1770(ji)(1:52),' = ', & 
    803                   & itypvar2mpp(ji+1) 
    804             ENDIF 
     780      ENDIF 
     781 
     782      IF (ldsatt) THEN 
     783         profdata%nvprot(:)    = ip3dt 
     784         profdata%nvprotmpp(:) = ip3dtmpp 
     785      ELSE 
     786         DO jvar = 1, kvars 
     787            profdata%nvprot(jvar)    = ivart(jvar) 
     788            profdata%nvprotmpp(jvar) = ivartmpp(jvar) 
    805789         END DO 
    806          WRITE(numout,'(1X,A)') & 
    807             & '---------------------------------------------------------------' 
    808          WRITE(numout,'(1X,A55,I8)') & 
    809             & 'Total profile data for variable '//TRIM( profdata%cvars(2) )// & 
    810             & '             = ', ivar2tmpp 
    811          WRITE(numout,'(1X,A)') & 
    812             & '---------------------------------------------------------------' 
    813          WRITE(numout,*)  
    814       ENDIF 
    815  
    816       IF (ldsatt) THEN 
    817          profdata%nvprot(1)    = ip3dt 
    818          profdata%nvprot(2)    = ip3dt 
    819          profdata%nvprotmpp(1) = ip3dtmpp 
    820          profdata%nvprotmpp(2) = ip3dtmpp 
    821       ELSE 
    822          profdata%nvprot(1)    = ivar1t 
    823          profdata%nvprot(2)    = ivar2t 
    824          profdata%nvprotmpp(1) = ivar1tmpp 
    825          profdata%nvprotmpp(2) = ivar2tmpp 
    826790      ENDIF 
    827791      profdata%nprof        = iprof 
     
    830794      ! Model level search 
    831795      !----------------------------------------------------------------------- 
    832       IF ( ldvar1 ) THEN 
    833          CALL obs_level_search( jpk, gdept_1d, & 
    834             & profdata%nvprot(1), profdata%var(1)%vdep, & 
    835             & profdata%var(1)%mvk ) 
    836       ENDIF 
    837       IF ( ldvar2 ) THEN 
    838          CALL obs_level_search( jpk, gdept_1d, & 
    839             & profdata%nvprot(2), profdata%var(2)%vdep, & 
    840             & profdata%var(2)%mvk ) 
    841       ENDIF 
     796      DO jvar = 1, kvars 
     797         IF ( ldvar(jvar) ) THEN 
     798            CALL obs_level_search( jpk, gdept_1d, & 
     799               & profdata%nvprot(jvar), profdata%var(jvar)%vdep, & 
     800               & profdata%var(jvar)%mvk ) 
     801         ENDIF 
     802      END DO 
    842803 
    843804      !----------------------------------------------------------------------- 
     
    852813      ! Deallocate temporary data 
    853814      !----------------------------------------------------------------------- 
    854       DEALLOCATE( ifileidx, iprofidx, zdat, clvars ) 
     815      DEALLOCATE( ifileidx, iprofidx, zdat, clvarsin ) 
    855816 
    856817      !----------------------------------------------------------------------- 
  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/OBS/obs_read_surf.F90

    r13226 r13758  
    4040   SUBROUTINE obs_rea_surf( surfdata, knumfiles, cdfilenames, & 
    4141      &                     kvars, kextr, kstp, ddobsini, ddobsend, & 
    42       &                     ldignmis, ldmod, ldnightav ) 
     42      &                     ldignmis, ldmod, ldnightav, cdvars ) 
    4343      !!--------------------------------------------------------------------- 
    4444      !! 
     
    7373      REAL(dp), INTENT(IN) :: ddobsini   ! Obs. ini time in YYYYMMDD.HHMMSS 
    7474      REAL(dp), INTENT(IN) :: ddobsend   ! Obs. end time in YYYYMMDD.HHMMSS 
     75      CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars 
    7576 
    7677      !! * Local declarations 
    7778      CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_surf' 
    7879      CHARACTER(len=8) :: clrefdate 
    79       CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 
     80      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvarsin 
    8081      INTEGER :: ji 
    8182      INTEGER :: jj 
     
    178179               &                ldgrid = .TRUE. ) 
    179180 
     181            IF ( inpfiles(jj)%nvar /= kvars ) THEN 
     182               CALL ctl_stop( 'Feedback format error: ', & 
     183                  &           ' unexpected number of vars in feedback file' ) 
     184            ENDIF 
     185 
    180186            IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN 
    181187               CALL ctl_stop( 'Model not in input data' ) 
     
    184190 
    185191            IF ( jj == 1 ) THEN 
    186                ALLOCATE( clvars( inpfiles(jj)%nvar ) ) 
     192               ALLOCATE( clvarsin( inpfiles(jj)%nvar ) ) 
    187193               DO ji = 1, inpfiles(jj)%nvar 
    188                  clvars(ji) = inpfiles(jj)%cname(ji) 
     194                 clvarsin(ji) = inpfiles(jj)%cname(ji) 
     195                 IF ( clvarsin(ji) /= cdvars(ji) ) THEN 
     196                    CALL ctl_stop( 'Feedback file variables do not match', & 
     197                        &           ' expected variable names for this type' ) 
     198                 ENDIF 
    189199               END DO 
    190200            ELSE 
    191201               DO ji = 1, inpfiles(jj)%nvar 
    192                   IF ( inpfiles(jj)%cname(ji) /= clvars(ji) ) THEN 
     202                  IF ( inpfiles(jj)%cname(ji) /= clvarsin(ji) ) THEN 
    193203                     CALL ctl_stop( 'Feedback file variables not consistent', & 
    194204                        &           ' with previous files for this type' ) 
     
    347357      iobs = 0 
    348358 
    349       surfdata%cvars(:)  = clvars(:) 
     359      surfdata%cvars(:)  = clvarsin(:) 
    350360 
    351361      ityp   (:) = 0 
     
    480490      ! Deallocate temporary data 
    481491      !----------------------------------------------------------------------- 
    482       DEALLOCATE( ifileidx, isurfidx, zdat, clvars ) 
     492      DEALLOCATE( ifileidx, isurfidx, zdat, clvarsin ) 
    483493 
    484494      !----------------------------------------------------------------------- 
  • NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OCE/OBS/obs_write.F90

    r12933 r13758  
    8686      CHARACTER(LEN=40) :: clfname 
    8787      CHARACTER(LEN=10) :: clfiletype 
     88      CHARACTER(LEN=ilenlong) :: cllongname  ! Long name of variable 
     89      CHARACTER(LEN=ilenunit) :: clunits     ! Units of variable 
     90      CHARACTER(LEN=ilengrid) :: clgrid      ! Grid of variable 
    8891      CHARACTER(LEN=12) :: clfmt            ! writing format 
    8992      INTEGER :: idg                        ! number of digits 
     
    115118      ! Find maximum level 
    116119      ilevel = 0 
    117       DO jvar = 1, 2 
     120      DO jvar = 1, profdata%nvar 
    118121         ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) ) 
    119122      END DO 
     
    180183 
    181184      END SELECT 
     185       
     186      IF ( ( TRIM(profdata%cvars(1)) /= 'POTM' ) .AND. & 
     187         & ( TRIM(profdata%cvars(1)) /= 'UVEL' ) ) THEN 
     188         CALL alloc_obfbdata( fbdata, 1, profdata%nprof, ilevel, & 
     189            &                 1 + iadd_clm + iadd, iext, .TRUE. ) 
     190         fbdata%cname(1)      = profdata%cvars(1) 
     191         fbdata%coblong(1)    = cllongname 
     192         fbdata%cobunit(1)    = clunits 
     193         fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(cllongname) 
     194         fbdata%caddunit(1,1) = clunits 
     195         IF ( profdata%lclim ) THEN 
     196            fbdata%caddlong(2,1) = 'Climatological interpolated ' // TRIM(cllongname) 
     197            fbdata%caddunit(2,1) = clunits 
     198         ENDIF          
     199         fbdata%cgrid(:)      = clgrid 
     200         DO je = 1, iext 
     201            fbdata%cextname(je) = pext%cdname(je) 
     202            fbdata%cextlong(je) = pext%cdlong(je,1) 
     203            fbdata%cextunit(je) = pext%cdunit(je,1) 
     204         END DO 
     205         DO ja = 1, iadd 
     206            fbdata%caddname(1+iadd_clm+ja) = padd%cdname(ja) 
     207            fbdata%caddlong(1+iadd_clm+ja,1) = padd%cdlong(ja,1) 
     208            fbdata%caddunit(1+iadd_clm+ja,1) = padd%cdunit(ja,1) 
     209         END DO 
     210      ENDIF 
    182211 
    183212      fbdata%caddname(1)   = 'Hx' 
     
    234263            &           krefdate = 19500101 ) 
    235264         ! Reform the profiles arrays for output 
    236          DO jvar = 1, 2 
     265         DO jvar = 1, profdata%nvar 
    237266            DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar) 
    238267               ik = profdata%var(jvar)%nvlidx(jk) 
     
    329358      CHARACTER(LEN=40) :: clfname         ! netCDF filename 
    330359      CHARACTER(LEN=10) :: clfiletype 
     360      CHARACTER(LEN=ilenlong), DIMENSION(surfdata%nvar) :: cllongname  ! Long name of variable 
     361      CHARACTER(LEN=ilenunit), DIMENSION(surfdata%nvar) :: clunits     ! Units of variable 
     362      CHARACTER(LEN=ilengrid), DIMENSION(surfdata%nvar) :: clgrid      ! Grid of variable 
    331363      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf' 
    332364      CHARACTER(LEN=12) :: clfmt           ! writing format 
     
    354386      SELECT CASE ( TRIM(surfdata%cvars(1)) ) 
    355387      CASE('SLA') 
     388          
     389         ! SLA needs special treatment because of MDT, so is all done here 
     390         ! Other variables are done more generically 
     391         ! No climatology for SLA, MDT is our best estimate of that and is already output. 
    356392 
    357393         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
     
    384420      CASE('SST') 
    385421 
    386          CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
    387             &                 1 + iadd, iext, .TRUE. ) 
    388  
    389          clfiletype = 'sstfb' 
    390          fbdata%cname(1)      = surfdata%cvars(1) 
    391          fbdata%coblong(1)    = 'Sea surface temperature' 
    392          fbdata%cobunit(1)    = 'Degree centigrade' 
     422         clfiletype    = 'sstfb' 
     423         cllongname(1) = 'Sea surface temperature' 
     424         clunits(1)    = 'Degree centigrade' 
     425         clgrid(1)     = 'T' 
     426          
     427      CASE('ICECONC') 
     428 
     429         clfiletype    = 'sicfb' 
     430         cllongname(1) = 'Sea ice concentration' 
     431         clunits(1)    = 'Fraction' 
     432         clgrid(1)     = 'T' 
     433 
     434      CASE('SSS') 
     435 
     436         clfiletype    = 'sssfb' 
     437         cllongname(1) = 'Sea surface salinity' 
     438         clunits(1)    = 'psu' 
     439         clgrid(1)     = 'T' 
     440         END DO 
     441 
     442      CASE DEFAULT 
     443 
     444         CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' ) 
     445 
     446      END SELECT 
     447 
     448      ! SLA needs special treatment because of MDT, so is done above 
     449      ! Remaining variables treated more generically 
     450 
     451      IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN 
     452       
     453         CALL alloc_obfbdata( fbdata, surfdata%nvar, surfdata%nsurf, 1, & 
     454            &                 1 + iadd_std + iadd_clm + iadd, iext, .TRUE. ) 
     455 
     456         DO jv = 1, surfdata%nvar 
     457            fbdata%cname(jv)      = surfdata%cvars(jv) 
     458            fbdata%coblong(jv)    = cllongname(jv) 
     459            fbdata%cobunit(jv)    = clunits(jv) 
     460         END DO 
    393461         DO je = 1, iext 
    394462            fbdata%cextname(je) = pext%cdname(je) 
     
    396464            fbdata%cextunit(je) = pext%cdunit(je,1) 
    397465         END DO 
    398          fbdata%caddlong(1,1) = 'Model interpolated SST' 
    399          fbdata%caddunit(1,1) = 'Degree centigrade' 
    400          fbdata%cgrid(1)      = 'T' 
     466         DO jv = 1, surfdata%nvar          
     467            IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN 
     468               fbdata%caddlong(1,jv) = 'Model interpolated ICE' 
     469            ELSE 
     470               fbdata%caddlong(1,jv) = 'Model interpolated ' // TRIM(surfdata%cvars(jv)) 
     471            ENDIF 
     472            fbdata%caddunit(1,jv) = clunits(jv) 
     473            fbdata%cgrid(jv)      = clgrid(jv) 
     474         END DO             
    401475         DO ja = 1, iadd 
    402             fbdata%caddname(1+ja) = padd%cdname(ja) 
    403             fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
    404             fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
    405          END DO 
    406  
    407       CASE('ICECONC') 
    408  
    409          CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
    410             &                 1 + iadd, iext, .TRUE. ) 
    411  
    412          clfiletype = 'sicfb' 
    413          fbdata%cname(1)      = surfdata%cvars(1) 
    414          fbdata%coblong(1)    = 'Sea ice' 
    415          fbdata%cobunit(1)    = 'Fraction' 
    416          DO je = 1, iext 
    417             fbdata%cextname(je) = pext%cdname(je) 
    418             fbdata%cextlong(je) = pext%cdlong(je,1) 
    419             fbdata%cextunit(je) = pext%cdunit(je,1) 
    420          END DO 
    421          fbdata%caddlong(1,1) = 'Model interpolated ICE' 
    422          fbdata%caddunit(1,1) = 'Fraction' 
    423          fbdata%cgrid(1)      = 'T' 
    424          DO ja = 1, iadd 
    425             fbdata%caddname(1+ja) = padd%cdname(ja) 
    426             fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
    427             fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
    428          END DO 
    429  
    430       CASE('SSS') 
    431  
    432          CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
    433             &                 1 + iadd, iext, .TRUE. ) 
    434  
    435          clfiletype = 'sssfb' 
    436          fbdata%cname(1)      = surfdata%cvars(1) 
    437          fbdata%coblong(1)    = 'Sea surface salinity' 
    438          fbdata%cobunit(1)    = 'psu' 
    439          DO je = 1, iext 
    440             fbdata%cextname(je) = pext%cdname(je) 
    441             fbdata%cextlong(je) = pext%cdlong(je,1) 
    442             fbdata%cextunit(je) = pext%cdunit(je,1) 
    443          END DO 
    444          fbdata%caddlong(1,1) = 'Model interpolated SSS' 
    445          fbdata%caddunit(1,1) = 'psu' 
    446          fbdata%cgrid(1)      = 'T' 
    447          DO ja = 1, iadd 
    448             fbdata%caddname(1+ja) = padd%cdname(ja) 
    449             fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
    450             fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
    451          END DO 
    452  
    453       CASE DEFAULT 
    454  
    455          CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' ) 
    456  
    457       END SELECT 
     476            fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm+ja) = padd%cdname(ja) 
     477            DO jv = 1, surfdata%nvar                      
     478               fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = padd%cdlong(ja,jv) 
     479               fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = padd%cdunit(ja,jv) 
     480            END DO 
     481         END DO 
     482      ENDIF 
    458483 
    459484      fbdata%caddname(1)   = 'Hx' 
Note: See TracChangeset for help on using the changeset viewer.