- Timestamp:
- 2017-12-13T18:08:50+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
r7646 r9023 23 23 USE obs_oper ! Observation operators 24 24 USE lib_mpp, ONLY : ctl_warn, ctl_stop 25 USE bdy_oce, ONLY : & ! Boundary information 26 idx_bdy, nb_bdy, ln_bdy 25 27 26 28 IMPLICIT NONE … … 40 42 CONTAINS 41 43 42 SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea ) 44 SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, & 45 kqc_cutoff ) 43 46 !!---------------------------------------------------------------------- 44 47 !! *** ROUTINE obs_pre_sla *** … … 57 60 !! ! 2015-02 (M. Martin) Combined routine for surface types. 58 61 !!---------------------------------------------------------------------- 62 !! * Modules used 59 63 USE par_oce ! Ocean parameters 60 64 USE dom_oce, ONLY : glamt, gphit, tmask, nproc ! Geographical information … … 63 67 TYPE(obs_surf), INTENT(INOUT) :: surfdataqc ! Subset of surface data not failing screening 64 68 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 65 ! 69 LOGICAL, INTENT(IN) :: ld_bound_reject ! Switch for rejecting obs near the boundary 70 INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value 71 !! * Local declarations 72 INTEGER :: iqc_cutoff = 255 ! cut off for QC value 66 73 INTEGER :: iyea0 ! Initial date 67 74 INTEGER :: imon0 ! - (year, month, day, hour, minute) … … 76 83 INTEGER :: inlasobs ! - close to land 77 84 INTEGER :: igrdobs ! - fail the grid search 85 INTEGER :: ibdysobs ! - close to open boundary 78 86 ! Global counters for observations that 79 87 INTEGER :: iotdobsmpp ! - outside time domain … … 82 90 INTEGER :: inlasobsmpp ! - close to land 83 91 INTEGER :: igrdobsmpp ! - fail the grid search 84 LOGICAL, DIMENSION(:), ALLOCATABLE :: llvalid ! SLA data selection 92 INTEGER :: ibdysobsmpp ! - close to open boundary 93 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 94 & llvalid ! SLA data selection 85 95 INTEGER :: jobs ! Obs. loop variable 86 96 INTEGER :: jstp ! Time loop variable … … 107 117 ilansobs = 0 108 118 inlasobs = 0 119 ibdysobs = 0 120 121 ! Set QC cutoff to optional value if provided 122 IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 109 123 110 124 ! ----------------------------------------------------------------------- … … 140 154 & tmask(:,:,1), surfdata%nqc, & 141 155 & iosdsobs, ilansobs, & 142 & inlasobs, ld_nea ) 156 & inlasobs, ld_nea, & 157 & ibdysobs, ld_bound_reject, & 158 & iqc_cutoff ) 143 159 144 160 CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 145 161 CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 146 162 CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 163 CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 147 164 148 165 ! ----------------------------------------------------------------------- … … 155 172 ALLOCATE( llvalid(surfdata%nsurf) ) 156 173 157 ! We want all data which has qc flags <= 10158 159 llvalid(:) = ( surfdata%nqc(:) <= 10)174 ! We want all data which has qc flags <= iqc_cutoff 175 176 llvalid(:) = ( surfdata%nqc(:) <= iqc_cutoff ) 160 177 161 178 ! The actual copying … … 190 207 & inlasobsmpp 191 208 ENDIF 209 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', & 210 & ibdysobsmpp 192 211 WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted = ', & 193 212 & surfdataqc%nsurfmpp … … 225 244 & kpi, kpj, kpk, & 226 245 & zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2, & 227 & ld_nea, kdailyavtypes)246 & ld_nea, ld_bound_reject, kdailyavtypes, kqc_cutoff ) 228 247 229 248 !!---------------------------------------------------------------------- … … 241 260 !! 242 261 !!---------------------------------------------------------------------- 243 USE par_oce ! Ocean parameters 244 USE dom_oce, ONLY : gdept_1d, nproc ! Geographical information 262 !! * Modules used 263 USE par_oce ! Ocean parameters 264 USE dom_oce, ONLY : & ! Geographical information 265 & gdept_1d, & 266 & nproc 245 267 246 268 !! * Arguments … … 250 272 LOGICAL, INTENT(IN) :: ld_var2 251 273 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 274 LOGICAL, INTENT(IN) :: ld_bound_reject ! Switch for rejecting observations near the boundary 252 275 INTEGER, INTENT(IN) :: kpi, kpj, kpk ! Local domain sizes 253 276 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & … … 261 284 & pgphi1, & 262 285 & pgphi2 286 INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value 263 287 264 288 !! * Local declarations 289 INTEGER :: iqc_cutoff = 255 ! cut off for QC value 265 290 INTEGER :: iyea0 ! Initial date 266 291 INTEGER :: imon0 ! - (year, month, day, hour, minute) … … 277 302 INTEGER :: inlav1obs ! - close to land (variable 1) 278 303 INTEGER :: inlav2obs ! - close to land (variable 2) 304 INTEGER :: ibdyv1obs ! - boundary (variable 1) 305 INTEGER :: ibdyv2obs ! - boundary (variable 2) 279 306 INTEGER :: igrdobs ! - fail the grid search 280 307 INTEGER :: iuvchku ! - reject u if v rejected and vice versa … … 288 315 INTEGER :: inlav1obsmpp ! - close to land (variable 1) 289 316 INTEGER :: inlav2obsmpp ! - close to land (variable 2) 317 INTEGER :: ibdyv1obsmpp ! - boundary (variable 1) 318 INTEGER :: ibdyv2obsmpp ! - boundary (variable 2) 290 319 INTEGER :: igrdobsmpp ! - fail the grid search 291 320 INTEGER :: iuvchkumpp ! - reject var1 if var2 rejected and vice versa … … 322 351 inlav1obs = 0 323 352 inlav2obs = 0 324 iuvchku = 0 325 iuvchkv = 0 353 ibdyv1obs = 0 354 ibdyv2obs = 0 355 iuvchku = 0 356 iuvchkv = 0 357 358 359 ! Set QC cutoff to optional value if provided 360 IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 326 361 327 362 ! ----------------------------------------------------------------------- … … 335 370 & profdata%nday, profdata%nhou, profdata%nmin, & 336 371 & profdata%ntyp, profdata%nqc, profdata%mstp, & 337 & iotdobs, kdailyavtypes = kdailyavtypes ) 372 & iotdobs, kdailyavtypes = kdailyavtypes, & 373 & kqc_cutoff = iqc_cutoff ) 338 374 ELSE 339 375 CALL obs_coo_tim_prof( icycle, & … … 342 378 & profdata%nday, profdata%nhou, profdata%nmin, & 343 379 & profdata%ntyp, profdata%nqc, profdata%mstp, & 344 & iotdobs )380 & iotdobs, kqc_cutoff = iqc_cutoff ) 345 381 ENDIF 346 382 … … 359 395 360 396 ! ----------------------------------------------------------------------- 361 ! Reject all observations for profiles with nqc > 10362 ! ----------------------------------------------------------------------- 363 364 CALL obs_pro_rej( profdata )397 ! Reject all observations for profiles with nqc > iqc_cutoff 398 ! ----------------------------------------------------------------------- 399 400 CALL obs_pro_rej( profdata, kqc_cutoff = iqc_cutoff ) 365 401 366 402 ! ----------------------------------------------------------------------- … … 381 417 & gdept_1d, zmask1, & 382 418 & profdata%nqc, profdata%var(1)%nvqc, & 383 & iosdv1obs, ilanv1obs, & 384 & inlav1obs, ld_nea ) 419 & iosdv1obs, ilanv1obs, & 420 & inlav1obs, ld_nea, & 421 & ibdyv1obs, ld_bound_reject, & 422 & iqc_cutoff ) 385 423 386 424 CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp ) 387 425 CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp ) 388 426 CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp ) 427 CALL obs_mpp_sum_integer( ibdyv1obs, ibdyv1obsmpp ) 389 428 390 429 ! Variable 2 … … 400 439 & gdept_1d, zmask2, & 401 440 & profdata%nqc, profdata%var(2)%nvqc, & 402 & iosdv2obs, ilanv2obs, & 403 & inlav2obs, ld_nea ) 441 & iosdv2obs, ilanv2obs, & 442 & inlav2obs, ld_nea, & 443 & ibdyv2obs, ld_bound_reject, & 444 & iqc_cutoff ) 404 445 405 446 CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp ) 406 447 CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp ) 407 448 CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp ) 449 CALL obs_mpp_sum_integer( ibdyv2obs, ibdyv2obsmpp ) 408 450 409 451 ! ----------------------------------------------------------------------- … … 412 454 413 455 IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 414 CALL obs_uv_rej( profdata, iuvchku, iuvchkv )456 CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff ) 415 457 CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp ) 416 458 CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp ) … … 429 471 END DO 430 472 431 ! We want all data which has qc flags = 0432 433 llvalid%luse(:) = ( profdata%nqc(:) <= 10)473 ! We want all data which has qc flags <= iqc_cutoff 474 475 llvalid%luse(:) = ( profdata%nqc(:) <= iqc_cutoff ) 434 476 DO jvar = 1,profdata%nvar 435 llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10)477 llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff ) 436 478 END DO 437 479 … … 475 517 & iuvchku 476 518 ENDIF 519 WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near open boundary (removed) = ',& 520 & ibdyv1obsmpp 477 521 WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted = ', & 478 522 & prodatqc%nvprotmpp(1) … … 492 536 & iuvchkv 493 537 ENDIF 538 WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near open boundary (removed) = ',& 539 & ibdyv2obsmpp 494 540 WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted = ', & 495 541 & prodatqc%nvprotmpp(2) … … 644 690 & .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN 645 691 kobsstp(jobs) = -1 646 kobsqc(jobs) = kobsqc(jobs) + 11692 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 647 693 kotdobs = kotdobs + 1 648 694 CYCLE … … 695 741 IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) & 696 742 & .OR.( kobsstp(jobs) > nitend ) ) THEN 697 kobsqc(jobs) = kobsqc(jobs) + 12743 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 698 744 kotdobs = kotdobs + 1 699 745 CYCLE … … 739 785 & kobsno, & 740 786 & kobsyea, kobsmon, kobsday, kobshou, kobsmin, & 741 & ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes ) 787 & ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes, & 788 & kqc_cutoff ) 742 789 !!---------------------------------------------------------------------- 743 790 !! *** ROUTINE obs_coo_tim *** … … 783 830 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 784 831 & kdailyavtypes ! Types for daily averages 832 INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff ! QC cutoff value 833 785 834 !! * Local declarations 786 835 INTEGER :: jobs 836 INTEGER :: iqc_cutoff=255 787 837 788 838 !----------------------------------------------------------------------- … … 803 853 DO jobs = 1, kobsno 804 854 805 IF ( kobsqc(jobs) <= 10) THEN855 IF ( kobsqc(jobs) <= iqc_cutoff ) THEN 806 856 807 857 IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.& 808 858 & ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN 809 kobsqc(jobs) = kobsqc(jobs) + 14859 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 810 860 kotdobs = kotdobs + 1 811 861 CYCLE … … 850 900 DO jobs = 1, kobsno 851 901 IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN 852 kobsqc(jobs) = kobsqc(jobs) + 18902 kobsqc(jobs) = IBSET(kobsqc(jobs),12) 853 903 kgrdobs = kgrdobs + 1 854 904 ENDIF … … 861 911 & plam, pphi, pmask, & 862 912 & kobsqc, kosdobs, klanobs, & 863 & knlaobs,ld_nea ) 913 & knlaobs,ld_nea, & 914 & kbdyobs,ld_bound_reject, & 915 & kqc_cutoff ) 864 916 !!---------------------------------------------------------------------- 865 917 !! *** ROUTINE obs_coo_spc_2d *** … … 894 946 INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: & 895 947 & kobsqc ! Observation quality control 896 INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain 897 INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell 898 INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land 899 LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land 948 INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain 949 INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell 950 INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land 951 INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary 952 LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land 953 LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary 954 INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value 955 900 956 !! * Local declarations 901 957 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 902 958 & zgmsk ! Grid mask 959 960 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 961 & zbmsk ! Boundary mask 962 REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 903 963 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 904 964 & zglam, & ! Model longitude at grid points … … 917 977 ! For invalid points use 2,2 918 978 919 IF ( kobsqc(jobs) >= 10) THEN979 IF ( kobsqc(jobs) >= kqc_cutoff ) THEN 920 980 921 981 igrdi(1,1,jobs) = 1 … … 942 1002 943 1003 END DO 1004 1005 IF (ln_bdy) THEN 1006 ! Create a mask grid points in boundary rim 1007 IF (ld_bound_reject) THEN 1008 zbdymask(:,:) = 1.0_wp 1009 DO ji = 1, nb_bdy 1010 DO jj = 1, idx_bdy(ji)%nblen(1) 1011 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 1012 ENDDO 1013 ENDDO 1014 1015 CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk ) 1016 ENDIF 1017 ENDIF 1018 944 1019 945 1020 CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk ) … … 950 1025 951 1026 ! Skip bad observations 952 IF ( kobsqc(jobs) >= 10) CYCLE1027 IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE 953 1028 954 1029 ! Flag if the observation falls outside the model spatial domain … … 957 1032 & .OR. ( pobsphi(jobs) < -90. ) & 958 1033 & .OR. ( pobsphi(jobs) > 90. ) ) THEN 959 kobsqc(jobs) = kobsqc(jobs) + 111034 kobsqc(jobs) = IBSET(kobsqc(jobs),11) 960 1035 kosdobs = kosdobs + 1 961 1036 CYCLE … … 964 1039 ! Flag if the observation falls with a model land cell 965 1040 IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 966 kobsqc(jobs) = kobsqc(jobs) + 121041 kobsqc(jobs) = IBSET(kobsqc(jobs),10) 967 1042 klanobs = klanobs + 1 968 1043 CYCLE … … 978 1053 IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) & 979 1054 & .AND. & 980 & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp )&981 & ) THEN1055 & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) & 1056 & < 1.0e-6_wp ) ) THEN 982 1057 lgridobs = .TRUE. 983 1058 iig = ji … … 992 1067 IF (lgridobs) THEN 993 1068 IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN 994 kobsqc(jobs) = kobsqc(jobs) + 121069 kobsqc(jobs) = IBSET(kobsqc(jobs),10) 995 1070 klanobs = klanobs + 1 996 1071 CYCLE … … 1000 1075 ! Flag if the observation falls is close to land 1001 1076 IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN 1002 IF (ld_nea) kobsqc(jobs) = kobsqc(jobs) + 141003 1077 knlaobs = knlaobs + 1 1004 CYCLE 1078 IF (ld_nea) THEN 1079 kobsqc(jobs) = IBSET(kobsqc(jobs),9) 1080 CYCLE 1081 ENDIF 1082 ENDIF 1083 1084 IF (ln_bdy) THEN 1085 ! Flag if the observation falls close to the boundary rim 1086 IF (ld_bound_reject) THEN 1087 IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 1088 kobsqc(jobs) = IBSET(kobsqc(jobs),8) 1089 kbdyobs = kbdyobs + 1 1090 CYCLE 1091 ENDIF 1092 ! for observations on the grid... 1093 IF (lgridobs) THEN 1094 IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 1095 kobsqc(jobs) = IBSET(kobsqc(jobs),8) 1096 kbdyobs = kbdyobs + 1 1097 CYCLE 1098 ENDIF 1099 ENDIF 1100 ENDIF 1005 1101 ENDIF 1006 1102 … … 1015 1111 & plam, pphi, pdep, pmask, & 1016 1112 & kpobsqc, kobsqc, kosdobs, & 1017 & klanobs, knlaobs, ld_nea ) 1113 & klanobs, knlaobs, ld_nea, & 1114 & kbdyobs, ld_bound_reject, & 1115 & kqc_cutoff ) 1018 1116 !!---------------------------------------------------------------------- 1019 1117 !! *** ROUTINE obs_coo_spc_3d *** … … 1077 1175 INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell 1078 1176 INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land 1177 INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary 1079 1178 LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land 1179 LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary 1180 INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value 1181 1080 1182 !! * Local declarations 1081 1183 REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 1082 1184 & zgmsk ! Grid mask 1185 REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 1186 & zbmsk ! Boundary mask 1187 REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 1083 1188 REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 1084 1189 & zgdepw … … 1100 1205 ! For invalid points use 2,2 1101 1206 1102 IF ( kpobsqc(jobs) >= 10) THEN1207 IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN 1103 1208 1104 1209 igrdi(1,1,jobs) = 1 … … 1125 1230 1126 1231 END DO 1232 1233 IF (ln_bdy) THEN 1234 ! Create a mask grid points in boundary rim 1235 IF (ld_bound_reject) THEN 1236 zbdymask(:,:) = 1.0_wp 1237 DO ji = 1, nb_bdy 1238 DO jj = 1, idx_bdy(ji)%nblen(1) 1239 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 1240 ENDDO 1241 ENDDO 1242 ENDIF 1243 1244 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk ) 1245 ENDIF 1127 1246 1128 1247 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk ) 1129 1248 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam ) 1130 1249 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 1131 IF ( .NOT.( ln_zps .OR. ln_zco ) ) THEN 1132 ! Need to know the bathy depth for each observation for sco 1133 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), & 1250 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), & 1134 1251 & zgdepw ) 1135 ENDIF1136 1252 1137 1253 DO jobs = 1, kprofno 1138 1254 1139 1255 ! Skip bad profiles 1140 IF ( kpobsqc(jobs) >= 10) CYCLE1256 IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE 1141 1257 1142 1258 ! Check if this observation is on a grid point … … 1149 1265 IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) & 1150 1266 & .AND. & 1151 & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &1267 & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) < 1.0e-6_wp ) & 1152 1268 & ) THEN 1153 1269 lgridobs = .TRUE. … … 1176 1292 & .OR. ( pobsdep(jobsp) < 0.0 ) & 1177 1293 & .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN 1178 kobsqc(jobsp) = kobsqc(jobsp) + 111294 kobsqc(jobsp) = IBSET(kobsqc(jobsp),11) 1179 1295 kosdobs = kosdobs + 1 1180 1296 CYCLE 1181 1297 ENDIF 1182 1298 1183 ! To check if an observations falls within land there are two cases: 1184 ! 1: z-coordibnates, where the check uses the mask 1185 ! 2: terrain following (eg s-coordinates), 1186 ! where we use the depth of the bottom cell to mask observations 1299 ! To check if an observations falls within land: 1187 1300 1188 IF( ln_zps .OR. ln_zco ) THEN !(CASE 1) 1301 ! Flag if the observation is deeper than the bathymetry 1302 ! Or if it is within the mask 1303 IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 1304 & .OR. & 1305 & ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1306 & == 0.0_wp) ) THEN 1307 kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 1308 klanobs = klanobs + 1 1309 CYCLE 1310 ENDIF 1189 1311 1190 ! Flag if the observation falls with a model land cell 1191 IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1192 & == 0.0_wp ) THEN 1193 kobsqc(jobsp) = kobsqc(jobsp) + 12 1194 klanobs = klanobs + 1 1195 CYCLE 1196 ENDIF 1197 1198 ! Flag if the observation is close to land 1199 IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 1200 & 0.0_wp) THEN 1201 knlaobs = knlaobs + 1 1202 IF (ld_nea) THEN 1203 kobsqc(jobsp) = kobsqc(jobsp) + 14 1204 ENDIF 1205 ENDIF 1206 1207 ELSE ! Case 2 1208 1209 ! Flag if the observation is deeper than the bathymetry 1210 ! Or if it is within the mask 1211 IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 1212 & .OR. & 1213 & ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1214 & == 0.0_wp) ) THEN 1215 kobsqc(jobsp) = kobsqc(jobsp) + 12 1216 klanobs = klanobs + 1 1217 CYCLE 1218 ENDIF 1219 1220 ! Flag if the observation is close to land 1221 IF ( ll_next_to_land ) THEN 1222 knlaobs = knlaobs + 1 1223 IF (ld_nea) THEN 1224 kobsqc(jobsp) = kobsqc(jobsp) + 14 1225 ENDIF 1226 ENDIF 1312 ! Flag if the observation is close to land 1313 IF ( ll_next_to_land ) THEN 1314 knlaobs = knlaobs + 1 1315 IF (ld_nea) THEN 1316 kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 1317 ENDIF 1227 1318 ENDIF 1228 1319 … … 1232 1323 IF (lgridobs) THEN 1233 1324 IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN 1234 kobsqc(jobsp) = kobsqc(jobsp) + 121325 kobsqc(jobsp) = IBSET(kobsqc(jobs),10) 1235 1326 klanobs = klanobs + 1 1236 1327 CYCLE … … 1250 1341 ENDIF 1251 1342 1343 IF (ln_bdy) THEN 1344 ! Flag if the observation falls close to the boundary rim 1345 IF (ld_bound_reject) THEN 1346 IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 1347 kobsqc(jobsp) = IBSET(kobsqc(jobs),8) 1348 kbdyobs = kbdyobs + 1 1349 CYCLE 1350 ENDIF 1351 ! for observations on the grid... 1352 IF (lgridobs) THEN 1353 IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 1354 kobsqc(jobsp) = IBSET(kobsqc(jobs),8) 1355 kbdyobs = kbdyobs + 1 1356 CYCLE 1357 ENDIF 1358 ENDIF 1359 ENDIF 1360 ENDIF 1361 1252 1362 END DO 1253 1363 END DO … … 1255 1365 END SUBROUTINE obs_coo_spc_3d 1256 1366 1257 SUBROUTINE obs_pro_rej( profdata )1367 SUBROUTINE obs_pro_rej( profdata, kqc_cutoff ) 1258 1368 !!---------------------------------------------------------------------- 1259 1369 !! *** ROUTINE obs_pro_rej *** … … 1273 1383 !! * Arguments 1274 1384 TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data 1385 INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value 1386 1275 1387 !! * Local declarations 1276 1388 INTEGER :: jprof … … 1282 1394 DO jprof = 1, profdata%nprof 1283 1395 1284 IF ( profdata%nqc(jprof) > 10) THEN1396 IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN 1285 1397 1286 1398 DO jvar = 1, profdata%nvar … … 1290 1402 1291 1403 profdata%var(jvar)%nvqc(jobs) = & 1292 & profdata%var(jvar)%nvqc(jobs) + 261404 & IBSET(profdata%var(jvar)%nvqc(jobs),14) 1293 1405 1294 1406 END DO … … 1302 1414 END SUBROUTINE obs_pro_rej 1303 1415 1304 SUBROUTINE obs_uv_rej( profdata, knumu, knumv )1416 SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff ) 1305 1417 !!---------------------------------------------------------------------- 1306 1418 !! *** ROUTINE obs_uv_rej *** … … 1322 1434 INTEGER, INTENT(INOUT) :: knumu ! Number of u rejected 1323 1435 INTEGER, INTENT(INOUT) :: knumv ! Number of v rejected 1436 INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value 1437 1324 1438 !! * Local declarations 1325 1439 INTEGER :: jprof … … 1341 1455 DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1) 1342 1456 1343 IF ( ( profdata%var(1)%nvqc(jobs) > 10) .AND. &1344 & ( profdata%var(2)%nvqc(jobs) <= 10) ) THEN1345 profdata%var(2)%nvqc(jobs) = profdata%var(2)%nvqc(jobs) + 421457 IF ( ( profdata%var(1)%nvqc(jobs) > kqc_cutoff ) .AND. & 1458 & ( profdata%var(2)%nvqc(jobs) <= kqc_cutoff) ) THEN 1459 profdata%var(2)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 1346 1460 knumv = knumv + 1 1347 1461 ENDIF 1348 IF ( ( profdata%var(2)%nvqc(jobs) > 10) .AND. &1349 & ( profdata%var(1)%nvqc(jobs) <= 10) ) THEN1350 profdata%var(1)%nvqc(jobs) = profdata%var(1)%nvqc(jobs) + 421462 IF ( ( profdata%var(2)%nvqc(jobs) > kqc_cutoff ) .AND. & 1463 & ( profdata%var(1)%nvqc(jobs) <= kqc_cutoff) ) THEN 1464 profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 1351 1465 knumu = knumu + 1 1352 1466 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
