- Timestamp:
- 2018-02-05T16:07:40+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
r8653 r9306 60 60 61 61 62 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, 63 & kit000, kdaystp, 64 & pvar 1, pvar2, pgdept, pgdepw,&65 & pmask 1, pmask2, &66 & plam 1, plam2, pphi1, pphi2,&62 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 63 & kit000, kdaystp, kvar, & 64 & pvar, pgdept, pgdepw, & 65 & pmask, & 66 & plam, pphi, & 67 67 & k1dint, k2dint, kdailyavtypes ) 68 68 … … 134 134 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 135 135 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 136 INTEGER, INTENT(IN) :: kvar ! Number of variable in prodatqc 136 137 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 137 & pvar1, & ! Model field 1 138 & pvar2, & ! Model field 2 139 & pmask1, & ! Land-sea mask 1 140 & pmask2 ! Land-sea mask 2 138 & pvar, & ! Model field for variable 139 & pmask ! Land-sea mask for variable 141 140 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 142 & plam1, & ! Model longitudes for variable 1 143 & plam2, & ! Model longitudes for variable 2 144 & pphi1, & ! Model latitudes for variable 1 145 & pphi2 ! Model latitudes for variable 2 141 & plam, & ! Model longitudes for variable 142 & pphi ! Model latitudes for variable 146 143 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 147 144 & pgdept, & ! Model array of depth T levels … … 166 163 & idailyavtypes 167 164 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 168 & igrdi1, & 169 & igrdi2, & 170 & igrdj1, & 171 & igrdj2 165 & igrdi, & 166 & igrdj 172 167 INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 173 168 … … 176 171 REAL(KIND=wp) :: zdaystp 177 172 REAL(KIND=wp), DIMENSION(kpk) :: & 178 & zobsmask1, &179 & zobsmask2, &180 173 & zobsk, & 181 174 & zobs2k 182 175 REAL(KIND=wp), DIMENSION(2,2,1) :: & 183 176 & zweig1, & 184 & zweig2, &185 177 & zweig 186 178 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 187 & zmask1, & 188 & zmask2, & 189 & zint1, & 190 & zint2, & 191 & zinm1, & 192 & zinm2, & 179 & zmask, & 180 & zint, & 181 & zinm, & 193 182 & zgdept, & 194 183 & zgdepw 195 184 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 196 & zglam1, & 197 & zglam2, & 198 & zgphi1, & 199 & zgphi2 200 REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2 185 & zglam, & 186 & zgphi 187 REAL(KIND=wp), DIMENSION(1) :: zmsk 201 188 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 202 189 … … 230 217 DO jj = 1, jpj 231 218 DO ji = 1, jpi 232 prodatqc%vdmean(ji,jj,jk,1) = 0.0 233 prodatqc%vdmean(ji,jj,jk,2) = 0.0 219 prodatqc%vdmean(ji,jj,jk,kvar) = 0.0 234 220 END DO 235 221 END DO … … 240 226 DO jj = 1, jpj 241 227 DO ji = 1, jpi 242 ! Increment field 1 for computing daily mean 243 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 244 & + pvar1(ji,jj,jk) 245 ! Increment field 2 for computing daily mean 246 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 247 & + pvar2(ji,jj,jk) 228 ! Increment field for computing daily mean 229 prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) & 230 & + pvar(ji,jj,jk) 248 231 END DO 249 232 END DO … … 258 241 DO jj = 1, jpj 259 242 DO ji = 1, jpi 260 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 261 & * zdaystp 262 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 263 & * zdaystp 243 prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) & 244 & * zdaystp 264 245 END DO 265 246 END DO … … 271 252 ! Get the data for interpolation 272 253 ALLOCATE( & 273 & igrdi1(2,2,ipro), & 274 & igrdi2(2,2,ipro), & 275 & igrdj1(2,2,ipro), & 276 & igrdj2(2,2,ipro), & 277 & zglam1(2,2,ipro), & 278 & zglam2(2,2,ipro), & 279 & zgphi1(2,2,ipro), & 280 & zgphi2(2,2,ipro), & 281 & zmask1(2,2,kpk,ipro), & 282 & zmask2(2,2,kpk,ipro), & 283 & zint1(2,2,kpk,ipro), & 284 & zint2(2,2,kpk,ipro), & 254 & igrdi(2,2,ipro), & 255 & igrdj(2,2,ipro), & 256 & zglam(2,2,ipro), & 257 & zgphi(2,2,ipro), & 258 & zmask(2,2,kpk,ipro), & 259 & zint(2,2,kpk,ipro), & 285 260 & zgdept(2,2,kpk,ipro), & 286 261 & zgdepw(2,2,kpk,ipro) & … … 289 264 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 290 265 iobs = jobs - prodatqc%nprofup 291 igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1 292 igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1 293 igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1 294 igrdj1(1,2,iobs) = prodatqc%mj(jobs,1) 295 igrdi1(2,1,iobs) = prodatqc%mi(jobs,1) 296 igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1 297 igrdi1(2,2,iobs) = prodatqc%mi(jobs,1) 298 igrdj1(2,2,iobs) = prodatqc%mj(jobs,1) 299 igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1 300 igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1 301 igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1 302 igrdj2(1,2,iobs) = prodatqc%mj(jobs,2) 303 igrdi2(2,1,iobs) = prodatqc%mi(jobs,2) 304 igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1 305 igrdi2(2,2,iobs) = prodatqc%mi(jobs,2) 306 igrdj2(2,2,iobs) = prodatqc%mj(jobs,2) 266 igrdi(1,1,iobs) = prodatqc%mi(jobs,kvar)-1 267 igrdj(1,1,iobs) = prodatqc%mj(jobs,kvar)-1 268 igrdi(1,2,iobs) = prodatqc%mi(jobs,kvar)-1 269 igrdj(1,2,iobs) = prodatqc%mj(jobs,kvar) 270 igrdi(2,1,iobs) = prodatqc%mi(jobs,kvar) 271 igrdj(2,1,iobs) = prodatqc%mj(jobs,kvar)-1 272 igrdi(2,2,iobs) = prodatqc%mi(jobs,kvar) 273 igrdj(2,2,iobs) = prodatqc%mj(jobs,kvar) 307 274 END DO 308 275 … … 311 278 zgdepw(:,:,:,:) = 0.0 312 279 313 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 ) 314 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 ) 315 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pmask1, zmask1 ) 316 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1, zint1 ) 317 318 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 ) 319 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 ) 320 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 ) 321 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2, zint2 ) 322 323 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept ) 324 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw ) 280 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, plam, zglam ) 281 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 282 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pmask, zmask ) 283 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pvar, zint ) 284 285 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept, zgdept ) 286 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw ) 325 287 326 288 ! At the end of the day also get interpolated means 327 289 IF ( ld_dailyav .AND. idayend == 0 ) THEN 328 290 329 ALLOCATE( & 330 & zinm1(2,2,kpk,ipro), & 331 & zinm2(2,2,kpk,ipro) & 332 & ) 333 334 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, & 335 & prodatqc%vdmean(:,:,:,1), zinm1 ) 336 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, & 337 & prodatqc%vdmean(:,:,:,2), zinm2 ) 291 ALLOCATE( zinm(2,2,kpk,ipro) ) 292 293 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 294 & prodatqc%vdmean(:,:,:,kvar), zinm ) 338 295 339 296 ENDIF … … 370 327 ! Horizontal weights 371 328 ! Masked values are calculated later. 372 IF ( prodatqc%npvend(jobs, 1) > 0 ) THEN329 IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN 373 330 374 331 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 375 & zglam1(:,:,iobs), zgphi1(:,:,iobs), & 376 & zmask1(:,:,1,iobs), zweig1, zmsk_1 ) 377 378 ENDIF 379 380 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 381 382 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 383 & zglam2(:,:,iobs), zgphi2(:,:,iobs), & 384 & zmask2(:,:,1,iobs), zweig2, zmsk_2 ) 385 386 ENDIF 387 388 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 332 & zglam(:,:,iobs), zgphi(:,:,iobs), & 333 & zmask(:,:,1,iobs), zweig1, zmsk ) 334 335 ENDIF 336 337 IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN 389 338 390 339 zobsk(:) = obfillflt … … 396 345 397 346 ! vertically interpolate all 4 corners 398 ista = prodatqc%npvsta(jobs, 1)399 iend = prodatqc%npvend(jobs, 1)347 ista = prodatqc%npvsta(jobs,kvar) 348 iend = prodatqc%npvend(jobs,kvar) 400 349 inum_obs = iend - ista + 1 401 350 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) … … 406 355 IF ( k1dint == 1 ) THEN 407 356 CALL obs_int_z1d_spl( kpk, & 408 & zinm 1(iin,ijn,:,iobs), &357 & zinm(iin,ijn,:,iobs), & 409 358 & zobs2k, zgdept(iin,ijn,:,iobs), & 410 & zmask 1(iin,ijn,:,iobs))359 & zmask(iin,ijn,:,iobs)) 411 360 ENDIF 412 361 413 362 CALL obs_level_search(kpk, & 414 363 & zgdept(iin,ijn,:,iobs), & 415 & inum_obs, prodatqc%var( 1)%vdep(ista:iend), &364 & inum_obs, prodatqc%var(kvar)%vdep(ista:iend), & 416 365 & iv_indic) 417 366 418 367 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 419 & prodatqc%var( 1)%vdep(ista:iend), &420 & zinm 1(iin,ijn,:,iobs), &368 & prodatqc%var(kvar)%vdep(ista:iend), & 369 & zinm(iin,ijn,:,iobs), & 421 370 & zobs2k, interp_corner(iin,ijn,:), & 422 371 & zgdept(iin,ijn,:,iobs), & 423 & zmask 1(iin,ijn,:,iobs))372 & zmask(iin,ijn,:,iobs)) 424 373 425 374 ENDDO … … 433 382 434 383 ! vertically interpolate all 4 corners 435 ista = prodatqc%npvsta(jobs, 1)436 iend = prodatqc%npvend(jobs, 1)384 ista = prodatqc%npvsta(jobs,kvar) 385 iend = prodatqc%npvend(jobs,kvar) 437 386 inum_obs = iend - ista + 1 438 387 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) … … 442 391 IF ( k1dint == 1 ) THEN 443 392 CALL obs_int_z1d_spl( kpk, & 444 & zint 1(iin,ijn,:,iobs),&393 & zint(iin,ijn,:,iobs),& 445 394 & zobs2k, zgdept(iin,ijn,:,iobs), & 446 & zmask 1(iin,ijn,:,iobs))395 & zmask(iin,ijn,:,iobs)) 447 396 448 397 ENDIF … … 450 399 CALL obs_level_search(kpk, & 451 400 & zgdept(iin,ijn,:,iobs),& 452 & inum_obs, prodatqc%var( 1)%vdep(ista:iend), &401 & inum_obs, prodatqc%var(kvar)%vdep(ista:iend), & 453 402 & iv_indic) 454 403 455 404 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 456 & prodatqc%var( 1)%vdep(ista:iend), &457 & zint 1(iin,ijn,:,iobs), &405 & prodatqc%var(kvar)%vdep(ista:iend), & 406 & zint(iin,ijn,:,iobs), & 458 407 & zobs2k,interp_corner(iin,ijn,:), & 459 408 & zgdept(iin,ijn,:,iobs), & 460 & zmask 1(iin,ijn,:,iobs) )409 & zmask(iin,ijn,:,iobs) ) 461 410 462 411 ENDDO … … 482 431 DO ijn=1,2 483 432 484 depth_loop 1: DO ik=kpk,2,-1485 IF(zmask 1(iin,ijn,ik-1,iobs ) > 0.9 )THEN433 depth_loop: DO ik=kpk,2,-1 434 IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN 486 435 487 436 zweig(iin,ijn,1) = & 488 437 & zweig1(iin,ijn,1) * & 489 438 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 490 & - prodatqc%var( 1)%vdep(iend)),0._wp)439 & - prodatqc%var(kvar)%vdep(iend)),0._wp) 491 440 492 EXIT depth_loop 1441 EXIT depth_loop 493 442 494 443 ENDIF 495 444 496 ENDDO depth_loop 1445 ENDDO depth_loop 497 446 498 447 ENDDO … … 500 449 501 450 CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 502 & prodatqc%var( 1)%vmod(iend:iend) )451 & prodatqc%var(kvar)%vmod(iend:iend) ) 503 452 504 453 ! Set QC flag for any observations found below the bottom 505 454 ! needed as the check here is more strict than that in obs_prep 506 IF (sum(zweig) == 0.0_wp) prodatqc%var( 1)%nvqc(iend:iend)=4455 IF (sum(zweig) == 0.0_wp) prodatqc%var(kvar)%nvqc(iend:iend)=4 507 456 508 457 ENDDO … … 510 459 DEALLOCATE(interp_corner,iv_indic) 511 460 512 ENDIF 513 514 ! For the second variable 515 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 516 517 zobsk(:) = obfillflt 518 519 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 520 521 IF ( idayend == 0 ) THEN 522 ! Daily averaged data 523 524 ! vertically interpolate all 4 corners 525 ista = prodatqc%npvsta(jobs,2) 526 iend = prodatqc%npvend(jobs,2) 527 inum_obs = iend - ista + 1 528 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 529 530 DO iin=1,2 531 DO ijn=1,2 532 533 IF ( k1dint == 1 ) THEN 534 CALL obs_int_z1d_spl( kpk, & 535 & zinm2(iin,ijn,:,iobs), & 536 & zobs2k, zgdept(iin,ijn,:,iobs), & 537 & zmask2(iin,ijn,:,iobs)) 538 ENDIF 539 540 CALL obs_level_search(kpk, & 541 & zgdept(iin,ijn,:,iobs), & 542 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 543 & iv_indic) 544 545 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 546 & prodatqc%var(2)%vdep(ista:iend), & 547 & zinm2(iin,ijn,:,iobs), & 548 & zobs2k, interp_corner(iin,ijn,:), & 549 & zgdept(iin,ijn,:,iobs), & 550 & zmask2(iin,ijn,:,iobs)) 551 552 ENDDO 553 ENDDO 554 555 ENDIF !idayend 556 557 ELSE 558 559 ! Point data 560 561 ! vertically interpolate all 4 corners 562 ista = prodatqc%npvsta(jobs,2) 563 iend = prodatqc%npvend(jobs,2) 564 inum_obs = iend - ista + 1 565 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 566 DO iin=1,2 567 DO ijn=1,2 568 569 IF ( k1dint == 1 ) THEN 570 CALL obs_int_z1d_spl( kpk, & 571 & zint2(iin,ijn,:,iobs),& 572 & zobs2k, zgdept(iin,ijn,:,iobs), & 573 & zmask2(iin,ijn,:,iobs)) 574 575 ENDIF 576 577 CALL obs_level_search(kpk, & 578 & zgdept(iin,ijn,:,iobs),& 579 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 580 & iv_indic) 581 582 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 583 & prodatqc%var(2)%vdep(ista:iend), & 584 & zint2(iin,ijn,:,iobs), & 585 & zobs2k,interp_corner(iin,ijn,:), & 586 & zgdept(iin,ijn,:,iobs), & 587 & zmask2(iin,ijn,:,iobs) ) 588 589 ENDDO 590 ENDDO 591 592 ENDIF 593 594 !------------------------------------------------------------- 595 ! Compute the horizontal interpolation for every profile level 596 !------------------------------------------------------------- 597 598 DO ikn=1,inum_obs 599 iend=ista+ikn-1 600 601 zweig(:,:,1) = 0._wp 602 603 ! This code forces the horizontal weights to be 604 ! zero IF the observation is below the bottom of the 605 ! corners of the interpolation nodes, Or if it is in 606 ! the mask. This is important for observations near 607 ! steep bathymetry 608 DO iin=1,2 609 DO ijn=1,2 610 611 depth_loop2: DO ik=kpk,2,-1 612 IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN 613 614 zweig(iin,ijn,1) = & 615 & zweig2(iin,ijn,1) * & 616 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 617 & - prodatqc%var(2)%vdep(iend)),0._wp) 618 619 EXIT depth_loop2 620 621 ENDIF 622 623 ENDDO depth_loop2 624 625 ENDDO 626 ENDDO 627 628 CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 629 & prodatqc%var(2)%vmod(iend:iend) ) 630 631 ! Set QC flag for any observations found below the bottom 632 ! needed as the check here is more strict than that in obs_prep 633 IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 634 635 ENDDO 636 637 DEALLOCATE(interp_corner,iv_indic) 638 639 ENDIF 461 ENDIF 640 462 641 463 ENDDO 642 464 643 465 ! Deallocate the data for interpolation 644 DEALLOCATE( & 645 & igrdi1, & 646 & igrdi2, & 647 & igrdj1, & 648 & igrdj2, & 649 & zglam1, & 650 & zglam2, & 651 & zgphi1, & 652 & zgphi2, & 653 & zmask1, & 654 & zmask2, & 655 & zint1, & 656 & zint2, & 466 DEALLOCATE( & 467 & igrdi, & 468 & igrdj, & 469 & zglam, & 470 & zgphi, & 471 & zmask, & 472 & zint, & 657 473 & zgdept, & 658 474 & zgdepw & … … 661 477 ! At the end of the day also get interpolated means 662 478 IF ( ld_dailyav .AND. idayend == 0 ) THEN 663 DEALLOCATE( & 664 & zinm1, & 665 & zinm2 & 666 & ) 479 DEALLOCATE( zinm ) 667 480 ENDIF 668 481 669 prodatqc%nprofup = prodatqc%nprofup + ipro 482 IF ( kvar == prodatqc%nvar ) THEN 483 prodatqc%nprofup = prodatqc%nprofup + ipro 484 ENDIF 670 485 671 486 END SUBROUTINE obs_prof_opt
Note: See TracChangeset
for help on using the changeset viewer.
