Changeset 9186


Ignore:
Timestamp:
2018-01-05T14:29:29+01:00 (3 years ago)
Author:
dford
Message:

Initial implementation of 3D biogeochemistry observation operator.

Location:
branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/CONFIG/SHARED/namelist_ref

    r9016 r9186  
    12121212   ln_sfco2    = .false.            ! Logical switch for surface fugacity         of carbon dioxide obs 
    12131213   ln_spco2    = .false.            ! Logical switch for surface partial pressure of carbon dioxide obs 
     1214   ln_plchltot = .false.            ! Logical switch for profile total log10(chlorophyll) obs 
     1215   ln_pchltot  = .false.            ! Logical switch for profile total chlorophyll obs 
     1216   ln_pno3     = .false.            ! Logical switch for profile nitrate obs 
     1217   ln_psi4     = .false.            ! Logical switch for profile silicate obs 
     1218   ln_ppo4     = .false.            ! Logical switch for profile phosphate obs 
     1219   ln_pdic     = .false.            ! Logical switch for profile dissolved inorganic carbon obs 
     1220   ln_palk     = .false.            ! Logical switch for profile alkalinity obs 
     1221   ln_pph      = .false.            ! Logical switch for profile pH obs 
     1222   ln_po2      = .false.            ! Logical switch for profile dissolved oxygen obs 
    12141223   ln_altbias = .false.             ! Logical switch for altimeter bias correction 
    12151224   ln_sstbias = .false.             ! Logical switch for SST bias correction 
     
    12261235   ln_sic_fp_indegs = .true. 
    12271236! All of the *files* variables below are arrays. Use namelist_cfg to add more files 
    1228    cn_profbfiles = 'profiles_01.nc'    ! Profile feedback input observation file names 
    1229    cn_slafbfiles = 'sla_01.nc'         ! SLA feedback input observation file names 
    1230    cn_sstfbfiles = 'sst_01.nc'         ! SST feedback input observation file names 
    1231    cn_sicfbfiles = 'sic_01.nc'         ! SIC feedback input observation file names 
    1232    cn_velfbfiles = 'vel_01.nc'         ! Velocity feedback input observation file names 
    1233    cn_sssfbfiles = 'sss_01.nc'         ! SSS feedback input observation file names 
    1234    cn_altbiasfile = 'altbias.nc'       ! Altimeter bias input file name 
    1235    cn_sstbiasfiles = 'sstbias.nc'      ! SST bias input file names 
    1236    cn_gridsearchfile = 'gridsearch.nc' ! Grid search file name 
    1237    rn_gridsearchres = 0.5              ! Grid search resolution 
    1238    rn_default_avglamscl = 0.           ! Default E/W diameter of observation footprint (metres/degrees) 
    1239    rn_default_avgphiscl = 0.           ! Default N/S diameter of observation footprint (metres/degrees) 
    1240    rn_sla_avglamscl = 0.               ! E/W diameter of SLA observation footprint (metres/degrees) 
    1241    rn_sla_avgphiscl = 0.               ! N/S diameter of SLA observation footprint (metres/degrees) 
    1242    rn_sst_avglamscl = 0.               ! E/W diameter of SST observation footprint (metres/degrees) 
    1243    rn_sst_avgphiscl = 0.               ! N/S diameter of SST observation footprint (metres/degrees) 
    1244    rn_sss_avglamscl = 0.               ! E/W diameter of SSS observation footprint (metres/degrees) 
    1245    rn_sss_avgphiscl = 0.               ! N/S diameter of SSS observation footprint (metres/degrees) 
    1246    rn_sic_avglamscl = 0.               ! E/W diameter of SIC observation footprint (metres/degrees) 
    1247    rn_sic_avgphiscl = 0.               ! N/S diameter of SIC observation footprint (metres/degrees) 
    1248    nn_1dint = 0                        ! Type of vertical interpolation method 
    1249    nn_2dint_default = 0                ! Default horizontal interpolation method 
    1250    nn_2dint_sla = -1                   ! Horizontal interpolation method for SLA 
    1251    nn_2dint_sst = -1                   ! Horizontal interpolation method for SST 
    1252    nn_2dint_sss = -1                   ! Horizontal interpolation method for SSS 
    1253    nn_2dint_sic = -1                   ! Horizontal interpolation method for SIC 
    1254    nn_msshc = 0                        ! MSSH correction scheme 
    1255    rn_mdtcorr = 1.61                   ! MDT  correction 
    1256    rn_mdtcutoff = 65.0                 ! MDT cutoff for computed correction 
    1257    nn_profdavtypes = -1                ! Profile daily average types - array 
     1237   cn_profbfiles = 'profiles_01.nc'      ! Profile feedback input observation file names 
     1238   cn_slafbfiles = 'sla_01.nc'           ! SLA feedback input observation file names 
     1239   cn_sstfbfiles = 'sst_01.nc'           ! SST feedback input observation file names 
     1240   cn_sicfbfiles = 'sic_01.nc'           ! SIC feedback input observation file names 
     1241   cn_velfbfiles = 'vel_01.nc'           ! Velocity feedback input observation file names 
     1242   cn_sssfbfiles = 'sss_01.nc'           ! SSS feedback input observation file names 
     1243   cn_slchltotfbfiles = 'slchltot_01.nc' ! Surface total              log10(chlorophyll) obs file names 
     1244   cn_slchldiafbfiles = 'slchldia_01.nc' ! Surface diatom             log10(chlorophyll) obs file names 
     1245   cn_slchlnonfbfiles = 'slchlnon_01.nc' ! Surface non-diatom         log10(chlorophyll) obs file names 
     1246   cn_slchldinfbfiles = 'slchldin_01.nc' ! Surface dinoflagellate     log10(chlorophyll) obs file names 
     1247   cn_slchlmicfbfiles = 'slchlmic_01.nc' ! Surface microphytoplankton log10(chlorophyll) obs file names 
     1248   cn_slchlnanfbfiles = 'slchlnan_01.nc' ! Surface nanophytoplankton  log10(chlorophyll) obs file names 
     1249   cn_slchlpicfbfiles = 'slchlpic_01.nc' ! Surface picophytoplankton  log10(chlorophyll) obs file names 
     1250   cn_schltotfbfiles  = 'schltot_01.nc'  ! Surface total              chlorophyll        obs file names 
     1251   cn_sspmfbfiles     = 'sspm_01.nc'     ! Surface suspended particulate matter obs file names 
     1252   cn_sfco2fbfiles    = 'sfco2_01.nc'    ! Surface fugacity         of carbon dioxide obs file names 
     1253   cn_spco2fbfiles    = 'spco2_01.nc'    ! Surface partial pressure of carbon dioxide obs file names 
     1254   cn_plchltotfbfiles = 'plchltot_01.nc' ! Profile total log10(chlorophyll) obs file names 
     1255   cn_pchltotfbfiles  = 'pchltot_01.nc'  ! Profile total chlorophyll obs file names 
     1256   cn_pno3fbfiles     = 'pno3_01.nc'     ! Profile nitrate obs file names 
     1257   cn_psi4fbfiles     = 'psi4_01.nc'     ! Profile silicate obs file names 
     1258   cn_ppo4fbfiles     = 'ppo4_01.nc'     ! Profile phosphate obs file names 
     1259   cn_pdicfbfiles     = 'pdic_01.nc'     ! Profile dissolved inorganic carbon obs file names 
     1260   cn_palkfbfiles     = 'palk_01.nc'     ! Profile alkalinity obs file names 
     1261   cn_pphfbfiles      = 'pph_01.nc'      ! Profile pH obs file names 
     1262   cn_po2fbfiles      = 'po2_01.nc'      ! Profile dissolved oxygen obs file names 
     1263   cn_altbiasfile = 'altbias.nc'         ! Altimeter bias input file name 
     1264   cn_sstbiasfiles = 'sstbias.nc'        ! SST bias input file names 
     1265   cn_gridsearchfile = 'gridsearch.nc'   ! Grid search file name 
     1266   rn_gridsearchres = 0.5                ! Grid search resolution 
     1267   rn_default_avglamscl = 0.             ! Default E/W diameter of observation footprint (metres/degrees) 
     1268   rn_default_avgphiscl = 0.             ! Default N/S diameter of observation footprint (metres/degrees) 
     1269   rn_sla_avglamscl = 0.                 ! E/W diameter of SLA observation footprint (metres/degrees) 
     1270   rn_sla_avgphiscl = 0.                 ! N/S diameter of SLA observation footprint (metres/degrees) 
     1271   rn_sst_avglamscl = 0.                 ! E/W diameter of SST observation footprint (metres/degrees) 
     1272   rn_sst_avgphiscl = 0.                 ! N/S diameter of SST observation footprint (metres/degrees) 
     1273   rn_sss_avglamscl = 0.                 ! E/W diameter of SSS observation footprint (metres/degrees) 
     1274   rn_sss_avgphiscl = 0.                 ! N/S diameter of SSS observation footprint (metres/degrees) 
     1275   rn_sic_avglamscl = 0.                 ! E/W diameter of SIC observation footprint (metres/degrees) 
     1276   rn_sic_avgphiscl = 0.                 ! N/S diameter of SIC observation footprint (metres/degrees) 
     1277   nn_1dint = 0                          ! Type of vertical interpolation method 
     1278   nn_2dint_default = 0                  ! Default horizontal interpolation method 
     1279   nn_2dint_sla = -1                     ! Horizontal interpolation method for SLA 
     1280   nn_2dint_sst = -1                     ! Horizontal interpolation method for SST 
     1281   nn_2dint_sss = -1                     ! Horizontal interpolation method for SSS 
     1282   nn_2dint_sic = -1                     ! Horizontal interpolation method for SIC 
     1283   nn_msshc = 0                          ! MSSH correction scheme 
     1284   rn_mdtcorr = 1.61                     ! MDT  correction 
     1285   rn_mdtcutoff = 65.0                   ! MDT cutoff for computed correction 
     1286   nn_profdavtypes = -1                  ! Profile daily average types - array 
    12581287    
    12591288/ 
  • branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

    r9016 r9186  
    164164         & cn_sfco2fbfiles,    & ! Surface fugacity         of carbon dioxide input filenames 
    165165         & cn_spco2fbfiles,    & ! Surface partial pressure of carbon dioxide input filenames 
     166         & cn_plchltotfbfiles, & ! Profile total log10(chlorophyll) input filenames 
     167         & cn_pchltotfbfiles,  & ! Profile total chlorophyll input filenames 
     168         & cn_pno3fbfiles,     & ! Profile nitrate input filenames 
     169         & cn_psi4fbfiles,     & ! Profile silicate input filenames 
     170         & cn_ppo4fbfiles,     & ! Profile phosphate input filenames 
     171         & cn_pdicfbfiles,     & ! Profile dissolved inorganic carbon input filenames 
     172         & cn_palkfbfiles,     & ! Profile alkalinity input filenames 
     173         & cn_pphfbfiles,      & ! Profile pH input filenames 
     174         & cn_po2fbfiles,      & ! Profile dissolved oxygen input filenames 
    166175         & cn_sstbiasfiles       ! SST bias input filenames 
    167176 
     
    188197      LOGICAL :: ln_sfco2        ! Logical switch for surface fugacity         of carbon dioxide obs 
    189198      LOGICAL :: ln_spco2        ! Logical switch for surface partial pressure of carbon dioxide obs 
     199      LOGICAL :: ln_plchltot     ! Logical switch for profile total log10(chlorophyll) obs 
     200      LOGICAL :: ln_pchltot      ! Logical switch for profile total chlorophyll obs 
     201      LOGICAL :: ln_pno3         ! Logical switch for profile nitrate obs 
     202      LOGICAL :: ln_psi4         ! Logical switch for profile silicate obs 
     203      LOGICAL :: ln_ppo4         ! Logical switch for profile phosphate obs 
     204      LOGICAL :: ln_pdic         ! Logical switch for profile dissolved inorganic carbon obs 
     205      LOGICAL :: ln_palk         ! Logical switch for profile alkalinity obs 
     206      LOGICAL :: ln_pph          ! Logical switch for profile pH obs 
     207      LOGICAL :: ln_po2          ! Logical switch for profile dissolved oxygen obs 
    190208      LOGICAL :: ln_nea          ! Logical switch to remove obs near land 
    191209      LOGICAL :: ln_altbias      ! Logical switch for altimeter bias 
     
    227245         &            ln_slchlpic, ln_schltot,                        & 
    228246         &            ln_sspm,     ln_sfco2,    ln_spco2,             & 
     247         &            ln_plchltot, ln_pchltot,  ln_pno3,              & 
     248         &            ln_psi4,     ln_ppo4,     ln_pdic,              & 
     249         &            ln_palk,     ln_pph,      ln_po2,               & 
    229250         &            ln_altbias, ln_sstbias, ln_nea,                 & 
    230251         &            ln_grid_global, ln_grid_search_lookup,          & 
     
    242263         &            cn_schltotfbfiles, cn_sspmfbfiles,              & 
    243264         &            cn_sfco2fbfiles, cn_spco2fbfiles,               & 
     265         &            cn_plchltotfbfiles, cn_pchltotfbfiles,          & 
     266         &            cn_pno3fbfiles, cn_psi4fbfiles, cn_ppo4fbfiles, & 
     267         &            cn_pdicfbfiles, cn_palkfbfiles, cn_pphfbfiles,  & 
     268         &            cn_po2fbfiles,                                  & 
    244269         &            cn_sstbiasfiles, cn_altbiasfile,                & 
    245270         &            cn_gridsearchfile, rn_gridsearchres,            & 
     
    285310      cn_sfco2fbfiles(:)    = '' 
    286311      cn_spco2fbfiles(:)    = '' 
     312      cn_plchltotfbfiles(:) = '' 
     313      cn_pchltotfbfiles(:)  = '' 
     314      cn_pno3fbfiles(:)     = '' 
     315      cn_psi4fbfiles(:)     = '' 
     316      cn_ppo4fbfiles(:)     = '' 
     317      cn_pdicfbfiles(:)     = '' 
     318      cn_palkfbfiles(:)     = '' 
     319      cn_pphfbfiles(:)      = '' 
     320      cn_po2fbfiles(:)      = '' 
    287321      cn_sstbiasfiles(:)    = '' 
    288322      nn_profdavtypes(:)    = -1 
     
    335369         WRITE(numout,*) '             Logical switch for surface fCO2 observations           ln_sfco2 = ', ln_sfco2 
    336370         WRITE(numout,*) '             Logical switch for surface pCO2 observations           ln_spco2 = ', ln_spco2 
     371         WRITE(numout,*) '             Logical switch for profile total logchl obs         ln_plchltot = ', ln_plchltot 
     372         WRITE(numout,*) '             Logical switch for profile total chl obs             ln_pchltot = ', ln_pchltot 
     373         WRITE(numout,*) '             Logical switch for profile nitrate obs                  ln_pno3 = ', ln_pno3 
     374         WRITE(numout,*) '             Logical switch for profile silicate obs                 ln_psi4 = ', ln_psi4 
     375         WRITE(numout,*) '             Logical switch for profile phosphate obs                ln_ppo4 = ', ln_ppo4 
     376         WRITE(numout,*) '             Logical switch for profile DIC obs                      ln_pdic = ', ln_pdic 
     377         WRITE(numout,*) '             Logical switch for profile alkalinity obs               ln_palk = ', ln_palk 
     378         WRITE(numout,*) '             Logical switch for profile pH obs                        ln_pph = ', ln_pph 
     379         WRITE(numout,*) '             Logical switch for profile oxygen obs                    ln_po2 = ', ln_po2 
    337380         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global 
    338381         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup 
     
    375418      !----------------------------------------------------------------------- 
    376419 
    377       nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 
     420      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d, ln_plchltot,          & 
     421         &                  ln_pchltot,  ln_pno3,     ln_psi4,     ln_ppo4,     & 
     422         &                  ln_pdic,     ln_palk,     ln_pph,      ln_po2 /) ) 
    378423      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss,                     & 
    379424         &                  ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, & 
    380425         &                  ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot,  & 
    381          &                  ln_sspm, ln_sfco2, ln_spco2 /) ) 
     426         &                  ln_sspm,     ln_sfco2,    ln_spco2 /) ) 
    382427 
    383428      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 
     
    400445         IF (ln_t3d .OR. ln_s3d) THEN 
    401446            jtype = jtype + 1 
    402             CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof  ', & 
     447            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof', & 
    403448               &                   cn_profbfiles, ifilesprof, cobstypesprof, clproffiles ) 
    404449         ENDIF 
    405450         IF (ln_vel3d) THEN 
    406451            jtype = jtype + 1 
    407             CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel   ', & 
     452            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel', & 
    408453               &                   cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles ) 
     454         ENDIF 
     455         IF (ln_plchltot) THEN 
     456            jtype = jtype + 1 
     457            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'plchltot', & 
     458               &                   cn_plchltotfbfiles, ifilesprof, cobstypesprof, clproffiles ) 
     459         ENDIF 
     460         IF (ln_pchltot) THEN 
     461            jtype = jtype + 1 
     462            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'pchltot', & 
     463               &                   cn_pchltotfbfiles, ifilesprof, cobstypesprof, clproffiles ) 
     464         ENDIF 
     465         IF (ln_pno3) THEN 
     466            jtype = jtype + 1 
     467            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'pno3', & 
     468               &                   cn_pno3fbfiles, ifilesprof, cobstypesprof, clproffiles ) 
     469         ENDIF 
     470         IF (ln_psi4) THEN 
     471            jtype = jtype + 1 
     472            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'psi4', & 
     473               &                   cn_psi4fbfiles, ifilesprof, cobstypesprof, clproffiles ) 
     474         ENDIF 
     475         IF (ln_ppo4) THEN 
     476            jtype = jtype + 1 
     477            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'ppo4', & 
     478               &                   cn_ppo4fbfiles, ifilesprof, cobstypesprof, clproffiles ) 
     479         ENDIF 
     480         IF (ln_pdic) THEN 
     481            jtype = jtype + 1 
     482            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'pdic', & 
     483               &                   cn_pdicfbfiles, ifilesprof, cobstypesprof, clproffiles ) 
     484         ENDIF 
     485         IF (ln_palk) THEN 
     486            jtype = jtype + 1 
     487            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'palk', & 
     488               &                   cn_palkfbfiles, ifilesprof, cobstypesprof, clproffiles ) 
     489         ENDIF 
     490         IF (ln_pph) THEN 
     491            jtype = jtype + 1 
     492            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'pph', & 
     493               &                   cn_pphfbfiles, ifilesprof, cobstypesprof, clproffiles ) 
     494         ENDIF 
     495         IF (ln_po2) THEN 
     496            jtype = jtype + 1 
     497            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'po2', & 
     498               &                   cn_po2fbfiles, ifilesprof, cobstypesprof, clproffiles ) 
    409499         ENDIF 
    410500 
     
    615705         DO jtype = 1, nproftypes 
    616706 
    617             nvarsprof(jtype) = 2 
    618707            IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN 
     708               nvarsprof(jtype) = 2 
    619709               nextrprof(jtype) = 1 
    620710               llvar1 = ln_t3d 
     
    626716               zgphi2 = gphit 
    627717               zmask2 = tmask 
    628             ENDIF 
    629             IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN 
     718            ELSE IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN 
     719               nvarsprof(jtype) = 2 
    630720               nextrprof(jtype) = 2 
    631721               llvar1 = ln_vel3d 
     
    637727               zgphi2 = gphiv 
    638728               zmask2 = vmask 
     729            ELSE 
     730               nvarsprof(jtype) = 1 
     731               nextrprof(jtype) = 0 
     732               llvar1 = .TRUE. 
     733               llvar2 = .FALSE. 
     734               zglam1 = glamt 
     735               zgphi1 = gphit 
     736               zmask1 = tmask 
     737               zglam2 = glamt 
     738               zgphi2 = gphit 
     739               zmask2 = tmask 
    639740            ENDIF 
    640741 
     
    763864#endif 
    764865#if defined key_hadocc 
    765       USE trc, ONLY :  &           ! HadOCC chlorophyll, fCO2 and pCO2 
     866      USE trc, ONLY :  &           ! HadOCC variables 
     867         & trn, & 
    766868         & HADOCC_CHL, & 
    767869         & HADOCC_FCO2, & 
    768870         & HADOCC_PCO2, & 
    769871         & HADOCC_FILL_FLT 
     872      USE par_hadocc 
    770873#elif defined key_medusa && defined key_foam_medusa 
    771       USE trc, ONLY :  &           ! MEDUSA chlorophyll, fCO2 and pCO2 
     874      USE trc, ONLY :  &           ! MEDUSA variables 
    772875         & trn 
    773       USE par_medusa, ONLY: & 
    774          & jpchn, & 
    775          & jpchd 
     876      USE par_medusa 
    776877#if defined key_roam 
    777878      USE sms_medusa, ONLY: & 
    778879         & f2_pco2w, & 
    779          & f2_fco2w 
     880         & f2_fco2w, & 
     881         & f3_pH 
    780882#endif 
    781883#elif defined key_fabm 
     
    864966               zgphi1(:,:) = gphiu(:,:) 
    865967               zgphi2(:,:) = gphiv(:,:) 
     968 
     969            CASE('plchltot') 
     970#if defined key_hadocc 
     971               ! Chlorophyll from HadOCC 
     972               zprofvar1(:,:,:) = HADOCC_CHL(:,:,:) 
     973#elif defined key_medusa && defined key_foam_medusa 
     974               ! Add non-diatom and diatom chlorophyll from MEDUSA 
     975               zprofvar1(:,:,:) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 
     976#elif defined key_fabm 
     977               ! Add all chlorophyll groups from ERSEM 
     978               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_chl1) + trn(:,:,:,jp_fabm_chl2) + & 
     979                  &               trn(:,:,:,jp_fabm_chl3) + trn(:,:,:,jp_fabm_chl4) 
     980#else 
     981               CALL ctl_stop( ' Trying to run plchltot observation operator', & 
     982                  &           ' but no biogeochemical model appears to have been defined' ) 
     983#endif 
     984               zprofmask1(:,:,:) = tmask(:,:,:) 
     985               ! Take the log10 where we can, otherwise exclude 
     986               tiny = 1.0e-20 
     987               WHERE(zprofvar1(:,:,:) > tiny .AND. zprofvar1(:,:,:) /= obfillflt ) 
     988                  zprofvar1(:,:,:)  = LOG10(zprofvar1(:,:,:)) 
     989               ELSEWHERE 
     990                  zprofvar1(:,:,:)  = obfillflt 
     991                  zprofmask1(:,:,:) = 0 
     992               END WHERE 
     993               zglam1(:,:) = glamt(:,:) 
     994               zgphi1(:,:) = gphit(:,:) 
     995 
     996            CASE('pchltot') 
     997#if defined key_hadocc 
     998               ! Chlorophyll from HadOCC 
     999               zprofvar1(:,:,:) = HADOCC_CHL(:,:,:) 
     1000#elif defined key_medusa && defined key_foam_medusa 
     1001               ! Add non-diatom and diatom chlorophyll from MEDUSA 
     1002               zprofvar1(:,:,:) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 
     1003#elif defined key_fabm 
     1004               ! Add all chlorophyll groups from ERSEM 
     1005               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_chl1) + trn(:,:,:,jp_fabm_chl2) + & 
     1006                  &               trn(:,:,:,jp_fabm_chl3) + trn(:,:,:,jp_fabm_chl4) 
     1007#else 
     1008               CALL ctl_stop( ' Trying to run pchltot observation operator', & 
     1009                  &           ' but no biogeochemical model appears to have been defined' ) 
     1010#endif 
     1011               zprofmask1(:,:,:) = tmask(:,:,:) 
     1012               zglam1(:,:) = glamt(:,:) 
     1013               zgphi1(:,:) = gphit(:,:) 
     1014 
     1015            CASE('pno3') 
     1016#if defined key_hadocc 
     1017               ! Dissolved inorganic nitrogen from HadOCC 
     1018               zprofvar1(:,:,:) = trn(:,:,:,jp_had_nut) 
     1019#elif defined key_medusa && defined key_foam_medusa 
     1020               ! Dissolved inorganic nitrogen from MEDUSA 
     1021               zprofvar1(:,:,:) = trn(:,:,:,jpdin) 
     1022#elif defined key_fabm 
     1023               ! Nitrate from ERSEM 
     1024               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_n3n) 
     1025#else 
     1026               CALL ctl_stop( ' Trying to run pno3 observation operator', & 
     1027                  &           ' but no biogeochemical model appears to have been defined' ) 
     1028#endif 
     1029               zprofmask1(:,:,:) = tmask(:,:,:) 
     1030               zglam1(:,:) = glamt(:,:) 
     1031               zgphi1(:,:) = gphit(:,:) 
     1032 
     1033            CASE('psi4') 
     1034#if defined key_hadocc 
     1035               CALL ctl_stop( ' Trying to run psi4 observation operator', & 
     1036                  &           ' but HadOCC does not simulate silicate' ) 
     1037#elif defined key_medusa && defined key_foam_medusa 
     1038               ! Silicate from MEDUSA 
     1039               zprofvar1(:,:,:) = trn(:,:,:,jpsil) 
     1040#elif defined key_fabm 
     1041               ! Silicate from ERSEM 
     1042               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_n5s) 
     1043#else 
     1044               CALL ctl_stop( ' Trying to run psi4 observation operator', & 
     1045                  &           ' but no biogeochemical model appears to have been defined' ) 
     1046#endif 
     1047               zprofmask1(:,:,:) = tmask(:,:,:) 
     1048               zglam1(:,:) = glamt(:,:) 
     1049               zgphi1(:,:) = gphit(:,:) 
     1050 
     1051            CASE('ppo4') 
     1052#if defined key_hadocc 
     1053               CALL ctl_stop( ' Trying to run ppo4 observation operator', & 
     1054                  &           ' but HadOCC does not simulate phosphate' ) 
     1055#elif defined key_medusa && defined key_foam_medusa 
     1056               CALL ctl_stop( ' Trying to run ppo4 observation operator', & 
     1057                  &           ' but MEDUSA does not simulate phosphate' ) 
     1058#elif defined key_fabm 
     1059               ! Phosphate from ERSEM 
     1060               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_n1p) 
     1061#else 
     1062               CALL ctl_stop( ' Trying to run ppo4 observation operator', & 
     1063                  &           ' but no biogeochemical model appears to have been defined' ) 
     1064#endif 
     1065               zprofmask1(:,:,:) = tmask(:,:,:) 
     1066               zglam1(:,:) = glamt(:,:) 
     1067               zgphi1(:,:) = gphit(:,:) 
     1068 
     1069            CASE('pdic') 
     1070#if defined key_hadocc 
     1071               ! Dissolved inorganic carbon from HadOCC 
     1072               zprofvar1(:,:,:) = trn(:,:,:,jp_had_dic) 
     1073#elif defined key_medusa && defined key_foam_medusa 
     1074               ! Dissolved inorganic carbon from MEDUSA 
     1075               zprofvar1(:,:,:) = trn(:,:,:,jpdic) 
     1076#elif defined key_fabm 
     1077               ! Dissolved inorganic carbon from ERSEM 
     1078               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_o3c) 
     1079#else 
     1080               CALL ctl_stop( ' Trying to run pdic observation operator', & 
     1081                  &           ' but no biogeochemical model appears to have been defined' ) 
     1082#endif 
     1083               zprofmask1(:,:,:) = tmask(:,:,:) 
     1084               zglam1(:,:) = glamt(:,:) 
     1085               zgphi1(:,:) = gphit(:,:) 
     1086 
     1087            CASE('palk') 
     1088#if defined key_hadocc 
     1089               ! Alkalinity from HadOCC 
     1090               zprofvar1(:,:,:) = trn(:,:,:,jp_had_alk) 
     1091#elif defined key_medusa && defined key_foam_medusa 
     1092               ! Alkalinity from MEDUSA 
     1093               zprofvar1(:,:,:) = trn(:,:,:,jpalk) 
     1094#elif defined key_fabm 
     1095               ! Alkalinity from ERSEM 
     1096               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_o3a) 
     1097#else 
     1098               CALL ctl_stop( ' Trying to run palk observation operator', & 
     1099                  &           ' but no biogeochemical model appears to have been defined' ) 
     1100#endif 
     1101               zprofmask1(:,:,:) = tmask(:,:,:) 
     1102               zglam1(:,:) = glamt(:,:) 
     1103               zgphi1(:,:) = gphit(:,:) 
     1104 
     1105            CASE('pph') 
     1106#if defined key_hadocc 
     1107               CALL ctl_stop( ' Trying to run pph observation operator', & 
     1108                  &           ' but HadOCC has no pH diagnostic defined' ) 
     1109#elif defined key_medusa && defined key_foam_medusa 
     1110               ! pH from MEDUSA 
     1111               zprofvar1(:,:,:) = f3_pH(:,:,:) 
     1112#elif defined key_fabm 
     1113               ! pH from ERSEM 
     1114               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_o3ph) 
     1115#else 
     1116               CALL ctl_stop( ' Trying to run pph observation operator', & 
     1117                  &           ' but no biogeochemical model appears to have been defined' ) 
     1118#endif 
     1119               zprofmask1(:,:,:) = tmask(:,:,:) 
     1120               zglam1(:,:) = glamt(:,:) 
     1121               zgphi1(:,:) = gphit(:,:) 
     1122 
     1123            CASE('po2') 
     1124#if defined key_hadocc 
     1125               CALL ctl_stop( ' Trying to run po2 observation operator', & 
     1126                  &           ' but HadOCC does not simulate oxygen' ) 
     1127#elif defined key_medusa && defined key_foam_medusa 
     1128               ! Oxygen from MEDUSA 
     1129               zprofvar1(:,:,:) = trn(:,:,:,jpoxy) 
     1130#elif defined key_fabm 
     1131               ! Oxygen from ERSEM 
     1132               zprofvar1(:,:,:) = trn(:,:,:,jp_fabm_o2o) 
     1133#else 
     1134               CALL ctl_stop( ' Trying to run po2 observation operator', & 
     1135                  &           ' but no biogeochemical model appears to have been defined' ) 
     1136#endif 
     1137               zprofmask1(:,:,:) = tmask(:,:,:) 
     1138               zglam1(:,:) = glamt(:,:) 
     1139               zgphi1(:,:) = gphit(:,:) 
     1140 
    8661141            CASE DEFAULT 
    8671142               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 
    8681143            END SELECT 
     1144             
     1145            IF ( ( TRIM(cobstypesprof(jtype)) /= 'prof' ) .AND. ( TRIM(cobstypesprof(jtype)) /= 'vel' ) ) THEN 
     1146               zprofvar2(:,:,:)  = zprofvar1(:,:,:) 
     1147               zprofmask2(:,:,:) = zprofmask1(:,:,:) 
     1148               zglam2(:,:)       = zglam1(:,:) 
     1149               zgphi2(:,:)       = zgphi1(:,:) 
     1150            ENDIF 
    8691151 
    8701152            CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  & 
  • branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r8653 r9186  
    231231                  DO ji = 1, jpi 
    232232                     prodatqc%vdmean(ji,jj,jk,1) = 0.0 
    233                      prodatqc%vdmean(ji,jj,jk,2) = 0.0 
     233                     IF ( prodatqc%nvar == 2 ) THEN 
     234                        prodatqc%vdmean(ji,jj,jk,2) = 0.0 
     235                     ENDIF 
    234236                  END DO 
    235237               END DO 
     
    243245                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
    244246                     &                        + pvar1(ji,jj,jk) 
    245                   ! Increment field 2 for computing daily mean 
    246                   prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
    247                      &                        + pvar2(ji,jj,jk) 
     247                  IF ( prodatqc%nvar == 2 ) THEN 
     248                     ! Increment field 2 for computing daily mean 
     249                     prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
     250                        &                        + pvar2(ji,jj,jk) 
     251                  ENDIF 
    248252               END DO 
    249253            END DO 
     
    260264                     prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
    261265                        &                        * zdaystp 
    262                      prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
    263                         &                        * zdaystp 
     266                     IF ( prodatqc%nvar == 2 ) THEN 
     267                        prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
     268                           &                        * zdaystp 
     269                     ENDIF 
    264270                  END DO 
    265271               END DO 
     
    297303         igrdi1(2,2,iobs) = prodatqc%mi(jobs,1) 
    298304         igrdj1(2,2,iobs) = prodatqc%mj(jobs,1) 
    299          igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1 
    300          igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1 
    301          igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1 
    302          igrdj2(1,2,iobs) = prodatqc%mj(jobs,2) 
    303          igrdi2(2,1,iobs) = prodatqc%mi(jobs,2) 
    304          igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1 
    305          igrdi2(2,2,iobs) = prodatqc%mi(jobs,2) 
    306          igrdj2(2,2,iobs) = prodatqc%mj(jobs,2) 
     305         IF ( prodatqc%nvar == 2 ) THEN 
     306            igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1 
     307            igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1 
     308            igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1 
     309            igrdj2(1,2,iobs) = prodatqc%mj(jobs,2) 
     310            igrdi2(2,1,iobs) = prodatqc%mi(jobs,2) 
     311            igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1 
     312            igrdi2(2,2,iobs) = prodatqc%mi(jobs,2) 
     313            igrdj2(2,2,iobs) = prodatqc%mj(jobs,2) 
     314         ENDIF 
    307315      END DO 
    308316 
     
    316324      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1,   zint1 ) 
    317325       
    318       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 ) 
    319       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 ) 
    320       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 ) 
    321       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2,   zint2 ) 
     326      IF ( prodatqc%nvar == 2 ) THEN 
     327         CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 ) 
     328         CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 ) 
     329         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 ) 
     330         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2,   zint2 ) 
     331      ENDIF 
    322332 
    323333      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept )  
     
    334344         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, & 
    335345            &                  prodatqc%vdmean(:,:,:,1), zinm1 ) 
    336          CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, & 
    337             &                  prodatqc%vdmean(:,:,:,2), zinm2 ) 
     346         IF ( prodatqc%nvar == 2 ) THEN 
     347            CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, & 
     348               &                  prodatqc%vdmean(:,:,:,2), zinm2 ) 
     349         ENDIF 
    338350 
    339351      ENDIF 
     
    378390         ENDIF 
    379391 
    380          IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    381  
    382             CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
    383                &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), & 
    384                &                   zmask2(:,:,1,iobs), zweig2, zmsk_2 ) 
    385   
     392         IF ( prodatqc%nvar == 2 ) THEN 
     393            IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
     394 
     395               CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
     396                  &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), & 
     397                  &                   zmask2(:,:,1,iobs), zweig2, zmsk_2 ) 
     398 
     399            ENDIF 
    386400         ENDIF 
    387401 
     
    512526         ENDIF  
    513527 
    514          ! For the second variable 
    515          IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    516  
    517             zobsk(:) = obfillflt 
    518  
    519             IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 
    520  
    521                IF ( idayend == 0 )  THEN 
    522                   ! Daily averaged data 
     528         IF ( prodatqc%nvar == 2 ) THEN 
     529            ! For the second variable 
     530            IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
     531 
     532               zobsk(:) = obfillflt 
     533 
     534               IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 
     535 
     536                  IF ( idayend == 0 )  THEN 
     537                     ! Daily averaged data 
     538 
     539                     ! vertically interpolate all 4 corners  
     540                     ista = prodatqc%npvsta(jobs,2)  
     541                     iend = prodatqc%npvend(jobs,2)  
     542                     inum_obs = iend - ista + 1  
     543                     ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
     544 
     545                     DO iin=1,2  
     546                        DO ijn=1,2  
     547 
     548                           IF ( k1dint == 1 ) THEN  
     549                              CALL obs_int_z1d_spl( kpk, &  
     550                                 &     zinm2(iin,ijn,:,iobs), &  
     551                                 &     zobs2k, zgdept(iin,ijn,:,iobs), &  
     552                                 &     zmask2(iin,ijn,:,iobs))  
     553                           ENDIF  
     554 
     555                           CALL obs_level_search(kpk, &  
     556                              &    zgdept(iin,ijn,:,iobs), &  
     557                              &    inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
     558                              &    iv_indic)  
     559 
     560                           CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
     561                              &    prodatqc%var(2)%vdep(ista:iend), &  
     562                              &    zinm2(iin,ijn,:,iobs), &  
     563                              &    zobs2k, interp_corner(iin,ijn,:), &  
     564                              &    zgdept(iin,ijn,:,iobs), &  
     565                              &    zmask2(iin,ijn,:,iobs))  
     566 
     567                        ENDDO  
     568                     ENDDO  
     569 
     570                  ENDIF !idayend 
     571 
     572               ELSE    
     573 
     574                  ! Point data  
    523575 
    524576                  ! vertically interpolate all 4 corners  
     
    526578                  iend = prodatqc%npvend(jobs,2)  
    527579                  inum_obs = iend - ista + 1  
    528                   ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
    529  
     580                  ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     581                  DO iin=1,2   
     582                     DO ijn=1,2  
     583 
     584                        IF ( k1dint == 1 ) THEN  
     585                           CALL obs_int_z1d_spl( kpk, &  
     586                              &    zint2(iin,ijn,:,iobs),&  
     587                              &    zobs2k, zgdept(iin,ijn,:,iobs), &  
     588                              &    zmask2(iin,ijn,:,iobs))  
     589 
     590                        ENDIF  
     591 
     592                        CALL obs_level_search(kpk, &  
     593                            &        zgdept(iin,ijn,:,iobs),&  
     594                            &        inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
     595                            &        iv_indic)  
     596 
     597                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
     598                            &          prodatqc%var(2)%vdep(ista:iend),     &  
     599                            &          zint2(iin,ijn,:,iobs),            &  
     600                            &          zobs2k,interp_corner(iin,ijn,:), &  
     601                            &          zgdept(iin,ijn,:,iobs),         &  
     602                            &          zmask2(iin,ijn,:,iobs) )       
     603 
     604                     ENDDO  
     605                  ENDDO  
     606 
     607               ENDIF  
     608 
     609               !-------------------------------------------------------------  
     610               ! Compute the horizontal interpolation for every profile level  
     611               !-------------------------------------------------------------  
     612 
     613               DO ikn=1,inum_obs  
     614                  iend=ista+ikn-1 
     615 
     616                  zweig(:,:,1) = 0._wp  
     617 
     618                  ! This code forces the horizontal weights to be   
     619                  ! zero IF the observation is below the bottom of the   
     620                  ! corners of the interpolation nodes, Or if it is in   
     621                  ! the mask. This is important for observations near   
     622                  ! steep bathymetry  
    530623                  DO iin=1,2  
    531624                     DO ijn=1,2  
    532625 
    533                         IF ( k1dint == 1 ) THEN  
    534                            CALL obs_int_z1d_spl( kpk, &  
    535                               &     zinm2(iin,ijn,:,iobs), &  
    536                               &     zobs2k, zgdept(iin,ijn,:,iobs), &  
    537                               &     zmask2(iin,ijn,:,iobs))  
    538                         ENDIF  
    539         
    540                         CALL obs_level_search(kpk, &  
    541                            &    zgdept(iin,ijn,:,iobs), &  
    542                            &    inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
    543                            &    iv_indic)  
    544  
    545                         CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
    546                            &    prodatqc%var(2)%vdep(ista:iend), &  
    547                            &    zinm2(iin,ijn,:,iobs), &  
    548                            &    zobs2k, interp_corner(iin,ijn,:), &  
    549                            &    zgdept(iin,ijn,:,iobs), &  
    550                            &    zmask2(iin,ijn,:,iobs))  
    551         
     626                        depth_loop2: DO ik=kpk,2,-1  
     627                           IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
     628 
     629                              zweig(iin,ijn,1) = &   
     630                                 & zweig2(iin,ijn,1) * &  
     631                                 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
     632                                 &  - prodatqc%var(2)%vdep(iend)),0._wp)  
     633 
     634                              EXIT depth_loop2  
     635 
     636                           ENDIF  
     637 
     638                        ENDDO depth_loop2  
     639 
    552640                     ENDDO  
    553641                  ENDDO  
    554642 
    555                ENDIF !idayend 
    556  
    557             ELSE    
    558  
    559                ! Point data  
    560       
    561                ! vertically interpolate all 4 corners  
    562                ista = prodatqc%npvsta(jobs,2)  
    563                iend = prodatqc%npvend(jobs,2)  
    564                inum_obs = iend - ista + 1  
    565                ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
    566                DO iin=1,2   
    567                   DO ijn=1,2  
    568                      
    569                      IF ( k1dint == 1 ) THEN  
    570                         CALL obs_int_z1d_spl( kpk, &  
    571                            &    zint2(iin,ijn,:,iobs),&  
    572                            &    zobs2k, zgdept(iin,ijn,:,iobs), &  
    573                            &    zmask2(iin,ijn,:,iobs))  
    574    
    575                      ENDIF  
    576         
    577                      CALL obs_level_search(kpk, &  
    578                          &        zgdept(iin,ijn,:,iobs),&  
    579                          &        inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
    580                          &        iv_indic)  
    581  
    582                      CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
    583                          &          prodatqc%var(2)%vdep(ista:iend),     &  
    584                          &          zint2(iin,ijn,:,iobs),            &  
    585                          &          zobs2k,interp_corner(iin,ijn,:), &  
    586                          &          zgdept(iin,ijn,:,iobs),         &  
    587                          &          zmask2(iin,ijn,:,iobs) )       
    588           
    589                   ENDDO  
     643                  CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
     644                     &              prodatqc%var(2)%vmod(iend:iend) )  
     645 
     646                     ! Set QC flag for any observations found below the bottom 
     647                     ! needed as the check here is more strict than that in obs_prep 
     648                  IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 
     649 
    590650               ENDDO  
    591               
     651 
     652               DEALLOCATE(interp_corner,iv_indic)  
     653 
    592654            ENDIF  
    593  
    594             !-------------------------------------------------------------  
    595             ! Compute the horizontal interpolation for every profile level  
    596             !-------------------------------------------------------------  
    597               
    598             DO ikn=1,inum_obs  
    599                iend=ista+ikn-1 
    600                    
    601                zweig(:,:,1) = 0._wp  
    602     
    603                ! This code forces the horizontal weights to be   
    604                ! zero IF the observation is below the bottom of the   
    605                ! corners of the interpolation nodes, Or if it is in   
    606                ! the mask. This is important for observations near   
    607                ! steep bathymetry  
    608                DO iin=1,2  
    609                   DO ijn=1,2  
    610       
    611                      depth_loop2: DO ik=kpk,2,-1  
    612                         IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
    613                              
    614                            zweig(iin,ijn,1) = &   
    615                               & zweig2(iin,ijn,1) * &  
    616                               & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
    617                               &  - prodatqc%var(2)%vdep(iend)),0._wp)  
    618                              
    619                            EXIT depth_loop2  
    620  
    621                         ENDIF  
    622  
    623                      ENDDO depth_loop2  
    624       
    625                   ENDDO  
    626                ENDDO  
    627     
    628                CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
    629                   &              prodatqc%var(2)%vmod(iend:iend) )  
    630  
    631                   ! Set QC flag for any observations found below the bottom 
    632                   ! needed as the check here is more strict than that in obs_prep 
    633                IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 
    634   
    635             ENDDO  
    636   
    637             DEALLOCATE(interp_corner,iv_indic)  
    638            
    639655         ENDIF  
    640656 
  • branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

    r8222 r9186  
    403403      CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,1), profdata%mj(:,1), & 
    404404         &              profdata%nqc,     igrdobs                         ) 
    405       CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,2), profdata%mj(:,2), & 
    406          &              profdata%nqc,     igrdobs                         ) 
     405      IF ( ld_var2 ) THEN 
     406         CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,2), profdata%mj(:,2), & 
     407            &              profdata%nqc,     igrdobs                         ) 
     408      ENDIF 
    407409 
    408410      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 
     
    441443      CALL obs_mpp_sum_integer( ibdyv1obs, ibdyv1obsmpp ) 
    442444 
    443       ! Variable 2 
    444       CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(2),   & 
    445          &                 profdata%npvsta(:,2),  profdata%npvend(:,2), & 
    446          &                 jpi,                   jpj,                  & 
    447          &                 jpk,                                         & 
    448          &                 profdata%mi,           profdata%mj,          &  
    449          &                 profdata%var(2)%mvk,                         & 
    450          &                 profdata%rlam,         profdata%rphi,        & 
    451          &                 profdata%var(2)%vdep,                        & 
    452          &                 pglam2,                pgphi2,               & 
    453          &                 gdept_1d,              zmask2,               & 
    454          &                 profdata%nqc,          profdata%var(2)%nvqc, & 
    455          &                 iosdv2obs,             ilanv2obs,            & 
    456          &                 inlav2obs,             ld_nea,               & 
    457          &                 ibdyv2obs,             ld_bound_reject,      & 
    458          &                 iqc_cutoff       ) 
    459  
    460       CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp ) 
    461       CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp ) 
    462       CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp ) 
    463       CALL obs_mpp_sum_integer( ibdyv2obs, ibdyv2obsmpp ) 
     445      IF ( ld_var2 ) THEN 
     446         ! Variable 2 
     447         CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(2),   & 
     448            &                 profdata%npvsta(:,2),  profdata%npvend(:,2), & 
     449            &                 jpi,                   jpj,                  & 
     450            &                 jpk,                                         & 
     451            &                 profdata%mi,           profdata%mj,          &  
     452            &                 profdata%var(2)%mvk,                         & 
     453            &                 profdata%rlam,         profdata%rphi,        & 
     454            &                 profdata%var(2)%vdep,                        & 
     455            &                 pglam2,                pgphi2,               & 
     456            &                 gdept_1d,              zmask2,               & 
     457            &                 profdata%nqc,          profdata%var(2)%nvqc, & 
     458            &                 iosdv2obs,             ilanv2obs,            & 
     459            &                 inlav2obs,             ld_nea,               & 
     460            &                 ibdyv2obs,             ld_bound_reject,      & 
     461            &                 iqc_cutoff       ) 
     462 
     463         CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp ) 
     464         CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp ) 
     465         CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp ) 
     466         CALL obs_mpp_sum_integer( ibdyv2obs, ibdyv2obsmpp ) 
     467      ENDIF 
    464468 
    465469      ! ----------------------------------------------------------------------- 
     
    535539         WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted                             = ', & 
    536540            &            prodatqc%nvprotmpp(1) 
    537          WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data outside space domain       = ', & 
    538             &            iosdv2obsmpp 
    539          WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data at land points             = ', & 
    540             &            ilanv2obsmpp 
    541          IF (ld_nea) THEN 
    542             WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (removed) = ',& 
    543                &            inlav2obsmpp 
    544          ELSE 
    545             WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (kept)    = ',& 
    546                &            inlav2obsmpp 
    547          ENDIF 
    548          IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 
    549             WRITE(numout,*) ' V observation rejected since U rejected     = ', & 
    550                &            iuvchkv 
    551          ENDIF 
    552          WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near open boundary (removed) = ',& 
    553                &            ibdyv2obsmpp 
    554          WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted                             = ', & 
    555             &            prodatqc%nvprotmpp(2) 
     541         IF ( ld_var2 ) THEN 
     542            WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data outside space domain       = ', & 
     543               &            iosdv2obsmpp 
     544            WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data at land points             = ', & 
     545               &            ilanv2obsmpp 
     546            IF (ld_nea) THEN 
     547               WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (removed) = ',& 
     548                  &            inlav2obsmpp 
     549            ELSE 
     550               WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (kept)    = ',& 
     551                  &            inlav2obsmpp 
     552            ENDIF 
     553            IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 
     554               WRITE(numout,*) ' V observation rejected since U rejected     = ', & 
     555                  &            iuvchkv 
     556            ENDIF 
     557            WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near open boundary (removed) = ',& 
     558                  &            ibdyv2obsmpp 
     559            WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted                             = ', & 
     560               &            prodatqc%nvprotmpp(2) 
     561         ENDIF 
    556562 
    557563         WRITE(numout,*) 
    558564         WRITE(numout,*) ' Number of observations per time step :' 
    559565         WRITE(numout,*) 
    560          WRITE(numout,'(10X,A,5X,A,5X,A,A)')'Time step','Profiles', & 
    561             &                               '     '//prodatqc%cvars(1)//'     ', & 
    562             &                               '     '//prodatqc%cvars(2)//'     ' 
     566         IF ( ld_var2 ) THEN 
     567            WRITE(numout,'(10X,A,5X,A,5X,A,A)')'Time step','Profiles', & 
     568               &                               '   '//prodatqc%cvars(1)//'     ', & 
     569               &                               '   '//prodatqc%cvars(2)//'     ' 
     570         ELSE 
     571            WRITE(numout,'(10X,A,5X,A,5X,A)')'Time step','Profiles', & 
     572               &                               '   '//prodatqc%cvars(1)//'     ' 
     573         ENDIF 
    563574         WRITE(numout,998) 
    564575      ENDIF 
     
    588599         DO jstp = nit000 - 1, nitend 
    589600            inrc = jstp - nit000 + 2 
    590             WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), & 
    591                &                    prodatqc%nvstpmpp(inrc,1), & 
    592                &                    prodatqc%nvstpmpp(inrc,2) 
     601            IF ( ld_var2 ) THEN 
     602               WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), & 
     603                  &                    prodatqc%nvstpmpp(inrc,1), & 
     604                  &                    prodatqc%nvstpmpp(inrc,2) 
     605            ELSE 
     606               WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), & 
     607                  &                    prodatqc%nvstpmpp(inrc,1) 
     608            ENDIF 
    593609         END DO 
    594610      ENDIF 
  • branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90

    r7992 r9186  
    219219               &                ldgrid = .TRUE. ) 
    220220 
    221             IF ( inpfiles(jj)%nvar < 2 ) THEN 
     221            IF ( inpfiles(jj)%nvar /= kvars ) THEN 
    222222               CALL ctl_stop( 'Feedback format error: ', & 
    223                   &           ' less than 2 vars in profile file' ) 
     223                  &           ' unexpected number of vars in profile file' ) 
    224224            ENDIF 
    225225 
     
    320320            ALLOCATE( iobsj1(inowin) ) 
    321321            ALLOCATE( iproc1(inowin) ) 
    322             ALLOCATE( iobsi2(inowin) ) 
    323             ALLOCATE( iobsj2(inowin) ) 
    324             ALLOCATE( iproc2(inowin) ) 
     322            IF ( kvars == 2 ) THEN 
     323               ALLOCATE( iobsi2(inowin) ) 
     324               ALLOCATE( iobsj2(inowin) ) 
     325               ALLOCATE( iproc2(inowin) ) 
     326            ENDIF 
    325327            inowin = 0 
    326328            DO ji = 1, inpfiles(jj)%nobs 
    327329               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 
     330               IF ( kvars == 2 ) THEN 
     331                  IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     332                     & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
     333               ELSE 
     334                  IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) ) CYCLE 
     335               ENDIF 
    330336               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    331337                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    347353               CALL obs_grid_search( inowin, zlam, zphi, iobsi2, iobsj2, & 
    348354                  &                  iproc2, 'V' ) 
     355            ELSE 
     356               CALL obs_grid_search( inowin, zlam, zphi, iobsi1, iobsj1, & 
     357                  &                  iproc1, 'T' ) 
    349358            ENDIF 
    350359 
     
    352361            DO ji = 1, inpfiles(jj)%nobs 
    353362               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 
     363               IF ( kvars == 2 ) THEN 
     364                  IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     365                     & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
     366               ELSE 
     367                  IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) ) CYCLE 
     368               ENDIF 
    356369               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    357370                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    360373                  inpfiles(jj)%iobsi(ji,1) = iobsi1(inowin) 
    361374                  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') 
     375                  IF ( kvars == 2 ) THEN 
     376                     inpfiles(jj)%iproc(ji,2) = iproc2(inowin) 
     377                     inpfiles(jj)%iobsi(ji,2) = iobsi2(inowin) 
     378                     inpfiles(jj)%iobsj(ji,2) = iobsj2(inowin) 
     379                     IF ( inpfiles(jj)%iproc(ji,1) /= & 
     380                        & inpfiles(jj)%iproc(ji,2) ) THEN 
     381                        CALL ctl_stop( 'Error in obs_read_prof:', & 
     382                           & 'var1 and var2 observation on different processors') 
     383                     ENDIF 
    369384                  ENDIF 
    370385               ENDIF 
    371386            END DO 
    372             DEALLOCATE( zlam, zphi, iobsi1, iobsj1, iproc1, iobsi2, iobsj2, iproc2 ) 
     387            DEALLOCATE( zlam, zphi, iobsi1, iobsj1, iproc1 ) 
     388            IF ( kvars == 2 ) THEN 
     389               DEALLOCATE( iobsi2, iobsj2, iproc2 ) 
     390            ENDIF 
    373391 
    374392            DO ji = 1, inpfiles(jj)%nobs 
    375393               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 
     394               IF ( kvars == 2 ) THEN 
     395                  IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     396                     & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
     397               ELSE 
     398                  IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) ) CYCLE 
     399               ENDIF 
    378400               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    379401                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    407429                     IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    408430                        & 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. 
     431                     IF ( kvars == 2 ) THEN 
     432                        IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
     433                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     434                           &    ldvar1 ) .OR. & 
     435                           & ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
     436                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     437                           &     ldvar2 ) ) THEN 
     438                           ip3dt = ip3dt + 1 
     439                           llvalprof = .TRUE. 
     440                        ENDIF 
     441                     ELSE 
     442                        IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
     443                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     444                           &    ldvar1 ) ) THEN 
     445                           ip3dt = ip3dt + 1 
     446                           llvalprof = .TRUE. 
     447                        ENDIF 
    417448                     ENDIF 
    418449                  END DO loop_p_count 
     
    438469         DO ji = 1, inpfiles(jj)%nobs 
    439470            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 
     471            IF ( kvars == 2 ) THEN 
     472               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     473                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
     474            ELSE 
     475               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) ) CYCLE 
     476            ENDIF 
    442477            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    443478               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    453488         DO ji = 1, inpfiles(jj)%nobs 
    454489            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 
     490            IF ( kvars == 2 ) THEN 
     491               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     492                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
     493            ELSE 
     494               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) ) CYCLE 
     495            ENDIF 
    457496            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    458497               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    501540         ji = iprofidx(iindx(jk)) 
    502541 
    503             IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     542         IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     543         IF ( kvars == 2 ) THEN 
    504544            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
    505545               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
     546         ELSE 
     547            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) ) CYCLE 
     548         ENDIF 
    506549 
    507550         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  & 
     
    519562 
    520563            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 
     564            IF ( kvars == 2 ) THEN 
     565               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     566                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
     567            ELSE 
     568               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) ) CYCLE 
     569            ENDIF 
    523570 
    524571            loop_prof : DO ij = 1, inpfiles(jj)%nlev 
     
    535582               ENDIF 
    536583 
    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  
     584               IF ( kvars == 2 ) THEN 
     585                  IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
     586                     & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
     587 
     588                     llvalprof = .TRUE.  
     589                     EXIT loop_prof 
     590 
     591                  ENDIF 
    543592               ENDIF 
    544593 
     
    575624               profdata%mi  (iprof,1) = inpfiles(jj)%iobsi(ji,1) 
    576625               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) 
     626               IF ( kvars == 2 ) THEN 
     627                  profdata%mi  (iprof,2) = inpfiles(jj)%iobsi(ji,2) 
     628                  profdata%mj  (iprof,2) = inpfiles(jj)%iobsj(ji,2) 
     629               ENDIF 
    579630 
    580631               ! Profile WMO number 
     
    616667                  IF (ldsatt) THEN 
    617668 
    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 
     669                     IF ( kvars == 2 ) THEN 
     670                        IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
     671                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     672                           &    ldvar1 ) .OR. & 
     673                           & ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
     674                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     675                           &   ldvar2 ) ) THEN 
     676                           ip3dt = ip3dt + 1 
     677                        ELSE 
     678                           CYCLE 
     679                        ENDIF 
    625680                     ELSE 
    626                         CYCLE 
     681                        IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
     682                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     683                           &    ldvar1 ) ) THEN 
     684                           ip3dt = ip3dt + 1 
     685                        ELSE 
     686                           CYCLE 
     687                        ENDIF 
    627688                     ENDIF 
    628689 
     
    693754                  ENDIF 
    694755 
    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 
     756                  IF ( kvars == 2 ) THEN 
     757                     IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
     758                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2)    .AND. & 
     759                        &   ldvar2 ) .OR. ldsatt ) THEN 
     760 
     761                        IF (ldsatt) THEN 
     762 
     763                           ivar2t = ip3dt 
     764 
     765                        ELSE 
     766 
     767                           ivar2t = ivar2t + 1 
     768 
     769                        ENDIF 
     770 
     771                        ! Depth of var2 observation 
     772                        profdata%var(2)%vdep(ivar2t) = & 
     773                           &                inpfiles(jj)%pdep(ij,ji) 
     774 
     775                        ! Depth of var2 observation QC 
     776                        profdata%var(2)%idqc(ivar2t) = & 
     777                           &                inpfiles(jj)%idqc(ij,ji) 
     778 
     779                        ! Depth of var2 observation QC flags 
     780                        profdata%var(2)%idqcf(:,ivar2t) = & 
     781                           &                inpfiles(jj)%idqcf(:,ij,ji) 
     782 
     783                        ! Profile index 
     784                        profdata%var(2)%nvpidx(ivar2t) = iprof 
     785 
     786                        ! Vertical index in original profile 
     787                        profdata%var(2)%nvlidx(ivar2t) = ij 
     788 
     789                        ! Profile var2 value 
     790                     IF (  ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) ) .AND. & 
     791                       &   ( .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2)    )  ) THEN 
     792                           profdata%var(2)%vobs(ivar2t) = & 
     793                              &                inpfiles(jj)%pob(ij,ji,2) 
     794                           IF ( ldmod ) THEN 
     795                              profdata%var(2)%vmod(ivar2t) = & 
     796                                 &                inpfiles(jj)%padd(ij,ji,1,2) 
     797                           ENDIF 
     798                           ! Count number of profile var2 data as function of type 
     799                           itypvar2( profdata%ntyp(iprof) + 1 ) = & 
     800                              & itypvar2( profdata%ntyp(iprof) + 1 ) + 1 
     801                        ELSE 
     802                           profdata%var(2)%vobs(ivar2t) = fbrmdi 
     803                        ENDIF 
     804 
     805                        ! Profile var2 qc 
     806                        profdata%var(2)%nvqc(ivar2t) = & 
     807                           & inpfiles(jj)%ivlqc(ij,ji,2) 
     808 
     809                        ! Profile var2 qc flags 
     810                        profdata%var(2)%nvqcf(:,ivar2t) = & 
     811                           & inpfiles(jj)%ivlqcf(:,ij,ji,2) 
    706812 
    707813                     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  
    751814                  ENDIF 
    752815 
     
    764827 
    765828      CALL obs_mpp_sum_integer ( ivar1t0, ivar1tmpp ) 
    766       CALL obs_mpp_sum_integer ( ivar2t0, ivar2tmpp ) 
     829      IF ( kvars == 2 ) THEN 
     830         CALL obs_mpp_sum_integer ( ivar2t0, ivar2tmpp ) 
     831      ENDIF 
    767832      CALL obs_mpp_sum_integer ( ip3dt,   ip3dtmpp  ) 
    768833 
    769834      CALL obs_mpp_sum_integers( itypvar1, itypvar1mpp, ntyp1770 + 1 ) 
    770       CALL obs_mpp_sum_integers( itypvar2, itypvar2mpp, ntyp1770 + 1 ) 
     835      IF ( kvars == 2 ) THEN 
     836         CALL obs_mpp_sum_integers( itypvar2, itypvar2mpp, ntyp1770 + 1 ) 
     837      ENDIF 
    771838 
    772839      !----------------------------------------------------------------------- 
     
    795862            & '---------------------------------------------------------------' 
    796863         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 
    805          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,*)  
     864         IF ( kvars == 2 ) THEN 
     865            WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(2) ) 
     866            WRITE(numout,'(1X,A)') '------------------------' 
     867            DO ji = 0, ntyp1770 
     868               IF ( itypvar2mpp(ji+1) > 0 ) THEN 
     869                  WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), & 
     870                     & cwmonam1770(ji)(1:52),' = ', & 
     871                     & itypvar2mpp(ji+1) 
     872               ENDIF 
     873            END DO 
     874            WRITE(numout,'(1X,A)') & 
     875               & '---------------------------------------------------------------' 
     876            WRITE(numout,'(1X,A55,I8)') & 
     877               & 'Total profile data for variable '//TRIM( profdata%cvars(2) )// & 
     878               & '             = ', ivar2tmpp 
     879            WRITE(numout,'(1X,A)') & 
     880               & '---------------------------------------------------------------' 
     881            WRITE(numout,*)  
     882         ENDIF 
    814883      ENDIF 
    815884 
    816885      IF (ldsatt) THEN 
    817886         profdata%nvprot(1)    = ip3dt 
    818          profdata%nvprot(2)    = ip3dt 
    819887         profdata%nvprotmpp(1) = ip3dtmpp 
    820          profdata%nvprotmpp(2) = ip3dtmpp 
     888         IF ( kvars == 2 ) THEN 
     889            profdata%nvprot(2)    = ip3dt 
     890            profdata%nvprotmpp(2) = ip3dtmpp 
     891         ENDIF 
    821892      ELSE 
    822893         profdata%nvprot(1)    = ivar1t 
    823          profdata%nvprot(2)    = ivar2t 
    824894         profdata%nvprotmpp(1) = ivar1tmpp 
    825          profdata%nvprotmpp(2) = ivar2tmpp 
     895         IF ( kvars == 2 ) THEN 
     896            profdata%nvprot(2)    = ivar2t 
     897            profdata%nvprotmpp(2) = ivar2tmpp 
     898         ENDIF 
    826899      ENDIF 
    827900      profdata%nprof        = iprof 
  • branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90

    r9016 r9186  
    8484      CHARACTER(LEN=40) :: clfname 
    8585      CHARACTER(LEN=10) :: clfiletype 
     86      CHARACTER(LEN=ilenlong) :: cllongname  ! Long name of variable 
     87      CHARACTER(LEN=ilenunit) :: clunits     ! Units of variable 
     88      CHARACTER(LEN=ilengrid) :: clgrid      ! Grid of variable 
    8689      INTEGER :: ilevel 
    8790      INTEGER :: jvar 
     
    111114      ! Find maximum level 
    112115      ilevel = 0 
    113       DO jvar = 1, 2 
     116      DO jvar = 1, profdata%nvar 
    114117         ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) ) 
    115118      END DO 
     
    175178         END DO 
    176179 
     180      CASE('PLCHLTOT') 
     181 
     182         clfiletype = 'plchltotfb' 
     183         cllongname = 'log10(chlorophyll concentration)' 
     184         clunits    = 'log10(mg/m3)' 
     185         clgrid     = 'T' 
     186 
     187      CASE('PCHLTOT') 
     188 
     189         clfiletype = 'pchltotfb' 
     190         cllongname = 'chlorophyll concentration' 
     191         clunits    = 'mg/m3' 
     192         clgrid     = 'T' 
     193 
     194      CASE('PNO3') 
     195 
     196         clfiletype = 'pno3fb' 
     197         cllongname = 'nitrate' 
     198         clunits    = 'mmol/m3' 
     199         clgrid     = 'T' 
     200 
     201      CASE('PSI4') 
     202 
     203         clfiletype = 'psi4fb' 
     204         cllongname = 'silicate' 
     205         clunits    = 'mmol/m3' 
     206         clgrid     = 'T' 
     207 
     208      CASE('PPO4') 
     209 
     210         clfiletype = 'ppo4fb' 
     211         cllongname = 'phosphate' 
     212         clunits    = 'mmol/m3' 
     213         clgrid     = 'T' 
     214 
     215      CASE('PDIC') 
     216 
     217         clfiletype = 'pdicfb' 
     218         cllongname = 'dissolved inorganic carbon' 
     219         clunits    = 'mmol/m3' 
     220         clgrid     = 'T' 
     221 
     222      CASE('PALK') 
     223 
     224         clfiletype = 'palkfb' 
     225         cllongname = 'alkalinity' 
     226         clunits    = 'meq/m3' 
     227         clgrid     = 'T' 
     228 
     229      CASE('PPH') 
     230 
     231         clfiletype = 'pphfb' 
     232         cllongname = 'pH' 
     233         clunits    = '-' 
     234         clgrid     = 'T' 
     235 
     236      CASE('PO2') 
     237 
     238         clfiletype = 'po2fb' 
     239         cllongname = 'dissolved oxygen' 
     240         clunits    = 'mmol/m3' 
     241         clgrid     = 'T' 
     242 
    177243      END SELECT 
     244       
     245      IF ( ( TRIM(profdata%cvars(1)) /= 'POTM' ) .AND. & 
     246         & ( TRIM(profdata%cvars(1)) /= 'UVEL' ) ) THEN 
     247         CALL alloc_obfbdata( fbdata, 1, profdata%nprof, ilevel, & 
     248            &                 1 + iadd, iext, .TRUE. ) 
     249         fbdata%cname(1)      = profdata%cvars(1) 
     250         fbdata%coblong(1)    = cllongname 
     251         fbdata%cobunit(1)    = clunits 
     252         fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(cllongname) 
     253         fbdata%caddunit(1,1) = clunits 
     254         fbdata%cgrid(:)      = clgrid 
     255         DO je = 1, iext 
     256            fbdata%cextname(je) = pext%cdname(je) 
     257            fbdata%cextlong(je) = pext%cdlong(je,1) 
     258            fbdata%cextunit(je) = pext%cdunit(je,1) 
     259         END DO 
     260         DO ja = 1, iadd 
     261            fbdata%caddname(1+ja) = padd%cdname(ja) 
     262            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
     263            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     264         END DO 
     265      ENDIF 
    178266 
    179267      fbdata%caddname(1)   = 'Hx' 
     
    228316            &           krefdate = 19500101 ) 
    229317         ! Reform the profiles arrays for output 
    230          DO jvar = 1, 2 
     318         DO jvar = 1, profdata%nvar 
    231319            DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar) 
    232320               ik = profdata%var(jvar)%nvlidx(jk) 
Note: See TracChangeset for help on using the changeset viewer.