Changeset 14056


Ignore:
Timestamp:
2020-12-03T15:08:29+01:00 (3 months ago)
Author:
ayoung
Message:

Adding branch for ticket #2567 to trunk.

Location:
NEMO/trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/cfgs/SHARED/namelist_ref

    r14053 r14056  
    14001400   ln_sstnight = .false.             ! Logical switch for calculating night-time average for SST obs 
    14011401   ln_bound_reject  = .false.        ! Logical to remove obs near boundaries in LAMs. 
     1402   ln_default_fp_indegs = .true.     ! Logical: T=> averaging footprint is in degrees, F=> in metres 
    14021403   ln_sla_fp_indegs = .true.         ! Logical for SLA: T=> averaging footprint is in degrees, F=> in metres 
    14031404   ln_sst_fp_indegs = .true.         ! Logical for SST: T=> averaging footprint is in degrees, F=> in metres 
     
    14151416   cn_gridsearchfile ='gridsearch.nc' ! Grid search file name 
    14161417   rn_gridsearchres = 0.5            ! Grid search resolution 
     1418   rn_default_avglamscl = 0.         ! Default E/W diameter of observation footprint (metres/degrees) 
     1419   rn_default_avgphiscl = 0.         ! Default N/S diameter of observation footprint (metres/degrees) 
    14171420   rn_mdtcorr  = 1.61                ! MDT  correction 
    14181421   rn_mdtcutoff = 65.0               ! MDT cutoff for computed correction 
     
    14281431   rn_sic_avgphiscl = 0.             ! N/S diameter of SIC observation footprint (metres/degrees) 
    14291432   nn_1dint = 0                      ! Type of vertical interpolation method 
    1430    nn_2dint = 0                      ! Default horizontal interpolation method 
     1433   nn_2dint_default = 0              ! Default horizontal interpolation method 
    14311434   nn_2dint_sla = 0                  ! Horizontal interpolation method for SLA 
    14321435   nn_2dint_sst = 0                  ! Horizontal interpolation method for SST 
  • NEMO/trunk/doc/namelists/namobs

    r11703 r14056  
    2020   ln_sstnight = .false.             ! Logical switch for calculating night-time average for SST obs 
    2121   ln_bound_reject  = .false.        ! Logical to remove obs near boundaries in LAMs. 
     22   ln_default_fp_indegs = .true.     ! Logical: T=> averaging footprint is in degrees, F=> in metres 
    2223   ln_sla_fp_indegs = .true.         ! Logical for SLA: T=> averaging footprint is in degrees, F=> in metres 
    2324   ln_sst_fp_indegs = .true.         ! Logical for SST: T=> averaging footprint is in degrees, F=> in metres 
     
    3940   rn_dobsini  = 00010101.000000     ! Initial date in window YYYYMMDD.HHMMSS 
    4041   rn_dobsend  = 00010102.000000     ! Final date in window YYYYMMDD.HHMMSS 
     42   rn_default_avglamscl = 0.         ! Default E/W diameter of observation footprint (metres/degrees) 
     43   rn_default_avgphiscl = 0.         ! Default N/S diameter of observation footprint (metres/degrees) 
    4144   rn_sla_avglamscl = 0.             ! E/W diameter of SLA observation footprint (metres/degrees) 
    4245   rn_sla_avgphiscl = 0.             ! N/S diameter of SLA observation footprint (metres/degrees) 
     
    4851   rn_sic_avgphiscl = 0.             ! N/S diameter of SIC observation footprint (metres/degrees) 
    4952   nn_1dint = 0                      ! Type of vertical interpolation method 
    50    nn_2dint = 0                      ! Default horizontal interpolation method 
     53   nn_2dint_default = 0              ! Default horizontal interpolation method 
    5154   nn_2dint_sla = 0                  ! Horizontal interpolation method for SLA 
    5255   nn_2dint_sst = 0                  ! Horizontal interpolation method for SST 
  • NEMO/trunk/src/OCE/OBS/diaobs.F90

    r13216 r14056  
    5757   PUBLIC calc_date        ! Compute the date of a timestep 
    5858 
    59    LOGICAL, PUBLIC :: ln_diaobs          !: Logical switch for the obs operator 
    60    LOGICAL         :: ln_sstnight        !  Logical switch for night mean SST obs 
    61    LOGICAL         :: ln_sla_fp_indegs   !  T=> SLA obs footprint size specified in degrees, F=> in metres 
    62    LOGICAL         :: ln_sst_fp_indegs   !  T=> SST obs footprint size specified in degrees, F=> in metres 
    63    LOGICAL         :: ln_sss_fp_indegs   !  T=> SSS obs footprint size specified in degrees, F=> in metres 
    64    LOGICAL         :: ln_sic_fp_indegs   !  T=> sea-ice obs footprint size specified in degrees, F=> in metres 
    65  
    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  
     59   LOGICAL, PUBLIC :: ln_diaobs            !: Logical switch for the obs operator 
     60   LOGICAL         :: ln_sstnight          !  Logical switch for night mean SST obs 
     61   LOGICAL         :: ln_default_fp_indegs !  T=> Default obs footprint size specified in degrees, F=> in metres 
     62   LOGICAL         :: ln_sla_fp_indegs     !  T=> SLA obs footprint size specified in degrees, F=> in metres 
     63   LOGICAL         :: ln_sst_fp_indegs     !  T=> SST obs footprint size specified in degrees, F=> in metres 
     64   LOGICAL         :: ln_sss_fp_indegs     !  T=> SSS obs footprint size specified in degrees, F=> in metres 
     65   LOGICAL         :: ln_sic_fp_indegs     !  T=> sea-ice obs footprint size specified in degrees, F=> in metres 
     66 
     67   REAL(wp) ::   rn_default_avglamscl      ! E/W diameter of SLA observation footprint (metres) 
     68   REAL(wp) ::   rn_default_avgphiscl      ! N/S diameter of SLA observation footprint (metre 
     69   REAL(wp) ::   rn_sla_avglamscl          ! E/W diameter of SLA observation footprint (metres) 
     70   REAL(wp) ::   rn_sla_avgphiscl          ! N/S diameter of SLA observation footprint (metres) 
     71   REAL(wp) ::   rn_sst_avglamscl          ! E/W diameter of SST observation footprint (metres) 
     72   REAL(wp) ::   rn_sst_avgphiscl          ! N/S diameter of SST observation footprint (metres) 
     73   REAL(wp) ::   rn_sss_avglamscl          ! E/W diameter of SSS observation footprint (metres) 
     74   REAL(wp) ::   rn_sss_avgphiscl          ! N/S diameter of SSS observation footprint (metres) 
     75   REAL(wp) ::   rn_sic_avglamscl          ! E/W diameter of sea-ice observation footprint (metres) 
     76   REAL(wp) ::   rn_sic_avgphiscl          ! N/S diameter of sea-ice observation footprint (metres) 
     77 
     78   INTEGER :: nn_1dint                     ! Vertical interpolation method 
     79   INTEGER :: nn_2dint_default             ! Default horizontal interpolation method 
     80   INTEGER :: nn_2dint_sla                 ! SLA horizontal interpolation method  
     81   INTEGER :: nn_2dint_sst                 ! SST horizontal interpolation method  
     82   INTEGER :: nn_2dint_sss                 ! SSS horizontal interpolation method  
     83   INTEGER :: nn_2dint_sic                 ! Seaice horizontal interpolation method  
    8184   INTEGER, DIMENSION(imaxavtypes) ::   nn_profdavtypes   ! Profile data types representing a daily average 
    8285   INTEGER :: nproftypes     ! Number of profile obs types 
     
    9497   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) ::   profdataqc   !: Profile data after quality control 
    9598 
    96    CHARACTER(len=lca), PUBLIC, DIMENSION(:), ALLOCATABLE ::   cobstypesprof, cobstypessurf   !: Profile & surface obs types 
     99   CHARACTER(len=8), PUBLIC, DIMENSION(:), ALLOCATABLE ::   cobstypesprof, cobstypessurf   !: Profile & surface obs types 
    97100 
    98101   !!---------------------------------------------------------------------- 
     
    121124      INTEGER :: jvar            ! Counter for variables 
    122125      INTEGER :: jfile           ! Counter for files 
    123       INTEGER :: jnumsstbias 
     126      INTEGER :: jnumsstbias     ! Number of SST bias files to read and apply 
     127      INTEGER :: n2dint_type     ! Local version of nn_2dint* 
    124128      ! 
    125129      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 
     
    130134         & cn_sicfbfiles, &      ! Seaice concentration input filenames 
    131135         & cn_velfbfiles, &      ! Velocity profile input filenames 
    132          & cn_sstbiasfiles      ! SST bias input filenames 
     136         & cn_sstbiasfiles       ! SST bias input filenames 
    133137      CHARACTER(LEN=128) :: & 
    134138         & cn_altbiasfile        ! Altimeter bias input filename 
     
    136140         & clproffiles, &        ! Profile filenames 
    137141         & clsurffiles           ! Surface filenames 
     142      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: & 
     143         & clvars                ! Expected variable names 
    138144         ! 
    139145      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles 
     
    150156      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs 
    151157      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 
     158      LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs 
     159      LOGICAL :: ltype_night     ! Local version of ln_sstnight (false for other variables) 
     160      LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar   ! Logical for profile variable read 
    154161      LOGICAL, DIMENSION(jpmaxnfiles) :: lmask ! Used for finding number of sstbias files 
    155162      ! 
    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 
     163      REAL(dp) :: rn_dobsini      ! Obs window start date YYYYMMDD.HHMMSS 
     164      REAL(dp) :: rn_dobsend      ! Obs window end date   YYYYMMDD.HHMMSS 
     165      REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl 
     166      REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl 
     167      REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: zglam   ! Model longitudes for profile variables 
     168      REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: zgphi   ! Model latitudes  for profile variables 
     169      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zmask   ! Model land/sea mask associated with variables 
    161170      !! 
    162171      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              & 
     
    165174         &            ln_grid_global, ln_grid_search_lookup,          & 
    166175         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          & 
    167          &            ln_sstnight,                                    & 
     176         &            ln_sstnight, ln_default_fp_indegs,              & 
    168177         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
    169178         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             & 
     
    174183         &            cn_gridsearchfile, rn_gridsearchres,            & 
    175184         &            rn_dobsini, rn_dobsend,                         & 
     185         &            rn_default_avglamscl, rn_default_avgphiscl,     & 
    176186         &            rn_sla_avglamscl, rn_sla_avgphiscl,             & 
    177187         &            rn_sst_avglamscl, rn_sst_avgphiscl,             & 
    178188         &            rn_sss_avglamscl, rn_sss_avgphiscl,             & 
    179189         &            rn_sic_avglamscl, rn_sic_avgphiscl,             & 
    180          &            nn_1dint, nn_2dint,                             & 
     190         &            nn_1dint, nn_2dint_default,                     & 
    181191         &            nn_2dint_sla, nn_2dint_sst,                     & 
    182192         &            nn_2dint_sss, nn_2dint_sic,                     & 
     
    234244         WRITE(numout,*) '      Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend 
    235245         WRITE(numout,*) '      Type of vertical interpolation method                  nn_1dint = ', nn_1dint 
    236          WRITE(numout,*) '      Type of horizontal interpolation method                nn_2dint = ', nn_2dint 
     246         WRITE(numout,*) '      Default horizontal interpolation method        nn_2dint_default = ', nn_2dint_default 
     247         WRITE(numout,*) '      Type of horizontal interpolation method for SLA    nn_2dint_sla = ', nn_2dint_sla 
     248         WRITE(numout,*) '      Type of horizontal interpolation method for SST    nn_2dint_sst = ', nn_2dint_sst 
     249         WRITE(numout,*) '      Type of horizontal interpolation method for SSS    nn_2dint_sss = ', nn_2dint_sss         
     250         WRITE(numout,*) '      Type of horizontal interpolation method for SIC    nn_2dint_sic = ', nn_2dint_sic 
     251         WRITE(numout,*) '      Default E/W diameter of obs footprint      rn_default_avglamscl = ', rn_default_avglamscl 
     252         WRITE(numout,*) '      Default N/S diameter of obs footprint      rn_default_avgphiscl = ', rn_default_avgphiscl 
     253         WRITE(numout,*) '      Default obs footprint in deg [T] or m [F]  ln_default_fp_indegs = ', ln_default_fp_indegs 
     254         WRITE(numout,*) '      SLA E/W diameter of obs footprint              rn_sla_avglamscl = ', rn_sla_avglamscl 
     255         WRITE(numout,*) '      SLA N/S diameter of obs footprint              rn_sla_avgphiscl = ', rn_sla_avgphiscl 
     256         WRITE(numout,*) '      SLA obs footprint in deg [T] or m [F]          ln_sla_fp_indegs = ', ln_sla_fp_indegs 
     257         WRITE(numout,*) '      SST E/W diameter of obs footprint              rn_sst_avglamscl = ', rn_sst_avglamscl 
     258         WRITE(numout,*) '      SST N/S diameter of obs footprint              rn_sst_avgphiscl = ', rn_sst_avgphiscl 
     259         WRITE(numout,*) '      SST obs footprint in deg [T] or m [F]          ln_sst_fp_indegs = ', ln_sst_fp_indegs 
     260         WRITE(numout,*) '      SIC E/W diameter of obs footprint              rn_sic_avglamscl = ', rn_sic_avglamscl 
     261         WRITE(numout,*) '      SIC N/S diameter of obs footprint              rn_sic_avgphiscl = ', rn_sic_avgphiscl 
     262         WRITE(numout,*) '      SIC obs footprint in deg [T] or m [F]          ln_sic_fp_indegs = ', ln_sic_fp_indegs 
    237263         WRITE(numout,*) '      Rejection of observations near land switch               ln_nea = ', ln_nea 
    238264         WRITE(numout,*) '      Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject 
     
    278304         IF( ln_t3d .OR. ln_s3d ) THEN 
    279305            jtype = jtype + 1 
    280             CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof  ', & 
    281                &                   cn_profbfiles, ifilesprof, cobstypesprof, clproffiles ) 
     306            cobstypesprof(jtype) = 'prof' 
     307            clproffiles(jtype,:) = cn_profbfiles 
    282308         ENDIF 
    283309         IF( ln_vel3d ) THEN 
    284310            jtype = jtype + 1 
    285             CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel   ', & 
    286                &                   cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles ) 
     311            cobstypesprof(jtype) = 'vel' 
     312            clproffiles(jtype,:) = cn_velfbfiles 
    287313         ENDIF 
     314         ! 
     315         CALL obs_settypefiles( nproftypes, jpmaxnfiles, ifilesprof, cobstypesprof, clproffiles ) 
    288316         ! 
    289317      ENDIF 
     
    303331         IF( ln_sla ) THEN 
    304332            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 ) 
     333            cobstypessurf(jtype) = 'sla' 
     334            clsurffiles(jtype,:) = cn_slafbfiles 
    313335         ENDIF 
    314336         IF( ln_sst ) THEN 
    315337            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 ) 
     338            cobstypessurf(jtype) = 'sst' 
     339            clsurffiles(jtype,:) = cn_sstfbfiles 
    324340         ENDIF 
    325341#if defined key_si3 || defined key_cice 
    326342         IF( ln_sic ) THEN 
    327343            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 ) 
     344            cobstypessurf(jtype) = 'sic' 
     345            clsurffiles(jtype,:) = cn_sicfbfiles 
    336346         ENDIF 
    337347#endif 
    338348         IF( ln_sss ) THEN 
    339349            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 ) 
     350            cobstypessurf(jtype) = 'sss' 
     351            clsurffiles(jtype,:) = cn_sssfbfiles 
    348352         ENDIF 
     353         ! 
     354         CALL obs_settypefiles( nsurftypes, jpmaxnfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     355 
     356         DO jtype = 1, nsurftypes 
     357 
     358            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
     359               IF ( nn_2dint_sla == -1 ) THEN 
     360                  n2dint_type  = nn_2dint_default 
     361               ELSE 
     362                  n2dint_type  = nn_2dint_sla 
     363               ENDIF 
     364               ztype_avglamscl = rn_sla_avglamscl 
     365               ztype_avgphiscl = rn_sla_avgphiscl 
     366               ltype_fp_indegs = ln_sla_fp_indegs 
     367               ltype_night     = .FALSE. 
     368            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN 
     369               IF ( nn_2dint_sst == -1 ) THEN 
     370                  n2dint_type  = nn_2dint_default 
     371               ELSE 
     372                  n2dint_type  = nn_2dint_sst 
     373               ENDIF 
     374               ztype_avglamscl = rn_sst_avglamscl 
     375               ztype_avgphiscl = rn_sst_avgphiscl 
     376               ltype_fp_indegs = ln_sst_fp_indegs 
     377               ltype_night     = ln_sstnight 
     378            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN 
     379               IF ( nn_2dint_sic == -1 ) THEN 
     380                  n2dint_type  = nn_2dint_default 
     381               ELSE 
     382                  n2dint_type  = nn_2dint_sic 
     383               ENDIF 
     384               ztype_avglamscl = rn_sic_avglamscl 
     385               ztype_avgphiscl = rn_sic_avgphiscl 
     386               ltype_fp_indegs = ln_sic_fp_indegs 
     387               ltype_night     = .FALSE. 
     388            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 
     389               IF ( nn_2dint_sss == -1 ) THEN 
     390                  n2dint_type  = nn_2dint_default 
     391               ELSE 
     392                  n2dint_type  = nn_2dint_sss 
     393               ENDIF 
     394               ztype_avglamscl = rn_sss_avglamscl 
     395               ztype_avgphiscl = rn_sss_avgphiscl 
     396               ltype_fp_indegs = ln_sss_fp_indegs 
     397               ltype_night     = .FALSE. 
     398            ELSE 
     399               n2dint_type     = nn_2dint_default 
     400               ztype_avglamscl = rn_default_avglamscl 
     401               ztype_avgphiscl = rn_default_avgphiscl 
     402               ltype_fp_indegs = ln_default_fp_indegs 
     403               ltype_night     = .FALSE. 
     404            ENDIF 
     405             
     406            CALL obs_setinterpopts( nsurftypes, jtype, TRIM(cobstypessurf(jtype)), & 
     407               &                    nn_2dint_default, n2dint_type,                 & 
     408               &                    ztype_avglamscl, ztype_avgphiscl,              & 
     409               &                    ltype_fp_indegs, ltype_night,                  & 
     410               &                    n2dintsurf, zavglamscl, zavgphiscl,            & 
     411               &                    lfpindegs, llnightav ) 
     412 
     413         END DO 
    349414         ! 
    350415      ENDIF 
     
    368433      ENDIF 
    369434      ! 
    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') 
     435      IF( nn_2dint_default < 0  .OR.  nn_2dint_default > 6  ) THEN 
     436         CALL ctl_stop('dia_obs_init: Choice of default horizontal (2D) interpolation method is not available') 
    372437      ENDIF 
    373438      ! 
     
    388453         DO jtype = 1, nproftypes 
    389454            ! 
    390             nvarsprof(jtype) = 2 
    391455            IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN 
    392                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 
     456               nvarsprof(jtype) = 2 
     457               nextrprof(jtype) = 1              
     458               ALLOCATE( llvar (nvarsprof(jtype)) ) 
     459               ALLOCATE( clvars(nvarsprof(jtype)) ) 
     460               ALLOCATE( zglam(jpi, jpj,      nvarsprof(jtype)) ) 
     461               ALLOCATE( zgphi(jpi, jpj,      nvarsprof(jtype)) ) 
     462               ALLOCATE( zmask(jpi, jpj, jpk, nvarsprof(jtype)) ) 
     463               llvar(1)       = ln_t3d 
     464               llvar(2)       = ln_s3d 
     465               clvars(1)      = 'POTM' 
     466               clvars(2)      = 'PSAL' 
     467               zglam(:,:,1)   = glamt(:,:) 
     468               zglam(:,:,2)   = glamt(:,:) 
     469               zgphi(:,:,1)   = gphit(:,:) 
     470               zgphi(:,:,2)   = gphit(:,:) 
     471               zmask(:,:,:,1) = tmask(:,:,:) 
     472               zmask(:,:,:,2) = tmask(:,:,:) 
     473            ELSE IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN 
     474               nvarsprof(jtype) = 2 
    403475               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 
     476               ALLOCATE( llvar (nvarsprof(jtype)) ) 
     477               ALLOCATE( clvars(nvarsprof(jtype)) ) 
     478               ALLOCATE( zglam(jpi, jpj,      nvarsprof(jtype)) ) 
     479               ALLOCATE( zgphi(jpi, jpj,      nvarsprof(jtype)) ) 
     480               ALLOCATE( zmask(jpi, jpj, jpk, nvarsprof(jtype)) ) 
     481               llvar(1)       = ln_vel3d 
     482               llvar(2)       = ln_vel3d 
     483               clvars(1)      = 'UVEL' 
     484               clvars(2)      = 'VVEL' 
     485               zglam(:,:,1)   = glamu(:,:) 
     486               zglam(:,:,2)   = glamv(:,:) 
     487               zgphi(:,:,1)   = gphiu(:,:) 
     488               zgphi(:,:,2)   = gphiv(:,:) 
     489               zmask(:,:,:,1) = umask(:,:,:) 
     490               zmask(:,:,:,2) = vmask(:,:,:) 
     491            ELSE 
     492               nvarsprof(jtype) = 1 
     493               nextrprof(jtype) = 0 
     494               ALLOCATE( llvar (nvarsprof(jtype)) ) 
     495               ALLOCATE( clvars(nvarsprof(jtype)) ) 
     496               ALLOCATE( zglam(jpi, jpj,      nvarsprof(jtype)) ) 
     497               ALLOCATE( zgphi(jpi, jpj,      nvarsprof(jtype)) ) 
     498               ALLOCATE( zmask(jpi, jpj, jpk, nvarsprof(jtype)) ) 
     499               llvar(1)       = .TRUE. 
     500               zglam(:,:,1)   = glamt(:,:) 
     501               zgphi(:,:,1)   = gphit(:,:) 
     502               zmask(:,:,:,1) = tmask(:,:,:) 
    412503            ENDIF 
    413504            ! 
     
    416507               &               clproffiles(jtype,1:ifilesprof(jtype)), & 
    417508               &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & 
    418                &               rn_dobsini, rn_dobsend, llvar1, llvar2, & 
    419                &               ln_ignmis, ln_s_at_t, .FALSE., & 
     509               &               rn_dobsini, rn_dobsend, llvar, & 
     510               &               ln_ignmis, ln_s_at_t, .FALSE., clvars, & 
    420511               &               kdailyavtypes = nn_profdavtypes ) 
    421512               ! 
     
    425516            ! 
    426517            CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), & 
    427                &               llvar1, llvar2, & 
     518               &               llvar, & 
    428519               &               jpi, jpj, jpk, & 
    429                &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2, & 
     520               &               zmask, zglam, zgphi, & 
    430521               &               ln_nea, ln_bound_reject, Kmm, & 
    431522               &               kdailyavtypes = nn_profdavtypes ) 
     523            ! 
     524            DEALLOCATE( llvar, clvars, zglam, zgphi, zmask ) 
     525            ! 
    432526         END DO 
    433527         ! 
     
    449543            IF( TRIM(cobstypessurf(jtype)) == 'sst' )   llnightav(jtype) = ln_sstnight 
    450544            ! 
     545            ALLOCATE( clvars( nvarssurf(jtype) ) ) 
     546            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
     547               clvars(1) = 'SLA' 
     548            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN 
     549               clvars(1) = 'SST' 
     550            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN 
     551               clvars(1) = 'ICECONC' 
     552            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 
     553               clvars(1) = 'SSS' 
     554            ENDIF 
     555            ! 
    451556            ! Read in surface obs types 
    452557            CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & 
    453558               &               clsurffiles(jtype,1:ifilessurf(jtype)), & 
    454559               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 
    455                &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) 
     560               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype), & 
     561               &               clvars ) 
    456562               ! 
    457563            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 
     
    473579                  &                  jnumsstbias      , cn_sstbiasfiles(1:jnumsstbias) )  
    474580            ENDIF 
     581            ! 
     582            DEALLOCATE( clvars ) 
    475583         END DO 
    476584         ! 
     
    516624      INTEGER :: jvar              ! Variable number 
    517625      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 
     626      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
     627         & zprofvar                ! Model values for variables in a prof ob 
     628      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
     629         & zprofmask               ! Mask associated with zprofvar 
    524630      REAL(wp), DIMENSION(jpi,jpj) :: & 
    525631         & zsurfvar, &             ! Model values equivalent to surface ob. 
    526632         & 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 
     633      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
     634         & zglam,    &             ! Model longitudes for prof variables 
     635         & zgphi                   ! Model latitudes for prof variables 
    532636 
    533637      !----------------------------------------------------------------------- 
     
    549653         DO jtype = 1, nproftypes 
    550654 
     655            ! Allocate local work arrays 
     656            ALLOCATE( zprofvar (jpi, jpj, jpk, profdataqc(jtype)%nvar) ) 
     657            ALLOCATE( zprofmask(jpi, jpj, jpk, profdataqc(jtype)%nvar) ) 
     658            ALLOCATE( zglam    (jpi, jpj,      profdataqc(jtype)%nvar) ) 
     659            ALLOCATE( zgphi    (jpi, jpj,      profdataqc(jtype)%nvar) )   
     660                               
     661            ! Defaults which might change 
     662            DO jvar = 1, profdataqc(jtype)%nvar 
     663               zprofmask(:,:,:,jvar) = tmask(:,:,:) 
     664               zglam(:,:,jvar)       = glamt(:,:) 
     665               zgphi(:,:,jvar)       = gphit(:,:) 
     666            END DO 
     667 
    551668            SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 
    552669            CASE('prof') 
    553                zprofvar1(:,:,:) = ts(:,:,:,jp_tem,Kmm) 
    554                zprofvar2(:,:,:) = ts(:,:,:,jp_sal,Kmm) 
    555                zprofmask1(:,:,:) = tmask(:,:,:) 
    556                zprofmask2(:,:,:) = tmask(:,:,:) 
    557                zglam1(:,:) = glamt(:,:) 
    558                zglam2(:,:) = glamt(:,:) 
    559                zgphi1(:,:) = gphit(:,:) 
    560                zgphi2(:,:) = gphit(:,:) 
     670               zprofvar(:,:,:,1) = ts(:,:,:,jp_tem,Kmm) 
     671               zprofvar(:,:,:,2) = ts(:,:,:,jp_sal,Kmm) 
    561672            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(:,:) 
     673               zprofvar(:,:,:,1) = uu(:,:,:,Kmm) 
     674               zprofvar(:,:,:,2) = vv(:,:,:,Kmm) 
     675               zprofmask(:,:,:,1) = umask(:,:,:) 
     676               zprofmask(:,:,:,2) = vmask(:,:,:) 
     677               zglam(:,:,1) = glamu(:,:) 
     678               zglam(:,:,2) = glamv(:,:) 
     679               zgphi(:,:,1) = gphiu(:,:) 
     680               zgphi(:,:,2) = gphiv(:,:) 
    570681            CASE DEFAULT 
    571682               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 
    572683            END SELECT 
    573684 
    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 ) 
     685            DO jvar = 1, profdataqc(jtype)%nvar 
     686               CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  & 
     687                  &               nit000, idaystp, jvar,                   & 
     688                  &               zprofvar(:,:,:,jvar),                    & 
     689                  &               gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm),      &  
     690                  &               zprofmask(:,:,:,jvar),                   & 
     691                  &               zglam(:,:,jvar), zgphi(:,:,jvar),        & 
     692                  &               nn_1dint, nn_2dint_default,              & 
     693                  &               kdailyavtypes = nn_profdavtypes ) 
     694            END DO 
     695             
     696            DEALLOCATE( zprofvar, zprofmask, zglam, zgphi ) 
    582697 
    583698         END DO 
     
    680795                  & ) 
    681796 
    682                CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv ) 
     797               CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv ) 
    683798 
    684799               DO jo = 1, profdataqc(jtype)%nprof 
     
    8961011   END SUBROUTINE fin_date 
    8971012    
    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 
     1013   SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, ifiles, cobstypes, cfiles ) 
     1014 
     1015      INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types 
     1016      INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type 
     1017      INTEGER, DIMENSION(ntypes), INTENT(OUT) :: & 
     1018         &                   ifiles      ! Out number of files for each type 
     1019      CHARACTER(len=8), DIMENSION(ntypes), INTENT(IN) :: & 
     1020         &                   cobstypes   ! List of obs types 
     1021      CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(IN) :: & 
     1022         &                   cfiles      ! List of files for all types 
     1023 
     1024      !Local variables 
     1025      INTEGER :: jfile 
     1026      INTEGER :: jtype 
     1027 
     1028      DO jtype = 1, ntypes 
     1029 
     1030         ifiles(jtype) = 0 
     1031         DO jfile = 1, jpmaxnfiles 
     1032            IF ( trim(cfiles(jtype,jfile)) /= '' ) & 
     1033                      ifiles(jtype) = ifiles(jtype) + 1 
     1034         END DO 
     1035 
     1036         IF ( ifiles(jtype) == 0 ) THEN 
     1037              CALL ctl_stop( 'Logical for observation type '//TRIM(cobstypes(jtype))//  & 
     1038                 &           ' set to true but no files available to read' ) 
     1039         ENDIF 
     1040 
     1041         IF(lwp) THEN     
     1042            WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:' 
     1043            DO jfile = 1, ifiles(jtype) 
     1044               WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile)) 
     1045            END DO 
     1046         ENDIF 
     1047 
     1048      END DO 
     1049 
     1050   END SUBROUTINE obs_settypefiles 
     1051 
     1052   SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             & 
     1053              &                  n2dint_default, n2dint_type,        & 
     1054              &                  ravglamscl_type, ravgphiscl_type,   & 
     1055              &                  lfp_indegs_type, lavnight_type,     & 
     1056              &                  n2dint, ravglamscl, ravgphiscl,     & 
     1057              &                  lfpindegs, lavnight ) 
     1058 
     1059      INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types 
     1060      INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs 
     1061      INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type 
     1062      INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type 
     1063      REAL(wp), INTENT(IN) :: & 
     1064         &                    ravglamscl_type, & !E/W diameter of obs footprint for this type 
     1065         &                    ravgphiscl_type    !N/S diameter of obs footprint for this type 
     1066      LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres 
     1067      LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average 
     1068      CHARACTER(len=8), INTENT(IN) :: ctypein  
     1069 
     1070      INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1071         &                    n2dint  
     1072      REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & 
     1073         &                    ravglamscl, ravgphiscl 
     1074      LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1075         &                    lfpindegs, lavnight 
     1076 
     1077      lavnight(jtype) = lavnight_type 
     1078 
     1079      IF ( (n2dint_type >= 0) .AND. (n2dint_type <= 6) ) THEN 
     1080         n2dint(jtype) = n2dint_type 
     1081      ELSE IF ( n2dint_type == -1 ) THEN 
     1082         n2dint(jtype) = n2dint_default 
     1083      ELSE 
     1084         CALL ctl_stop(' Choice of '//TRIM(ctypein)//' horizontal (2D) interpolation method', & 
     1085           &                    ' is not available') 
     1086      ENDIF 
     1087 
     1088      ! For averaging observation footprints set options for size of footprint  
     1089      IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN 
     1090         IF ( ravglamscl_type > 0._wp ) THEN 
     1091            ravglamscl(jtype) = ravglamscl_type 
     1092         ELSE 
     1093            CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
     1094                           'scale (ravglamscl) for observation type '//TRIM(ctypein) )       
     1095         ENDIF 
     1096 
     1097         IF ( ravgphiscl_type > 0._wp ) THEN 
     1098            ravgphiscl(jtype) = ravgphiscl_type 
     1099         ELSE 
     1100            CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
     1101                           'scale (ravgphiscl) for observation type '//TRIM(ctypein) )       
     1102         ENDIF 
     1103 
     1104         lfpindegs(jtype) = lfp_indegs_type  
     1105 
     1106      ENDIF 
     1107 
     1108      ! Write out info  
     1109      IF(lwp) THEN 
     1110         IF ( n2dint(jtype) <= 4 ) THEN 
     1111            WRITE(numout,*) '             '//TRIM(ctypein)// & 
     1112               &            ' model counterparts will be interpolated horizontally' 
     1113         ELSE IF ( n2dint(jtype) <= 6 ) THEN 
     1114            WRITE(numout,*) '             '//TRIM(ctypein)// & 
     1115               &            ' model counterparts will be averaged horizontally' 
     1116            WRITE(numout,*) '             '//'    with E/W scale: ',ravglamscl(jtype) 
     1117            WRITE(numout,*) '             '//'    with N/S scale: ',ravgphiscl(jtype) 
     1118            IF ( lfpindegs(jtype) ) THEN 
     1119                WRITE(numout,*) '             '//'    (in degrees)' 
     1120            ELSE 
     1121                WRITE(numout,*) '             '//'    (in metres)' 
     1122            ENDIF 
     1123         ENDIF 
     1124      ENDIF 
     1125 
     1126   END SUBROUTINE obs_setinterpopts 
    10121127 
    10131128END MODULE diaobs 
  • NEMO/trunk/src/OCE/OBS/obs_oper.F90

    r13295 r14056  
    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,kvar)-1 
     221         igrdj(1,1,iobs) = prodatqc%mj(jobs,kvar)-1 
     222         igrdi(1,2,iobs) = prodatqc%mi(jobs,kvar)-1 
     223         igrdj(1,2,iobs) = prodatqc%mj(jobs,kvar) 
     224         igrdi(2,1,iobs) = prodatqc%mi(jobs,kvar) 
     225         igrdj(2,1,iobs) = prodatqc%mj(jobs,kvar)-1 
     226         igrdi(2,2,iobs) = prodatqc%mi(jobs,kvar) 
     227         igrdj(2,2,iobs) = prodatqc%mj(jobs,kvar) 
    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 ) 
     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  
     
    450403    
    451404               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
    452                   &              prodatqc%var(1)%vmod(iend:iend) )  
     405                  &              prodatqc%var(kvar)%vmod(iend:iend) )  
    453406 
    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/trunk/src/OCE/OBS/obs_prep.F90

    r12489 r14056  
    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,              Kmm       ) 
     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/trunk/src/OCE/OBS/obs_read_prof.F90

    r13226 r14056  
    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                           ! Count number of profile var1 data as function of type 
     705                           itypvar( profdata%ntyp(iprof) + 1, jvar ) = & 
     706                              & itypvar( profdata%ntyp(iprof) + 1, jvar ) + 1 
     707                        ELSE 
     708                           profdata%var(jvar)%vobs(ivart(jvar)) = fbrmdi 
     709                        ENDIF 
     710 
     711                        ! Profile jvar qc 
     712                        profdata%var(jvar)%nvqc(ivart(jvar)) = & 
     713                           & inpfiles(jj)%ivlqc(ij,ji,jvar) 
     714 
     715                        ! Profile jvar qc flags 
     716                        profdata%var(jvar)%nvqcf(:,ivart(jvar)) = & 
     717                           & inpfiles(jj)%ivlqcf(:,ij,ji,jvar) 
     718 
     719                        ! Profile insitu T value 
     720                        IF ( TRIM( inpfiles(jj)%cname(jvar) ) == 'POTM' ) THEN 
     721                           profdata%var(jvar)%vext(ivart(jvar),1) = & 
     722                              &                inpfiles(jj)%pext(ij,ji,1) 
     723                        ENDIF 
     724 
    627725                     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 
     726                   
     727                  END DO 
    752728 
    753729               END DO loop_p 
     
    763739      !----------------------------------------------------------------------- 
    764740 
    765       CALL obs_mpp_sum_integer ( ivar1t0, ivar1tmpp ) 
    766       CALL obs_mpp_sum_integer ( ivar2t0, ivar2tmpp ) 
     741      DO jvar = 1, kvars 
     742         CALL obs_mpp_sum_integer ( ivart0(jvar), ivartmpp(jvar) ) 
     743      END DO 
    767744      CALL obs_mpp_sum_integer ( ip3dt,   ip3dtmpp  ) 
    768745 
    769       CALL obs_mpp_sum_integers( itypvar1, itypvar1mpp, ntyp1770 + 1 ) 
    770       CALL obs_mpp_sum_integers( itypvar2, itypvar2mpp, ntyp1770 + 1 ) 
     746      DO jvar = 1, kvars 
     747         CALL obs_mpp_sum_integers( itypvar(:,jvar), itypvarmpp(:,jvar), ntyp1770 + 1 ) 
     748      END DO 
    771749 
    772750      !----------------------------------------------------------------------- 
     
    778756         WRITE(numout,'(1X,A)') '------------' 
    779757         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 
     758         DO jvar = 1, kvars 
     759            WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(jvar) ) 
     760            WRITE(numout,'(1X,A)') '------------------------' 
     761            DO ji = 0, ntyp1770 
     762               IF ( itypvarmpp(ji+1,jvar) > 0 ) THEN 
     763                  WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), & 
     764                     & cwmonam1770(ji)(1:52),' = ', & 
     765                     & itypvarmpp(ji+1,jvar) 
     766               ENDIF 
     767            END DO 
     768            WRITE(numout,'(1X,A)') & 
     769               & '---------------------------------------------------------------' 
     770            WRITE(numout,'(1X,A55,I8)') & 
     771               & 'Total profile data for variable '//TRIM( profdata%cvars(jvar) )// & 
     772               & '             = ', ivartmpp(jvar) 
     773            WRITE(numout,'(1X,A)') & 
     774               & '---------------------------------------------------------------' 
     775            WRITE(numout,*)  
    788776         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 
     777      ENDIF 
     778 
     779      IF (ldsatt) THEN 
     780         profdata%nvprot(:)    = ip3dt 
     781         profdata%nvprotmpp(:) = ip3dtmpp 
     782      ELSE 
     783         DO jvar = 1, kvars 
     784            profdata%nvprot(jvar)    = ivart(jvar) 
     785            profdata%nvprotmpp(jvar) = ivartmpp(jvar) 
    805786         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 
    826787      ENDIF 
    827788      profdata%nprof        = iprof 
     
    830791      ! Model level search 
    831792      !----------------------------------------------------------------------- 
    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 
     793      DO jvar = 1, kvars 
     794         IF ( ldvar(jvar) ) THEN 
     795            CALL obs_level_search( jpk, gdept_1d, & 
     796               & profdata%nvprot(jvar), profdata%var(jvar)%vdep, & 
     797               & profdata%var(jvar)%mvk ) 
     798         ENDIF 
     799      END DO 
    842800 
    843801      !----------------------------------------------------------------------- 
     
    852810      ! Deallocate temporary data 
    853811      !----------------------------------------------------------------------- 
    854       DEALLOCATE( ifileidx, iprofidx, zdat, clvars ) 
     812      DEALLOCATE( ifileidx, iprofidx, zdat, clvarsin ) 
    855813 
    856814      !----------------------------------------------------------------------- 
  • NEMO/trunk/src/OCE/OBS/obs_read_surf.F90

    r13226 r14056  
    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/trunk/src/OCE/OBS/obs_write.F90

    r12933 r14056  
    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, 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         fbdata%cgrid(:)      = clgrid 
     196         DO je = 1, iext 
     197            fbdata%cextname(je) = pext%cdname(je) 
     198            fbdata%cextlong(je) = pext%cdlong(je,1) 
     199            fbdata%cextunit(je) = pext%cdunit(je,1) 
     200         END DO 
     201         DO ja = 1, iadd 
     202            fbdata%caddname(1+ja) = padd%cdname(ja) 
     203            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
     204            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     205         END DO 
     206      ENDIF 
    182207 
    183208      fbdata%caddname(1)   = 'Hx' 
     
    234259            &           krefdate = 19500101 ) 
    235260         ! Reform the profiles arrays for output 
    236          DO jvar = 1, 2 
     261         DO jvar = 1, profdata%nvar 
    237262            DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar) 
    238263               ik = profdata%var(jvar)%nvlidx(jk) 
     
    329354      CHARACTER(LEN=40) :: clfname         ! netCDF filename 
    330355      CHARACTER(LEN=10) :: clfiletype 
     356      CHARACTER(LEN=ilenlong) :: cllongname  ! Long name of variable 
     357      CHARACTER(LEN=ilenunit) :: clunits     ! Units of variable 
     358      CHARACTER(LEN=ilengrid) :: clgrid      ! Grid of variable 
    331359      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf' 
    332360      CHARACTER(LEN=12) :: clfmt           ! writing format 
     
    354382      SELECT CASE ( TRIM(surfdata%cvars(1)) ) 
    355383      CASE('SLA') 
     384          
     385         ! SLA needs special treatment because of MDT, so is all done here 
     386         ! Other variables are done more generically 
     387         ! No climatology for SLA, MDT is our best estimate of that and is already output. 
    356388 
    357389         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
     
    384416      CASE('SST') 
    385417 
     418         clfiletype = 'sstfb' 
     419         cllongname = 'Sea surface temperature' 
     420         clunits    = 'Degree centigrade' 
     421         clgrid     = 'T' 
     422          
     423      CASE('ICECONC') 
     424 
     425         clfiletype = 'sicfb' 
     426         cllongname = 'Sea ice concentration' 
     427         clunits    = 'Fraction' 
     428         clgrid     = 'T' 
     429 
     430      CASE('SSS') 
     431 
     432         clfiletype = 'sssfb' 
     433         cllongname = 'Sea surface salinity' 
     434         clunits    = 'psu' 
     435         clgrid     = 'T' 
     436 
     437      CASE DEFAULT 
     438 
     439         CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' ) 
     440 
     441      END SELECT 
     442 
     443      ! SLA needs special treatment because of MDT, so is done above 
     444      ! Remaining variables treated more generically 
     445 
     446      IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN 
     447       
    386448         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
    387449            &                 1 + iadd, iext, .TRUE. ) 
    388450 
    389          clfiletype = 'sstfb' 
    390451         fbdata%cname(1)      = surfdata%cvars(1) 
    391          fbdata%coblong(1)    = 'Sea surface temperature' 
    392          fbdata%cobunit(1)    = 'Degree centigrade' 
     452         fbdata%coblong(1)    = cllongname 
     453         fbdata%cobunit(1)    = clunits 
    393454         DO je = 1, iext 
    394455            fbdata%cextname(je) = pext%cdname(je) 
    395456            fbdata%cextlong(je) = pext%cdlong(je,1) 
    396457            fbdata%cextunit(je) = pext%cdunit(je,1) 
    397          END DO 
    398          fbdata%caddlong(1,1) = 'Model interpolated SST' 
    399          fbdata%caddunit(1,1) = 'Degree centigrade' 
    400          fbdata%cgrid(1)      = 'T' 
     458         END DO         
     459         IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN 
     460            fbdata%caddlong(1,1) = 'Model interpolated ICE' 
     461         ELSE 
     462            fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(surfdata%cvars(1)) 
     463         ENDIF 
     464         fbdata%caddunit(1,1) = clunits 
     465         fbdata%cgrid(1)      = clgrid 
    401466         DO ja = 1, iadd 
    402467            fbdata%caddname(1+ja) = padd%cdname(ja) 
     
    404469            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
    405470         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 
     471      ENDIF 
    458472 
    459473      fbdata%caddname(1)   = 'Hx' 
Note: See TracChangeset for help on using the changeset viewer.