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