- Timestamp:
- 2017-05-02T13:21:57+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
r7960 r7992 7 7 8 8 !!--------------------------------------------------------------------- 9 !! obs_pre_pro : First level check and screening of T/S profiles 10 !! obs_pre_sla : First level check and screening of SLA observations 11 !! obs_pre_sst : First level check and screening of SLA observations 12 !! obs_pre_seaice : First level check and screening of sea ice observations 13 !! obs_pre_vel : First level check and screening of velocity obs. 14 !! obs_scr : Basic screening of the observations 15 !! obs_coo_tim : Compute number of time steps to the observation time 16 !! obs_sor : Sort the observation arrays 9 !! obs_pre_prof : First level check and screening of profile observations 10 !! obs_pre_surf : First level check and screening of surface observations 11 !! obs_scr : Basic screening of the observations 12 !! obs_coo_tim : Compute number of time steps to the observation time 13 !! obs_sor : Sort the observation arrays 17 14 !!--------------------------------------------------------------------- 18 15 !! * Modules used … … 27 24 USE obs_inter_sup ! Interpolation support 28 25 USE obs_oper ! Observation operators 26 #if defined key_bdy 27 USE bdy_oce, ONLY : & ! Boundary information 28 idx_bdy, nb_bdy 29 #endif 29 30 USE lib_mpp, ONLY : & 30 31 & ctl_warn, ctl_stop … … 36 37 37 38 PUBLIC & 38 & obs_pre_pro, & ! First level check and screening of profiles 39 & obs_pre_sla, & ! First level check and screening of SLA data 40 & obs_pre_sst, & ! First level check and screening of SLA data 41 & obs_pre_seaice, & ! First level check and screening of sea ice data 42 & obs_pre_vel, & ! First level check and screening of velocity profiles 43 & calc_month_len ! Calculate the number of days in the months of a year 39 & obs_pre_prof, & ! First level check and screening of profile obs 40 & obs_pre_surf, & ! First level check and screening of surface obs 41 & calc_month_len ! Calculate the number of days in the months of a year 44 42 45 43 !!---------------------------------------------------------------------- … … 49 47 !!---------------------------------------------------------------------- 50 48 49 !! * Substitutions 50 # include "domzgr_substitute.h90" 51 51 52 CONTAINS 52 53 53 SUBROUTINE obs_pre_pro( profdata, prodatqc, ld_t3d, ld_s3d, ld_nea, & 54 & kdailyavtypes ) 55 !!---------------------------------------------------------------------- 56 !! *** ROUTINE obs_pre_pro *** 57 !! 58 !! ** Purpose : First level check and screening of T and S profiles 59 !! 60 !! ** Method : First level check and screening of T and S profiles 61 !! 62 !! ** Action : 63 !! 64 !! References : 65 !! 66 !! History : 67 !! ! 2007-01 (K. Mogensen) Merge of obs_pre_t3d and obs_pre_s3d 68 !! ! 2007-03 (K. Mogensen) General handling of profiles 69 !! ! 2007-06 (K. Mogensen et al) Reject obs. near land. 70 !!---------------------------------------------------------------------- 71 !! * Modules used 72 USE domstp ! Domain: set the time-step 73 USE par_oce ! Ocean parameters 74 USE dom_oce, ONLY : & ! Geographical information 75 & glamt, & 76 & gphit, & 77 & gdept_1d,& 78 & tmask, & 79 & nproc 80 !! * Arguments 81 TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data 82 TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening 83 LOGICAL, INTENT(IN) :: ld_t3d ! Switch for temperature 84 LOGICAL, INTENT(IN) :: ld_s3d ! Switch for salinity 85 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 86 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 87 & kdailyavtypes! Types for daily averages 88 !! * Local declarations 89 INTEGER :: iyea0 ! Initial date 90 INTEGER :: imon0 ! - (year, month, day, hour, minute) 91 INTEGER :: iday0 92 INTEGER :: ihou0 93 INTEGER :: imin0 94 INTEGER :: icycle ! Current assimilation cycle 95 ! Counters for observations that 96 INTEGER :: iotdobs ! - outside time domain 97 INTEGER :: iosdtobs ! - outside space domain (temperature) 98 INTEGER :: iosdsobs ! - outside space domain (salinity) 99 INTEGER :: ilantobs ! - within a model land cell (temperature) 100 INTEGER :: ilansobs ! - within a model land cell (salinity) 101 INTEGER :: inlatobs ! - close to land (temperature) 102 INTEGER :: inlasobs ! - close to land (salinity) 103 INTEGER :: igrdobs ! - fail the grid search 104 ! Global counters for observations that 105 INTEGER :: iotdobsmpp ! - outside time domain 106 INTEGER :: iosdtobsmpp ! - outside space domain (temperature) 107 INTEGER :: iosdsobsmpp ! - outside space domain (salinity) 108 INTEGER :: ilantobsmpp ! - within a model land cell (temperature) 109 INTEGER :: ilansobsmpp ! - within a model land cell (salinity) 110 INTEGER :: inlatobsmpp ! - close to land (temperature) 111 INTEGER :: inlasobsmpp ! - close to land (salinity) 112 INTEGER :: igrdobsmpp ! - fail the grid search 113 TYPE(obs_prof_valid) :: llvalid ! Profile selection 114 TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: & 115 & llvvalid ! T,S selection 116 INTEGER :: jvar ! Variable loop variable 117 INTEGER :: jobs ! Obs. loop variable 118 INTEGER :: jstp ! Time loop variable 119 INTEGER :: inrc ! Time index variable 120 121 IF(lwp) WRITE(numout,*)'obs_pre_pro : Preparing the profile observations...' 122 123 ! Initial date initialization (year, month, day, hour, minute) 124 iyea0 = ndate0 / 10000 125 imon0 = ( ndate0 - iyea0 * 10000 ) / 100 126 iday0 = ndate0 - iyea0 * 10000 - imon0 * 100 127 ihou0 = 0 128 imin0 = 0 129 130 icycle = no ! Assimilation cycle 131 132 ! Diagnotics counters for various failures. 133 134 iotdobs = 0 135 igrdobs = 0 136 iosdtobs = 0 137 iosdsobs = 0 138 ilantobs = 0 139 ilansobs = 0 140 inlatobs = 0 141 inlasobs = 0 142 143 ! ----------------------------------------------------------------------- 144 ! Find time coordinate for profiles 145 ! ----------------------------------------------------------------------- 146 147 IF ( PRESENT(kdailyavtypes) ) THEN 148 CALL obs_coo_tim_prof( icycle, & 149 & iyea0, imon0, iday0, ihou0, imin0, & 150 & profdata%nprof, profdata%nyea, profdata%nmon, & 151 & profdata%nday, profdata%nhou, profdata%nmin, & 152 & profdata%ntyp, profdata%nqc, profdata%mstp, & 153 & iotdobs, kdailyavtypes = kdailyavtypes ) 154 ELSE 155 CALL obs_coo_tim_prof( icycle, & 156 & iyea0, imon0, iday0, ihou0, imin0, & 157 & profdata%nprof, profdata%nyea, profdata%nmon, & 158 & profdata%nday, profdata%nhou, profdata%nmin, & 159 & profdata%ntyp, profdata%nqc, profdata%mstp, & 160 & iotdobs ) 161 ENDIF 162 CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 163 164 ! ----------------------------------------------------------------------- 165 ! Check for profiles failing the grid search 166 ! ----------------------------------------------------------------------- 167 168 CALL obs_coo_grd( profdata%nprof, profdata%mi, profdata%mj, & 169 & profdata%nqc, igrdobs ) 170 171 CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 172 173 ! ----------------------------------------------------------------------- 174 ! Reject all observations for profiles with nqc > 10 175 ! ----------------------------------------------------------------------- 176 177 CALL obs_pro_rej( profdata ) 178 179 ! ----------------------------------------------------------------------- 180 ! Check for land points. This includes points below the model 181 ! bathymetry so this is done for every point in the profile 182 ! ----------------------------------------------------------------------- 183 184 ! Temperature 185 186 CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(1), & 187 & profdata%npvsta(:,1), profdata%npvend(:,1), & 188 & jpi, jpj, & 189 & jpk, & 190 & profdata%mi, profdata%mj, & 191 & profdata%var(1)%mvk, & 192 & profdata%rlam, profdata%rphi, & 193 & profdata%var(1)%vdep, & 194 & glamt, gphit, & 195 & gdept_1d, tmask, & 196 & profdata%nqc, profdata%var(1)%nvqc, & 197 & iosdtobs, ilantobs, & 198 & inlatobs, ld_nea ) 199 200 CALL obs_mpp_sum_integer( iosdtobs, iosdtobsmpp ) 201 CALL obs_mpp_sum_integer( ilantobs, ilantobsmpp ) 202 CALL obs_mpp_sum_integer( inlatobs, inlatobsmpp ) 203 204 ! Salinity 205 206 CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(2), & 207 & profdata%npvsta(:,2), profdata%npvend(:,2), & 208 & jpi, jpj, & 209 & jpk, & 210 & profdata%mi, profdata%mj, & 211 & profdata%var(2)%mvk, & 212 & profdata%rlam, profdata%rphi, & 213 & profdata%var(2)%vdep, & 214 & glamt, gphit, & 215 & gdept_1d, tmask, & 216 & profdata%nqc, profdata%var(2)%nvqc, & 217 & iosdsobs, ilansobs, & 218 & inlasobs, ld_nea ) 219 220 CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 221 CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 222 CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 223 224 ! ----------------------------------------------------------------------- 225 ! Copy useful data from the profdata data structure to 226 ! the prodatqc data structure 227 ! ----------------------------------------------------------------------- 228 229 ! Allocate the selection arrays 230 231 ALLOCATE( llvalid%luse(profdata%nprof) ) 232 DO jvar = 1,profdata%nvar 233 ALLOCATE( llvvalid(jvar)%luse(profdata%nvprot(jvar)) ) 234 END DO 235 236 ! We want all data which has qc flags <= 10 237 238 llvalid%luse(:) = ( profdata%nqc(:) <= 10 ) 239 DO jvar = 1,profdata%nvar 240 llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 ) 241 END DO 242 243 ! The actual copying 244 245 CALL obs_prof_compress( profdata, prodatqc, .TRUE., numout, & 246 & lvalid=llvalid, lvvalid=llvvalid ) 247 248 ! Dellocate the selection arrays 249 DEALLOCATE( llvalid%luse ) 250 DO jvar = 1,profdata%nvar 251 DEALLOCATE( llvvalid(jvar)%luse ) 252 END DO 253 254 ! ----------------------------------------------------------------------- 255 ! Print information about what observations are left after qc 256 ! ----------------------------------------------------------------------- 257 258 ! Update the total observation counter array 259 260 IF(lwp) THEN 261 WRITE(numout,*) 262 WRITE(numout,*) 'obs_pre_pro :' 263 WRITE(numout,*) '~~~~~~~~~~~' 264 WRITE(numout,*) 265 WRITE(numout,*) ' Profiles outside time domain = ', & 266 & iotdobsmpp 267 WRITE(numout,*) ' Remaining profiles that failed grid search = ', & 268 & igrdobsmpp 269 WRITE(numout,*) ' Remaining T data outside space domain = ', & 270 & iosdtobsmpp 271 WRITE(numout,*) ' Remaining T data at land points = ', & 272 & ilantobsmpp 273 IF (ld_nea) THEN 274 WRITE(numout,*) ' Remaining T data near land points (removed) = ',& 275 & inlatobsmpp 276 ELSE 277 WRITE(numout,*) ' Remaining T data near land points (kept) = ',& 278 & inlatobsmpp 279 ENDIF 280 WRITE(numout,*) ' T data accepted = ', & 281 & prodatqc%nvprotmpp(1) 282 WRITE(numout,*) ' Remaining S data outside space domain = ', & 283 & iosdsobsmpp 284 WRITE(numout,*) ' Remaining S data at land points = ', & 285 & ilansobsmpp 286 IF (ld_nea) THEN 287 WRITE(numout,*) ' Remaining S data near land points (removed) = ',& 288 & inlasobsmpp 289 ELSE 290 WRITE(numout,*) ' Remaining S data near land points (kept) = ',& 291 & inlasobsmpp 292 ENDIF 293 WRITE(numout,*) ' S data accepted = ', & 294 & prodatqc%nvprotmpp(2) 295 296 WRITE(numout,*) 297 WRITE(numout,*) ' Number of observations per time step :' 298 WRITE(numout,*) 299 WRITE(numout,997) 300 WRITE(numout,998) 301 ENDIF 302 303 DO jobs = 1, prodatqc%nprof 304 inrc = prodatqc%mstp(jobs) + 2 - nit000 305 prodatqc%npstp(inrc) = prodatqc%npstp(inrc) + 1 306 DO jvar = 1, prodatqc%nvar 307 IF ( prodatqc%npvend(jobs,jvar) > 0 ) THEN 308 prodatqc%nvstp(inrc,jvar) = prodatqc%nvstp(inrc,jvar) + & 309 & ( prodatqc%npvend(jobs,jvar) - & 310 & prodatqc%npvsta(jobs,jvar) + 1 ) 311 ENDIF 312 END DO 313 END DO 314 315 316 CALL obs_mpp_sum_integers( prodatqc%npstp, prodatqc%npstpmpp, & 317 & nitend - nit000 + 2 ) 318 DO jvar = 1, prodatqc%nvar 319 CALL obs_mpp_sum_integers( prodatqc%nvstp(:,jvar), & 320 & prodatqc%nvstpmpp(:,jvar), & 321 & nitend - nit000 + 2 ) 322 END DO 323 324 IF ( lwp ) THEN 325 DO jstp = nit000 - 1, nitend 326 inrc = jstp - nit000 + 2 327 WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), & 328 & prodatqc%nvstpmpp(inrc,1), & 329 & prodatqc%nvstpmpp(inrc,2) 330 END DO 331 ENDIF 332 333 997 FORMAT(10X,'Time step',5X,'Profiles',5X,'Temperature',5X,'Salinity') 334 998 FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'--------') 335 999 FORMAT(10X,I9,5X,I8,5X,I11,5X,I8) 336 337 END SUBROUTINE obs_pre_pro 338 339 SUBROUTINE obs_pre_sla( sladata, sladatqc, ld_sla, ld_nea ) 54 SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, & 55 kqc_cutoff ) 340 56 !!---------------------------------------------------------------------- 341 57 !! *** ROUTINE obs_pre_sla *** 342 58 !! 343 !! ** Purpose : First level check and screening of SLAobservations344 !! 345 !! ** Method : First level check and screening of SLAobservations59 !! ** Purpose : First level check and screening of surface observations 60 !! 61 !! ** Method : First level check and screening of surface observations 346 62 !! 347 63 !! ** Action : … … 352 68 !! ! 2007-03 (A. Weaver, K. Mogensen) Original 353 69 !! ! 2007-06 (K. Mogensen et al) Reject obs. near land. 70 !! ! 2015-02 (M. Martin) Combined routine for surface types. 354 71 !!---------------------------------------------------------------------- 355 72 !! * Modules used … … 362 79 & nproc 363 80 !! * Arguments 364 TYPE(obs_surf), INTENT(INOUT) :: sladata ! Full set of SLA data 365 TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of SLA data not failing screening 366 LOGICAL, INTENT(IN) :: ld_sla ! Switch for SLA data 367 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 81 TYPE(obs_surf), INTENT(INOUT) :: surfdata ! Full set of surface data 82 TYPE(obs_surf), INTENT(INOUT) :: surfdataqc ! Subset of surface data not failing screening 83 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 84 LOGICAL, INTENT(IN) :: ld_bound_reject ! Switch for rejecting obs near the boundary 85 INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value 368 86 !! * Local declarations 87 INTEGER :: iqc_cutoff = 255 ! cut off for QC value 369 88 INTEGER :: iyea0 ! Initial date 370 89 INTEGER :: imon0 ! - (year, month, day, hour, minute) … … 379 98 INTEGER :: inlasobs ! - close to land 380 99 INTEGER :: igrdobs ! - fail the grid search 100 INTEGER :: ibdysobs ! - close to open boundary 381 101 ! Global counters for observations that 382 102 INTEGER :: iotdobsmpp ! - outside time domain … … 385 105 INTEGER :: inlasobsmpp ! - close to land 386 106 INTEGER :: igrdobsmpp ! - fail the grid search 107 INTEGER :: ibdysobsmpp ! - close to open boundary 387 108 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 388 109 & llvalid ! SLA data selection … … 391 112 INTEGER :: inrc ! Time index variable 392 113 393 IF(lwp) WRITE(numout,*)'obs_pre_sla : Preparing the SLA observations...' 394 114 IF(lwp) WRITE(numout,*)'obs_pre_surf : Preparing the surface observations...' 115 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 116 395 117 ! Initial date initialization (year, month, day, hour, minute) 396 118 iyea0 = ndate0 / 10000 … … 409 131 ilansobs = 0 410 132 inlasobs = 0 411 412 ! ----------------------------------------------------------------------- 413 ! Find time coordinate for SLA data 133 ibdysobs = 0 134 135 ! Set QC cutoff to optional value if provided 136 IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 137 138 ! ----------------------------------------------------------------------- 139 ! Find time coordinate for surface data 414 140 ! ----------------------------------------------------------------------- 415 141 416 142 CALL obs_coo_tim( icycle, & 417 143 & iyea0, imon0, iday0, ihou0, imin0, & 418 & s ladata%nsurf, sladata%nyea, sladata%nmon, &419 & s ladata%nday, sladata%nhou, sladata%nmin, &420 & s ladata%nqc, sladata%mstp, iotdobs )144 & surfdata%nsurf, surfdata%nyea, surfdata%nmon, & 145 & surfdata%nday, surfdata%nhou, surfdata%nmin, & 146 & surfdata%nqc, surfdata%mstp, iotdobs ) 421 147 422 148 CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 423 149 424 150 ! ----------------------------------------------------------------------- 425 ! Check for SLAdata failing the grid search426 ! ----------------------------------------------------------------------- 427 428 CALL obs_coo_grd( s ladata%nsurf, sladata%mi, sladata%mj, &429 & s ladata%nqc, igrdobs )151 ! Check for surface data failing the grid search 152 ! ----------------------------------------------------------------------- 153 154 CALL obs_coo_grd( surfdata%nsurf, surfdata%mi, surfdata%mj, & 155 & surfdata%nqc, igrdobs ) 430 156 431 157 CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) … … 435 161 ! ----------------------------------------------------------------------- 436 162 437 CALL obs_coo_spc_2d( s ladata%nsurf, &163 CALL obs_coo_spc_2d( surfdata%nsurf, & 438 164 & jpi, jpj, & 439 & s ladata%mi, sladata%mj, &440 & s ladata%rlam, sladata%rphi, &165 & surfdata%mi, surfdata%mj, & 166 & surfdata%rlam, surfdata%rphi, & 441 167 & glamt, gphit, & 442 & tmask(:,:,1), s ladata%nqc, &168 & tmask(:,:,1), surfdata%nqc, & 443 169 & iosdsobs, ilansobs, & 444 & inlasobs, ld_nea ) 170 & inlasobs, ld_nea, & 171 & ibdysobs, ld_bound_reject, & 172 & iqc_cutoff ) 445 173 446 174 CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 447 175 CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 448 176 CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 449 450 ! ----------------------------------------------------------------------- 451 ! Copy useful data from the sladata data structure to 452 ! the sladatqc data structure 177 CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 178 179 ! ----------------------------------------------------------------------- 180 ! Copy useful data from the surfdata data structure to 181 ! the surfdataqc data structure 453 182 ! ----------------------------------------------------------------------- 454 183 455 184 ! Allocate the selection arrays 456 185 457 ALLOCATE( llvalid(s ladata%nsurf) )458 459 ! We want all data which has qc flags <= 10460 461 llvalid(:) = ( s ladata%nqc(:) <= 10)186 ALLOCATE( llvalid(surfdata%nsurf) ) 187 188 ! We want all data which has qc flags <= iqc_cutoff 189 190 llvalid(:) = ( surfdata%nqc(:) <= iqc_cutoff ) 462 191 463 192 ! The actual copying 464 193 465 CALL obs_surf_compress( s ladata, sladatqc, .TRUE., numout, &194 CALL obs_surf_compress( surfdata, surfdataqc, .TRUE., numout, & 466 195 & lvalid=llvalid ) 467 196 … … 477 206 IF(lwp) THEN 478 207 WRITE(numout,*) 479 WRITE(numout,*) 'obs_pre_sla :' 480 WRITE(numout,*) '~~~~~~~~~~~' 481 WRITE(numout,*) 482 WRITE(numout,*) ' SLA data outside time domain = ', & 208 WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data outside time domain = ', & 483 209 & iotdobsmpp 484 WRITE(numout,*) ' Remaining SLAdata that failed grid search = ', &210 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data that failed grid search = ', & 485 211 & igrdobsmpp 486 WRITE(numout,*) ' Remaining SLAdata outside space domain = ', &212 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data outside space domain = ', & 487 213 & iosdsobsmpp 488 WRITE(numout,*) ' Remaining SLAdata at land points = ', &214 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data at land points = ', & 489 215 & ilansobsmpp 490 216 IF (ld_nea) THEN 491 WRITE(numout,*) ' Remaining SLAdata near land points (removed) = ', &217 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (removed) = ', & 492 218 & inlasobsmpp 493 219 ELSE 494 WRITE(numout,*) ' Remaining SLAdata near land points (kept) = ', &220 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (kept) = ', & 495 221 & inlasobsmpp 496 222 ENDIF 497 WRITE(numout,*) ' SLA data accepted = ', & 498 & sladatqc%nsurfmpp 223 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', & 224 & ibdysobsmpp 225 WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted = ', & 226 & surfdataqc%nsurfmpp 499 227 500 228 WRITE(numout,*) 501 229 WRITE(numout,*) ' Number of observations per time step :' 502 230 WRITE(numout,*) 503 WRITE(numout,1997) 504 WRITE(numout,1998) 231 WRITE(numout,'(10X,A,10X,A)')'Time step',surfdataqc%cvars(1) 232 WRITE(numout,'(10X,A,5X,A)')'---------','-----------------' 233 CALL FLUSH(numout) 505 234 ENDIF 506 235 507 DO jobs = 1, s ladatqc%nsurf508 inrc = s ladatqc%mstp(jobs) + 2 - nit000509 s ladatqc%nsstp(inrc) = sladatqc%nsstp(inrc) + 1236 DO jobs = 1, surfdataqc%nsurf 237 inrc = surfdataqc%mstp(jobs) + 2 - nit000 238 surfdataqc%nsstp(inrc) = surfdataqc%nsstp(inrc) + 1 510 239 END DO 511 240 512 CALL obs_mpp_sum_integers( s ladatqc%nsstp, sladatqc%nsstpmpp, &241 CALL obs_mpp_sum_integers( surfdataqc%nsstp, surfdataqc%nsstpmpp, & 513 242 & nitend - nit000 + 2 ) 514 243 … … 516 245 DO jstp = nit000 - 1, nitend 517 246 inrc = jstp - nit000 + 2 518 WRITE(numout,1999) jstp, sladatqc%nsstpmpp(inrc) 247 WRITE(numout,1999) jstp, surfdataqc%nsstpmpp(inrc) 248 CALL FLUSH(numout) 519 249 END DO 520 250 ENDIF 521 251 522 1997 FORMAT(10X,'Time step',5X,'Sea level anomaly')523 1998 FORMAT(10X,'---------',5X,'-----------------')524 252 1999 FORMAT(10X,I9,5X,I17) 525 253 526 END SUBROUTINE obs_pre_sla 527 528 SUBROUTINE obs_pre_sst( sstdata, sstdatqc, ld_sst, ld_nea ) 529 !!---------------------------------------------------------------------- 530 !! *** ROUTINE obs_pre_sst *** 531 !! 532 !! ** Purpose : First level check and screening of SST observations 533 !! 534 !! ** Method : First level check and screening of SST observations 535 !! 536 !! ** Action : 537 !! 538 !! References : 539 !! 254 END SUBROUTINE obs_pre_surf 255 256 257 SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_var1, ld_var2, & 258 & kpi, kpj, kpk, & 259 & zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2, & 260 & ld_nea, ld_bound_reject, kdailyavtypes, kqc_cutoff ) 261 262 !!---------------------------------------------------------------------- 263 !! *** ROUTINE obs_pre_prof *** 264 !! 265 !! ** Purpose : First level check and screening of profiles 266 !! 267 !! ** Method : First level check and screening of profiles 268 !! 540 269 !! History : 541 !! ! 2007-03 (S. Ricci) SST data preparation 270 !! ! 2007-06 (K. Mogensen) original : T and S profile data 271 !! ! 2008-09 (M. Valdivieso) : TAO velocity data 272 !! ! 2009-01 (K. Mogensen) : New feedback stricture 273 !! ! 2015-02 (M. Martin) : Combined profile routine. 274 !! 542 275 !!---------------------------------------------------------------------- 543 276 !! * Modules used … … 545 278 USE par_oce ! Ocean parameters 546 279 USE dom_oce, ONLY : & ! Geographical information 547 & glamt, & 548 & gphit, & 549 & tmask, & 280 & gdept_1d, & 550 281 & nproc 282 551 283 !! * Arguments 552 TYPE(obs_surf), INTENT(INOUT) :: sstdata ! Full set of SST data 553 TYPE(obs_surf), INTENT(INOUT) :: sstdatqc ! Subset of SST data not failing screening 554 LOGICAL :: ld_sst ! Switch for SST data 555 LOGICAL :: ld_nea ! Switch for rejecting observation near land 284 TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data 285 TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening 286 LOGICAL, INTENT(IN) :: ld_var1 ! Observed variables switches 287 LOGICAL, INTENT(IN) :: ld_var2 288 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 289 LOGICAL, INTENT(IN) :: ld_bound_reject ! Switch for rejecting observations near the boundary 290 INTEGER, INTENT(IN) :: kpi, kpj, kpk ! Local domain sizes 291 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 292 & kdailyavtypes ! Types for daily averages 293 REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 294 & zmask1, & 295 & zmask2 296 REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 297 & pglam1, & 298 & pglam2, & 299 & pgphi1, & 300 & pgphi2 301 INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value 302 556 303 !! * Local declarations 304 INTEGER :: iqc_cutoff = 255 ! cut off for QC value 557 305 INTEGER :: iyea0 ! Initial date 558 306 INTEGER :: imon0 ! - (year, month, day, hour, minute) 559 INTEGER :: iday0 307 INTEGER :: iday0 560 308 INTEGER :: ihou0 561 309 INTEGER :: imin0 562 310 INTEGER :: icycle ! Current assimilation cycle 563 ! Counters for observations that 311 ! Counters for observations that are 564 312 INTEGER :: iotdobs ! - outside time domain 565 INTEGER :: iosdsobs ! - outside space domain 566 INTEGER :: ilansobs ! - within a model land cell 567 INTEGER :: inlasobs ! - close to land 313 INTEGER :: iosdv1obs ! - outside space domain (variable 1) 314 INTEGER :: iosdv2obs ! - outside space domain (variable 2) 315 INTEGER :: ilanv1obs ! - within a model land cell (variable 1) 316 INTEGER :: ilanv2obs ! - within a model land cell (variable 2) 317 INTEGER :: inlav1obs ! - close to land (variable 1) 318 INTEGER :: inlav2obs ! - close to land (variable 2) 319 INTEGER :: ibdyv1obs ! - boundary (variable 1) 320 INTEGER :: ibdyv2obs ! - boundary (variable 2) 568 321 INTEGER :: igrdobs ! - fail the grid search 569 ! Global counters for observations that 322 INTEGER :: iuvchku ! - reject u if v rejected and vice versa 323 INTEGER :: iuvchkv ! 324 ! Global counters for observations that are 570 325 INTEGER :: iotdobsmpp ! - outside time domain 571 INTEGER :: iosdsobsmpp ! - outside space domain 572 INTEGER :: ilansobsmpp ! - within a model land cell 573 INTEGER :: inlasobsmpp ! - close to land 326 INTEGER :: iosdv1obsmpp ! - outside space domain (variable 1) 327 INTEGER :: iosdv2obsmpp ! - outside space domain (variable 2) 328 INTEGER :: ilanv1obsmpp ! - within a model land cell (variable 1) 329 INTEGER :: ilanv2obsmpp ! - within a model land cell (variable 2) 330 INTEGER :: inlav1obsmpp ! - close to land (variable 1) 331 INTEGER :: inlav2obsmpp ! - close to land (variable 2) 332 INTEGER :: ibdyv1obsmpp ! - boundary (variable 1) 333 INTEGER :: ibdyv2obsmpp ! - boundary (variable 2) 574 334 INTEGER :: igrdobsmpp ! - fail the grid search 575 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 576 & llvalid ! SST data selection 335 INTEGER :: iuvchkumpp ! - reject var1 if var2 rejected and vice versa 336 INTEGER :: iuvchkvmpp ! 337 TYPE(obs_prof_valid) :: llvalid ! Profile selection 338 TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: & 339 & llvvalid ! var1,var2 selection 340 INTEGER :: jvar ! Variable loop variable 577 341 INTEGER :: jobs ! Obs. loop variable 578 342 INTEGER :: jstp ! Time loop variable 579 343 INTEGER :: inrc ! Time index variable 580 344 581 IF(lwp) WRITE(numout,*)'obs_pre_sst : Preparing the SST observations...' 345 IF(lwp) WRITE(numout,*)'obs_pre_prof: Preparing the profile data...' 346 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 582 347 583 348 ! Initial date initialization (year, month, day, hour, minute) … … 592 357 ! Diagnotics counters for various failures. 593 358 594 iotdobs = 0 595 igrdobs = 0 596 iosdsobs = 0 597 ilansobs = 0 598 inlasobs = 0 599 600 ! ----------------------------------------------------------------------- 601 ! Find time coordinate for SST data 602 ! ----------------------------------------------------------------------- 603 604 CALL obs_coo_tim( icycle, & 605 & iyea0, imon0, iday0, ihou0, imin0, & 606 & sstdata%nsurf, sstdata%nyea, sstdata%nmon, & 607 & sstdata%nday, sstdata%nhou, sstdata%nmin, & 608 & sstdata%nqc, sstdata%mstp, iotdobs ) 609 CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 610 ! ----------------------------------------------------------------------- 611 ! Check for SST data failing the grid search 612 ! ----------------------------------------------------------------------- 613 614 CALL obs_coo_grd( sstdata%nsurf, sstdata%mi, sstdata%mj, & 615 & sstdata%nqc, igrdobs ) 616 CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 617 618 ! ----------------------------------------------------------------------- 619 ! Check for land points. 620 ! ----------------------------------------------------------------------- 621 622 CALL obs_coo_spc_2d( sstdata%nsurf, & 623 & jpi, jpj, & 624 & sstdata%mi, sstdata%mj, & 625 & sstdata%rlam, sstdata%rphi, & 626 & glamt, gphit, & 627 & tmask(:,:,1), sstdata%nqc, & 628 & iosdsobs, ilansobs, & 629 & inlasobs, ld_nea ) 630 631 CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 632 CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 633 CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 634 635 ! ----------------------------------------------------------------------- 636 ! Copy useful data from the sstdata data structure to 637 ! the sstdatqc data structure 638 ! ----------------------------------------------------------------------- 639 640 ! Allocate the selection arrays 641 642 ALLOCATE( llvalid(sstdata%nsurf) ) 643 644 ! We want all data which has qc flags <= 0 645 646 llvalid(:) = ( sstdata%nqc(:) <= 10 ) 647 648 ! The actual copying 649 650 CALL obs_surf_compress( sstdata, sstdatqc, .TRUE., numout, & 651 & lvalid=llvalid ) 652 653 ! Dellocate the selection arrays 654 DEALLOCATE( llvalid ) 655 656 ! ----------------------------------------------------------------------- 657 ! Print information about what observations are left after qc 658 ! ----------------------------------------------------------------------- 659 660 ! Update the total observation counter array 661 662 IF(lwp) THEN 663 WRITE(numout,*) 664 WRITE(numout,*) 'obs_pre_sst :' 665 WRITE(numout,*) '~~~~~~~~~~~' 666 WRITE(numout,*) 667 WRITE(numout,*) ' SST data outside time domain = ', & 668 & iotdobsmpp 669 WRITE(numout,*) ' Remaining SST data that failed grid search = ', & 670 & igrdobsmpp 671 WRITE(numout,*) ' Remaining SST data outside space domain = ', & 672 & iosdsobsmpp 673 WRITE(numout,*) ' Remaining SST data at land points = ', & 674 & ilansobsmpp 675 IF (ld_nea) THEN 676 WRITE(numout,*) ' Remaining SST data near land points (removed) = ', & 677 & inlasobsmpp 678 ELSE 679 WRITE(numout,*) ' Remaining SST data near land points (kept) = ', & 680 & inlasobsmpp 681 ENDIF 682 WRITE(numout,*) ' SST data accepted = ', & 683 & sstdatqc%nsurfmpp 684 685 WRITE(numout,*) 686 WRITE(numout,*) ' Number of observations per time step :' 687 WRITE(numout,*) 688 WRITE(numout,1997) 689 WRITE(numout,1998) 359 iotdobs = 0 360 igrdobs = 0 361 iosdv1obs = 0 362 iosdv2obs = 0 363 ilanv1obs = 0 364 ilanv2obs = 0 365 inlav1obs = 0 366 inlav2obs = 0 367 ibdyv1obs = 0 368 ibdyv2obs = 0 369 iuvchku = 0 370 iuvchkv = 0 371 372 373 ! Set QC cutoff to optional value if provided 374 IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 375 376 ! ----------------------------------------------------------------------- 377 ! Find time coordinate for profiles 378 ! ----------------------------------------------------------------------- 379 380 IF ( PRESENT(kdailyavtypes) ) THEN 381 CALL obs_coo_tim_prof( icycle, & 382 & iyea0, imon0, iday0, ihou0, imin0, & 383 & profdata%nprof, profdata%nyea, profdata%nmon, & 384 & profdata%nday, profdata%nhou, profdata%nmin, & 385 & profdata%ntyp, profdata%nqc, profdata%mstp, & 386 & iotdobs, kdailyavtypes = kdailyavtypes, & 387 & kqc_cutoff = iqc_cutoff ) 388 ELSE 389 CALL obs_coo_tim_prof( icycle, & 390 & iyea0, imon0, iday0, ihou0, imin0, & 391 & profdata%nprof, profdata%nyea, profdata%nmon, & 392 & profdata%nday, profdata%nhou, profdata%nmin, & 393 & profdata%ntyp, profdata%nqc, profdata%mstp, & 394 & iotdobs, kqc_cutoff = iqc_cutoff ) 690 395 ENDIF 691 692 DO jobs = 1, sstdatqc%nsurf 693 inrc = sstdatqc%mstp(jobs) + 2 - nit000 694 sstdatqc%nsstp(inrc) = sstdatqc%nsstp(inrc) + 1 695 END DO 696 697 CALL obs_mpp_sum_integers( sstdatqc%nsstp, sstdatqc%nsstpmpp, & 698 & nitend - nit000 + 2 ) 699 700 IF ( lwp ) THEN 701 DO jstp = nit000 - 1, nitend 702 inrc = jstp - nit000 + 2 703 WRITE(numout,1999) jstp, sstdatqc%nsstpmpp(inrc) 704 END DO 705 ENDIF 706 707 1997 FORMAT(10X,'Time step',5X,'Sea surface temperature') 708 1998 FORMAT(10X,'---------',5X,'-----------------') 709 1999 FORMAT(10X,I9,5X,I17) 710 711 END SUBROUTINE obs_pre_sst 712 713 SUBROUTINE obs_pre_seaice( seaicedata, seaicedatqc, ld_seaice, ld_nea ) 714 !!---------------------------------------------------------------------- 715 !! *** ROUTINE obs_pre_seaice *** 716 !! 717 !! ** Purpose : First level check and screening of Sea Ice observations 718 !! 719 !! ** Method : First level check and screening of Sea Ice observations 720 !! 721 !! ** Action : 722 !! 723 !! References : 724 !! 725 !! History : 726 !! ! 2007-11 (D. Lea) based on obs_pre_sst 727 !!---------------------------------------------------------------------- 728 !! * Modules used 729 USE domstp ! Domain: set the time-step 730 USE par_oce ! Ocean parameters 731 USE dom_oce, ONLY : & ! Geographical information 732 & glamt, & 733 & gphit, & 734 & tmask, & 735 & nproc 736 !! * Arguments 737 TYPE(obs_surf), INTENT(INOUT) :: seaicedata ! Full set of Sea Ice data 738 TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc ! Subset of sea ice data not failing screening 739 LOGICAL :: ld_seaice ! Switch for sea ice data 740 LOGICAL :: ld_nea ! Switch for rejecting observation near land 741 !! * Local declarations 742 INTEGER :: iyea0 ! Initial date 743 INTEGER :: imon0 ! - (year, month, day, hour, minute) 744 INTEGER :: iday0 745 INTEGER :: ihou0 746 INTEGER :: imin0 747 INTEGER :: icycle ! Current assimilation cycle 748 ! Counters for observations that 749 INTEGER :: iotdobs ! - outside time domain 750 INTEGER :: iosdsobs ! - outside space domain 751 INTEGER :: ilansobs ! - within a model land cell 752 INTEGER :: inlasobs ! - close to land 753 INTEGER :: igrdobs ! - fail the grid search 754 ! Global counters for observations that 755 INTEGER :: iotdobsmpp ! - outside time domain 756 INTEGER :: iosdsobsmpp ! - outside space domain 757 INTEGER :: ilansobsmpp ! - within a model land cell 758 INTEGER :: inlasobsmpp ! - close to land 759 INTEGER :: igrdobsmpp ! - fail the grid search 760 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 761 & llvalid ! data selection 762 INTEGER :: jobs ! Obs. loop variable 763 INTEGER :: jstp ! Time loop variable 764 INTEGER :: inrc ! Time index variable 765 766 IF (lwp) WRITE(numout,*)'obs_pre_seaice : Preparing the sea ice observations...' 767 768 ! Initial date initialization (year, month, day, hour, minute) 769 iyea0 = ndate0 / 10000 770 imon0 = ( ndate0 - iyea0 * 10000 ) / 100 771 iday0 = ndate0 - iyea0 * 10000 - imon0 * 100 772 ihou0 = 0 773 imin0 = 0 774 775 icycle = no ! Assimilation cycle 776 777 ! Diagnotics counters for various failures. 778 779 iotdobs = 0 780 igrdobs = 0 781 iosdsobs = 0 782 ilansobs = 0 783 inlasobs = 0 784 785 ! ----------------------------------------------------------------------- 786 ! Find time coordinate for sea ice data 787 ! ----------------------------------------------------------------------- 788 789 CALL obs_coo_tim( icycle, & 790 & iyea0, imon0, iday0, ihou0, imin0, & 791 & seaicedata%nsurf, seaicedata%nyea, seaicedata%nmon, & 792 & seaicedata%nday, seaicedata%nhou, seaicedata%nmin, & 793 & seaicedata%nqc, seaicedata%mstp, iotdobs ) 794 CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 795 ! ----------------------------------------------------------------------- 796 ! Check for sea ice data failing the grid search 797 ! ----------------------------------------------------------------------- 798 799 CALL obs_coo_grd( seaicedata%nsurf, seaicedata%mi, seaicedata%mj, & 800 & seaicedata%nqc, igrdobs ) 801 CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 802 803 ! ----------------------------------------------------------------------- 804 ! Check for land points. 805 ! ----------------------------------------------------------------------- 806 807 CALL obs_coo_spc_2d( seaicedata%nsurf, & 808 & jpi, jpj, & 809 & seaicedata%mi, seaicedata%mj, & 810 & seaicedata%rlam, seaicedata%rphi, & 811 & glamt, gphit, & 812 & tmask(:,:,1), seaicedata%nqc, & 813 & iosdsobs, ilansobs, & 814 & inlasobs, ld_nea ) 815 816 CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 817 CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 818 CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 819 820 ! ----------------------------------------------------------------------- 821 ! Copy useful data from the seaicedata data structure to 822 ! the seaicedatqc data structure 823 ! ----------------------------------------------------------------------- 824 825 ! Allocate the selection arrays 826 827 ALLOCATE( llvalid(seaicedata%nsurf) ) 828 829 ! We want all data which has qc flags <= 0 830 831 llvalid(:) = ( seaicedata%nqc(:) <= 10 ) 832 833 ! The actual copying 834 835 CALL obs_surf_compress( seaicedata, seaicedatqc, .TRUE., numout, & 836 & lvalid=llvalid ) 837 838 ! Dellocate the selection arrays 839 DEALLOCATE( llvalid ) 840 841 ! ----------------------------------------------------------------------- 842 ! Print information about what observations are left after qc 843 ! ----------------------------------------------------------------------- 844 845 ! Update the total observation counter array 846 847 IF(lwp) THEN 848 WRITE(numout,*) 849 WRITE(numout,*) 'obs_pre_seaice :' 850 WRITE(numout,*) '~~~~~~~~~~~' 851 WRITE(numout,*) 852 WRITE(numout,*) ' Sea ice data outside time domain = ', & 853 & iotdobsmpp 854 WRITE(numout,*) ' Remaining sea ice data that failed grid search = ', & 855 & igrdobsmpp 856 WRITE(numout,*) ' Remaining sea ice data outside space domain = ', & 857 & iosdsobsmpp 858 WRITE(numout,*) ' Remaining sea ice data at land points = ', & 859 & ilansobsmpp 860 IF (ld_nea) THEN 861 WRITE(numout,*) ' Remaining sea ice data near land points (removed) = ', & 862 & inlasobsmpp 863 ELSE 864 WRITE(numout,*) ' Remaining sea ice data near land points (kept) = ', & 865 & inlasobsmpp 866 ENDIF 867 WRITE(numout,*) ' Sea ice data accepted = ', & 868 & seaicedatqc%nsurfmpp 869 870 WRITE(numout,*) 871 WRITE(numout,*) ' Number of observations per time step :' 872 WRITE(numout,*) 873 WRITE(numout,1997) 874 WRITE(numout,1998) 875 ENDIF 876 877 DO jobs = 1, seaicedatqc%nsurf 878 inrc = seaicedatqc%mstp(jobs) + 2 - nit000 879 seaicedatqc%nsstp(inrc) = seaicedatqc%nsstp(inrc) + 1 880 END DO 881 882 CALL obs_mpp_sum_integers( seaicedatqc%nsstp, seaicedatqc%nsstpmpp, & 883 & nitend - nit000 + 2 ) 884 885 IF ( lwp ) THEN 886 DO jstp = nit000 - 1, nitend 887 inrc = jstp - nit000 + 2 888 WRITE(numout,1999) jstp, seaicedatqc%nsstpmpp(inrc) 889 END DO 890 ENDIF 891 892 1997 FORMAT(10X,'Time step',5X,'Sea ice data ') 893 1998 FORMAT(10X,'---------',5X,'-----------------') 894 1999 FORMAT(10X,I9,5X,I17) 895 896 END SUBROUTINE obs_pre_seaice 897 898 SUBROUTINE obs_pre_vel( profdata, prodatqc, ld_vel3d, ld_nea, ld_dailyav ) 899 !!---------------------------------------------------------------------- 900 !! *** ROUTINE obs_pre_taovel *** 901 !! 902 !! ** Purpose : First level check and screening of U and V profiles 903 !! 904 !! ** Method : First level check and screening of U and V profiles 905 !! 906 !! History : 907 !! ! 2007-06 (K. Mogensen) original : T and S profile data 908 !! ! 2008-09 (M. Valdivieso) : TAO velocity data 909 !! ! 2009-01 (K. Mogensen) : New feedback strictuer 910 !! 911 !!---------------------------------------------------------------------- 912 !! * Modules used 913 USE domstp ! Domain: set the time-step 914 USE par_oce ! Ocean parameters 915 USE dom_oce, ONLY : & ! Geographical information 916 & glamt, glamu, glamv, & 917 & gphit, gphiu, gphiv, & 918 & gdept_1d, & 919 & tmask, umask, vmask, & 920 & nproc 921 !! * Arguments 922 TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data 923 TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening 924 LOGICAL, INTENT(IN) :: ld_vel3d ! Switch for zonal and meridional velocity components 925 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 926 LOGICAL, INTENT(IN) :: ld_dailyav ! Switch for daily average data 927 !! * Local declarations 928 INTEGER :: iyea0 ! Initial date 929 INTEGER :: imon0 ! - (year, month, day, hour, minute) 930 INTEGER :: iday0 931 INTEGER :: ihou0 932 INTEGER :: imin0 933 INTEGER :: icycle ! Current assimilation cycle 934 ! Counters for observations that 935 INTEGER :: iotdobs ! - outside time domain 936 INTEGER :: iosduobs ! - outside space domain (zonal velocity component) 937 INTEGER :: iosdvobs ! - outside space domain (meridional velocity component) 938 INTEGER :: ilanuobs ! - within a model land cell (zonal velocity component) 939 INTEGER :: ilanvobs ! - within a model land cell (meridional velocity component) 940 INTEGER :: inlauobs ! - close to land (zonal velocity component) 941 INTEGER :: inlavobs ! - close to land (meridional velocity component) 942 INTEGER :: igrdobs ! - fail the grid search 943 INTEGER :: iuvchku ! - reject u if v rejected and vice versa 944 INTEGER :: iuvchkv ! 945 ! Global counters for observations that 946 INTEGER :: iotdobsmpp ! - outside time domain 947 INTEGER :: iosduobsmpp ! - outside space domain (zonal velocity component) 948 INTEGER :: iosdvobsmpp ! - outside space domain (meridional velocity component) 949 INTEGER :: ilanuobsmpp ! - within a model land cell (zonal velocity component) 950 INTEGER :: ilanvobsmpp ! - within a model land cell (meridional velocity component) 951 INTEGER :: inlauobsmpp ! - close to land (zonal velocity component) 952 INTEGER :: inlavobsmpp ! - close to land (meridional velocity component) 953 INTEGER :: igrdobsmpp ! - fail the grid search 954 INTEGER :: iuvchkumpp ! - reject u if v rejected and vice versa 955 INTEGER :: iuvchkvmpp ! 956 TYPE(obs_prof_valid) :: llvalid ! Profile selection 957 TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: & 958 & llvvalid ! U,V selection 959 INTEGER :: jvar ! Variable loop variable 960 INTEGER :: jobs ! Obs. loop variable 961 INTEGER :: jstp ! Time loop variable 962 INTEGER :: inrc ! Time index variable 963 964 IF(lwp) WRITE(numout,*)'obs_pre_vel: Preparing the velocity profile data' 965 966 ! Initial date initialization (year, month, day, hour, minute) 967 iyea0 = ndate0 / 10000 968 imon0 = ( ndate0 - iyea0 * 10000 ) / 100 969 iday0 = ndate0 - iyea0 * 10000 - imon0 * 100 970 ihou0 = 0 971 imin0 = 0 972 973 icycle = no ! Assimilation cycle 974 975 ! Diagnotics counters for various failures. 976 977 iotdobs = 0 978 igrdobs = 0 979 iosduobs = 0 980 iosdvobs = 0 981 ilanuobs = 0 982 ilanvobs = 0 983 inlauobs = 0 984 inlavobs = 0 985 iuvchku = 0 986 iuvchkv = 0 987 988 ! ----------------------------------------------------------------------- 989 ! Find time coordinate for profiles 990 ! ----------------------------------------------------------------------- 991 992 CALL obs_coo_tim_prof( icycle, & 993 & iyea0, imon0, iday0, ihou0, imin0, & 994 & profdata%nprof, profdata%nyea, profdata%nmon, & 995 & profdata%nday, profdata%nhou, profdata%nmin, & 996 & profdata%ntyp, profdata%nqc, profdata%mstp, & 997 & iotdobs, ld_dailyav = ld_dailyav ) 998 396 999 397 CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 1000 398 … … 1011 409 1012 410 ! ----------------------------------------------------------------------- 1013 ! Reject all observations for profiles with nqc > 101014 ! ----------------------------------------------------------------------- 1015 1016 CALL obs_pro_rej( profdata )411 ! Reject all observations for profiles with nqc > iqc_cutoff 412 ! ----------------------------------------------------------------------- 413 414 CALL obs_pro_rej( profdata, kqc_cutoff = iqc_cutoff ) 1017 415 1018 416 ! ----------------------------------------------------------------------- … … 1021 419 ! ----------------------------------------------------------------------- 1022 420 1023 ! Zonal Velocity Component 1024 421 ! Variable 1 1025 422 CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(1), & 1026 423 & profdata%npvsta(:,1), profdata%npvend(:,1), & 1027 424 & jpi, jpj, & 1028 425 & jpk, & 1029 & profdata%mi, profdata%mj, & 426 & profdata%mi, profdata%mj, & 1030 427 & profdata%var(1)%mvk, & 1031 428 & profdata%rlam, profdata%rphi, & 1032 429 & profdata%var(1)%vdep, & 1033 & glamu, gphiu,&1034 & gdept_1d, umask,&430 & pglam1, pgphi1, & 431 & gdept_1d, zmask1, & 1035 432 & profdata%nqc, profdata%var(1)%nvqc, & 1036 & iosduobs, ilanuobs, & 1037 & inlauobs, ld_nea ) 1038 1039 CALL obs_mpp_sum_integer( iosduobs, iosduobsmpp ) 1040 CALL obs_mpp_sum_integer( ilanuobs, ilanuobsmpp ) 1041 CALL obs_mpp_sum_integer( inlauobs, inlauobsmpp ) 1042 1043 ! Meridional Velocity Component 1044 433 & iosdv1obs, ilanv1obs, & 434 & inlav1obs, ld_nea, & 435 & ibdyv1obs, ld_bound_reject, & 436 & iqc_cutoff ) 437 438 CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp ) 439 CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp ) 440 CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp ) 441 CALL obs_mpp_sum_integer( ibdyv1obs, ibdyv1obsmpp ) 442 443 ! Variable 2 1045 444 CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(2), & 1046 445 & profdata%npvsta(:,2), profdata%npvend(:,2), & … … 1051 450 & profdata%rlam, profdata%rphi, & 1052 451 & profdata%var(2)%vdep, & 1053 & glamv, gphiv,&1054 & gdept_1d, vmask,&452 & pglam2, pgphi2, & 453 & gdept_1d, zmask2, & 1055 454 & profdata%nqc, profdata%var(2)%nvqc, & 1056 & iosdvobs, ilanvobs, & 1057 & inlavobs, ld_nea ) 1058 1059 CALL obs_mpp_sum_integer( iosdvobs, iosdvobsmpp ) 1060 CALL obs_mpp_sum_integer( ilanvobs, ilanvobsmpp ) 1061 CALL obs_mpp_sum_integer( inlavobs, inlavobsmpp ) 455 & iosdv2obs, ilanv2obs, & 456 & inlav2obs, ld_nea, & 457 & ibdyv2obs, ld_bound_reject, & 458 & iqc_cutoff ) 459 460 CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp ) 461 CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp ) 462 CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp ) 463 CALL obs_mpp_sum_integer( ibdyv2obs, ibdyv2obsmpp ) 1062 464 1063 465 ! ----------------------------------------------------------------------- … … 1065 467 ! ----------------------------------------------------------------------- 1066 468 1067 CALL obs_uv_rej( profdata, iuvchku, iuvchkv ) 1068 CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp ) 1069 CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp ) 469 IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 470 CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff ) 471 CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp ) 472 CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp ) 473 ENDIF 1070 474 1071 475 ! ----------------------------------------------------------------------- … … 1081 485 END DO 1082 486 1083 ! We want all data which has qc flags = 01084 1085 llvalid%luse(:) = ( profdata%nqc(:) <= 10)487 ! We want all data which has qc flags <= iqc_cutoff 488 489 llvalid%luse(:) = ( profdata%nqc(:) <= iqc_cutoff ) 1086 490 DO jvar = 1,profdata%nvar 1087 llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10)491 llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff ) 1088 492 END DO 1089 493 … … 1106 510 1107 511 IF(lwp) THEN 512 1108 513 WRITE(numout,*) 1109 WRITE(numout,*) 'obs_pre_vel :' 1110 WRITE(numout,*) '~~~~~~~~~~~' 1111 WRITE(numout,*) 1112 WRITE(numout,*) ' Profiles outside time domain = ', & 514 WRITE(numout,*) ' Profiles outside time domain = ', & 1113 515 & iotdobsmpp 1114 WRITE(numout,*) ' Remaining profiles that failed grid search = ', &516 WRITE(numout,*) ' Remaining profiles that failed grid search = ', & 1115 517 & igrdobsmpp 1116 WRITE(numout,*) ' Remaining Udata outside space domain = ', &1117 & iosd uobsmpp1118 WRITE(numout,*) ' Remaining Udata at land points = ', &1119 & ilan uobsmpp518 WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data outside space domain = ', & 519 & iosdv1obsmpp 520 WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data at land points = ', & 521 & ilanv1obsmpp 1120 522 IF (ld_nea) THEN 1121 WRITE(numout,*) ' Remaining Udata near land points (removed) = ',&1122 & inla uobsmpp523 WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (removed) = ',& 524 & inlav1obsmpp 1123 525 ELSE 1124 WRITE(numout,*) ' Remaining U data near land points (kept) = ',& 1125 & inlauobsmpp 1126 ENDIF 1127 WRITE(numout,*) ' U observation rejected since V rejected = ', & 1128 & iuvchku 1129 WRITE(numout,*) ' U data accepted = ', & 526 WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (kept) = ',& 527 & inlav1obsmpp 528 ENDIF 529 IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 530 WRITE(numout,*) ' U observation rejected since V rejected = ', & 531 & iuvchku 532 ENDIF 533 WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near open boundary (removed) = ',& 534 & ibdyv1obsmpp 535 WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted = ', & 1130 536 & prodatqc%nvprotmpp(1) 1131 WRITE(numout,*) ' Remaining Vdata outside space domain = ', &1132 & iosdv obsmpp1133 WRITE(numout,*) ' Remaining Vdata at land points = ', &1134 & ilanv obsmpp537 WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data outside space domain = ', & 538 & iosdv2obsmpp 539 WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data at land points = ', & 540 & ilanv2obsmpp 1135 541 IF (ld_nea) THEN 1136 WRITE(numout,*) ' Remaining Vdata near land points (removed) = ',&1137 & inlav obsmpp542 WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (removed) = ',& 543 & inlav2obsmpp 1138 544 ELSE 1139 WRITE(numout,*) ' Remaining V data near land points (kept) = ',& 1140 & inlavobsmpp 1141 ENDIF 1142 WRITE(numout,*) ' V observation rejected since U rejected = ', & 1143 & iuvchkv 1144 WRITE(numout,*) ' V data accepted = ', & 545 WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (kept) = ',& 546 & inlav2obsmpp 547 ENDIF 548 IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 549 WRITE(numout,*) ' V observation rejected since U rejected = ', & 550 & iuvchkv 551 ENDIF 552 WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near open boundary (removed) = ',& 553 & ibdyv2obsmpp 554 WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted = ', & 1145 555 & prodatqc%nvprotmpp(2) 1146 556 … … 1148 558 WRITE(numout,*) ' Number of observations per time step :' 1149 559 WRITE(numout,*) 1150 WRITE(numout,997) 560 WRITE(numout,'(10X,A,5X,A,5X,A,A)')'Time step','Profiles', & 561 & ' '//prodatqc%cvars(1)//' ', & 562 & ' '//prodatqc%cvars(2)//' ' 1151 563 WRITE(numout,998) 1152 564 ENDIF … … 1182 594 ENDIF 1183 595 1184 997 FORMAT(10X,'Time step',5X,'Profiles',5X,'Zonal Comp.',5X,'Meridional Comp.')1185 596 998 FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'----------------') 1186 597 999 FORMAT(10X,I9,5X,I8,5X,I11,5X,I8) 1187 598 1188 END SUBROUTINE obs_pre_ vel599 END SUBROUTINE obs_pre_prof 1189 600 1190 601 SUBROUTINE obs_coo_tim( kcycle, & … … 1293 704 & .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN 1294 705 kobsstp(jobs) = -1 1295 kobsqc(jobs) = kobsqc(jobs) + 11706 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 1296 707 kotdobs = kotdobs + 1 1297 708 CYCLE … … 1344 755 IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) & 1345 756 & .OR.( kobsstp(jobs) > nitend ) ) THEN 1346 kobsqc(jobs) = kobsqc(jobs) + 12757 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 1347 758 kotdobs = kotdobs + 1 1348 759 CYCLE … … 1389 800 & kobsyea, kobsmon, kobsday, kobshou, kobsmin, & 1390 801 & ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes, & 1391 & ld_dailyav)802 & kqc_cutoff ) 1392 803 !!---------------------------------------------------------------------- 1393 804 !! *** ROUTINE obs_coo_tim *** … … 1433 844 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 1434 845 & kdailyavtypes ! Types for daily averages 1435 LOGICAL, OPTIONAL :: ld_dailyav ! All types are daily averages 846 INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff ! QC cutoff value 847 1436 848 !! * Local declarations 1437 849 INTEGER :: jobs 850 INTEGER :: iqc_cutoff=255 1438 851 1439 852 !----------------------------------------------------------------------- … … 1454 867 DO jobs = 1, kobsno 1455 868 1456 IF ( kobsqc(jobs) <= 10) THEN869 IF ( kobsqc(jobs) <= iqc_cutoff ) THEN 1457 870 1458 871 IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.& 1459 872 & ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN 1460 kobsqc(jobs) = kobsqc(jobs) + 14873 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 1461 874 kotdobs = kotdobs + 1 1462 875 CYCLE … … 1467 880 ENDIF 1468 881 1469 !------------------------------------------------------------------------1470 ! If ld_dailyav is set then all data assumed to be daily averaged1471 !------------------------------------------------------------------------1472 1473 IF ( PRESENT( ld_dailyav) ) THEN1474 IF (ld_dailyav) THEN1475 DO jobs = 1, kobsno1476 1477 IF ( kobsqc(jobs) <= 10 ) THEN1478 1479 IF ( kobsstp(jobs) == (nit000 - 1) ) THEN1480 kobsqc(jobs) = kobsqc(jobs) + 141481 kotdobs = kotdobs + 11482 CYCLE1483 ENDIF1484 1485 ENDIF1486 END DO1487 ENDIF1488 ENDIF1489 882 1490 883 END SUBROUTINE obs_coo_tim_prof … … 1521 914 DO jobs = 1, kobsno 1522 915 IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN 1523 kobsqc(jobs) = kobsqc(jobs) + 18916 kobsqc(jobs) = IBSET(kobsqc(jobs),12) 1524 917 kgrdobs = kgrdobs + 1 1525 918 ENDIF … … 1532 925 & plam, pphi, pmask, & 1533 926 & kobsqc, kosdobs, klanobs, & 1534 & knlaobs,ld_nea ) 927 & knlaobs,ld_nea, & 928 & kbdyobs,ld_bound_reject, & 929 & kqc_cutoff ) 1535 930 !!---------------------------------------------------------------------- 1536 931 !! *** ROUTINE obs_coo_spc_2d *** … … 1565 960 INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: & 1566 961 & kobsqc ! Observation quality control 1567 INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain 1568 INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell 1569 INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land 1570 LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land 962 INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain 963 INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell 964 INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land 965 INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary 966 LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land 967 LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary 968 INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value 969 1571 970 !! * Local declarations 1572 971 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 1573 972 & zgmsk ! Grid mask 973 #if defined key_bdy 974 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 975 & zbmsk ! Boundary mask 976 REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 977 #endif 1574 978 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 1575 979 & zglam, & ! Model longitude at grid points … … 1588 992 ! For invalid points use 2,2 1589 993 1590 IF ( kobsqc(jobs) >= 10) THEN994 IF ( kobsqc(jobs) >= kqc_cutoff ) THEN 1591 995 1592 996 igrdi(1,1,jobs) = 1 … … 1613 1017 1614 1018 END DO 1615 1616 CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, pmask, zgmsk ) 1617 CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, plam, zglam ) 1618 CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, pphi, zgphi ) 1019 1020 #if defined key_bdy 1021 ! Create a mask grid points in boundary rim 1022 IF (ld_bound_reject) THEN 1023 zbdymask(:,:) = 1.0_wp 1024 DO ji = 1, nb_bdy 1025 DO jj = 1, idx_bdy(ji)%nblen(1) 1026 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 1027 ENDDO 1028 ENDDO 1029 1030 CALL obs_int_comm_2d( 2, 2, kobsno, igrdi, igrdj, zbdymask, zbmsk ) 1031 ENDIF 1032 #endif 1033 1034 CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk ) 1035 CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, plam, zglam ) 1036 CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 1619 1037 1620 1038 DO jobs = 1, kobsno 1621 1039 1622 1040 ! Skip bad observations 1623 IF ( kobsqc(jobs) >= 10) CYCLE1041 IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE 1624 1042 1625 1043 ! Flag if the observation falls outside the model spatial domain … … 1628 1046 & .OR. ( pobsphi(jobs) < -90. ) & 1629 1047 & .OR. ( pobsphi(jobs) > 90. ) ) THEN 1630 kobsqc(jobs) = kobsqc(jobs) + 111048 kobsqc(jobs) = IBSET(kobsqc(jobs),11) 1631 1049 kosdobs = kosdobs + 1 1632 1050 CYCLE … … 1635 1053 ! Flag if the observation falls with a model land cell 1636 1054 IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 1637 kobsqc(jobs) = kobsqc(jobs) + 121055 kobsqc(jobs) = IBSET(kobsqc(jobs),10) 1638 1056 klanobs = klanobs + 1 1639 1057 CYCLE … … 1657 1075 END DO 1658 1076 END DO 1659 1660 ! For observations on the grid reject them if their are at 1661 ! a masked point 1662 1077 1663 1078 IF (lgridobs) THEN 1664 1079 IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN 1665 kobsqc(jobs) = kobsqc(jobs) + 121080 kobsqc(jobs) = IBSET(kobsqc(jobs),10) 1666 1081 klanobs = klanobs + 1 1667 1082 CYCLE 1668 1083 ENDIF 1669 1084 ENDIF 1670 1085 1086 1671 1087 ! Flag if the observation falls is close to land 1672 1088 IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN 1673 IF (ld_nea) kobsqc(jobs) = kobsqc(jobs) + 141674 1089 knlaobs = knlaobs + 1 1675 CYCLE 1676 ENDIF 1090 IF (ld_nea) THEN 1091 kobsqc(jobs) = IBSET(kobsqc(jobs),9) 1092 CYCLE 1093 ENDIF 1094 ENDIF 1095 1096 #if defined key_bdy 1097 ! Flag if the observation falls close to the boundary rim 1098 IF (ld_bound_reject) THEN 1099 IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 1100 kobsqc(jobs) = IBSET(kobsqc(jobs),8) 1101 kbdyobs = kbdyobs + 1 1102 CYCLE 1103 ENDIF 1104 ! for observations on the grid... 1105 IF (lgridobs) THEN 1106 IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 1107 kobsqc(jobs) = IBSET(kobsqc(jobs),8) 1108 kbdyobs = kbdyobs + 1 1109 CYCLE 1110 ENDIF 1111 ENDIF 1112 ENDIF 1113 #endif 1677 1114 1678 1115 END DO … … 1686 1123 & plam, pphi, pdep, pmask, & 1687 1124 & kpobsqc, kobsqc, kosdobs, & 1688 & klanobs, knlaobs, ld_nea ) 1125 & klanobs, knlaobs, ld_nea, & 1126 & kbdyobs, ld_bound_reject, & 1127 & kqc_cutoff ) 1689 1128 !!---------------------------------------------------------------------- 1690 1129 !! *** ROUTINE obs_coo_spc_3d *** … … 1709 1148 !! * Modules used 1710 1149 USE dom_oce, ONLY : & ! Geographical information 1711 & gdepw_1d 1150 & gdepw_1d, & 1151 & gdepw_0, & 1152 #if defined key_vvl 1153 & gdepw_n, & 1154 & gdept_n, & 1155 #endif 1156 & ln_zco, & 1157 & ln_zps, & 1158 & lk_vvl 1712 1159 1713 1160 !! * Arguments … … 1743 1190 INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell 1744 1191 INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land 1192 INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary 1745 1193 LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land 1194 LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary 1195 INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value 1196 1746 1197 !! * Local declarations 1747 1198 REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 1748 1199 & zgmsk ! Grid mask 1200 #if defined key_bdy 1201 REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 1202 & zbmsk ! Boundary mask 1203 REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 1204 #endif 1205 REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 1206 & zgdepw 1749 1207 REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 1750 1208 & zglam, & ! Model longitude at grid points … … 1754 1212 & igrdj 1755 1213 LOGICAL :: lgridobs ! Is observation on a model grid point. 1214 LOGICAL :: ll_next_to_land ! Is a profile next to land 1756 1215 INTEGER :: iig, ijg ! i,j of observation on model grid point. 1757 1216 INTEGER :: jobs, jobsp, jk, ji, jj … … 1763 1222 ! For invalid points use 2,2 1764 1223 1765 IF ( kpobsqc(jobs) >= 10) THEN1224 IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN 1766 1225 1767 1226 igrdi(1,1,jobs) = 1 … … 1788 1247 1789 1248 END DO 1790 1791 CALL obs_int_comm_3d( 2, 2, kprofno, kpk, igrdi, igrdj, pmask, zgmsk ) 1792 CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, plam, zglam ) 1793 CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, pphi, zgphi ) 1249 1250 #if defined key_bdy 1251 ! Create a mask grid points in boundary rim 1252 IF (ld_bound_reject) THEN 1253 zbdymask(:,:) = 1.0_wp 1254 DO ji = 1, nb_bdy 1255 DO jj = 1, idx_bdy(ji)%nblen(1) 1256 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 1257 ENDDO 1258 ENDDO 1259 ENDIF 1260 1261 CALL obs_int_comm_2d( 2, 2, kprofno, igrdi, igrdj, zbdymask, zbmsk ) 1262 #endif 1263 1264 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk ) 1265 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam ) 1266 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 1267 ! Need to know the bathy depth for each observation for sco 1268 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, fsdepw(:,:,:), & 1269 & zgdepw ) 1794 1270 1795 1271 DO jobs = 1, kprofno 1796 1272 1797 1273 ! Skip bad profiles 1798 IF ( kpobsqc(jobs) >= 10) CYCLE1274 IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE 1799 1275 1800 1276 ! Check if this observation is on a grid point … … 1816 1292 END DO 1817 1293 1294 ! Check if next to land 1295 IF ( ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN 1296 ll_next_to_land=.TRUE. 1297 ELSE 1298 ll_next_to_land=.FALSE. 1299 ENDIF 1300 1818 1301 ! Reject observations 1819 1302 … … 1827 1310 & .OR. ( pobsdep(jobsp) < 0.0 ) & 1828 1311 & .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN 1829 kobsqc(jobsp) = kobsqc(jobsp) + 111312 kobsqc(jobsp) = IBSET(kobsqc(jobsp),11) 1830 1313 kosdobs = kosdobs + 1 1831 1314 CYCLE 1832 1315 ENDIF 1833 1316 1834 ! Flag if the observation falls with a model land cell 1835 IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1836 & == 0.0_wp ) THEN 1837 kobsqc(jobsp) = kobsqc(jobsp) + 12 1838 klanobs = klanobs + 1 1839 CYCLE 1317 ! To check if an observations falls within land there are two cases: 1318 ! 1: z-coordibnates, where the check uses the mask 1319 ! 2: terrain following (eg s-coordinates), 1320 ! where we use the depth of the bottom cell to mask observations 1321 1322 IF( (.NOT. lk_vvl) .AND. ( ln_zps .OR. ln_zco ) ) THEN !(CASE 1) 1323 1324 ! Flag if the observation falls with a model land cell 1325 IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1326 & == 0.0_wp ) THEN 1327 kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 1328 klanobs = klanobs + 1 1329 CYCLE 1330 ENDIF 1331 1332 ! Flag if the observation is close to land 1333 IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 1334 & 0.0_wp) THEN 1335 knlaobs = knlaobs + 1 1336 IF (ld_nea) THEN 1337 kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 1338 ENDIF 1339 ENDIF 1340 1341 ELSE ! Case 2 1342 ! Flag if the observation is deeper than the bathymetry 1343 ! Or if it is within the mask 1344 IF ( ANY( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 1345 & .OR. & 1346 & ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1347 & == 0.0_wp) ) THEN 1348 kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 1349 klanobs = klanobs + 1 1350 CYCLE 1351 ENDIF 1352 1353 ! Flag if the observation is close to land 1354 IF ( ll_next_to_land ) THEN 1355 knlaobs = knlaobs + 1 1356 IF (ld_nea) THEN 1357 kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 1358 ENDIF 1359 ENDIF 1360 1840 1361 ENDIF 1841 1362 … … 1845 1366 IF (lgridobs) THEN 1846 1367 IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN 1847 kobsqc(jobsp) = kobsqc(jobsp) + 121368 kobsqc(jobsp) = IBSET(kobsqc(jobs),10) 1848 1369 klanobs = klanobs + 1 1849 1370 CYCLE … … 1851 1372 ENDIF 1852 1373 1853 ! Flag if the observation falls is close to land1854 IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &1855 & 0.0_wp) THEN1856 IF (ld_nea) kobsqc(jobsp) = kobsqc(jobsp) + 141857 knlaobs = knlaobs + 11858 ENDIF1859 1860 1374 ! Set observation depth equal to that of the first model depth 1861 1375 IF ( pobsdep(jobsp) <= pdep(1) ) THEN … … 1863 1377 ENDIF 1864 1378 1379 #if defined key_bdy 1380 ! Flag if the observation falls close to the boundary rim 1381 IF (ld_bound_reject) THEN 1382 IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 1383 kobsqc(jobsp) = IBSET(kobsqc(jobs),8) 1384 kbdyobs = kbdyobs + 1 1385 CYCLE 1386 ENDIF 1387 ! for observations on the grid... 1388 IF (lgridobs) THEN 1389 IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 1390 kobsqc(jobsp) = IBSET(kobsqc(jobs),8) 1391 kbdyobs = kbdyobs + 1 1392 CYCLE 1393 ENDIF 1394 ENDIF 1395 ENDIF 1396 #endif 1397 1865 1398 END DO 1866 1399 END DO … … 1868 1401 END SUBROUTINE obs_coo_spc_3d 1869 1402 1870 SUBROUTINE obs_pro_rej( profdata )1403 SUBROUTINE obs_pro_rej( profdata, kqc_cutoff ) 1871 1404 !!---------------------------------------------------------------------- 1872 1405 !! *** ROUTINE obs_pro_rej *** … … 1886 1419 !! * Arguments 1887 1420 TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data 1421 INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value 1422 1888 1423 !! * Local declarations 1889 1424 INTEGER :: jprof … … 1895 1430 DO jprof = 1, profdata%nprof 1896 1431 1897 IF ( profdata%nqc(jprof) > 10) THEN1432 IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN 1898 1433 1899 1434 DO jvar = 1, profdata%nvar … … 1903 1438 1904 1439 profdata%var(jvar)%nvqc(jobs) = & 1905 & profdata%var(jvar)%nvqc(jobs) + 261440 & IBSET(profdata%var(jvar)%nvqc(jobs),14) 1906 1441 1907 1442 END DO … … 1915 1450 END SUBROUTINE obs_pro_rej 1916 1451 1917 SUBROUTINE obs_uv_rej( profdata, knumu, knumv )1452 SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff ) 1918 1453 !!---------------------------------------------------------------------- 1919 1454 !! *** ROUTINE obs_uv_rej *** … … 1935 1470 INTEGER, INTENT(INOUT) :: knumu ! Number of u rejected 1936 1471 INTEGER, INTENT(INOUT) :: knumv ! Number of v rejected 1472 INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value 1473 1937 1474 !! * Local declarations 1938 1475 INTEGER :: jprof … … 1954 1491 DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1) 1955 1492 1956 IF ( ( profdata%var(1)%nvqc(jobs) > 10) .AND. &1957 & ( profdata%var(2)%nvqc(jobs) <= 10) ) THEN1958 profdata%var(2)%nvqc(jobs) = profdata%var(2)%nvqc(jobs) + 421493 IF ( ( profdata%var(1)%nvqc(jobs) > kqc_cutoff ) .AND. & 1494 & ( profdata%var(2)%nvqc(jobs) <= kqc_cutoff) ) THEN 1495 profdata%var(2)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 1959 1496 knumv = knumv + 1 1960 1497 ENDIF 1961 IF ( ( profdata%var(2)%nvqc(jobs) > 10) .AND. &1962 & ( profdata%var(1)%nvqc(jobs) <= 10) ) THEN1963 profdata%var(1)%nvqc(jobs) = profdata%var(1)%nvqc(jobs) + 421498 IF ( ( profdata%var(2)%nvqc(jobs) > kqc_cutoff ) .AND. & 1499 & ( profdata%var(1)%nvqc(jobs) <= kqc_cutoff) ) THEN 1500 profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 1964 1501 knumu = knumu + 1 1965 1502 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
