Changeset 9023 for branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS
- Timestamp:
- 2017-12-13T18:08:50+01:00 (6 years ago)
- Location:
- branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS
- Files:
-
- 12 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
r9019 r9023 46 46 LOGICAL, PUBLIC :: ln_diaobs !: Logical switch for the obs operator 47 47 LOGICAL :: ln_sstnight !: Logical switch for night mean SST obs 48 48 LOGICAL :: ln_sla_fp_indegs !: T=> SLA obs footprint size specified in degrees, F=> in metres 49 LOGICAL :: ln_sst_fp_indegs !: T=> SST obs footprint size specified in degrees, F=> in metres 50 LOGICAL :: ln_sss_fp_indegs !: T=> SSS obs footprint size specified in degrees, F=> in metres 51 LOGICAL :: ln_sic_fp_indegs !: T=> sea-ice obs footprint size specified in degrees, F=> in metres 52 53 REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint (metres) 54 REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint (metres) 55 REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint (metres) 56 REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint (metres) 57 REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint (metres) 58 REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint (metres) 59 REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint (metres) 60 REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint (metres) 61 49 62 INTEGER :: nn_1dint !: Vertical interpolation method 50 INTEGER :: nn_2dint !: Horizontal interpolation method 63 INTEGER :: nn_2dint !: Default horizontal interpolation method 64 INTEGER :: nn_2dint_sla !: SLA horizontal interpolation method 65 INTEGER :: nn_2dint_sst !: SST horizontal interpolation method 66 INTEGER :: nn_2dint_sss !: SSS horizontal interpolation method 67 INTEGER :: nn_2dint_sic !: Seaice horizontal interpolation method 51 68 INTEGER, DIMENSION(imaxavtypes) :: nn_profdavtypes !: Profile data types representing a daily average 52 69 INTEGER :: nproftypes !: Number of profile obs types 53 70 INTEGER :: nsurftypes !: Number of surface obs types 54 INTEGER, DIMENSION(:), ALLOCATABLE :: nvarsprof, nvarssurf !: Number of profile & surface variables 55 INTEGER, DIMENSION(:), ALLOCATABLE :: nextrprof, nextrsurf !: Number of profile & surface extra variables 56 INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: sstbias_type !: SST bias type 57 TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: surfdata, surfdataqc !: Initial surface data before & after quality control 58 TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: profdata, profdataqc !: Initial profile data before & after quality control 71 INTEGER, DIMENSION(:), ALLOCATABLE :: & 72 & nvarsprof, & !: Number of profile variables 73 & nvarssurf !: Number of surface variables 74 INTEGER, DIMENSION(:), ALLOCATABLE :: & 75 & nextrprof, & !: Number of profile extra variables 76 & nextrsurf !: Number of surface extra variables 77 INTEGER, DIMENSION(:), ALLOCATABLE :: & 78 & n2dintsurf !: Interpolation option for surface variables 79 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 80 & zavglamscl, & !: E/W diameter of averaging footprint for surface variables 81 & zavgphiscl !: N/S diameter of averaging footprint for surface variables 82 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 83 & lfpindegs, & !: T=> surface obs footprint size specified in degrees, F=> in metres 84 & llnightav !: Logical for calculating night-time averages 85 86 TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: & 87 & surfdata, & !: Initial surface data 88 & surfdataqc !: Surface data after quality control 89 TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: & 90 & profdata, & !: Initial profile data 91 & profdataqc !: Profile data after quality control 59 92 60 93 CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: cobstypesprof, cobstypessurf !: Profile & surface obs types … … 95 128 & cn_profbfiles, & ! T/S profile input filenames 96 129 & cn_sstfbfiles, & ! Sea surface temperature input filenames 130 & cn_sssfbfiles, & ! Sea surface salinity input filenames 97 131 & cn_slafbfiles, & ! Sea level anomaly input filenames 98 132 & cn_sicfbfiles, & ! Seaice concentration input filenames 99 133 & cn_velfbfiles, & ! Velocity profile input filenames 100 & cn_sstbias _files ! SST bias input filenames134 & cn_sstbiasfiles ! SST bias input filenames 101 135 CHARACTER(LEN=128) :: & 102 136 & cn_altbiasfile ! Altimeter bias input filename … … 109 143 LOGICAL :: ln_sla ! Logical switch for sea level anomalies 110 144 LOGICAL :: ln_sst ! Logical switch for sea surface temperature 145 LOGICAL :: ln_sss ! Logical switch for sea surface salinity 111 146 LOGICAL :: ln_sic ! Logical switch for sea ice concentration 112 147 LOGICAL :: ln_vel3d ! Logical switch for velocity (u,v) obs … … 116 151 LOGICAL :: ln_ignmis ! Logical switch for ignoring missing files 117 152 LOGICAL :: ln_s_at_t ! Logical switch to compute model S at T obs 153 LOGICAL :: ln_bound_reject ! Logical to remove obs near boundaries in LAMs. 118 154 LOGICAL :: llvar1 ! Logical for profile variable 1 119 155 LOGICAL :: llvar2 ! Logical for profile variable 1 120 LOGICAL :: llnightav ! Logical for calculating night-time averages121 156 LOGICAL, DIMENSION(jpmaxnfiles) :: lmask ! Used for finding number of sstbias files 122 157 … … 134 169 135 170 NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, & 136 & ln_sst, ln_sic, ln_vel3d, & 137 & ln_altbias, ln_nea, ln_grid_global, & 138 & ln_grid_search_lookup, & 139 & ln_ignmis, ln_s_at_t, ln_sstnight, & 171 & ln_sst, ln_sic, ln_sss, ln_vel3d, & 172 & ln_altbias, ln_sstbias, ln_nea, & 173 & ln_grid_global, ln_grid_search_lookup, & 174 & ln_ignmis, ln_s_at_t, ln_bound_reject, & 175 & ln_sstnight, & 176 & ln_sla_fp_indegs, ln_sst_fp_indegs, & 177 & ln_sss_fp_indegs, ln_sic_fp_indegs, & 140 178 & cn_profbfiles, cn_slafbfiles, & 141 179 & cn_sstfbfiles, cn_sicfbfiles, & 142 & cn_velfbfiles, cn_altbiasfile, & 180 & cn_velfbfiles, cn_sssfbfiles, & 181 & cn_sstbiasfiles, cn_altbiasfile, & 143 182 & cn_gridsearchfile, rn_gridsearchres, & 144 & rn_dobsini, rn_dobsend, nn_1dint, nn_2dint, & 183 & rn_dobsini, rn_dobsend, & 184 & rn_sla_avglamscl, rn_sla_avgphiscl, & 185 & rn_sst_avglamscl, rn_sst_avgphiscl, & 186 & rn_sss_avglamscl, rn_sss_avgphiscl, & 187 & rn_sic_avglamscl, rn_sic_avgphiscl, & 188 & nn_1dint, nn_2dint, & 189 & nn_2dint_sla, nn_2dint_sst, & 190 & nn_2dint_sss, nn_2dint_sic, & 145 191 & nn_msshc, rn_mdtcorr, rn_mdtcutoff, & 146 & nn_profdavtypes , ln_sstbias, cn_sstbias_files192 & nn_profdavtypes 147 193 148 194 INTEGER :: jnumsstbias … … 157 203 ! Read namelist parameters 158 204 !----------------------------------------------------------------------- 159 160 !Initalise all values in namelist arrays161 ALLOCATE(sstbias_type(jpmaxnfiles))162 205 ! Some namelist arrays need initialising 163 206 cn_profbfiles(:) = '' … … 166 209 cn_sicfbfiles(:) = '' 167 210 cn_velfbfiles(:) = '' 168 cn_sstbias_files(:) = '' 211 cn_sssfbfiles(:) = '' 212 cn_sstbiasfiles(:) = '' 169 213 nn_profdavtypes(:) = -1 170 214 … … 187 231 RETURN 188 232 ENDIF 189 190 !----------------------------------------------------------------------- 191 ! Set up list of observation types to be used 192 ! and the files associated with each type 193 !----------------------------------------------------------------------- 194 195 nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 196 nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic /) ) 197 198 IF (ln_sstbias) THEN 199 lmask(:) = .FALSE. 200 WHERE (cn_sstbias_files(:) /= '') lmask(:) = .TRUE. 201 jnumsstbias = COUNT(lmask) 202 lmask(:) = .FALSE. 203 ENDIF 204 205 IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 206 IF(lwp) WRITE(numout,cform_war) 207 IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 208 & ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 209 & ' are set to .FALSE. so turning off calls to dia_obs' 210 nwarn = nwarn + 1 211 ln_diaobs = .FALSE. 212 RETURN 213 ENDIF 214 215 IF ( nproftypes > 0 ) THEN 216 217 ALLOCATE( cobstypesprof(nproftypes) ) 218 ALLOCATE( ifilesprof(nproftypes) ) 219 ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 220 221 jtype = 0 222 IF (ln_t3d .OR. ln_s3d) THEN 223 jtype = jtype + 1 224 clproffiles(jtype,:) = cn_profbfiles(:) 225 cobstypesprof(jtype) = 'prof ' 226 ifilesprof(jtype) = 0 227 DO jfile = 1, jpmaxnfiles 228 IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 229 ifilesprof(jtype) = ifilesprof(jtype) + 1 230 END DO 231 ENDIF 232 IF (ln_vel3d) THEN 233 jtype = jtype + 1 234 clproffiles(jtype,:) = cn_velfbfiles(:) 235 cobstypesprof(jtype) = 'vel ' 236 ifilesprof(jtype) = 0 237 DO jfile = 1, jpmaxnfiles 238 IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 239 ifilesprof(jtype) = ifilesprof(jtype) + 1 240 END DO 241 ENDIF 242 243 ENDIF 244 245 IF ( nsurftypes > 0 ) THEN 246 247 ALLOCATE( cobstypessurf(nsurftypes) ) 248 ALLOCATE( ifilessurf(nsurftypes) ) 249 ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 250 251 jtype = 0 252 IF (ln_sla) THEN 253 jtype = jtype + 1 254 clsurffiles(jtype,:) = cn_slafbfiles(:) 255 cobstypessurf(jtype) = 'sla ' 256 ifilessurf(jtype) = 0 257 DO jfile = 1, jpmaxnfiles 258 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 259 ifilessurf(jtype) = ifilessurf(jtype) + 1 260 END DO 261 ENDIF 262 IF (ln_sst) THEN 263 jtype = jtype + 1 264 clsurffiles(jtype,:) = cn_sstfbfiles(:) 265 cobstypessurf(jtype) = 'sst ' 266 ifilessurf(jtype) = 0 267 DO jfile = 1, jpmaxnfiles 268 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 269 ifilessurf(jtype) = ifilessurf(jtype) + 1 270 END DO 271 ENDIF 272 #if defined key_lim3 273 IF (ln_sic) THEN 274 jtype = jtype + 1 275 clsurffiles(jtype,:) = cn_sicfbfiles(:) 276 cobstypessurf(jtype) = 'sic ' 277 ifilessurf(jtype) = 0 278 DO jfile = 1, jpmaxnfiles 279 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 280 ifilessurf(jtype) = ifilessurf(jtype) + 1 281 END DO 282 ENDIF 283 #endif 284 285 ENDIF 286 287 !Write namelist settings to stdout 233 288 234 IF(lwp) THEN 289 235 WRITE(numout,*) … … 297 243 WRITE(numout,*) ' Logical switch for Sea Ice observations ln_sic = ', ln_sic 298 244 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 299 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ',ln_grid_global300 WRITE(numout,*) ' Logical switch for SST bias correction ln_sstbias = ', ln_sstbias301 WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup245 WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss 246 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ', ln_grid_global 247 WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup 302 248 IF (ln_grid_search_lookup) & 303 249 WRITE(numout,*) ' Grid search lookup file header cn_gridsearchfile = ', cn_gridsearchfile … … 307 253 WRITE(numout,*) ' Type of horizontal interpolation method nn_2dint = ', nn_2dint 308 254 WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea 255 WRITE(numout,*) ' Rejection of obs near open bdys ln_bound_reject = ', ln_bound_reject 309 256 WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc 310 257 WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr 311 258 WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff 312 259 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 260 WRITE(numout,*) ' Logical switch for sst bias ln_sstbias = ', ln_sstbias 313 261 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 314 262 WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes 315 263 WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight 316 WRITE(numout,*) ' Number of profile obs types: ',nproftypes 317 318 IF ( nproftypes > 0 ) THEN 319 DO jtype = 1, nproftypes 320 DO jfile = 1, ifilesprof(jtype) 321 WRITE(numout,'(1X,2A)') ' '//cobstypesprof(jtype)//' input observation file names = ', & 322 TRIM(clproffiles(jtype,jfile)) 323 END DO 324 END DO 264 ENDIF 265 !----------------------------------------------------------------------- 266 ! Set up list of observation types to be used 267 ! and the files associated with each type 268 !----------------------------------------------------------------------- 269 270 nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 271 nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss /) ) 272 273 IF (ln_sstbias) THEN 274 lmask(:) = .FALSE. 275 WHERE (cn_sstbiasfiles(:) /= '') lmask(:) = .TRUE. 276 jnumsstbias = COUNT(lmask) 277 lmask(:) = .FALSE. 278 ENDIF 279 280 IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 281 IF(lwp) WRITE(numout,cform_war) 282 IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 283 & ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 284 & ' are set to .FALSE. so turning off calls to dia_obs' 285 nwarn = nwarn + 1 286 ln_diaobs = .FALSE. 287 RETURN 288 ENDIF 289 290 IF ( nproftypes > 0 ) THEN 291 292 ALLOCATE( cobstypesprof(nproftypes) ) 293 ALLOCATE( ifilesprof(nproftypes) ) 294 ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 295 296 jtype = 0 297 IF (ln_t3d .OR. ln_s3d) THEN 298 jtype = jtype + 1 299 CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof ', & 300 & cn_profbfiles, ifilesprof, cobstypesprof, clproffiles ) 325 301 ENDIF 326 327 WRITE(numout,*)' Number of surface obs types: ',nsurftypes 328 IF ( nsurftypes > 0 ) THEN 329 DO jtype = 1, nsurftypes 330 DO jfile = 1, ifilessurf(jtype) 331 WRITE(numout,'(1X,2A)') ' '//cobstypessurf(jtype)//' input observation file names = ', & 332 TRIM(clsurffiles(jtype,jfile)) 333 END DO 334 END DO 302 IF (ln_vel3d) THEN 303 jtype = jtype + 1 304 CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel ', & 305 & cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles ) 335 306 ENDIF 336 WRITE(numout,*) '~~~~~~~~~~~~' 337 338 ENDIF 307 308 ENDIF 309 310 IF ( nsurftypes > 0 ) THEN 311 312 ALLOCATE( cobstypessurf(nsurftypes) ) 313 ALLOCATE( ifilessurf(nsurftypes) ) 314 ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 315 ALLOCATE(n2dintsurf(nsurftypes)) 316 ALLOCATE(zavglamscl(nsurftypes)) 317 ALLOCATE(zavgphiscl(nsurftypes)) 318 ALLOCATE(lfpindegs(nsurftypes)) 319 ALLOCATE(llnightav(nsurftypes)) 320 321 jtype = 0 322 IF (ln_sla) THEN 323 jtype = jtype + 1 324 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla ', & 325 & cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles ) 326 CALL obs_setinterpopts( nsurftypes, jtype, 'sla ', & 327 & nn_2dint, nn_2dint_sla, & 328 & rn_sla_avglamscl, rn_sla_avgphiscl, & 329 & ln_sla_fp_indegs, .FALSE., & 330 & n2dintsurf, zavglamscl, zavgphiscl, & 331 & lfpindegs, llnightav ) 332 ENDIF 333 IF (ln_sst) THEN 334 jtype = jtype + 1 335 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst ', & 336 & cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 337 CALL obs_setinterpopts( nsurftypes, jtype, 'sst ', & 338 & nn_2dint, nn_2dint_sst, & 339 & rn_sst_avglamscl, rn_sst_avgphiscl, & 340 & ln_sst_fp_indegs, ln_sstnight, & 341 & n2dintsurf, zavglamscl, zavgphiscl, & 342 & lfpindegs, llnightav ) 343 ENDIF 344 #if defined key_lim3 || defined key_cice 345 IF (ln_sic) THEN 346 jtype = jtype + 1 347 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic ', & 348 & cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 349 CALL obs_setinterpopts( nsurftypes, jtype, 'sic ', & 350 & nn_2dint, nn_2dint_sic, & 351 & rn_sic_avglamscl, rn_sic_avgphiscl, & 352 & ln_sic_fp_indegs, .FALSE., & 353 & n2dintsurf, zavglamscl, zavgphiscl, & 354 & lfpindegs, llnightav ) 355 ENDIF 356 #endif 357 IF (ln_sss) THEN 358 jtype = jtype + 1 359 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss ', & 360 & cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 361 CALL obs_setinterpopts( nsurftypes, jtype, 'sss ', & 362 & nn_2dint, nn_2dint_sss, & 363 & rn_sss_avglamscl, rn_sss_avgphiscl, & 364 & ln_sss_fp_indegs, .FALSE., & 365 & n2dintsurf, zavglamscl, zavgphiscl, & 366 & lfpindegs, llnightav ) 367 ENDIF 368 369 ENDIF 370 371 339 372 340 373 !----------------------------------------------------------------------- … … 356 389 ENDIF 357 390 358 IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 4) ) THEN391 IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 6 ) ) THEN 359 392 CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & 360 393 & ' is not available') … … 421 454 & jpi, jpj, jpk, & 422 455 & zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2, & 423 & ln_nea, kdailyavtypes = nn_profdavtypes ) 456 & ln_nea, ln_bound_reject, & 457 & kdailyavtypes = nn_profdavtypes ) 424 458 425 459 END DO … … 440 474 nvarssurf(jtype) = 1 441 475 nextrsurf(jtype) = 0 442 llnightav = .FALSE.476 llnightav(jtype) = .FALSE. 443 477 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 444 IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav = ln_sstnight478 IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav(jtype) = ln_sstnight 445 479 446 480 !Read in surface obs types … … 448 482 & clsurffiles(jtype,1:ifilessurf(jtype)), & 449 483 & nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 450 & rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav ) 451 452 453 CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea ) 484 & rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) 485 486 CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 454 487 455 488 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 456 CALL obs_rea_mdt( surfdataqc(jtype), nn_2dint ) 457 IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), nn_2dint, cn_altbiasfile ) 489 CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 490 IF ( ln_altbias ) & 491 & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 458 492 ENDIF 493 459 494 IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 460 !Read in bias field and correct SST. 461 IF ( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set,"// & 462 " but no bias"// & 463 " files to read in") 464 CALL obs_app_sstbias( surfdataqc(jtype), nn_2dint, & 465 jnumsstbias, cn_sstbias_files(1:jnumsstbias) ) 495 jnumsstbias = 0 496 DO jfile = 1, jpmaxnfiles 497 IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 498 & jnumsstbias = jnumsstbias + 1 499 END DO 500 IF ( jnumsstbias == 0 ) THEN 501 CALL ctl_stop("ln_sstbias set but no bias files to read in") 502 ENDIF 503 504 CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & 505 & jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) 506 466 507 ENDIF 467 508 END DO … … 512 553 USE ice , ONLY : at_i ! LIM3 Ice model variables 513 554 #endif 555 #if defined key_cice 556 USE sbc_oce, ONLY : fr_i ! ice fraction 557 #endif 558 514 559 IMPLICIT NONE 515 560 … … 528 573 & zprofmask2 ! Mask associated with zprofvar2 529 574 REAL(wp), POINTER, DIMENSION(:,:) :: & 530 & zsurfvar ! Model values equivalent to surface ob. 575 & zsurfvar, & ! Model values equivalent to surface ob. 576 & zsurfmask ! Mask associated with surface variable 531 577 REAL(wp), POINTER, DIMENSION(:,:) :: & 532 578 & zglam1, & ! Model longitudes for prof variable 1 … … 534 580 & zgphi1, & ! Model latitudes for prof variable 1 535 581 & zgphi2 ! Model latitudes for prof variable 2 536 #if ! defined key_lim3537 REAL(wp), POINTER, DIMENSION(:,:) :: at_i538 #endif539 LOGICAL :: llnightav ! Logical for calculating night-time average540 582 541 583 !Allocate local work arrays … … 545 587 CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 ) 546 588 CALL wrk_alloc( jpi, jpj, zsurfvar ) 589 CALL wrk_alloc( jpi, jpj, zsurfmask ) 547 590 CALL wrk_alloc( jpi, jpj, zglam1 ) 548 591 CALL wrk_alloc( jpi, jpj, zglam2 ) 549 592 CALL wrk_alloc( jpi, jpj, zgphi1 ) 550 593 CALL wrk_alloc( jpi, jpj, zgphi2 ) 551 #if ! defined key_lim3552 CALL wrk_alloc(jpi,jpj,at_i)553 #endif554 594 !----------------------------------------------------------------------- 555 595 … … 562 602 idaystp = NINT( rday / rdt ) 563 603 564 !-----------------------------------------------------------------------565 ! No LIM => at_i == 0.0_wp566 !-----------------------------------------------------------------------567 #if ! defined key_lim3568 at_i(:,:) = 0.0_wp569 #endif570 604 !----------------------------------------------------------------------- 571 605 ! Call the profile and surface observation operators … … 595 629 zgphi1(:,:) = gphiu(:,:) 596 630 zgphi2(:,:) = gphiv(:,:) 631 CASE DEFAULT 632 CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 597 633 END SELECT 598 634 599 IF( ln_zco .OR. ln_zps ) THEN 600 CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, & 601 & nit000, idaystp, & 602 & zprofvar1, zprofvar2, & 603 & gdept_1d, zprofmask1, zprofmask2, & 604 & zglam1, zglam2, zgphi1, zgphi2, & 605 & nn_1dint, nn_2dint, & 606 & kdailyavtypes = nn_profdavtypes ) 607 ELSE IF(TRIM(cobstypesprof(jtype)) == 'prof') THEN 608 !TG - THIS NEEDS MODIFICATION TO MATCH SIMPLIFICATION 609 CALL obs_pro_sco_opt( profdataqc(jtype), & 610 & kstp, jpi, jpj, jpk, nit000, idaystp, & 611 & zprofvar1, zprofvar2, & 612 & gdept_n(:,:,:), gdepw_n(:,:,:), & 613 & tmask, nn_1dint, nn_2dint, & 614 & kdailyavtypes = nn_profdavtypes ) 615 ELSE 616 CALL ctl_stop('DIA_OBS: Generalised vertical interpolation not'// & 617 'yet working for velocity data (turn off velocity observations') 618 ENDIF 635 CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, & 636 & nit000, idaystp, & 637 & zprofvar1, zprofvar2, & 638 & gdept_n(:,:,:), gdepw_n(:,:,:), & 639 & zprofmask1, zprofmask2, & 640 & zglam1, zglam2, zgphi1, zgphi2, & 641 & nn_1dint, nn_2dint, & 642 & kdailyavtypes = nn_profdavtypes ) 619 643 620 644 END DO … … 625 649 626 650 DO jtype = 1, nsurftypes 651 652 !Defaults which might be changed 653 zsurfmask(:,:) = tmask(:,:,1) 627 654 628 655 SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 629 656 CASE('sst') 630 657 zsurfvar(:,:) = tsn(:,:,1,jp_tem) 631 llnightav = ln_sstnight632 658 CASE('sla') 633 659 zsurfvar(:,:) = sshn(:,:) 634 llnightav = .FALSE.635 #if defined key_lim3 660 CASE('sss') 661 zsurfvar(:,:) = tsn(:,:,1,jp_sal) 636 662 CASE('sic') 637 663 IF ( kstp == 0 ) THEN … … 646 672 CYCLE 647 673 ELSE 648 zsurfvar(:,:) = at_i(:,:) 674 #if defined key_cice 675 zsurfvar(:,:) = fr_i(:,:) 676 #elif defined key_lim2 || defined key_lim3 677 zsurfvar(:,:) = 1._wp - frld(:,:) 678 #else 679 CALL ctl_stop( ' Trying to run sea-ice observation operator', & 680 & ' but no sea-ice model appears to have been defined' ) 681 #endif 649 682 ENDIF 650 683 651 llnightav = .FALSE.652 #endif653 684 END SELECT 654 685 655 686 CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, & 656 & nit000, idaystp, zsurfvar, tmask(:,:,1), & 657 & nn_2dint, llnightav ) 687 & nit000, idaystp, zsurfvar, zsurfmask, & 688 & n2dintsurf(jtype), llnightav(jtype), & 689 & zavglamscl(jtype), zavgphiscl(jtype), & 690 & lfpindegs(jtype) ) 658 691 659 692 END DO … … 666 699 CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 ) 667 700 CALL wrk_dealloc( jpi, jpj, zsurfvar ) 701 CALL wrk_dealloc( jpi, jpj, zsurfmask ) 668 702 CALL wrk_dealloc( jpi, jpj, zglam1 ) 669 703 CALL wrk_dealloc( jpi, jpj, zglam2 ) 670 704 CALL wrk_dealloc( jpi, jpj, zgphi1 ) 671 705 CALL wrk_dealloc( jpi, jpj, zgphi2 ) 672 #if ! defined key_lim3673 CALL wrk_dealloc(jpi,jpj,at_i)674 #endif675 706 676 707 END SUBROUTINE dia_obs … … 789 820 790 821 IF ( nsurftypes > 0 ) & 791 & DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf ) 822 & DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, & 823 & n2dintsurf, zavglamscl, zavgphiscl, lfpindegs, llnightav ) 792 824 793 825 END SUBROUTINE dia_obs_dealloc … … 938 970 END SUBROUTINE fin_date 939 971 972 SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, jtype, ctypein, & 973 & cfilestype, ifiles, cobstypes, cfiles ) 974 975 INTEGER, INTENT(IN) :: ntypes ! Total number of obs types 976 INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type 977 INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs 978 INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 979 & ifiles ! Out appended number of files for this type 980 981 CHARACTER(len=6), INTENT(IN) :: ctypein 982 CHARACTER(len=128), DIMENSION(jpmaxnfiles), INTENT(IN) :: & 983 & cfilestype ! In list of files for this obs type 984 CHARACTER(len=6), DIMENSION(ntypes), INTENT(INOUT) :: & 985 & cobstypes ! Out appended list of obs types 986 CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(INOUT) :: & 987 & cfiles ! Out appended list of files for all types 988 989 !Local variables 990 INTEGER :: jfile 991 992 cfiles(jtype,:) = cfilestype(:) 993 cobstypes(jtype) = ctypein 994 ifiles(jtype) = 0 995 DO jfile = 1, jpmaxnfiles 996 IF ( trim(cfiles(jtype,jfile)) /= '' ) & 997 ifiles(jtype) = ifiles(jtype) + 1 998 END DO 999 1000 IF ( ifiles(jtype) == 0 ) THEN 1001 CALL ctl_stop( 'Logical for observation type '//TRIM(ctypein)// & 1002 & ' set to true but no files available to read' ) 1003 ENDIF 1004 1005 IF(lwp) THEN 1006 WRITE(numout,*) ' '//cobstypes(jtype)//' input observation file names:' 1007 DO jfile = 1, ifiles(jtype) 1008 WRITE(numout,*) ' '//TRIM(cfiles(jtype,jfile)) 1009 END DO 1010 ENDIF 1011 1012 END SUBROUTINE obs_settypefiles 1013 1014 SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein, & 1015 & n2dint_default, n2dint_type, & 1016 & zavglamscl_type, zavgphiscl_type, & 1017 & lfp_indegs_type, lavnight_type, & 1018 & n2dint, zavglamscl, zavgphiscl, & 1019 & lfpindegs, lavnight ) 1020 1021 INTEGER, INTENT(IN) :: ntypes ! Total number of obs types 1022 INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs 1023 INTEGER, INTENT(IN) :: n2dint_default ! Default option for interpolation type 1024 INTEGER, INTENT(IN) :: n2dint_type ! Option for interpolation type 1025 REAL(wp), INTENT(IN) :: & 1026 & zavglamscl_type, & !E/W diameter of obs footprint for this type 1027 & zavgphiscl_type !N/S diameter of obs footprint for this type 1028 LOGICAL, INTENT(IN) :: lfp_indegs_type !T=> footprint in degrees, F=> in metres 1029 LOGICAL, INTENT(IN) :: lavnight_type !T=> obs represent night time average 1030 CHARACTER(len=6), INTENT(IN) :: ctypein 1031 1032 INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 1033 & n2dint 1034 REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & 1035 & zavglamscl, zavgphiscl 1036 LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & 1037 & lfpindegs, lavnight 1038 1039 lavnight(jtype) = lavnight_type 1040 1041 IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN 1042 n2dint(jtype) = n2dint_type 1043 ELSE 1044 n2dint(jtype) = n2dint_default 1045 ENDIF 1046 1047 ! For averaging observation footprints set options for size of footprint 1048 IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN 1049 IF ( zavglamscl_type > 0._wp ) THEN 1050 zavglamscl(jtype) = zavglamscl_type 1051 ELSE 1052 CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 1053 'scale (zavglamscl) for observation type '//TRIM(ctypein) ) 1054 ENDIF 1055 1056 IF ( zavgphiscl_type > 0._wp ) THEN 1057 zavgphiscl(jtype) = zavgphiscl_type 1058 ELSE 1059 CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 1060 'scale (zavgphiscl) for observation type '//TRIM(ctypein) ) 1061 ENDIF 1062 1063 lfpindegs(jtype) = lfp_indegs_type 1064 1065 ENDIF 1066 1067 ! Write out info 1068 IF(lwp) THEN 1069 IF ( n2dint(jtype) <= 4 ) THEN 1070 WRITE(numout,*) ' '//TRIM(ctypein)// & 1071 & ' model counterparts will be interpolated horizontally' 1072 ELSE IF ( n2dint(jtype) <= 6 ) THEN 1073 WRITE(numout,*) ' '//TRIM(ctypein)// & 1074 & ' model counterparts will be averaged horizontally' 1075 WRITE(numout,*) ' '//' with E/W scale: ',zavglamscl(jtype) 1076 WRITE(numout,*) ' '//' with N/S scale: ',zavgphiscl(jtype) 1077 IF ( lfpindegs(jtype) ) THEN 1078 WRITE(numout,*) ' '//' (in degrees)' 1079 ELSE 1080 WRITE(numout,*) ' '//' (in metres)' 1081 ENDIF 1082 ENDIF 1083 ENDIF 1084 1085 END SUBROUTINE obs_setinterpopts 1086 940 1087 END MODULE diaobs -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_mpp.F90
r6140 r9023 142 142 !! 143 143 144 iobsp =kobsp144 iobsp(:)=kobsp(:) 145 145 146 146 WHERE( iobsp(:) == -1 ) … … 148 148 END WHERE 149 149 150 iobsp =-1*iobsp150 iobsp(:)=-1*iobsp(:) 151 151 152 152 CALL obs_mpp_max_integer( iobsp, kno ) 153 153 154 kobsp =-1*iobsp154 kobsp(:)=-1*iobsp(:) 155 155 156 156 isum=0 -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
r7646 r9023 9 9 !! obs_prof_opt : Compute the model counterpart of profile data 10 10 !! obs_surf_opt : Compute the model counterpart of surface data 11 !! obs_pro_sco_opt: Compute the model counterpart of temperature and12 !! salinity observations from profiles in generalised13 !! vertical coordinates14 11 !!---------------------------------------------------------------------- 15 12 … … 22 19 & obs_int_h2d, & 23 20 & obs_int_h2d_init 21 USE obs_averg_h2d, ONLY : & ! Horizontal averaging to the obs footprint 22 & obs_avg_h2d, & 23 & obs_avg_h2d_init, & 24 & obs_max_fpsize 24 25 USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the obs pt 25 26 & obs_int_z1d, & 26 27 & obs_int_z1d_spl 27 USE obs_const, ONLY : &28 & obfillflt ! Fillvalue28 USE obs_const, ONLY : & ! Obs fill value 29 & obfillflt 29 30 USE dom_oce, ONLY : & 30 & glamt, glamu, glamv, & 31 & gphit, gphiu, gphiv, & 32 & gdept_n, gdept_0 33 USE lib_mpp, ONLY : & 31 & glamt, glamf, & 32 & gphit, gphif 33 USE lib_mpp, ONLY : & ! Warning and stopping routines 34 34 & ctl_warn, ctl_stop 35 USE sbcdcy, ONLY : & ! For calculation of where it is night-time 36 & sbc_dcy, nday_qsr 35 37 USE obs_grid, ONLY : & 36 38 & obs_level_search 37 USE sbcdcy, ONLY : & ! For calculation of where it is night-time38 & sbc_dcy, nday_qsr39 39 40 40 IMPLICIT NONE … … 44 44 45 45 PUBLIC obs_prof_opt, & ! Compute the model counterpart of profile obs 46 & obs_pro_sco_opt, & ! Compute the model counterpart of profile observations47 46 & obs_surf_opt ! Compute the model counterpart of surface obs 48 47 … … 58 57 CONTAINS 59 58 59 60 60 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 61 61 & kit000, kdaystp, & 62 & pvar1, pvar2, pgdept, pmask1, pmask2, & 62 & pvar1, pvar2, pgdept, pgdepw, & 63 & pmask1, pmask2, & 63 64 & plam1, plam2, pphi1, pphi2, & 64 65 & k1dint, k2dint, kdailyavtypes ) … … 111 112 !! ! 07-03 (K. Mogensen) General handling of profiles 112 113 !! ! 15-02 (M. Martin) Combined routine for all profile types 114 !! ! 17-02 (M. Martin) Include generalised vertical coordinate changes 113 115 !!----------------------------------------------------------------------- 114 116 … … 140 142 & pphi1, & ! Model latitudes for variable 1 141 143 & pphi2 ! Model latitudes for variable 2 142 REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 143 & pgdept ! Model array of depth levels 144 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 145 & pgdept, & ! Model array of depth T levels 146 & pgdepw ! Model array of depth W levels 144 147 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 145 148 & kdailyavtypes ! Types for daily averages … … 156 159 INTEGER :: iend 157 160 INTEGER :: iobs 161 INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes 162 INTEGER :: inum_obs 158 163 INTEGER, DIMENSION(imaxavtypes) :: & 159 164 & idailyavtypes … … 163 168 & igrdj1, & 164 169 & igrdj2 170 INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 171 165 172 REAL(KIND=wp) :: zlam 166 173 REAL(KIND=wp) :: zphi … … 171 178 & zobsk, & 172 179 & zobs2k 173 REAL(KIND=wp), DIMENSION(2,2, kpk) :: &180 REAL(KIND=wp), DIMENSION(2,2,1) :: & 174 181 & zweig1, & 175 & zweig2 182 & zweig2, & 183 & zweig 176 184 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 177 185 & zmask1, & 178 186 & zmask2, & 179 & zint1, & 180 & zint2, & 181 & zinm1, & 182 & zinm2 187 & zint1, & 188 & zint2, & 189 & zinm1, & 190 & zinm2, & 191 & zgdept, & 192 & zgdepw 183 193 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 184 194 & zglam1, & … … 186 196 & zgphi1, & 187 197 & zgphi2 198 REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2 199 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 200 188 201 LOGICAL :: ld_dailyav 189 202 … … 266 279 & zmask1(2,2,kpk,ipro), & 267 280 & zmask2(2,2,kpk,ipro), & 268 & zint1(2,2,kpk,ipro), & 269 & zint2(2,2,kpk,ipro) & 281 & zint1(2,2,kpk,ipro), & 282 & zint2(2,2,kpk,ipro), & 283 & zgdept(2,2,kpk,ipro), & 284 & zgdepw(2,2,kpk,ipro) & 270 285 & ) 271 286 … … 290 305 END DO 291 306 307 ! Initialise depth arrays 308 zgdept(:,:,:,:) = 0.0 309 zgdepw(:,:,:,:) = 0.0 310 292 311 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 ) 293 312 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 ) … … 300 319 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2, zint2 ) 301 320 321 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept ) 322 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw ) 323 302 324 ! At the end of the day also get interpolated means 303 325 IF ( ld_dailyav .AND. idayend == 0 ) THEN … … 314 336 315 337 ENDIF 338 339 ! Return if no observations to process 340 ! Has to be done after comm commands to ensure processors 341 ! stay in sync 342 IF ( ipro == 0 ) RETURN 316 343 317 344 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro … … 339 366 zphi = prodatqc%rphi(jobs) 340 367 341 ! Horizontal weights and vertical mask342 368 ! Horizontal weights 369 ! Masked values are calculated later. 343 370 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 344 371 345 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &372 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 346 373 & zglam1(:,:,iobs), zgphi1(:,:,iobs), & 347 & zmask1(:,:, :,iobs), zweig1, zobsmask1 )374 & zmask1(:,:,1,iobs), zweig1, zmsk_1 ) 348 375 349 376 ENDIF … … 351 378 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 352 379 353 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &380 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 354 381 & zglam2(:,:,iobs), zgphi2(:,:,iobs), & 355 & zmask2(:,:, :,iobs), zweig2, zobsmask2 )382 & zmask2(:,:,1,iobs), zweig2, zmsk_2 ) 356 383 357 384 ENDIF … … 365 392 IF ( idayend == 0 ) THEN 366 393 ! Daily averaged data 367 CALL obs_int_h2d( kpk, kpk, & 368 & zweig1, zinm1(:,:,:,iobs), zobsk ) 369 370 ENDIF 371 372 ELSE 373 374 ! Point data 375 CALL obs_int_h2d( kpk, kpk, & 376 & zweig1, zint1(:,:,:,iobs), zobsk ) 377 378 ENDIF 379 380 !------------------------------------------------------------- 381 ! Compute vertical second-derivative of the interpolating 382 ! polynomial at obs points 383 !------------------------------------------------------------- 384 385 IF ( k1dint == 1 ) THEN 386 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 387 & pgdept, zobsmask1 ) 388 ENDIF 389 390 !----------------------------------------------------------------- 391 ! Vertical interpolation to the observation point 392 !----------------------------------------------------------------- 393 ista = prodatqc%npvsta(jobs,1) 394 iend = prodatqc%npvend(jobs,1) 395 CALL obs_int_z1d( kpk, & 396 & prodatqc%var(1)%mvk(ista:iend), & 397 & k1dint, iend - ista + 1, & 398 & prodatqc%var(1)%vdep(ista:iend), & 399 & zobsk, zobs2k, & 400 & prodatqc%var(1)%vmod(ista:iend), & 401 & pgdept, zobsmask1 ) 402 403 ENDIF 404 394 395 ! vertically interpolate all 4 corners 396 ista = prodatqc%npvsta(jobs,1) 397 iend = prodatqc%npvend(jobs,1) 398 inum_obs = iend - ista + 1 399 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 400 401 DO iin=1,2 402 DO ijn=1,2 403 404 IF ( k1dint == 1 ) THEN 405 CALL obs_int_z1d_spl( kpk, & 406 & zinm1(iin,ijn,:,iobs), & 407 & zobs2k, zgdept(iin,ijn,:,iobs), & 408 & zmask1(iin,ijn,:,iobs)) 409 ENDIF 410 411 CALL obs_level_search(kpk, & 412 & zgdept(iin,ijn,:,iobs), & 413 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 414 & iv_indic) 415 416 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 417 & prodatqc%var(1)%vdep(ista:iend), & 418 & zinm1(iin,ijn,:,iobs), & 419 & zobs2k, interp_corner(iin,ijn,:), & 420 & zgdept(iin,ijn,:,iobs), & 421 & zmask1(iin,ijn,:,iobs)) 422 423 ENDDO 424 ENDDO 425 426 ENDIF !idayend 427 428 ELSE 429 430 ! Point data 431 432 ! vertically interpolate all 4 corners 433 ista = prodatqc%npvsta(jobs,1) 434 iend = prodatqc%npvend(jobs,1) 435 inum_obs = iend - ista + 1 436 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 437 DO iin=1,2 438 DO ijn=1,2 439 440 IF ( k1dint == 1 ) THEN 441 CALL obs_int_z1d_spl( kpk, & 442 & zint1(iin,ijn,:,iobs),& 443 & zobs2k, zgdept(iin,ijn,:,iobs), & 444 & zmask1(iin,ijn,:,iobs)) 445 446 ENDIF 447 448 CALL obs_level_search(kpk, & 449 & zgdept(iin,ijn,:,iobs),& 450 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 451 & iv_indic) 452 453 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 454 & prodatqc%var(1)%vdep(ista:iend), & 455 & zint1(iin,ijn,:,iobs), & 456 & zobs2k,interp_corner(iin,ijn,:), & 457 & zgdept(iin,ijn,:,iobs), & 458 & zmask1(iin,ijn,:,iobs) ) 459 460 ENDDO 461 ENDDO 462 463 ENDIF 464 465 !------------------------------------------------------------- 466 ! Compute the horizontal interpolation for every profile level 467 !------------------------------------------------------------- 468 469 DO ikn=1,inum_obs 470 iend=ista+ikn-1 471 472 zweig(:,:,1) = 0._wp 473 474 ! This code forces the horizontal weights to be 475 ! zero IF the observation is below the bottom of the 476 ! corners of the interpolation nodes, Or if it is in 477 ! the mask. This is important for observations near 478 ! steep bathymetry 479 DO iin=1,2 480 DO ijn=1,2 481 482 depth_loop1: DO ik=kpk,2,-1 483 IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN 484 485 zweig(iin,ijn,1) = & 486 & zweig1(iin,ijn,1) * & 487 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 488 & - prodatqc%var(1)%vdep(iend)),0._wp) 489 490 EXIT depth_loop1 491 492 ENDIF 493 494 ENDDO depth_loop1 495 496 ENDDO 497 ENDDO 498 499 CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 500 & prodatqc%var(1)%vmod(iend:iend) ) 501 502 ! Set QC flag for any observations found below the bottom 503 ! needed as the check here is more strict than that in obs_prep 504 IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4 505 506 ENDDO 507 508 DEALLOCATE(interp_corner,iv_indic) 509 510 ENDIF 511 512 ! For the second variable 405 513 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 406 514 … … 410 518 411 519 IF ( idayend == 0 ) THEN 412 413 520 ! Daily averaged data 414 CALL obs_int_h2d( kpk, kpk, & 415 & zweig2, zinm2(:,:,:,iobs), zobsk ) 416 417 ENDIF 418 419 ELSE 420 421 ! Point data 422 CALL obs_int_h2d( kpk, kpk, & 423 & zweig2, zint2(:,:,:,iobs), zobsk ) 424 425 ENDIF 426 427 428 !------------------------------------------------------------- 429 ! Compute vertical second-derivative of the interpolating 430 ! polynomial at obs points 431 !------------------------------------------------------------- 432 433 IF ( k1dint == 1 ) THEN 434 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 435 & pgdept, zobsmask2 ) 436 ENDIF 437 438 !---------------------------------------------------------------- 439 ! Vertical interpolation to the observation point 440 !---------------------------------------------------------------- 441 ista = prodatqc%npvsta(jobs,2) 442 iend = prodatqc%npvend(jobs,2) 443 CALL obs_int_z1d( kpk, & 444 & prodatqc%var(2)%mvk(ista:iend),& 445 & k1dint, iend - ista + 1, & 446 & prodatqc%var(2)%vdep(ista:iend),& 447 & zobsk, zobs2k, & 448 & prodatqc%var(2)%vmod(ista:iend),& 449 & pgdept, zobsmask2 ) 450 451 ENDIF 452 453 END DO 521 522 ! vertically interpolate all 4 corners 523 ista = prodatqc%npvsta(jobs,2) 524 iend = prodatqc%npvend(jobs,2) 525 inum_obs = iend - ista + 1 526 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 527 528 DO iin=1,2 529 DO ijn=1,2 530 531 IF ( k1dint == 1 ) THEN 532 CALL obs_int_z1d_spl( kpk, & 533 & zinm2(iin,ijn,:,iobs), & 534 & zobs2k, zgdept(iin,ijn,:,iobs), & 535 & zmask2(iin,ijn,:,iobs)) 536 ENDIF 537 538 CALL obs_level_search(kpk, & 539 & zgdept(iin,ijn,:,iobs), & 540 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 541 & iv_indic) 542 543 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 544 & prodatqc%var(2)%vdep(ista:iend), & 545 & zinm2(iin,ijn,:,iobs), & 546 & zobs2k, interp_corner(iin,ijn,:), & 547 & zgdept(iin,ijn,:,iobs), & 548 & zmask2(iin,ijn,:,iobs)) 549 550 ENDDO 551 ENDDO 552 553 ENDIF !idayend 554 555 ELSE 556 557 ! Point data 558 559 ! vertically interpolate all 4 corners 560 ista = prodatqc%npvsta(jobs,2) 561 iend = prodatqc%npvend(jobs,2) 562 inum_obs = iend - ista + 1 563 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 564 DO iin=1,2 565 DO ijn=1,2 566 567 IF ( k1dint == 1 ) THEN 568 CALL obs_int_z1d_spl( kpk, & 569 & zint2(iin,ijn,:,iobs),& 570 & zobs2k, zgdept(iin,ijn,:,iobs), & 571 & zmask2(iin,ijn,:,iobs)) 572 573 ENDIF 574 575 CALL obs_level_search(kpk, & 576 & zgdept(iin,ijn,:,iobs),& 577 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 578 & iv_indic) 579 580 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 581 & prodatqc%var(2)%vdep(ista:iend), & 582 & zint2(iin,ijn,:,iobs), & 583 & zobs2k,interp_corner(iin,ijn,:), & 584 & zgdept(iin,ijn,:,iobs), & 585 & zmask2(iin,ijn,:,iobs) ) 586 587 ENDDO 588 ENDDO 589 590 ENDIF 591 592 !------------------------------------------------------------- 593 ! Compute the horizontal interpolation for every profile level 594 !------------------------------------------------------------- 595 596 DO ikn=1,inum_obs 597 iend=ista+ikn-1 598 599 zweig(:,:,1) = 0._wp 600 601 ! This code forces the horizontal weights to be 602 ! zero IF the observation is below the bottom of the 603 ! corners of the interpolation nodes, Or if it is in 604 ! the mask. This is important for observations near 605 ! steep bathymetry 606 DO iin=1,2 607 DO ijn=1,2 608 609 depth_loop2: DO ik=kpk,2,-1 610 IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN 611 612 zweig(iin,ijn,1) = & 613 & zweig2(iin,ijn,1) * & 614 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 615 & - prodatqc%var(2)%vdep(iend)),0._wp) 616 617 EXIT depth_loop2 618 619 ENDIF 620 621 ENDDO depth_loop2 622 623 ENDDO 624 ENDDO 625 626 CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 627 & prodatqc%var(2)%vmod(iend:iend) ) 628 629 ! Set QC flag for any observations found below the bottom 630 ! needed as the check here is more strict than that in obs_prep 631 IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 632 633 ENDDO 634 635 DEALLOCATE(interp_corner,iv_indic) 636 637 ENDIF 638 639 ENDDO 454 640 455 641 ! Deallocate the data for interpolation … … 466 652 & zmask2, & 467 653 & zint1, & 468 & zint2 & 654 & zint2, & 655 & zgdept, & 656 & zgdepw & 469 657 & ) 470 658 … … 481 669 END SUBROUTINE obs_prof_opt 482 670 483 SUBROUTINE obs_pro_sco_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 484 & ptn, psn, pgdept, pgdepw, ptmask, k1dint, k2dint, & 485 & kdailyavtypes ) 486 !!----------------------------------------------------------------------- 487 !! 488 !! *** ROUTINE obs_pro_opt *** 489 !! 490 !! ** Purpose : Compute the model counterpart of profiles 491 !! data by interpolating from the model grid to the 492 !! observation point. Generalised vertical coordinate version 493 !! 494 !! ** Method : Linearly interpolate to each observation point using 495 !! the model values at the corners of the surrounding grid box. 496 !! 497 !! First, model values on the model grid are interpolated vertically to the 498 !! Depths of the profile observations. Two vertical interpolation schemes are 499 !! available: 500 !! - linear (k1dint = 0) 501 !! - Cubic spline (k1dint = 1) 502 !! 503 !! 504 !! Secondly the interpolated values are interpolated horizontally to the 505 !! obs (lon, lat) point. 506 !! Several horizontal interpolation schemes are available: 507 !! - distance-weighted (great circle) (k2dint = 0) 508 !! - distance-weighted (small angle) (k2dint = 1) 509 !! - bilinear (geographical grid) (k2dint = 2) 510 !! - bilinear (quadrilateral grid) (k2dint = 3) 511 !! - polynomial (quadrilateral grid) (k2dint = 4) 512 !! 513 !! For the cubic spline the 2nd derivative of the interpolating 514 !! polynomial is computed before entering the vertical interpolation 515 !! routine. 516 !! 517 !! For ENACT moored buoy data (e.g., TAO), the model equivalent is 518 !! a daily mean model temperature field. So, we first compute 519 !! the mean, then interpolate only at the end of the day. 520 !! 521 !! This is the procedure to be used with generalised vertical model 522 !! coordinates (ie s-coordinates. It is ~4x slower than the equivalent 523 !! horizontal then vertical interpolation algorithm, but can deal with situations 524 !! where the model levels are not flat. 525 !! ONLY PERFORMED if ln_sco=.TRUE. 526 !! 527 !! Note: the in situ temperature observations must be converted 528 !! to potential temperature (the model variable) prior to 529 !! assimilation. 530 !!?????????????????????????????????????????????????????????????? 531 !! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR??? 532 !!?????????????????????????????????????????????????????????????? 533 !! 534 !! ** Action : 535 !! 536 !! History : 537 !! ! 2014-08 (J. While) Adapted from obs_pro_opt to handel generalised 538 !! vertical coordinates 539 !!----------------------------------------------------------------------- 540 541 !! * Modules used 542 USE obs_profiles_def ! Definition of storage space for profile obs. 543 544 IMPLICIT NONE 545 546 !! * Arguments 547 TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening 548 INTEGER, INTENT(IN) :: kt ! Time step 549 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 550 INTEGER, INTENT(IN) :: kpj 551 INTEGER, INTENT(IN) :: kpk 552 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 553 ! (kit000-1 = restart time) 554 INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header) 555 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 556 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 557 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 558 & ptn, & ! Model temperature field 559 & psn, & ! Model salinity field 560 & ptmask ! Land-sea mask 561 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 562 & pgdept, & ! Model array of depth T levels 563 & pgdepw ! Model array of depth W levels 564 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 565 & kdailyavtypes ! Types for daily averages 566 567 !! * Local declarations 568 INTEGER :: ji 569 INTEGER :: jj 570 INTEGER :: jk 571 INTEGER :: iico, ijco 572 INTEGER :: jobs 573 INTEGER :: inrc 574 INTEGER :: ipro 575 INTEGER :: idayend 576 INTEGER :: ista 577 INTEGER :: iend 578 INTEGER :: iobs 579 INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes 580 INTEGER, DIMENSION(imaxavtypes) :: & 581 & idailyavtypes 582 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 583 & igrdi, & 584 & igrdj 585 INTEGER :: & 586 & inum_obs 587 INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 588 REAL(KIND=wp) :: zlam 589 REAL(KIND=wp) :: zphi 590 REAL(KIND=wp) :: zdaystp 591 REAL(KIND=wp), DIMENSION(kpk) :: & 592 & zobsmask, & 593 & zobsk, & 594 & zobs2k 595 REAL(KIND=wp), DIMENSION(2,2,1) :: & 596 & zweig, & 597 & l_zweig 598 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 599 & zmask, & 600 & zintt, & 601 & zints, & 602 & zinmt, & 603 & zgdept,& 604 & zgdepw,& 605 & zinms 606 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 607 & zglam, & 608 & zgphi 609 REAL(KIND=wp), DIMENSION(1) :: zmsk_1 610 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 611 612 !------------------------------------------------------------------------ 613 ! Local initialization 614 !------------------------------------------------------------------------ 615 ! ... Record and data counters 616 inrc = kt - kit000 + 2 617 ipro = prodatqc%npstp(inrc) 618 619 ! Daily average types 620 IF ( PRESENT(kdailyavtypes) ) THEN 621 idailyavtypes(:) = kdailyavtypes(:) 622 ELSE 623 idailyavtypes(:) = -1 624 ENDIF 625 626 ! Initialize daily mean for first time-step 627 idayend = MOD( kt - kit000 + 1, kdaystp ) 628 629 ! Added kt == 0 test to catch restart case 630 IF ( idayend == 1 .OR. kt == 0) THEN 631 632 IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 633 DO jk = 1, jpk 634 DO jj = 1, jpj 635 DO ji = 1, jpi 636 prodatqc%vdmean(ji,jj,jk,1) = 0.0 637 prodatqc%vdmean(ji,jj,jk,2) = 0.0 638 END DO 639 END DO 640 END DO 641 642 ENDIF 643 644 DO jk = 1, jpk 645 DO jj = 1, jpj 646 DO ji = 1, jpi 647 ! Increment the temperature field for computing daily mean 648 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 649 & + ptn(ji,jj,jk) 650 ! Increment the salinity field for computing daily mean 651 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 652 & + psn(ji,jj,jk) 653 END DO 654 END DO 655 END DO 656 657 ! Compute the daily mean at the end of day 658 zdaystp = 1.0 / REAL( kdaystp ) 659 IF ( idayend == 0 ) THEN 660 DO jk = 1, jpk 661 DO jj = 1, jpj 662 DO ji = 1, jpi 663 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 664 & * zdaystp 665 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 666 & * zdaystp 667 END DO 668 END DO 669 END DO 670 ENDIF 671 672 ! Get the data for interpolation 673 ALLOCATE( & 674 & igrdi(2,2,ipro), & 675 & igrdj(2,2,ipro), & 676 & zglam(2,2,ipro), & 677 & zgphi(2,2,ipro), & 678 & zmask(2,2,kpk,ipro), & 679 & zintt(2,2,kpk,ipro), & 680 & zints(2,2,kpk,ipro), & 681 & zgdept(2,2,kpk,ipro), & 682 & zgdepw(2,2,kpk,ipro) & 683 & ) 684 685 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 686 iobs = jobs - prodatqc%nprofup 687 igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1 688 igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1 689 igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1 690 igrdj(1,2,iobs) = prodatqc%mj(jobs,1) 691 igrdi(2,1,iobs) = prodatqc%mi(jobs,1) 692 igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1 693 igrdi(2,2,iobs) = prodatqc%mi(jobs,1) 694 igrdj(2,2,iobs) = prodatqc%mj(jobs,1) 695 END DO 696 697 ! Initialise depth arrays 698 zgdept = 0.0 699 zgdepw = 0.0 700 701 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, glamt, zglam ) 702 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, gphit, zgphi ) 703 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, ptmask,zmask ) 704 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, ptn, zintt ) 705 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, psn, zints ) 706 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept(:,:,:), & 707 & zgdept ) 708 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw(:,:,:), & 709 & zgdepw ) 710 711 ! At the end of the day also get interpolated means 712 IF ( idayend == 0 ) THEN 713 714 ALLOCATE( & 715 & zinmt(2,2,kpk,ipro), & 716 & zinms(2,2,kpk,ipro) & 717 & ) 718 719 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 720 & prodatqc%vdmean(:,:,:,1), zinmt ) 721 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 722 & prodatqc%vdmean(:,:,:,2), zinms ) 723 724 ENDIF 725 726 ! Return if no observations to process 727 ! Has to be done after comm commands to ensure processors 728 ! stay in sync 729 IF ( ipro == 0 ) RETURN 730 731 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 732 733 iobs = jobs - prodatqc%nprofup 734 735 IF ( kt /= prodatqc%mstp(jobs) ) THEN 736 737 IF(lwp) THEN 738 WRITE(numout,*) 739 WRITE(numout,*) ' E R R O R : Observation', & 740 & ' time step is not consistent with the', & 741 & ' model time step' 742 WRITE(numout,*) ' =========' 743 WRITE(numout,*) 744 WRITE(numout,*) ' Record = ', jobs, & 745 & ' kt = ', kt, & 746 & ' mstp = ', prodatqc%mstp(jobs), & 747 & ' ntyp = ', prodatqc%ntyp(jobs) 748 ENDIF 749 CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 750 ENDIF 751 752 zlam = prodatqc%rlam(jobs) 753 zphi = prodatqc%rphi(jobs) 754 755 ! Horizontal weights 756 ! Only calculated once, for both T and S. 757 ! Masked values are calculated later. 758 759 IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 760 & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 761 762 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 763 & zglam(:,:,iobs), zgphi(:,:,iobs), & 764 & zmask(:,:,1,iobs), zweig, zmsk_1 ) 765 766 ENDIF 767 768 ! IF zmsk_1 = 0; then ob is on land 769 IF (zmsk_1(1) < 0.1) THEN 770 WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask' 771 772 ELSE 773 774 ! Temperature 775 776 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 777 778 zobsk(:) = obfillflt 779 780 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 781 782 IF ( idayend == 0 ) THEN 783 784 ! Daily averaged moored buoy (MRB) data 785 786 ! vertically interpolate all 4 corners 787 ista = prodatqc%npvsta(jobs,1) 788 iend = prodatqc%npvend(jobs,1) 789 inum_obs = iend - ista + 1 790 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 791 792 DO iin=1,2 793 DO ijn=1,2 794 795 796 797 IF ( k1dint == 1 ) THEN 798 CALL obs_int_z1d_spl( kpk, & 799 & zinmt(iin,ijn,:,iobs), & 800 & zobs2k, zgdept(iin,ijn,:,iobs), & 801 & zmask(iin,ijn,:,iobs)) 802 ENDIF 803 804 CALL obs_level_search(kpk, & 805 & zgdept(iin,ijn,:,iobs), & 806 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 807 & iv_indic) 808 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 809 & prodatqc%var(1)%vdep(ista:iend), & 810 & zinmt(iin,ijn,:,iobs), & 811 & zobs2k, interp_corner(iin,ijn,:), & 812 & zgdept(iin,ijn,:,iobs), & 813 & zmask(iin,ijn,:,iobs)) 814 815 ENDDO 816 ENDDO 817 818 819 ELSE 820 821 CALL ctl_stop( ' A nonzero' // & 822 & ' number of profile T BUOY data should' // & 823 & ' only occur at the end of a given day' ) 824 825 ENDIF 826 827 ELSE 828 829 ! Point data 830 831 ! vertically interpolate all 4 corners 832 ista = prodatqc%npvsta(jobs,1) 833 iend = prodatqc%npvend(jobs,1) 834 inum_obs = iend - ista + 1 835 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 836 DO iin=1,2 837 DO ijn=1,2 838 839 840 IF ( k1dint == 1 ) THEN 841 CALL obs_int_z1d_spl( kpk, & 842 & zintt(iin,ijn,:,iobs),& 843 & zobs2k, zgdept(iin,ijn,:,iobs), & 844 & zmask(iin,ijn,:,iobs)) 845 846 ENDIF 847 848 CALL obs_level_search(kpk, & 849 & zgdept(iin,ijn,:,iobs),& 850 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 851 & iv_indic) 852 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 853 & prodatqc%var(1)%vdep(ista:iend), & 854 & zintt(iin,ijn,:,iobs), & 855 & zobs2k,interp_corner(iin,ijn,:), & 856 & zgdept(iin,ijn,:,iobs), & 857 & zmask(iin,ijn,:,iobs) ) 858 859 ENDDO 860 ENDDO 861 862 ENDIF 863 864 !------------------------------------------------------------- 865 ! Compute the horizontal interpolation for every profile level 866 !------------------------------------------------------------- 867 868 DO ikn=1,inum_obs 869 iend=ista+ikn-1 870 871 l_zweig(:,:,1) = 0._wp 872 873 ! This code forces the horizontal weights to be 874 ! zero IF the observation is below the bottom of the 875 ! corners of the interpolation nodes, Or if it is in 876 ! the mask. This is important for observations are near 877 ! steep bathymetry 878 DO iin=1,2 879 DO ijn=1,2 880 881 depth_loop1: DO ik=kpk,2,-1 882 IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN 883 884 l_zweig(iin,ijn,1) = & 885 & zweig(iin,ijn,1) * & 886 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 887 & - prodatqc%var(1)%vdep(iend)),0._wp) 888 889 EXIT depth_loop1 890 ENDIF 891 ENDDO depth_loop1 892 893 ENDDO 894 ENDDO 895 896 CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 897 & prodatqc%var(1)%vmod(iend:iend) ) 898 899 ENDDO 900 901 902 DEALLOCATE(interp_corner,iv_indic) 903 904 ENDIF 905 906 907 ! Salinity 908 909 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 910 911 zobsk(:) = obfillflt 912 913 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 914 915 IF ( idayend == 0 ) THEN 916 917 ! Daily averaged moored buoy (MRB) data 918 919 ! vertically interpolate all 4 corners 920 ista = prodatqc%npvsta(iobs,2) 921 iend = prodatqc%npvend(iobs,2) 922 inum_obs = iend - ista + 1 923 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 924 925 DO iin=1,2 926 DO ijn=1,2 927 928 929 930 IF ( k1dint == 1 ) THEN 931 CALL obs_int_z1d_spl( kpk, & 932 & zinms(iin,ijn,:,iobs), & 933 & zobs2k, zgdept(iin,ijn,:,iobs), & 934 & zmask(iin,ijn,:,iobs)) 935 ENDIF 936 937 CALL obs_level_search(kpk, & 938 & zgdept(iin,ijn,:,iobs), & 939 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 940 & iv_indic) 941 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 942 & prodatqc%var(2)%vdep(ista:iend), & 943 & zinms(iin,ijn,:,iobs), & 944 & zobs2k, interp_corner(iin,ijn,:), & 945 & zgdept(iin,ijn,:,iobs), & 946 & zmask(iin,ijn,:,iobs)) 947 948 ENDDO 949 ENDDO 950 951 952 ELSE 953 954 CALL ctl_stop( ' A nonzero' // & 955 & ' number of profile T BUOY data should' // & 956 & ' only occur at the end of a given day' ) 957 958 ENDIF 959 960 ELSE 961 962 ! Point data 963 964 ! vertically interpolate all 4 corners 965 ista = prodatqc%npvsta(jobs,2) 966 iend = prodatqc%npvend(jobs,2) 967 inum_obs = iend - ista + 1 968 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 969 970 DO iin=1,2 971 DO ijn=1,2 972 973 974 IF ( k1dint == 1 ) THEN 975 CALL obs_int_z1d_spl( kpk, & 976 & zints(iin,ijn,:,iobs),& 977 & zobs2k, zgdept(iin,ijn,:,iobs), & 978 & zmask(iin,ijn,:,iobs)) 979 980 ENDIF 981 982 CALL obs_level_search(kpk, & 983 & zgdept(iin,ijn,:,iobs),& 984 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 985 & iv_indic) 986 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 987 & prodatqc%var(2)%vdep(ista:iend), & 988 & zints(iin,ijn,:,iobs), & 989 & zobs2k,interp_corner(iin,ijn,:), & 990 & zgdept(iin,ijn,:,iobs), & 991 & zmask(iin,ijn,:,iobs) ) 992 993 ENDDO 994 ENDDO 995 996 ENDIF 997 998 !------------------------------------------------------------- 999 ! Compute the horizontal interpolation for every profile level 1000 !------------------------------------------------------------- 1001 1002 DO ikn=1,inum_obs 1003 iend=ista+ikn-1 1004 1005 l_zweig(:,:,1) = 0._wp 1006 1007 ! This code forces the horizontal weights to be 1008 ! zero IF the observation is below the bottom of the 1009 ! corners of the interpolation nodes, Or if it is in 1010 ! the mask. This is important for observations are near 1011 ! steep bathymetry 1012 DO iin=1,2 1013 DO ijn=1,2 1014 1015 depth_loop2: DO ik=kpk,2,-1 1016 IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN 1017 1018 l_zweig(iin,ijn,1) = & 1019 & zweig(iin,ijn,1) * & 1020 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 1021 & - prodatqc%var(2)%vdep(iend)),0._wp) 1022 1023 EXIT depth_loop2 1024 ENDIF 1025 ENDDO depth_loop2 1026 1027 ENDDO 1028 ENDDO 1029 1030 CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 1031 & prodatqc%var(2)%vmod(iend:iend) ) 1032 1033 ENDDO 1034 1035 1036 DEALLOCATE(interp_corner,iv_indic) 1037 1038 ENDIF 1039 1040 ENDIF 1041 1042 END DO 1043 1044 ! Deallocate the data for interpolation 1045 DEALLOCATE( & 1046 & igrdi, & 1047 & igrdj, & 1048 & zglam, & 1049 & zgphi, & 1050 & zmask, & 1051 & zintt, & 1052 & zints, & 1053 & zgdept,& 1054 & zgdepw & 1055 & ) 1056 ! At the end of the day also get interpolated means 1057 IF ( idayend == 0 ) THEN 1058 DEALLOCATE( & 1059 & zinmt, & 1060 & zinms & 1061 & ) 1062 ENDIF 1063 1064 prodatqc%nprofup = prodatqc%nprofup + ipro 1065 1066 END SUBROUTINE obs_pro_sco_opt 1067 1068 SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, & 1069 & kit000, kdaystp, psurf, psurfmask, & 1070 & k2dint, ldnightav ) 671 SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, & 672 & kit000, kdaystp, psurf, psurfmask, & 673 & k2dint, ldnightav, plamscl, pphiscl, & 674 & lindegrees ) 1071 675 1072 676 !!----------------------------------------------------------------------- … … 1090 694 !! - polynomial (quadrilateral grid) (k2dint = 4) 1091 695 !! 696 !! Two horizontal averaging schemes are also available: 697 !! - weighted radial footprint (k2dint = 5) 698 !! - weighted rectangular footprint (k2dint = 6) 699 !! 1092 700 !! 1093 701 !! ** Action : … … 1096 704 !! ! 07-03 (A. Weaver) 1097 705 !! ! 15-02 (M. Martin) Combined routine for surface types 706 !! ! 17-03 (M. Martin) Added horizontal averaging options 1098 707 !!----------------------------------------------------------------------- 1099 708 … … 1117 726 & psurfmask ! Land-sea mask 1118 727 LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 728 REAL(KIND=wp), INTENT(IN) :: & 729 & plamscl, & ! Diameter in metres of obs footprint in E/W, N/S directions 730 & pphiscl ! This is the full width (rather than half-width) 731 LOGICAL, INTENT(IN) :: & 732 & lindegrees ! T=> plamscl and pphiscl are specified in degrees, F=> in metres 1119 733 1120 734 !! * Local declarations … … 1125 739 INTEGER :: isurf 1126 740 INTEGER :: iobs 741 INTEGER :: imaxifp, imaxjfp 742 INTEGER :: imodi, imodj 1127 743 INTEGER :: idayend 1128 744 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 1129 & igrdi, & 1130 & igrdj 745 & igrdi, & 746 & igrdj, & 747 & igrdip1, & 748 & igrdjp1 1131 749 INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 1132 750 & icount_night, & … … 1136 754 REAL(wp), DIMENSION(1) :: zext, zobsmask 1137 755 REAL(wp) :: zdaystp 1138 REAL(wp), DIMENSION(2,2,1) :: &1139 & zweig1140 756 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 757 & zweig, & 1141 758 & zmask, & 1142 759 & zsurf, & 1143 760 & zsurfm, & 761 & zsurftmp, & 1144 762 & zglam, & 1145 & zgphi 763 & zgphi, & 764 & zglamf, & 765 & zgphif 766 1146 767 REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 1147 768 & zintmp, & … … 1155 776 inrc = kt - kit000 + 2 1156 777 isurf = surfdataqc%nsstp(inrc) 778 779 ! Work out the maximum footprint size for the 780 ! interpolation/averaging in model grid-points - has to be even. 781 782 CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 783 1157 784 1158 785 IF ( ldnightav ) THEN … … 1221 848 1222 849 ALLOCATE( & 1223 & igrdi(2,2,isurf), & 1224 & igrdj(2,2,isurf), & 1225 & zglam(2,2,isurf), & 1226 & zgphi(2,2,isurf), & 1227 & zmask(2,2,isurf), & 1228 & zsurf(2,2,isurf) & 850 & zweig(imaxifp,imaxjfp,1), & 851 & igrdi(imaxifp,imaxjfp,isurf), & 852 & igrdj(imaxifp,imaxjfp,isurf), & 853 & zglam(imaxifp,imaxjfp,isurf), & 854 & zgphi(imaxifp,imaxjfp,isurf), & 855 & zmask(imaxifp,imaxjfp,isurf), & 856 & zsurf(imaxifp,imaxjfp,isurf), & 857 & zsurftmp(imaxifp,imaxjfp,isurf), & 858 & zglamf(imaxifp+1,imaxjfp+1,isurf), & 859 & zgphif(imaxifp+1,imaxjfp+1,isurf), & 860 & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 861 & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 1229 862 & ) 1230 863 1231 864 DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 1232 865 iobs = jobs - surfdataqc%nsurfup 1233 igrdi(1,1,iobs) = surfdataqc%mi(jobs)-1 1234 igrdj(1,1,iobs) = surfdataqc%mj(jobs)-1 1235 igrdi(1,2,iobs) = surfdataqc%mi(jobs)-1 1236 igrdj(1,2,iobs) = surfdataqc%mj(jobs) 1237 igrdi(2,1,iobs) = surfdataqc%mi(jobs) 1238 igrdj(2,1,iobs) = surfdataqc%mj(jobs)-1 1239 igrdi(2,2,iobs) = surfdataqc%mi(jobs) 1240 igrdj(2,2,iobs) = surfdataqc%mj(jobs) 866 DO ji = 0, imaxifp 867 imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1 868 869 !Deal with wrap around in longitude 870 IF ( imodi < 1 ) imodi = imodi + jpiglo 871 IF ( imodi > jpiglo ) imodi = imodi - jpiglo 872 873 DO jj = 0, imaxjfp 874 imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1 875 !If model values are out of the domain to the north/south then 876 !set them to be the edge of the domain 877 IF ( imodj < 1 ) imodj = 1 878 IF ( imodj > jpjglo ) imodj = jpjglo 879 880 igrdip1(ji+1,jj+1,iobs) = imodi 881 igrdjp1(ji+1,jj+1,iobs) = imodj 882 883 IF ( ji >= 1 .AND. jj >= 1 ) THEN 884 igrdi(ji,jj,iobs) = imodi 885 igrdj(ji,jj,iobs) = imodj 886 ENDIF 887 888 END DO 889 END DO 1241 890 END DO 1242 891 1243 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &892 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 1244 893 & igrdi, igrdj, glamt, zglam ) 1245 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &894 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 1246 895 & igrdi, igrdj, gphit, zgphi ) 1247 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &896 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 1248 897 & igrdi, igrdj, psurfmask, zmask ) 1249 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &898 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 1250 899 & igrdi, igrdj, psurf, zsurf ) 900 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 901 & igrdip1, igrdjp1, glamf, zglamf ) 902 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 903 & igrdip1, igrdjp1, gphif, zgphif ) 1251 904 1252 905 ! At the end of the day get interpolated means 1253 IF (ldnightav ) THEN 1254 IF ( idayend == 0 ) THEN 1255 1256 ALLOCATE( & 1257 & zsurfm(2,2,isurf) & 1258 & ) 1259 1260 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, igrdi, igrdj, & 1261 & surfdataqc%vdmean(:,:), zsurfm ) 1262 1263 ENDIF 906 IF ( idayend == 0 .AND. ldnightav ) THEN 907 908 ALLOCATE( & 909 & zsurfm(imaxifp,imaxjfp,isurf) & 910 & ) 911 912 CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, & 913 & surfdataqc%vdmean(:,:), zsurfm ) 914 1264 915 ENDIF 1265 916 … … 1290 941 zphi = surfdataqc%rphi(jobs) 1291 942 1292 ! Get weights to interpolate the model value to the observation point1293 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &1294 & zglam(:,:,iobs), zgphi(:,:,iobs), &1295 & zmask(:,:,iobs), zweig, zobsmask )1296 1297 ! Interpolate the model field to the observation point1298 943 IF ( ldnightav .AND. idayend == 0 ) THEN 1299 944 ! Night-time averaged data 1300 CALL obs_int_h2d( 1, 1, zweig, zsurfm(:,:,iobs), zext)945 zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 1301 946 ELSE 1302 CALL obs_int_h2d( 1, 1, zweig, zsurf(:,:,iobs), zext ) 947 zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 948 ENDIF 949 950 IF ( k2dint <= 4 ) THEN 951 952 ! Get weights to interpolate the model value to the observation point 953 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 954 & zglam(:,:,iobs), zgphi(:,:,iobs), & 955 & zmask(:,:,iobs), zweig, zobsmask ) 956 957 ! Interpolate the model value to the observation point 958 CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 959 960 ELSE 961 962 ! Get weights to average the model SLA to the observation footprint 963 CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam, zphi, & 964 & zglam(:,:,iobs), zgphi(:,:,iobs), & 965 & zglamf(:,:,iobs), zgphif(:,:,iobs), & 966 & zmask(:,:,iobs), plamscl, pphiscl, & 967 & lindegrees, zweig, zobsmask ) 968 969 ! Average the model SST to the observation footprint 970 CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 971 & zweig, zsurftmp(:,:,iobs), zext ) 972 1303 973 ENDIF 1304 974 … … 1310 980 surfdataqc%rmod(jobs,1) = zext(1) 1311 981 ENDIF 982 983 IF ( zext(1) == obfillflt ) THEN 984 ! If the observation value is a fill value, set QC flag to bad 985 surfdataqc%nqc(jobs) = 4 986 ENDIF 1312 987 1313 988 END DO … … 1315 990 ! Deallocate the data for interpolation 1316 991 DEALLOCATE( & 992 & zweig, & 1317 993 & igrdi, & 1318 994 & igrdj, & … … 1320 996 & zgphi, & 1321 997 & zmask, & 1322 & zsurf & 998 & zsurf, & 999 & zsurftmp, & 1000 & zglamf, & 1001 & zgphif, & 1002 & igrdip1,& 1003 & igrdjp1 & 1323 1004 & ) 1324 1005 1325 1006 ! At the end of the day also deallocate night-time mean array 1326 IF ( ldnightav ) THEN 1327 IF ( idayend == 0 ) THEN 1328 DEALLOCATE( & 1329 & zsurfm & 1330 & ) 1331 ENDIF 1007 IF ( idayend == 0 .AND. ldnightav ) THEN 1008 DEALLOCATE( & 1009 & zsurfm & 1010 & ) 1332 1011 ENDIF 1333 1012 -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
r7646 r9023 23 23 USE obs_oper ! Observation operators 24 24 USE lib_mpp, ONLY : ctl_warn, ctl_stop 25 USE bdy_oce, ONLY : & ! Boundary information 26 idx_bdy, nb_bdy, ln_bdy 25 27 26 28 IMPLICIT NONE … … 40 42 CONTAINS 41 43 42 SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea ) 44 SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, & 45 kqc_cutoff ) 43 46 !!---------------------------------------------------------------------- 44 47 !! *** ROUTINE obs_pre_sla *** … … 57 60 !! ! 2015-02 (M. Martin) Combined routine for surface types. 58 61 !!---------------------------------------------------------------------- 62 !! * Modules used 59 63 USE par_oce ! Ocean parameters 60 64 USE dom_oce, ONLY : glamt, gphit, tmask, nproc ! Geographical information … … 63 67 TYPE(obs_surf), INTENT(INOUT) :: surfdataqc ! Subset of surface data not failing screening 64 68 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 65 ! 69 LOGICAL, INTENT(IN) :: ld_bound_reject ! Switch for rejecting obs near the boundary 70 INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value 71 !! * Local declarations 72 INTEGER :: iqc_cutoff = 255 ! cut off for QC value 66 73 INTEGER :: iyea0 ! Initial date 67 74 INTEGER :: imon0 ! - (year, month, day, hour, minute) … … 76 83 INTEGER :: inlasobs ! - close to land 77 84 INTEGER :: igrdobs ! - fail the grid search 85 INTEGER :: ibdysobs ! - close to open boundary 78 86 ! Global counters for observations that 79 87 INTEGER :: iotdobsmpp ! - outside time domain … … 82 90 INTEGER :: inlasobsmpp ! - close to land 83 91 INTEGER :: igrdobsmpp ! - fail the grid search 84 LOGICAL, DIMENSION(:), ALLOCATABLE :: llvalid ! SLA data selection 92 INTEGER :: ibdysobsmpp ! - close to open boundary 93 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 94 & llvalid ! SLA data selection 85 95 INTEGER :: jobs ! Obs. loop variable 86 96 INTEGER :: jstp ! Time loop variable … … 107 117 ilansobs = 0 108 118 inlasobs = 0 119 ibdysobs = 0 120 121 ! Set QC cutoff to optional value if provided 122 IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 109 123 110 124 ! ----------------------------------------------------------------------- … … 140 154 & tmask(:,:,1), surfdata%nqc, & 141 155 & iosdsobs, ilansobs, & 142 & inlasobs, ld_nea ) 156 & inlasobs, ld_nea, & 157 & ibdysobs, ld_bound_reject, & 158 & iqc_cutoff ) 143 159 144 160 CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 145 161 CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 146 162 CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 163 CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 147 164 148 165 ! ----------------------------------------------------------------------- … … 155 172 ALLOCATE( llvalid(surfdata%nsurf) ) 156 173 157 ! We want all data which has qc flags <= 10158 159 llvalid(:) = ( surfdata%nqc(:) <= 10)174 ! We want all data which has qc flags <= iqc_cutoff 175 176 llvalid(:) = ( surfdata%nqc(:) <= iqc_cutoff ) 160 177 161 178 ! The actual copying … … 190 207 & inlasobsmpp 191 208 ENDIF 209 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', & 210 & ibdysobsmpp 192 211 WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted = ', & 193 212 & surfdataqc%nsurfmpp … … 225 244 & kpi, kpj, kpk, & 226 245 & zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2, & 227 & ld_nea, kdailyavtypes)246 & ld_nea, ld_bound_reject, kdailyavtypes, kqc_cutoff ) 228 247 229 248 !!---------------------------------------------------------------------- … … 241 260 !! 242 261 !!---------------------------------------------------------------------- 243 USE par_oce ! Ocean parameters 244 USE dom_oce, ONLY : gdept_1d, nproc ! Geographical information 262 !! * Modules used 263 USE par_oce ! Ocean parameters 264 USE dom_oce, ONLY : & ! Geographical information 265 & gdept_1d, & 266 & nproc 245 267 246 268 !! * Arguments … … 250 272 LOGICAL, INTENT(IN) :: ld_var2 251 273 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 274 LOGICAL, INTENT(IN) :: ld_bound_reject ! Switch for rejecting observations near the boundary 252 275 INTEGER, INTENT(IN) :: kpi, kpj, kpk ! Local domain sizes 253 276 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & … … 261 284 & pgphi1, & 262 285 & pgphi2 286 INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value 263 287 264 288 !! * Local declarations 289 INTEGER :: iqc_cutoff = 255 ! cut off for QC value 265 290 INTEGER :: iyea0 ! Initial date 266 291 INTEGER :: imon0 ! - (year, month, day, hour, minute) … … 277 302 INTEGER :: inlav1obs ! - close to land (variable 1) 278 303 INTEGER :: inlav2obs ! - close to land (variable 2) 304 INTEGER :: ibdyv1obs ! - boundary (variable 1) 305 INTEGER :: ibdyv2obs ! - boundary (variable 2) 279 306 INTEGER :: igrdobs ! - fail the grid search 280 307 INTEGER :: iuvchku ! - reject u if v rejected and vice versa … … 288 315 INTEGER :: inlav1obsmpp ! - close to land (variable 1) 289 316 INTEGER :: inlav2obsmpp ! - close to land (variable 2) 317 INTEGER :: ibdyv1obsmpp ! - boundary (variable 1) 318 INTEGER :: ibdyv2obsmpp ! - boundary (variable 2) 290 319 INTEGER :: igrdobsmpp ! - fail the grid search 291 320 INTEGER :: iuvchkumpp ! - reject var1 if var2 rejected and vice versa … … 322 351 inlav1obs = 0 323 352 inlav2obs = 0 324 iuvchku = 0 325 iuvchkv = 0 353 ibdyv1obs = 0 354 ibdyv2obs = 0 355 iuvchku = 0 356 iuvchkv = 0 357 358 359 ! Set QC cutoff to optional value if provided 360 IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 326 361 327 362 ! ----------------------------------------------------------------------- … … 335 370 & profdata%nday, profdata%nhou, profdata%nmin, & 336 371 & profdata%ntyp, profdata%nqc, profdata%mstp, & 337 & iotdobs, kdailyavtypes = kdailyavtypes ) 372 & iotdobs, kdailyavtypes = kdailyavtypes, & 373 & kqc_cutoff = iqc_cutoff ) 338 374 ELSE 339 375 CALL obs_coo_tim_prof( icycle, & … … 342 378 & profdata%nday, profdata%nhou, profdata%nmin, & 343 379 & profdata%ntyp, profdata%nqc, profdata%mstp, & 344 & iotdobs )380 & iotdobs, kqc_cutoff = iqc_cutoff ) 345 381 ENDIF 346 382 … … 359 395 360 396 ! ----------------------------------------------------------------------- 361 ! Reject all observations for profiles with nqc > 10362 ! ----------------------------------------------------------------------- 363 364 CALL obs_pro_rej( profdata )397 ! Reject all observations for profiles with nqc > iqc_cutoff 398 ! ----------------------------------------------------------------------- 399 400 CALL obs_pro_rej( profdata, kqc_cutoff = iqc_cutoff ) 365 401 366 402 ! ----------------------------------------------------------------------- … … 381 417 & gdept_1d, zmask1, & 382 418 & profdata%nqc, profdata%var(1)%nvqc, & 383 & iosdv1obs, ilanv1obs, & 384 & inlav1obs, ld_nea ) 419 & iosdv1obs, ilanv1obs, & 420 & inlav1obs, ld_nea, & 421 & ibdyv1obs, ld_bound_reject, & 422 & iqc_cutoff ) 385 423 386 424 CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp ) 387 425 CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp ) 388 426 CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp ) 427 CALL obs_mpp_sum_integer( ibdyv1obs, ibdyv1obsmpp ) 389 428 390 429 ! Variable 2 … … 400 439 & gdept_1d, zmask2, & 401 440 & profdata%nqc, profdata%var(2)%nvqc, & 402 & iosdv2obs, ilanv2obs, & 403 & inlav2obs, ld_nea ) 441 & iosdv2obs, ilanv2obs, & 442 & inlav2obs, ld_nea, & 443 & ibdyv2obs, ld_bound_reject, & 444 & iqc_cutoff ) 404 445 405 446 CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp ) 406 447 CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp ) 407 448 CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp ) 449 CALL obs_mpp_sum_integer( ibdyv2obs, ibdyv2obsmpp ) 408 450 409 451 ! ----------------------------------------------------------------------- … … 412 454 413 455 IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 414 CALL obs_uv_rej( profdata, iuvchku, iuvchkv )456 CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff ) 415 457 CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp ) 416 458 CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp ) … … 429 471 END DO 430 472 431 ! We want all data which has qc flags = 0432 433 llvalid%luse(:) = ( profdata%nqc(:) <= 10)473 ! We want all data which has qc flags <= iqc_cutoff 474 475 llvalid%luse(:) = ( profdata%nqc(:) <= iqc_cutoff ) 434 476 DO jvar = 1,profdata%nvar 435 llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10)477 llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff ) 436 478 END DO 437 479 … … 475 517 & iuvchku 476 518 ENDIF 519 WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near open boundary (removed) = ',& 520 & ibdyv1obsmpp 477 521 WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted = ', & 478 522 & prodatqc%nvprotmpp(1) … … 492 536 & iuvchkv 493 537 ENDIF 538 WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near open boundary (removed) = ',& 539 & ibdyv2obsmpp 494 540 WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted = ', & 495 541 & prodatqc%nvprotmpp(2) … … 644 690 & .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN 645 691 kobsstp(jobs) = -1 646 kobsqc(jobs) = kobsqc(jobs) + 11692 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 647 693 kotdobs = kotdobs + 1 648 694 CYCLE … … 695 741 IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) & 696 742 & .OR.( kobsstp(jobs) > nitend ) ) THEN 697 kobsqc(jobs) = kobsqc(jobs) + 12743 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 698 744 kotdobs = kotdobs + 1 699 745 CYCLE … … 739 785 & kobsno, & 740 786 & kobsyea, kobsmon, kobsday, kobshou, kobsmin, & 741 & ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes ) 787 & ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes, & 788 & kqc_cutoff ) 742 789 !!---------------------------------------------------------------------- 743 790 !! *** ROUTINE obs_coo_tim *** … … 783 830 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 784 831 & kdailyavtypes ! Types for daily averages 832 INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff ! QC cutoff value 833 785 834 !! * Local declarations 786 835 INTEGER :: jobs 836 INTEGER :: iqc_cutoff=255 787 837 788 838 !----------------------------------------------------------------------- … … 803 853 DO jobs = 1, kobsno 804 854 805 IF ( kobsqc(jobs) <= 10) THEN855 IF ( kobsqc(jobs) <= iqc_cutoff ) THEN 806 856 807 857 IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.& 808 858 & ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN 809 kobsqc(jobs) = kobsqc(jobs) + 14859 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 810 860 kotdobs = kotdobs + 1 811 861 CYCLE … … 850 900 DO jobs = 1, kobsno 851 901 IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN 852 kobsqc(jobs) = kobsqc(jobs) + 18902 kobsqc(jobs) = IBSET(kobsqc(jobs),12) 853 903 kgrdobs = kgrdobs + 1 854 904 ENDIF … … 861 911 & plam, pphi, pmask, & 862 912 & kobsqc, kosdobs, klanobs, & 863 & knlaobs,ld_nea ) 913 & knlaobs,ld_nea, & 914 & kbdyobs,ld_bound_reject, & 915 & kqc_cutoff ) 864 916 !!---------------------------------------------------------------------- 865 917 !! *** ROUTINE obs_coo_spc_2d *** … … 894 946 INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: & 895 947 & kobsqc ! Observation quality control 896 INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain 897 INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell 898 INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land 899 LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land 948 INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain 949 INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell 950 INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land 951 INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary 952 LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land 953 LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary 954 INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value 955 900 956 !! * Local declarations 901 957 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 902 958 & zgmsk ! Grid mask 959 960 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 961 & zbmsk ! Boundary mask 962 REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 903 963 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 904 964 & zglam, & ! Model longitude at grid points … … 917 977 ! For invalid points use 2,2 918 978 919 IF ( kobsqc(jobs) >= 10) THEN979 IF ( kobsqc(jobs) >= kqc_cutoff ) THEN 920 980 921 981 igrdi(1,1,jobs) = 1 … … 942 1002 943 1003 END DO 1004 1005 IF (ln_bdy) THEN 1006 ! Create a mask grid points in boundary rim 1007 IF (ld_bound_reject) THEN 1008 zbdymask(:,:) = 1.0_wp 1009 DO ji = 1, nb_bdy 1010 DO jj = 1, idx_bdy(ji)%nblen(1) 1011 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 1012 ENDDO 1013 ENDDO 1014 1015 CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk ) 1016 ENDIF 1017 ENDIF 1018 944 1019 945 1020 CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk ) … … 950 1025 951 1026 ! Skip bad observations 952 IF ( kobsqc(jobs) >= 10) CYCLE1027 IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE 953 1028 954 1029 ! Flag if the observation falls outside the model spatial domain … … 957 1032 & .OR. ( pobsphi(jobs) < -90. ) & 958 1033 & .OR. ( pobsphi(jobs) > 90. ) ) THEN 959 kobsqc(jobs) = kobsqc(jobs) + 111034 kobsqc(jobs) = IBSET(kobsqc(jobs),11) 960 1035 kosdobs = kosdobs + 1 961 1036 CYCLE … … 964 1039 ! Flag if the observation falls with a model land cell 965 1040 IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 966 kobsqc(jobs) = kobsqc(jobs) + 121041 kobsqc(jobs) = IBSET(kobsqc(jobs),10) 967 1042 klanobs = klanobs + 1 968 1043 CYCLE … … 978 1053 IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) & 979 1054 & .AND. & 980 & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp )&981 & ) THEN1055 & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) & 1056 & < 1.0e-6_wp ) ) THEN 982 1057 lgridobs = .TRUE. 983 1058 iig = ji … … 992 1067 IF (lgridobs) THEN 993 1068 IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN 994 kobsqc(jobs) = kobsqc(jobs) + 121069 kobsqc(jobs) = IBSET(kobsqc(jobs),10) 995 1070 klanobs = klanobs + 1 996 1071 CYCLE … … 1000 1075 ! Flag if the observation falls is close to land 1001 1076 IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN 1002 IF (ld_nea) kobsqc(jobs) = kobsqc(jobs) + 141003 1077 knlaobs = knlaobs + 1 1004 CYCLE 1078 IF (ld_nea) THEN 1079 kobsqc(jobs) = IBSET(kobsqc(jobs),9) 1080 CYCLE 1081 ENDIF 1082 ENDIF 1083 1084 IF (ln_bdy) THEN 1085 ! Flag if the observation falls close to the boundary rim 1086 IF (ld_bound_reject) THEN 1087 IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 1088 kobsqc(jobs) = IBSET(kobsqc(jobs),8) 1089 kbdyobs = kbdyobs + 1 1090 CYCLE 1091 ENDIF 1092 ! for observations on the grid... 1093 IF (lgridobs) THEN 1094 IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 1095 kobsqc(jobs) = IBSET(kobsqc(jobs),8) 1096 kbdyobs = kbdyobs + 1 1097 CYCLE 1098 ENDIF 1099 ENDIF 1100 ENDIF 1005 1101 ENDIF 1006 1102 … … 1015 1111 & plam, pphi, pdep, pmask, & 1016 1112 & kpobsqc, kobsqc, kosdobs, & 1017 & klanobs, knlaobs, ld_nea ) 1113 & klanobs, knlaobs, ld_nea, & 1114 & kbdyobs, ld_bound_reject, & 1115 & kqc_cutoff ) 1018 1116 !!---------------------------------------------------------------------- 1019 1117 !! *** ROUTINE obs_coo_spc_3d *** … … 1077 1175 INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell 1078 1176 INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land 1177 INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary 1079 1178 LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land 1179 LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary 1180 INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value 1181 1080 1182 !! * Local declarations 1081 1183 REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 1082 1184 & zgmsk ! Grid mask 1185 REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 1186 & zbmsk ! Boundary mask 1187 REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 1083 1188 REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 1084 1189 & zgdepw … … 1100 1205 ! For invalid points use 2,2 1101 1206 1102 IF ( kpobsqc(jobs) >= 10) THEN1207 IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN 1103 1208 1104 1209 igrdi(1,1,jobs) = 1 … … 1125 1230 1126 1231 END DO 1232 1233 IF (ln_bdy) THEN 1234 ! Create a mask grid points in boundary rim 1235 IF (ld_bound_reject) THEN 1236 zbdymask(:,:) = 1.0_wp 1237 DO ji = 1, nb_bdy 1238 DO jj = 1, idx_bdy(ji)%nblen(1) 1239 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 1240 ENDDO 1241 ENDDO 1242 ENDIF 1243 1244 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk ) 1245 ENDIF 1127 1246 1128 1247 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk ) 1129 1248 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam ) 1130 1249 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 1131 IF ( .NOT.( ln_zps .OR. ln_zco ) ) THEN 1132 ! Need to know the bathy depth for each observation for sco 1133 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), & 1250 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), & 1134 1251 & zgdepw ) 1135 ENDIF1136 1252 1137 1253 DO jobs = 1, kprofno 1138 1254 1139 1255 ! Skip bad profiles 1140 IF ( kpobsqc(jobs) >= 10) CYCLE1256 IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE 1141 1257 1142 1258 ! Check if this observation is on a grid point … … 1149 1265 IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) & 1150 1266 & .AND. & 1151 & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &1267 & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) < 1.0e-6_wp ) & 1152 1268 & ) THEN 1153 1269 lgridobs = .TRUE. … … 1176 1292 & .OR. ( pobsdep(jobsp) < 0.0 ) & 1177 1293 & .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN 1178 kobsqc(jobsp) = kobsqc(jobsp) + 111294 kobsqc(jobsp) = IBSET(kobsqc(jobsp),11) 1179 1295 kosdobs = kosdobs + 1 1180 1296 CYCLE 1181 1297 ENDIF 1182 1298 1183 ! To check if an observations falls within land there are two cases: 1184 ! 1: z-coordibnates, where the check uses the mask 1185 ! 2: terrain following (eg s-coordinates), 1186 ! where we use the depth of the bottom cell to mask observations 1299 ! To check if an observations falls within land: 1187 1300 1188 IF( ln_zps .OR. ln_zco ) THEN !(CASE 1) 1301 ! Flag if the observation is deeper than the bathymetry 1302 ! Or if it is within the mask 1303 IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 1304 & .OR. & 1305 & ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1306 & == 0.0_wp) ) THEN 1307 kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 1308 klanobs = klanobs + 1 1309 CYCLE 1310 ENDIF 1189 1311 1190 ! Flag if the observation falls with a model land cell 1191 IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1192 & == 0.0_wp ) THEN 1193 kobsqc(jobsp) = kobsqc(jobsp) + 12 1194 klanobs = klanobs + 1 1195 CYCLE 1196 ENDIF 1197 1198 ! Flag if the observation is close to land 1199 IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 1200 & 0.0_wp) THEN 1201 knlaobs = knlaobs + 1 1202 IF (ld_nea) THEN 1203 kobsqc(jobsp) = kobsqc(jobsp) + 14 1204 ENDIF 1205 ENDIF 1206 1207 ELSE ! Case 2 1208 1209 ! Flag if the observation is deeper than the bathymetry 1210 ! Or if it is within the mask 1211 IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 1212 & .OR. & 1213 & ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1214 & == 0.0_wp) ) THEN 1215 kobsqc(jobsp) = kobsqc(jobsp) + 12 1216 klanobs = klanobs + 1 1217 CYCLE 1218 ENDIF 1219 1220 ! Flag if the observation is close to land 1221 IF ( ll_next_to_land ) THEN 1222 knlaobs = knlaobs + 1 1223 IF (ld_nea) THEN 1224 kobsqc(jobsp) = kobsqc(jobsp) + 14 1225 ENDIF 1226 ENDIF 1312 ! Flag if the observation is close to land 1313 IF ( ll_next_to_land ) THEN 1314 knlaobs = knlaobs + 1 1315 IF (ld_nea) THEN 1316 kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 1317 ENDIF 1227 1318 ENDIF 1228 1319 … … 1232 1323 IF (lgridobs) THEN 1233 1324 IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN 1234 kobsqc(jobsp) = kobsqc(jobsp) + 121325 kobsqc(jobsp) = IBSET(kobsqc(jobs),10) 1235 1326 klanobs = klanobs + 1 1236 1327 CYCLE … … 1250 1341 ENDIF 1251 1342 1343 IF (ln_bdy) THEN 1344 ! Flag if the observation falls close to the boundary rim 1345 IF (ld_bound_reject) THEN 1346 IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 1347 kobsqc(jobsp) = IBSET(kobsqc(jobs),8) 1348 kbdyobs = kbdyobs + 1 1349 CYCLE 1350 ENDIF 1351 ! for observations on the grid... 1352 IF (lgridobs) THEN 1353 IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 1354 kobsqc(jobsp) = IBSET(kobsqc(jobs),8) 1355 kbdyobs = kbdyobs + 1 1356 CYCLE 1357 ENDIF 1358 ENDIF 1359 ENDIF 1360 ENDIF 1361 1252 1362 END DO 1253 1363 END DO … … 1255 1365 END SUBROUTINE obs_coo_spc_3d 1256 1366 1257 SUBROUTINE obs_pro_rej( profdata )1367 SUBROUTINE obs_pro_rej( profdata, kqc_cutoff ) 1258 1368 !!---------------------------------------------------------------------- 1259 1369 !! *** ROUTINE obs_pro_rej *** … … 1273 1383 !! * Arguments 1274 1384 TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data 1385 INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value 1386 1275 1387 !! * Local declarations 1276 1388 INTEGER :: jprof … … 1282 1394 DO jprof = 1, profdata%nprof 1283 1395 1284 IF ( profdata%nqc(jprof) > 10) THEN1396 IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN 1285 1397 1286 1398 DO jvar = 1, profdata%nvar … … 1290 1402 1291 1403 profdata%var(jvar)%nvqc(jobs) = & 1292 & profdata%var(jvar)%nvqc(jobs) + 261404 & IBSET(profdata%var(jvar)%nvqc(jobs),14) 1293 1405 1294 1406 END DO … … 1302 1414 END SUBROUTINE obs_pro_rej 1303 1415 1304 SUBROUTINE obs_uv_rej( profdata, knumu, knumv )1416 SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff ) 1305 1417 !!---------------------------------------------------------------------- 1306 1418 !! *** ROUTINE obs_uv_rej *** … … 1322 1434 INTEGER, INTENT(INOUT) :: knumu ! Number of u rejected 1323 1435 INTEGER, INTENT(INOUT) :: knumv ! Number of v rejected 1436 INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value 1437 1324 1438 !! * Local declarations 1325 1439 INTEGER :: jprof … … 1341 1455 DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1) 1342 1456 1343 IF ( ( profdata%var(1)%nvqc(jobs) > 10) .AND. &1344 & ( profdata%var(2)%nvqc(jobs) <= 10) ) THEN1345 profdata%var(2)%nvqc(jobs) = profdata%var(2)%nvqc(jobs) + 421457 IF ( ( profdata%var(1)%nvqc(jobs) > kqc_cutoff ) .AND. & 1458 & ( profdata%var(2)%nvqc(jobs) <= kqc_cutoff) ) THEN 1459 profdata%var(2)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 1346 1460 knumv = knumv + 1 1347 1461 ENDIF 1348 IF ( ( profdata%var(2)%nvqc(jobs) > 10) .AND. &1349 & ( profdata%var(1)%nvqc(jobs) <= 10) ) THEN1350 profdata%var(1)%nvqc(jobs) = profdata%var(1)%nvqc(jobs) + 421462 IF ( ( profdata%var(2)%nvqc(jobs) > kqc_cutoff ) .AND. & 1463 & ( profdata%var(1)%nvqc(jobs) <= kqc_cutoff) ) THEN 1464 profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 1351 1465 knumu = knumu + 1 1352 1466 ENDIF -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90
r6140 r9023 104 104 ! Bookkeeping arrays with sizes equal to number of variables 105 105 106 CHARACTER(len= 6), POINTER, DIMENSION(:) :: &106 CHARACTER(len=8), POINTER, DIMENSION(:) :: & 107 107 & cvars !: Variable names 108 108 -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90
r6140 r9023 87 87 CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof' 88 88 CHARACTER(len=8) :: clrefdate 89 CHARACTER(len= 6), DIMENSION(:), ALLOCATABLE :: clvars89 CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 90 90 INTEGER :: jvar 91 91 INTEGER :: ji … … 307 307 inowin = 0 308 308 DO ji = 1, inpfiles(jj)%nobs 309 IF ( inpfiles(jj)%ioqc(ji) > 2) CYCLE310 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2) .AND. &311 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE309 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 312 312 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 313 313 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 325 325 inowin = 0 326 326 DO ji = 1, inpfiles(jj)%nobs 327 IF ( inpfiles(jj)%ioqc(ji) > 2) CYCLE328 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2) .AND. &329 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE327 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 330 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 331 331 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 351 351 inowin = 0 352 352 DO ji = 1, inpfiles(jj)%nobs 353 IF ( inpfiles(jj)%ioqc(ji) > 2) CYCLE354 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2) .AND. &355 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE353 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 356 356 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 357 357 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 373 373 374 374 DO ji = 1, inpfiles(jj)%nobs 375 IF ( inpfiles(jj)%ioqc(ji) > 2) CYCLE376 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2) .AND. &377 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE375 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 378 378 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 379 379 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 388 388 IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 389 389 & CYCLE 390 IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2) .AND. &391 & ( inpfiles(jj)%idqc(ij,ji) <= 2) ) THEN390 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 391 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 392 392 ivar1t0 = ivar1t0 + 1 393 393 ENDIF … … 398 398 IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 399 399 & CYCLE 400 IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2) .AND. &401 & ( inpfiles(jj)%idqc(ij,ji) <= 2) ) THEN400 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 401 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 402 402 ivar2t0 = ivar2t0 + 1 403 403 ENDIF … … 407 407 IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 408 408 & CYCLE 409 IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2) .AND. &410 & ( inpfiles(jj)%idqc(ij,ji) <= 2 )) .AND. &411 & 412 & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2) .AND. &413 & ( inpfiles(jj)%idqc(ij,ji) <= 2 )) .AND. &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 414 & ldvar2 ) ) THEN 415 415 ip3dt = ip3dt + 1 … … 437 437 DO jj = 1, inobf 438 438 DO ji = 1, inpfiles(jj)%nobs 439 IF ( inpfiles(jj)%ioqc(ji) > 2) CYCLE440 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2) .AND. &441 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE439 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 442 442 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 443 443 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 452 452 DO jj = 1, inobf 453 453 DO ji = 1, inpfiles(jj)%nobs 454 IF ( inpfiles(jj)%ioqc(ji) > 2) CYCLE455 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2) .AND. &456 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE454 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 457 457 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 458 458 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 501 501 ji = iprofidx(iindx(jk)) 502 502 503 IF ( inpfiles(jj)%ioqc(ji) > 2) CYCLE504 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2) .AND. &505 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE503 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 506 506 507 507 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & … … 518 518 IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 519 519 520 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 521 & ( inpfiles(jj)%ivqc(ji,2) > 2 ) ) CYCLE 520 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 522 523 523 524 loop_prof : DO ij = 1, inpfiles(jj)%nlev … … 526 527 & CYCLE 527 528 528 IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2) .AND. &529 & ( inpfiles(jj)%idqc(ij,ji) <= 2) ) THEN529 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 530 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 530 531 531 532 llvalprof = .TRUE. … … 534 535 ENDIF 535 536 536 IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2) .AND. &537 & ( inpfiles(jj)%idqc(ij,ji) <= 2) ) THEN537 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 538 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 538 539 539 540 llvalprof = .TRUE. … … 615 616 IF (ldsatt) THEN 616 617 617 IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2) .AND. &618 & ( inpfiles(jj)%idqc(ij,ji) <= 2 )) .AND. &619 & 620 & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2) .AND. &621 & ( inpfiles(jj)%idqc(ij,ji) <= 2 )) .AND. &622 & 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 623 624 ip3dt = ip3dt + 1 624 625 ELSE … … 628 629 ENDIF 629 630 630 IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2) .AND. &631 & ( inpfiles(jj)%idqc(ij,ji) <= 2 )) .AND. &632 &ldvar1 ) .OR. ldsatt ) THEN631 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 633 634 634 635 IF (ldsatt) THEN … … 661 662 662 663 ! Profile var1 value 663 IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2) .AND. &664 & ( inpfiles(jj)%idqc(ij,ji) <= 2) ) THEN664 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 665 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 665 666 profdata%var(1)%vobs(ivar1t) = & 666 667 & inpfiles(jj)%pob(ij,ji,1) … … 692 693 ENDIF 693 694 694 IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2) .AND. &695 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ).AND. &695 IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 696 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 696 697 & ldvar2 ) .OR. ldsatt ) THEN 697 698 … … 725 726 726 727 ! Profile var2 value 727 IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2) .AND. &728 & ( inpfiles(jj)%idqc(ij,ji) <= 2 )) THEN728 IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) ) .AND. & 729 & ( .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) ) THEN 729 730 profdata%var(2)%vobs(ivar2t) = & 730 731 & inpfiles(jj)%pob(ij,ji,2) -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90
r6140 r9023 77 77 CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_surf' 78 78 CHARACTER(len=8) :: clrefdate 79 CHARACTER(len= 6), DIMENSION(:), ALLOCATABLE :: clvars79 CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 80 80 INTEGER :: ji 81 81 INTEGER :: jj … … 172 172 173 173 !------------------------------------------------------------------ 174 ! Read the profile file into inpfiles174 ! Read the surface file into inpfiles 175 175 !------------------------------------------------------------------ 176 176 CALL init_obfbdata( inpfiles(jj) ) … … 296 296 ENDIF 297 297 llvalprof = .FALSE. 298 IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. & 299 & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN 298 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN 300 299 iobs = iobs + 1 301 300 ENDIF … … 370 369 ! Set observation information 371 370 372 IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. & 373 & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN 371 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN 374 372 375 373 iobs = iobs + 1 -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90
r6140 r9023 154 154 155 155 ! mark any masked data with a QC flag 156 IF( zobsmask(1) == 0 ) sladata%nqc(jobs) = 11156 IF( zobsmask(1) == 0 ) sladata%nqc(jobs) = IBSET(sladata%nqc(jobs),15) 157 157 158 158 END DO -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_sstbias.F90
r6140 r9023 1 1 MODULE obs_sstbias 2 2 !!====================================================================== 3 !! *** MODULE obs_ readsstbias ***4 !! Observation diagnostics: Read the bias for S LAdata3 !! *** MODULE obs_sstbias *** 4 !! Observation diagnostics: Read the bias for SST data 5 5 !!====================================================================== 6 6 !!---------------------------------------------------------------------- 7 !! obs_ rea_sstbias : Driver for reading altimeterbias7 !! obs_app_sstbias : Driver for reading and applying the SST bias 8 8 !!---------------------------------------------------------------------- 9 9 !! * Modules used … … 139 139 cl_bias_files(jtype) ) 140 140 ! Get the SST bias data 141 CALL iom_get( numsstbias, jpdom_data, ' sst', z_sstbias_2d(:,:), 1 )141 CALL iom_get( numsstbias, jpdom_data, 'tn', z_sstbias_2d(:,:), 1 ) 142 142 z_sstbias(:,:,jtype) = z_sstbias_2d(:,:) 143 143 ! Close the file … … 190 190 igrdj_tmp(:,:,jt) = igrdj(:,:,jobs) 191 191 zglam_tmp(:,:,jt) = zglam(:,:,jobs) 192 zgphi_tmp(:,:,jt) = zgphi(:,:,jobs)193 192 zgphi_tmp(:,:,jt) = zgphi(:,:,jobs) 194 193 zmask_tmp(:,:,jt) = zmask(:,:,jobs) -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90
r6140 r9023 50 50 INTEGER :: npj 51 51 INTEGER :: nsurfup !: Observation counter used in obs_oper 52 INTEGER :: nrec !: Number of surface observation records in window 52 53 53 54 ! Arrays with size equal to the number of surface observations … … 56 57 & mi, & !: i-th grid coord. for interpolating to surface observation 57 58 & mj, & !: j-th grid coord. for interpolating to surface observation 59 & mt, & !: time record number for gridded data 58 60 & nsidx,& !: Surface observation number 59 61 & nsfil,& !: Surface observation number in file … … 67 69 & ntyp !: Type of surface observation product 68 70 69 CHARACTER(len= 6), POINTER, DIMENSION(:) :: &71 CHARACTER(len=8), POINTER, DIMENSION(:) :: & 70 72 & cvars !: Variable names 71 73 … … 93 95 & nsstpmpp !: Global number of surface observations per time step 94 96 97 ! Arrays with size equal to the number of observation records in the window 98 INTEGER, POINTER, DIMENSION(:) :: & 99 & mrecstp ! Time step of the records 100 95 101 ! Arrays used to store source indices when 96 102 ! compressing obs_surf derived types … … 100 106 INTEGER, POINTER, DIMENSION(:) :: & 101 107 & nsind !: Source indices of surface data in compressed data 108 109 ! Is this a gridded product? 110 111 LOGICAL :: lgrid 102 112 103 113 END TYPE obs_surf … … 160 170 & surf%mi(ksurf), & 161 171 & surf%mj(ksurf), & 172 & surf%mt(ksurf), & 162 173 & surf%nsidx(ksurf), & 163 174 & surf%nsfil(ksurf), & … … 176 187 & ) 177 188 189 surf%mt(:) = -1 190 178 191 179 192 ! Allocate arrays of number of surface data size * number of variables … … 190 203 & ) 191 204 205 surf%rext(:,:) = 0.0_wp 206 192 207 ! Allocate arrays of number of time step size 193 208 … … 217 232 218 233 surf%nsurfup = 0 234 235 ! Not gridded by default 236 237 surf%lgrid = .FALSE. 219 238 220 239 END SUBROUTINE obs_surf_alloc … … 242 261 & surf%mi, & 243 262 & surf%mj, & 263 & surf%mt, & 244 264 & surf%nsidx, & 245 265 & surf%nsfil, & … … 370 390 newsurf%mi(insurf) = surf%mi(ji) 371 391 newsurf%mj(insurf) = surf%mj(ji) 392 newsurf%mt(insurf) = surf%mt(ji) 372 393 newsurf%nsidx(insurf) = surf%nsidx(ji) 373 394 newsurf%nsfil(insurf) = surf%nsfil(ji) … … 414 435 newsurf%nstp = surf%nstp 415 436 newsurf%cvars(:) = surf%cvars(:) 437 438 ! Set gridded stuff 439 440 newsurf%mt(insurf) = surf%mt(ji) 416 441 417 442 ! Deallocate temporary data … … 454 479 oldsurf%mi(jj) = surf%mi(ji) 455 480 oldsurf%mj(jj) = surf%mj(ji) 481 oldsurf%mt(jj) = surf%mt(ji) 456 482 oldsurf%nsidx(jj) = surf%nsidx(ji) 457 483 oldsurf%nsfil(jj) = surf%nsfil(ji) -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90
r6140 r9023 8 8 !! obs_wri_prof : Write profile observations in feedback format 9 9 !! obs_wri_surf : Write surface observations in feedback format 10 !! obs_wri_stats : Print basic statistics on the data being written out10 !! obs_wri_stats : Print basic statistics on the data being written out 11 11 !!---------------------------------------------------------------------- 12 12 … … 83 83 TYPE(obfbdata) :: fbdata 84 84 CHARACTER(LEN=40) :: clfname 85 CHARACTER(LEN= 6) :: clfiletype85 CHARACTER(LEN=10) :: clfiletype 86 86 INTEGER :: ilevel 87 87 INTEGER :: jvar … … 196 196 fbdata%ivqc(jo,:) = profdata%ivqc(jo,:) 197 197 fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:) 198 IF ( profdata%nqc(jo) > 10) THEN199 fbdata%ioqc(jo) = 4198 IF ( profdata%nqc(jo) > 255 ) THEN 199 fbdata%ioqc(jo) = IBSET(profdata%nqc(jo),2) 200 200 fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo) 201 fbdata%ioqcf(2,jo) = profdata%nqc(jo) - 10201 fbdata%ioqcf(2,jo) = profdata%nqc(jo) 202 202 ELSE 203 203 fbdata%ioqc(jo) = profdata%nqc(jo) … … 236 236 fbdata%idqc(ik,jo) = profdata%var(jvar)%idqc(jk) 237 237 fbdata%idqcf(:,ik,jo) = profdata%var(jvar)%idqcf(:,jk) 238 IF ( profdata%var(jvar)%nvqc(jk) > 10) THEN239 fbdata%ivlqc(ik,jo,jvar) = 4238 IF ( profdata%var(jvar)%nvqc(jk) > 255 ) THEN 239 fbdata%ivlqc(ik,jo,jvar) = IBSET(profdata%var(jvar)%nvqc(jk),2) 240 240 fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk) 241 fbdata%ivlqcf(2,ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) - 10 241 !$AGRIF_DO_NOT_TREAT 242 fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000000011111111') 243 !$AGRIF_END_DO_NOT_TREAT 242 244 ELSE 243 245 fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) … … 320 322 TYPE(obfbdata) :: fbdata 321 323 CHARACTER(LEN=40) :: clfname ! netCDF filename 322 CHARACTER(LEN= 6):: clfiletype324 CHARACTER(LEN=10) :: clfiletype 323 325 CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf' 324 326 INTEGER :: jo … … 395 397 END DO 396 398 397 CASE('ICECON ')399 CASE('ICECONC') 398 400 399 401 CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & … … 418 420 END DO 419 421 422 CASE('SSS') 423 424 CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 425 & 1 + iadd, iext, .TRUE. ) 426 427 clfiletype = 'sssfb' 428 fbdata%cname(1) = surfdata%cvars(1) 429 fbdata%coblong(1) = 'Sea surface salinity' 430 fbdata%cobunit(1) = 'psu' 431 DO je = 1, iext 432 fbdata%cextname(je) = pext%cdname(je) 433 fbdata%cextlong(je) = pext%cdlong(je,1) 434 fbdata%cextunit(je) = pext%cdunit(je,1) 435 END DO 436 fbdata%caddlong(1,1) = 'Model interpolated SSS' 437 fbdata%caddunit(1,1) = 'psu' 438 fbdata%cgrid(1) = 'T' 439 DO ja = 1, iadd 440 fbdata%caddname(1+ja) = padd%cdname(ja) 441 fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 442 fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 443 END DO 444 445 CASE DEFAULT 446 447 CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' ) 448 420 449 END SELECT 421 450 … … 439 468 fbdata%ivqc(jo,:) = 0 440 469 fbdata%ivqcf(:,jo,:) = 0 441 IF ( surfdata%nqc(jo) > 10) THEN470 IF ( surfdata%nqc(jo) > 255 ) THEN 442 471 fbdata%ioqc(jo) = 4 443 472 fbdata%ioqcf(1,jo) = 0 444 fbdata%ioqcf(2,jo) = surfdata%nqc(jo) - 10 473 !$AGRIF_DO_NOT_TREAT 474 fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000000011111111') 475 !$AGRIF_END_DO_NOT_TREAT 445 476 ELSE 446 477 fbdata%ioqc(jo) = surfdata%nqc(jo) … … 474 505 fbdata%idqc(1,jo) = 0 475 506 fbdata%idqcf(:,1,jo) = 0 476 IF ( surfdata%nqc(jo) > 10) THEN507 IF ( surfdata%nqc(jo) > 255 ) THEN 477 508 fbdata%ivqc(jo,1) = 4 478 509 fbdata%ivlqc(1,jo,1) = 4 479 510 fbdata%ivlqcf(1,1,jo,1) = 0 480 fbdata%ivlqcf(2,1,jo,1) = surfdata%nqc(jo) - 10 511 !$AGRIF_DO_NOT_TREAT 512 fbdata%ivlqcf(2,1,jo,1) = IAND(surfdata%nqc(jo),b'0000000011111111') 513 !$AGRIF_END_DO_NOT_TREAT 481 514 ELSE 482 515 fbdata%ivqc(jo,1) = surfdata%nqc(jo) -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_h2d.h90
r2474 r9023 1240 1240 & zdum, & 1241 1241 & zaamax 1242 1242 1243 imax = -1 1243 1244 ! Main computation 1244 1245 pflt = 1.0_wp
Note: See TracChangeset
for help on using the changeset viewer.