- Timestamp:
- 2016-11-28T17:04:10+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
r4292 r7351 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 … … 36 33 37 34 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 35 & obs_pre_prof, & ! First level check and screening of profile obs 36 & obs_pre_surf, & ! First level check and screening of surface obs 37 & calc_month_len ! Calculate the number of days in the months of a year 44 38 45 39 !!---------------------------------------------------------------------- … … 49 43 !!---------------------------------------------------------------------- 50 44 45 51 46 CONTAINS 52 47 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 ) 48 SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea ) 340 49 !!---------------------------------------------------------------------- 341 50 !! *** ROUTINE obs_pre_sla *** 342 51 !! 343 !! ** Purpose : First level check and screening of SLAobservations344 !! 345 !! ** Method : First level check and screening of SLAobservations52 !! ** Purpose : First level check and screening of surface observations 53 !! 54 !! ** Method : First level check and screening of surface observations 346 55 !! 347 56 !! ** Action : … … 352 61 !! ! 2007-03 (A. Weaver, K. Mogensen) Original 353 62 !! ! 2007-06 (K. Mogensen et al) Reject obs. near land. 63 !! ! 2015-02 (M. Martin) Combined routine for surface types. 354 64 !!---------------------------------------------------------------------- 355 65 !! * Modules used … … 362 72 & nproc 363 73 !! * 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 74 TYPE(obs_surf), INTENT(INOUT) :: surfdata ! Full set of surface data 75 TYPE(obs_surf), INTENT(INOUT) :: surfdataqc ! Subset of surface data not failing screening 367 76 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 368 77 !! * Local declarations … … 391 100 INTEGER :: inrc ! Time index variable 392 101 393 IF(lwp) WRITE(numout,*)'obs_pre_sla : Preparing the SLA observations...' 394 102 IF(lwp) WRITE(numout,*)'obs_pre_surf : Preparing the surface observations...' 103 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 104 395 105 ! Initial date initialization (year, month, day, hour, minute) 396 106 iyea0 = ndate0 / 10000 397 107 imon0 = ( ndate0 - iyea0 * 10000 ) / 100 398 108 iday0 = ndate0 - iyea0 * 10000 - imon0 * 100 399 ihou0 = 0400 imin0 = 0109 ihou0 = nn_time0 / 100 110 imin0 = ( nn_time0 - ihou0 * 100 ) 401 111 402 112 icycle = no ! Assimilation cycle … … 411 121 412 122 ! ----------------------------------------------------------------------- 413 ! Find time coordinate for SLAdata123 ! Find time coordinate for surface data 414 124 ! ----------------------------------------------------------------------- 415 125 416 126 CALL obs_coo_tim( icycle, & 417 127 & 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 )128 & surfdata%nsurf, surfdata%nyea, surfdata%nmon, & 129 & surfdata%nday, surfdata%nhou, surfdata%nmin, & 130 & surfdata%nqc, surfdata%mstp, iotdobs ) 421 131 422 132 CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 423 133 424 134 ! ----------------------------------------------------------------------- 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 )135 ! Check for surface data failing the grid search 136 ! ----------------------------------------------------------------------- 137 138 CALL obs_coo_grd( surfdata%nsurf, surfdata%mi, surfdata%mj, & 139 & surfdata%nqc, igrdobs ) 430 140 431 141 CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) … … 435 145 ! ----------------------------------------------------------------------- 436 146 437 CALL obs_coo_spc_2d( s ladata%nsurf, &147 CALL obs_coo_spc_2d( surfdata%nsurf, & 438 148 & jpi, jpj, & 439 & s ladata%mi, sladata%mj, &440 & s ladata%rlam, sladata%rphi, &149 & surfdata%mi, surfdata%mj, & 150 & surfdata%rlam, surfdata%rphi, & 441 151 & glamt, gphit, & 442 & tmask(:,:,1), s ladata%nqc, &152 & tmask(:,:,1), surfdata%nqc, & 443 153 & iosdsobs, ilansobs, & 444 154 & inlasobs, ld_nea ) … … 449 159 450 160 ! ----------------------------------------------------------------------- 451 ! Copy useful data from the s ladata data structure to452 ! the s ladatqc data structure161 ! Copy useful data from the surfdata data structure to 162 ! the surfdataqc data structure 453 163 ! ----------------------------------------------------------------------- 454 164 455 165 ! Allocate the selection arrays 456 166 457 ALLOCATE( llvalid(s ladata%nsurf) )167 ALLOCATE( llvalid(surfdata%nsurf) ) 458 168 459 169 ! We want all data which has qc flags <= 10 460 170 461 llvalid(:) = ( s ladata%nqc(:) <= 10 )171 llvalid(:) = ( surfdata%nqc(:) <= 10 ) 462 172 463 173 ! The actual copying 464 174 465 CALL obs_surf_compress( s ladata, sladatqc, .TRUE., numout, &175 CALL obs_surf_compress( surfdata, surfdataqc, .TRUE., numout, & 466 176 & lvalid=llvalid ) 467 177 … … 477 187 IF(lwp) THEN 478 188 WRITE(numout,*) 479 WRITE(numout,*) 'obs_pre_sla :' 480 WRITE(numout,*) '~~~~~~~~~~~' 481 WRITE(numout,*) 482 WRITE(numout,*) ' SLA data outside time domain = ', & 189 WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data outside time domain = ', & 483 190 & iotdobsmpp 484 WRITE(numout,*) ' Remaining SLAdata that failed grid search = ', &191 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data that failed grid search = ', & 485 192 & igrdobsmpp 486 WRITE(numout,*) ' Remaining SLAdata outside space domain = ', &193 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data outside space domain = ', & 487 194 & iosdsobsmpp 488 WRITE(numout,*) ' Remaining SLAdata at land points = ', &195 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data at land points = ', & 489 196 & ilansobsmpp 490 197 IF (ld_nea) THEN 491 WRITE(numout,*) ' Remaining SLAdata near land points (removed) = ', &198 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (removed) = ', & 492 199 & inlasobsmpp 493 200 ELSE 494 WRITE(numout,*) ' Remaining SLAdata near land points (kept) = ', &201 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (kept) = ', & 495 202 & inlasobsmpp 496 203 ENDIF 497 WRITE(numout,*) ' SLAdata accepted = ', &498 & s ladatqc%nsurfmpp204 WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted = ', & 205 & surfdataqc%nsurfmpp 499 206 500 207 WRITE(numout,*) 501 208 WRITE(numout,*) ' Number of observations per time step :' 502 209 WRITE(numout,*) 503 WRITE(numout,1997) 504 WRITE(numout,1998) 210 WRITE(numout,'(10X,A,10X,A)')'Time step',surfdataqc%cvars(1) 211 WRITE(numout,'(10X,A,5X,A)')'---------','-----------------' 212 CALL FLUSH(numout) 505 213 ENDIF 506 214 507 DO jobs = 1, s ladatqc%nsurf508 inrc = s ladatqc%mstp(jobs) + 2 - nit000509 s ladatqc%nsstp(inrc) = sladatqc%nsstp(inrc) + 1215 DO jobs = 1, surfdataqc%nsurf 216 inrc = surfdataqc%mstp(jobs) + 2 - nit000 217 surfdataqc%nsstp(inrc) = surfdataqc%nsstp(inrc) + 1 510 218 END DO 511 219 512 CALL obs_mpp_sum_integers( s ladatqc%nsstp, sladatqc%nsstpmpp, &220 CALL obs_mpp_sum_integers( surfdataqc%nsstp, surfdataqc%nsstpmpp, & 513 221 & nitend - nit000 + 2 ) 514 222 … … 516 224 DO jstp = nit000 - 1, nitend 517 225 inrc = jstp - nit000 + 2 518 WRITE(numout,1999) jstp, sladatqc%nsstpmpp(inrc) 226 WRITE(numout,1999) jstp, surfdataqc%nsstpmpp(inrc) 227 CALL FLUSH(numout) 519 228 END DO 520 229 ENDIF 521 230 522 1997 FORMAT(10X,'Time step',5X,'Sea level anomaly')523 1998 FORMAT(10X,'---------',5X,'-----------------')524 231 1999 FORMAT(10X,I9,5X,I17) 525 232 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 !! 233 END SUBROUTINE obs_pre_surf 234 235 236 SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_var1, ld_var2, & 237 & kpi, kpj, kpk, & 238 & zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2, & 239 & ld_nea, kdailyavtypes ) 240 241 !!---------------------------------------------------------------------- 242 !! *** ROUTINE obs_pre_prof *** 243 !! 244 !! ** Purpose : First level check and screening of profiles 245 !! 246 !! ** Method : First level check and screening of profiles 247 !! 540 248 !! History : 541 !! ! 2007-03 (S. Ricci) SST data preparation 249 !! ! 2007-06 (K. Mogensen) original : T and S profile data 250 !! ! 2008-09 (M. Valdivieso) : TAO velocity data 251 !! ! 2009-01 (K. Mogensen) : New feedback stricture 252 !! ! 2015-02 (M. Martin) : Combined profile routine. 253 !! 542 254 !!---------------------------------------------------------------------- 543 255 !! * Modules used … … 545 257 USE par_oce ! Ocean parameters 546 258 USE dom_oce, ONLY : & ! Geographical information 547 & glamt, & 548 & gphit, & 549 & tmask, & 259 & gdept_1d, & 550 260 & nproc 551 !! * 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 556 !! * Local declarations 557 INTEGER :: iyea0 ! Initial date 558 INTEGER :: imon0 ! - (year, month, day, hour, minute) 559 INTEGER :: iday0 560 INTEGER :: ihou0 561 INTEGER :: imin0 562 INTEGER :: icycle ! Current assimilation cycle 563 ! Counters for observations that 564 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 568 INTEGER :: igrdobs ! - fail the grid search 569 ! Global counters for observations that 570 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 574 INTEGER :: igrdobsmpp ! - fail the grid search 575 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 576 & llvalid ! SST data selection 577 INTEGER :: jobs ! Obs. loop variable 578 INTEGER :: jstp ! Time loop variable 579 INTEGER :: inrc ! Time index variable 580 581 IF(lwp) WRITE(numout,*)'obs_pre_sst : Preparing the SST observations...' 582 583 ! Initial date initialization (year, month, day, hour, minute) 584 iyea0 = ndate0 / 10000 585 imon0 = ( ndate0 - iyea0 * 10000 ) / 100 586 iday0 = ndate0 - iyea0 * 10000 - imon0 * 100 587 ihou0 = 0 588 imin0 = 0 589 590 icycle = no ! Assimilation cycle 591 592 ! Diagnotics counters for various failures. 593 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) 690 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 261 921 262 !! * Arguments 922 263 TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data 923 264 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 265 LOGICAL, INTENT(IN) :: ld_var1 ! Observed variables switches 266 LOGICAL, INTENT(IN) :: ld_var2 267 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 268 INTEGER, INTENT(IN) :: kpi, kpj, kpk ! Local domain sizes 269 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 270 & kdailyavtypes ! Types for daily averages 271 REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 272 & zmask1, & 273 & zmask2 274 REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 275 & pglam1, & 276 & pglam2, & 277 & pgphi1, & 278 & pgphi2 279 927 280 !! * Local declarations 928 281 INTEGER :: iyea0 ! Initial date … … 932 285 INTEGER :: imin0 933 286 INTEGER :: icycle ! Current assimilation cycle 934 ! Counters for observations that 287 ! Counters for observations that are 935 288 INTEGER :: iotdobs ! - outside time domain 936 INTEGER :: iosd uobs ! - outside space domain (zonal velocity component)937 INTEGER :: iosdv obs ! - outside space domain (meridional velocity component)938 INTEGER :: ilan uobs ! - within a model land cell (zonal velocity component)939 INTEGER :: ilanv obs ! - within a model land cell (meridional velocity component)940 INTEGER :: inla uobs ! - close to land (zonal velocity component)941 INTEGER :: inlav obs ! - close to land (meridional velocity component)289 INTEGER :: iosdv1obs ! - outside space domain (variable 1) 290 INTEGER :: iosdv2obs ! - outside space domain (variable 2) 291 INTEGER :: ilanv1obs ! - within a model land cell (variable 1) 292 INTEGER :: ilanv2obs ! - within a model land cell (variable 2) 293 INTEGER :: inlav1obs ! - close to land (variable 1) 294 INTEGER :: inlav2obs ! - close to land (variable 2) 942 295 INTEGER :: igrdobs ! - fail the grid search 943 296 INTEGER :: iuvchku ! - reject u if v rejected and vice versa 944 297 INTEGER :: iuvchkv ! 945 ! Global counters for observations that 298 ! Global counters for observations that are 946 299 INTEGER :: iotdobsmpp ! - outside time domain 947 INTEGER :: iosd uobsmpp ! - outside space domain (zonal velocity component)948 INTEGER :: iosdv obsmpp ! - outside space domain (meridional velocity component)949 INTEGER :: ilan uobsmpp ! - within a model land cell (zonal velocity component)950 INTEGER :: ilanv obsmpp ! - within a model land cell (meridional velocity component)951 INTEGER :: inla uobsmpp ! - close to land (zonal velocity component)952 INTEGER :: inlav obsmpp ! - close to land (meridional velocity component)300 INTEGER :: iosdv1obsmpp ! - outside space domain (variable 1) 301 INTEGER :: iosdv2obsmpp ! - outside space domain (variable 2) 302 INTEGER :: ilanv1obsmpp ! - within a model land cell (variable 1) 303 INTEGER :: ilanv2obsmpp ! - within a model land cell (variable 2) 304 INTEGER :: inlav1obsmpp ! - close to land (variable 1) 305 INTEGER :: inlav2obsmpp ! - close to land (variable 2) 953 306 INTEGER :: igrdobsmpp ! - fail the grid search 954 INTEGER :: iuvchkumpp ! - reject u if vrejected and vice versa307 INTEGER :: iuvchkumpp ! - reject var1 if var2 rejected and vice versa 955 308 INTEGER :: iuvchkvmpp ! 956 309 TYPE(obs_prof_valid) :: llvalid ! Profile selection 957 310 TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: & 958 & llvvalid ! U,Vselection311 & llvvalid ! var1,var2 selection 959 312 INTEGER :: jvar ! Variable loop variable 960 313 INTEGER :: jobs ! Obs. loop variable … … 962 315 INTEGER :: inrc ! Time index variable 963 316 964 IF(lwp) WRITE(numout,*)'obs_pre_vel: Preparing the velocity profile data' 317 IF(lwp) WRITE(numout,*)'obs_pre_prof: Preparing the profile data...' 318 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 965 319 966 320 ! Initial date initialization (year, month, day, hour, minute) … … 968 322 imon0 = ( ndate0 - iyea0 * 10000 ) / 100 969 323 iday0 = ndate0 - iyea0 * 10000 - imon0 * 100 970 ihou0 = 0971 imin0 = 0324 ihou0 = nn_time0 / 100 325 imin0 = ( nn_time0 - ihou0 * 100 ) 972 326 973 327 icycle = no ! Assimilation cycle … … 977 331 iotdobs = 0 978 332 igrdobs = 0 979 iosd uobs = 0980 iosdv obs = 0981 ilan uobs = 0982 ilanv obs = 0983 inla uobs = 0984 inlav obs = 0333 iosdv1obs = 0 334 iosdv2obs = 0 335 ilanv1obs = 0 336 ilanv2obs = 0 337 inlav1obs = 0 338 inlav2obs = 0 985 339 iuvchku = 0 986 340 iuvchkv = 0 … … 990 344 ! ----------------------------------------------------------------------- 991 345 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 346 IF ( PRESENT(kdailyavtypes) ) THEN 347 CALL obs_coo_tim_prof( icycle, & 348 & iyea0, imon0, iday0, ihou0, imin0, & 349 & profdata%nprof, profdata%nyea, profdata%nmon, & 350 & profdata%nday, profdata%nhou, profdata%nmin, & 351 & profdata%ntyp, profdata%nqc, profdata%mstp, & 352 & iotdobs, kdailyavtypes = kdailyavtypes ) 353 ELSE 354 CALL obs_coo_tim_prof( icycle, & 355 & iyea0, imon0, iday0, ihou0, imin0, & 356 & profdata%nprof, profdata%nyea, profdata%nmon, & 357 & profdata%nday, profdata%nhou, profdata%nmin, & 358 & profdata%ntyp, profdata%nqc, profdata%mstp, & 359 & iotdobs ) 360 ENDIF 361 999 362 CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 1000 363 … … 1021 384 ! ----------------------------------------------------------------------- 1022 385 1023 ! Zonal Velocity Component 1024 386 ! Variable 1 1025 387 CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(1), & 1026 388 & profdata%npvsta(:,1), profdata%npvend(:,1), & 1027 389 & jpi, jpj, & 1028 390 & jpk, & 1029 & profdata%mi, profdata%mj, & 391 & profdata%mi, profdata%mj, & 1030 392 & profdata%var(1)%mvk, & 1031 393 & profdata%rlam, profdata%rphi, & 1032 394 & profdata%var(1)%vdep, & 1033 & glamu, gphiu,&1034 & gdept_1d, umask,&395 & pglam1, pgphi1, & 396 & gdept_1d, zmask1, & 1035 397 & 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 398 & iosdv1obs, ilanv1obs, & 399 & inlav1obs, ld_nea ) 400 401 CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp ) 402 CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp ) 403 CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp ) 404 405 ! Variable 2 1045 406 CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(2), & 1046 407 & profdata%npvsta(:,2), profdata%npvend(:,2), & … … 1051 412 & profdata%rlam, profdata%rphi, & 1052 413 & profdata%var(2)%vdep, & 1053 & glamv, gphiv,&1054 & gdept_1d, vmask,&414 & pglam2, pgphi2, & 415 & gdept_1d, zmask2, & 1055 416 & profdata%nqc, profdata%var(2)%nvqc, & 1056 & iosdv obs, ilanvobs,&1057 & inlav obs, ld_nea )1058 1059 CALL obs_mpp_sum_integer( iosdv obs, iosdvobsmpp )1060 CALL obs_mpp_sum_integer( ilanv obs, ilanvobsmpp )1061 CALL obs_mpp_sum_integer( inlav obs, inlavobsmpp )417 & iosdv2obs, ilanv2obs, & 418 & inlav2obs, ld_nea ) 419 420 CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp ) 421 CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp ) 422 CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp ) 1062 423 1063 424 ! ----------------------------------------------------------------------- … … 1065 426 ! ----------------------------------------------------------------------- 1066 427 1067 CALL obs_uv_rej( profdata, iuvchku, iuvchkv ) 1068 CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp ) 1069 CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp ) 428 IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 429 CALL obs_uv_rej( profdata, iuvchku, iuvchkv ) 430 CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp ) 431 CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp ) 432 ENDIF 1070 433 1071 434 ! ----------------------------------------------------------------------- … … 1106 469 1107 470 IF(lwp) THEN 471 1108 472 WRITE(numout,*) 1109 WRITE(numout,*) 'obs_pre_vel :' 1110 WRITE(numout,*) '~~~~~~~~~~~' 1111 WRITE(numout,*) 1112 WRITE(numout,*) ' Profiles outside time domain = ', & 473 WRITE(numout,*) ' Profiles outside time domain = ', & 1113 474 & iotdobsmpp 1114 WRITE(numout,*) ' Remaining profiles that failed grid search = ', &475 WRITE(numout,*) ' Remaining profiles that failed grid search = ', & 1115 476 & igrdobsmpp 1116 WRITE(numout,*) ' Remaining Udata outside space domain = ', &1117 & iosd uobsmpp1118 WRITE(numout,*) ' Remaining Udata at land points = ', &1119 & ilan uobsmpp477 WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data outside space domain = ', & 478 & iosdv1obsmpp 479 WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data at land points = ', & 480 & ilanv1obsmpp 1120 481 IF (ld_nea) THEN 1121 WRITE(numout,*) ' Remaining Udata near land points (removed) = ',&1122 & inla uobsmpp482 WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (removed) = ',& 483 & inlav1obsmpp 1123 484 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 = ', & 485 WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (kept) = ',& 486 & inlav1obsmpp 487 ENDIF 488 IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 489 WRITE(numout,*) ' U observation rejected since V rejected = ', & 490 & iuvchku 491 ENDIF 492 WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted = ', & 1130 493 & prodatqc%nvprotmpp(1) 1131 WRITE(numout,*) ' Remaining Vdata outside space domain = ', &1132 & iosdv obsmpp1133 WRITE(numout,*) ' Remaining Vdata at land points = ', &1134 & ilanv obsmpp494 WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data outside space domain = ', & 495 & iosdv2obsmpp 496 WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data at land points = ', & 497 & ilanv2obsmpp 1135 498 IF (ld_nea) THEN 1136 WRITE(numout,*) ' Remaining Vdata near land points (removed) = ',&1137 & inlav obsmpp499 WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (removed) = ',& 500 & inlav2obsmpp 1138 501 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 = ', & 502 WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (kept) = ',& 503 & inlav2obsmpp 504 ENDIF 505 IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 506 WRITE(numout,*) ' V observation rejected since U rejected = ', & 507 & iuvchkv 508 ENDIF 509 WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted = ', & 1145 510 & prodatqc%nvprotmpp(2) 1146 511 … … 1148 513 WRITE(numout,*) ' Number of observations per time step :' 1149 514 WRITE(numout,*) 1150 WRITE(numout,997) 515 WRITE(numout,'(10X,A,5X,A,5X,A,A)')'Time step','Profiles', & 516 & ' '//prodatqc%cvars(1)//' ', & 517 & ' '//prodatqc%cvars(2)//' ' 1151 518 WRITE(numout,998) 1152 519 ENDIF … … 1182 549 ENDIF 1183 550 1184 997 FORMAT(10X,'Time step',5X,'Profiles',5X,'Zonal Comp.',5X,'Meridional Comp.')1185 551 998 FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'----------------') 1186 552 999 FORMAT(10X,I9,5X,I8,5X,I11,5X,I8) 1187 553 1188 END SUBROUTINE obs_pre_ vel554 END SUBROUTINE obs_pre_prof 1189 555 1190 556 SUBROUTINE obs_coo_tim( kcycle, & … … 1388 754 & kobsno, & 1389 755 & kobsyea, kobsmon, kobsday, kobshou, kobsmin, & 1390 & ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes, & 1391 & ld_dailyav ) 756 & ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes ) 1392 757 !!---------------------------------------------------------------------- 1393 758 !! *** ROUTINE obs_coo_tim *** … … 1433 798 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 1434 799 & kdailyavtypes ! Types for daily averages 1435 LOGICAL, OPTIONAL :: ld_dailyav ! All types are daily averages1436 800 !! * Local declarations 1437 801 INTEGER :: jobs … … 1467 831 ENDIF 1468 832 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 833 1490 834 END SUBROUTINE obs_coo_tim_prof … … 1614 958 END DO 1615 959 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 )960 CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk ) 961 CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, plam, zglam ) 962 CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 1619 963 1620 964 DO jobs = 1, kobsno … … 1709 1053 !! * Modules used 1710 1054 USE dom_oce, ONLY : & ! Geographical information 1711 & gdepw_1d 1055 & gdepw_1d, & 1056 & gdepw_0, & 1057 & gdepw_n, & 1058 & gdept_n, & 1059 & ln_zco, & 1060 & ln_zps 1712 1061 1713 1062 !! * Arguments … … 1747 1096 REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 1748 1097 & zgmsk ! Grid mask 1098 REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 1099 & zgdepw 1749 1100 REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 1750 1101 & zglam, & ! Model longitude at grid points … … 1754 1105 & igrdj 1755 1106 LOGICAL :: lgridobs ! Is observation on a model grid point. 1107 LOGICAL :: ll_next_to_land ! Is a profile next to land 1756 1108 INTEGER :: iig, ijg ! i,j of observation on model grid point. 1757 1109 INTEGER :: jobs, jobsp, jk, ji, jj … … 1789 1141 END DO 1790 1142 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 ) 1143 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk ) 1144 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam ) 1145 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 1146 IF ( .NOT.( ln_zps .OR. ln_zco ) ) THEN 1147 ! Need to know the bathy depth for each observation for sco 1148 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), & 1149 & zgdepw ) 1150 ENDIF 1794 1151 1795 1152 DO jobs = 1, kprofno … … 1816 1173 END DO 1817 1174 1175 ! Check if next to land 1176 IF ( ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN 1177 ll_next_to_land=.TRUE. 1178 ELSE 1179 ll_next_to_land=.FALSE. 1180 ENDIF 1181 1818 1182 ! Reject observations 1819 1183 … … 1832 1196 ENDIF 1833 1197 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 1198 ! To check if an observations falls within land there are two cases: 1199 ! 1: z-coordibnates, where the check uses the mask 1200 ! 2: terrain following (eg s-coordinates), 1201 ! where we use the depth of the bottom cell to mask observations 1202 1203 IF( ln_zps .OR. ln_zco ) THEN !(CASE 1) 1204 1205 ! Flag if the observation falls with a model land cell 1206 IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1207 & == 0.0_wp ) THEN 1208 kobsqc(jobsp) = kobsqc(jobsp) + 12 1209 klanobs = klanobs + 1 1210 CYCLE 1211 ENDIF 1212 1213 ! Flag if the observation is close to land 1214 IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 1215 & 0.0_wp) THEN 1216 knlaobs = knlaobs + 1 1217 IF (ld_nea) THEN 1218 kobsqc(jobsp) = kobsqc(jobsp) + 14 1219 ENDIF 1220 ENDIF 1221 1222 ELSE ! Case 2 1223 1224 ! Flag if the observation is deeper than the bathymetry 1225 ! Or if it is within the mask 1226 IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 1227 & .OR. & 1228 & ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1229 & == 0.0_wp) ) THEN 1230 kobsqc(jobsp) = kobsqc(jobsp) + 12 1231 klanobs = klanobs + 1 1232 CYCLE 1233 ENDIF 1234 1235 ! Flag if the observation is close to land 1236 IF ( ll_next_to_land ) THEN 1237 knlaobs = knlaobs + 1 1238 IF (ld_nea) THEN 1239 kobsqc(jobsp) = kobsqc(jobsp) + 14 1240 ENDIF 1241 ENDIF 1840 1242 ENDIF 1841 1243 1842 1244 ! For observations on the grid reject them if their are at 1843 1245 ! a masked point
Note: See TracChangeset
for help on using the changeset viewer.
