- 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_oper.F90
r7646 r9023 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, & 26 27 & obs_int_z1d_spl 27 USE obs_const, ONLY : &28 & obfillflt ! Fillvalue28 USE obs_const, ONLY : & ! Obs fill value 29 & obfillflt 29 30 USE dom_oce, ONLY : & 30 & glamt, glamu, glamv, & 31 & gphit, gphiu, gphiv, & 32 & gdept_n, gdept_0 33 USE lib_mpp, ONLY : & 31 & glamt, glamf, & 32 & gphit, gphif 33 USE lib_mpp, ONLY : & ! Warning and stopping routines 34 34 & ctl_warn, ctl_stop 35 USE sbcdcy, ONLY : & ! For calculation of where it is night-time 36 & sbc_dcy, nday_qsr 35 37 USE obs_grid, ONLY : & 36 38 & obs_level_search 37 USE sbcdcy, ONLY : & ! For calculation of where it is night-time38 & sbc_dcy, nday_qsr39 39 40 40 IMPLICIT NONE … … 44 44 45 45 PUBLIC obs_prof_opt, & ! Compute the model counterpart of profile obs 46 & obs_pro_sco_opt, & ! Compute the model counterpart of profile observations47 46 & obs_surf_opt ! Compute the model counterpart of surface obs 48 47 … … 58 57 CONTAINS 59 58 59 60 60 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 61 61 & kit000, kdaystp, & 62 & pvar1, pvar2, pgdept, pmask1, pmask2, & 62 & pvar1, pvar2, pgdept, pgdepw, & 63 & pmask1, pmask2, & 63 64 & plam1, plam2, pphi1, pphi2, & 64 65 & k1dint, k2dint, kdailyavtypes ) … … 111 112 !! ! 07-03 (K. Mogensen) General handling of profiles 112 113 !! ! 15-02 (M. Martin) Combined routine for all profile types 114 !! ! 17-02 (M. Martin) Include generalised vertical coordinate changes 113 115 !!----------------------------------------------------------------------- 114 116 … … 140 142 & pphi1, & ! Model latitudes for variable 1 141 143 & pphi2 ! Model latitudes for variable 2 142 REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 143 & pgdept ! Model array of depth levels 144 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 145 & pgdept, & ! Model array of depth T levels 146 & pgdepw ! Model array of depth W levels 144 147 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 145 148 & kdailyavtypes ! Types for daily averages … … 156 159 INTEGER :: iend 157 160 INTEGER :: iobs 161 INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes 162 INTEGER :: inum_obs 158 163 INTEGER, DIMENSION(imaxavtypes) :: & 159 164 & idailyavtypes … … 163 168 & igrdj1, & 164 169 & igrdj2 170 INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 171 165 172 REAL(KIND=wp) :: zlam 166 173 REAL(KIND=wp) :: zphi … … 171 178 & zobsk, & 172 179 & zobs2k 173 REAL(KIND=wp), DIMENSION(2,2, kpk) :: &180 REAL(KIND=wp), DIMENSION(2,2,1) :: & 174 181 & zweig1, & 175 & zweig2 182 & zweig2, & 183 & zweig 176 184 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 177 185 & zmask1, & 178 186 & zmask2, & 179 & zint1, & 180 & zint2, & 181 & zinm1, & 182 & zinm2 187 & zint1, & 188 & zint2, & 189 & zinm1, & 190 & zinm2, & 191 & zgdept, & 192 & zgdepw 183 193 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 184 194 & zglam1, & … … 186 196 & zgphi1, & 187 197 & zgphi2 198 REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2 199 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 200 188 201 LOGICAL :: ld_dailyav 189 202 … … 266 279 & zmask1(2,2,kpk,ipro), & 267 280 & zmask2(2,2,kpk,ipro), & 268 & zint1(2,2,kpk,ipro), & 269 & zint2(2,2,kpk,ipro) & 281 & zint1(2,2,kpk,ipro), & 282 & zint2(2,2,kpk,ipro), & 283 & zgdept(2,2,kpk,ipro), & 284 & zgdepw(2,2,kpk,ipro) & 270 285 & ) 271 286 … … 290 305 END DO 291 306 307 ! Initialise depth arrays 308 zgdept(:,:,:,:) = 0.0 309 zgdepw(:,:,:,:) = 0.0 310 292 311 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 ) 293 312 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 ) … … 300 319 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2, zint2 ) 301 320 321 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept ) 322 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw ) 323 302 324 ! At the end of the day also get interpolated means 303 325 IF ( ld_dailyav .AND. idayend == 0 ) THEN … … 314 336 315 337 ENDIF 338 339 ! Return if no observations to process 340 ! Has to be done after comm commands to ensure processors 341 ! stay in sync 342 IF ( ipro == 0 ) RETURN 316 343 317 344 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro … … 339 366 zphi = prodatqc%rphi(jobs) 340 367 341 ! Horizontal weights and vertical mask342 368 ! Horizontal weights 369 ! Masked values are calculated later. 343 370 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 344 371 345 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &372 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 346 373 & zglam1(:,:,iobs), zgphi1(:,:,iobs), & 347 & zmask1(:,:, :,iobs), zweig1, zobsmask1 )374 & zmask1(:,:,1,iobs), zweig1, zmsk_1 ) 348 375 349 376 ENDIF … … 351 378 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 352 379 353 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &380 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 354 381 & zglam2(:,:,iobs), zgphi2(:,:,iobs), & 355 & zmask2(:,:, :,iobs), zweig2, zobsmask2 )382 & zmask2(:,:,1,iobs), zweig2, zmsk_2 ) 356 383 357 384 ENDIF … … 365 392 IF ( idayend == 0 ) THEN 366 393 ! Daily averaged data 367 CALL obs_int_h2d( kpk, kpk, & 368 & zweig1, zinm1(:,:,:,iobs), zobsk ) 369 370 ENDIF 371 372 ELSE 373 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 394 395 ! vertically interpolate all 4 corners 396 ista = prodatqc%npvsta(jobs,1) 397 iend = prodatqc%npvend(jobs,1) 398 inum_obs = iend - ista + 1 399 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 400 401 DO iin=1,2 402 DO ijn=1,2 403 404 IF ( k1dint == 1 ) THEN 405 CALL obs_int_z1d_spl( kpk, & 406 & zinm1(iin,ijn,:,iobs), & 407 & zobs2k, zgdept(iin,ijn,:,iobs), & 408 & zmask1(iin,ijn,:,iobs)) 409 ENDIF 410 411 CALL obs_level_search(kpk, & 412 & zgdept(iin,ijn,:,iobs), & 413 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 414 & iv_indic) 415 416 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 417 & prodatqc%var(1)%vdep(ista:iend), & 418 & zinm1(iin,ijn,:,iobs), & 419 & zobs2k, interp_corner(iin,ijn,:), & 420 & zgdept(iin,ijn,:,iobs), & 421 & zmask1(iin,ijn,:,iobs)) 422 423 ENDDO 424 ENDDO 425 426 ENDIF !idayend 427 428 ELSE 429 430 ! Point data 431 432 ! vertically interpolate all 4 corners 433 ista = prodatqc%npvsta(jobs,1) 434 iend = prodatqc%npvend(jobs,1) 435 inum_obs = iend - ista + 1 436 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 437 DO iin=1,2 438 DO ijn=1,2 439 440 IF ( k1dint == 1 ) THEN 441 CALL obs_int_z1d_spl( kpk, & 442 & zint1(iin,ijn,:,iobs),& 443 & zobs2k, zgdept(iin,ijn,:,iobs), & 444 & zmask1(iin,ijn,:,iobs)) 445 446 ENDIF 447 448 CALL obs_level_search(kpk, & 449 & zgdept(iin,ijn,:,iobs),& 450 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 451 & iv_indic) 452 453 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 454 & prodatqc%var(1)%vdep(ista:iend), & 455 & zint1(iin,ijn,:,iobs), & 456 & zobs2k,interp_corner(iin,ijn,:), & 457 & zgdept(iin,ijn,:,iobs), & 458 & zmask1(iin,ijn,:,iobs) ) 459 460 ENDDO 461 ENDDO 462 463 ENDIF 464 465 !------------------------------------------------------------- 466 ! Compute the horizontal interpolation for every profile level 467 !------------------------------------------------------------- 468 469 DO ikn=1,inum_obs 470 iend=ista+ikn-1 471 472 zweig(:,:,1) = 0._wp 473 474 ! This code forces the horizontal weights to be 475 ! zero IF the observation is below the bottom of the 476 ! corners of the interpolation nodes, Or if it is in 477 ! the mask. This is important for observations near 478 ! steep bathymetry 479 DO iin=1,2 480 DO ijn=1,2 481 482 depth_loop1: DO ik=kpk,2,-1 483 IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN 484 485 zweig(iin,ijn,1) = & 486 & zweig1(iin,ijn,1) * & 487 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 488 & - prodatqc%var(1)%vdep(iend)),0._wp) 489 490 EXIT depth_loop1 491 492 ENDIF 493 494 ENDDO depth_loop1 495 496 ENDDO 497 ENDDO 498 499 CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 500 & prodatqc%var(1)%vmod(iend:iend) ) 501 502 ! Set QC flag for any observations found below the bottom 503 ! needed as the check here is more strict than that in obs_prep 504 IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4 505 506 ENDDO 507 508 DEALLOCATE(interp_corner,iv_indic) 509 510 ENDIF 511 512 ! For the second variable 405 513 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 406 514 … … 410 518 411 519 IF ( idayend == 0 ) THEN 412 413 520 ! 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 521 522 ! vertically interpolate all 4 corners 523 ista = prodatqc%npvsta(jobs,2) 524 iend = prodatqc%npvend(jobs,2) 525 inum_obs = iend - ista + 1 526 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 527 528 DO iin=1,2 529 DO ijn=1,2 530 531 IF ( k1dint == 1 ) THEN 532 CALL obs_int_z1d_spl( kpk, & 533 & zinm2(iin,ijn,:,iobs), & 534 & zobs2k, zgdept(iin,ijn,:,iobs), & 535 & zmask2(iin,ijn,:,iobs)) 536 ENDIF 537 538 CALL obs_level_search(kpk, & 539 & zgdept(iin,ijn,:,iobs), & 540 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 541 & iv_indic) 542 543 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 544 & prodatqc%var(2)%vdep(ista:iend), & 545 & zinm2(iin,ijn,:,iobs), & 546 & zobs2k, interp_corner(iin,ijn,:), & 547 & zgdept(iin,ijn,:,iobs), & 548 & zmask2(iin,ijn,:,iobs)) 549 550 ENDDO 551 ENDDO 552 553 ENDIF !idayend 554 555 ELSE 556 557 ! Point data 558 559 ! vertically interpolate all 4 corners 560 ista = prodatqc%npvsta(jobs,2) 561 iend = prodatqc%npvend(jobs,2) 562 inum_obs = iend - ista + 1 563 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 564 DO iin=1,2 565 DO ijn=1,2 566 567 IF ( k1dint == 1 ) THEN 568 CALL obs_int_z1d_spl( kpk, & 569 & zint2(iin,ijn,:,iobs),& 570 & zobs2k, zgdept(iin,ijn,:,iobs), & 571 & zmask2(iin,ijn,:,iobs)) 572 573 ENDIF 574 575 CALL obs_level_search(kpk, & 576 & zgdept(iin,ijn,:,iobs),& 577 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 578 & iv_indic) 579 580 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 581 & prodatqc%var(2)%vdep(ista:iend), & 582 & zint2(iin,ijn,:,iobs), & 583 & zobs2k,interp_corner(iin,ijn,:), & 584 & zgdept(iin,ijn,:,iobs), & 585 & zmask2(iin,ijn,:,iobs) ) 586 587 ENDDO 588 ENDDO 589 590 ENDIF 591 592 !------------------------------------------------------------- 593 ! Compute the horizontal interpolation for every profile level 594 !------------------------------------------------------------- 595 596 DO ikn=1,inum_obs 597 iend=ista+ikn-1 598 599 zweig(:,:,1) = 0._wp 600 601 ! This code forces the horizontal weights to be 602 ! zero IF the observation is below the bottom of the 603 ! corners of the interpolation nodes, Or if it is in 604 ! the mask. This is important for observations near 605 ! steep bathymetry 606 DO iin=1,2 607 DO ijn=1,2 608 609 depth_loop2: DO ik=kpk,2,-1 610 IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN 611 612 zweig(iin,ijn,1) = & 613 & zweig2(iin,ijn,1) * & 614 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 615 & - prodatqc%var(2)%vdep(iend)),0._wp) 616 617 EXIT depth_loop2 618 619 ENDIF 620 621 ENDDO depth_loop2 622 623 ENDDO 624 ENDDO 625 626 CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 627 & prodatqc%var(2)%vmod(iend:iend) ) 628 629 ! Set QC flag for any observations found below the bottom 630 ! needed as the check here is more strict than that in obs_prep 631 IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 632 633 ENDDO 634 635 DEALLOCATE(interp_corner,iv_indic) 636 637 ENDIF 638 639 ENDDO 454 640 455 641 ! Deallocate the data for interpolation … … 466 652 & zmask2, & 467 653 & zint1, & 468 & zint2 & 654 & zint2, & 655 & zgdept, & 656 & zgdepw & 469 657 & ) 470 658 … … 481 669 END SUBROUTINE obs_prof_opt 482 670 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 ) 671 SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, & 672 & kit000, kdaystp, psurf, psurfmask, & 673 & k2dint, ldnightav, plamscl, pphiscl, & 674 & lindegrees ) 1071 675 1072 676 !!----------------------------------------------------------------------- … … 1090 694 !! - polynomial (quadrilateral grid) (k2dint = 4) 1091 695 !! 696 !! Two horizontal averaging schemes are also available: 697 !! - weighted radial footprint (k2dint = 5) 698 !! - weighted rectangular footprint (k2dint = 6) 699 !! 1092 700 !! 1093 701 !! ** Action : … … 1096 704 !! ! 07-03 (A. Weaver) 1097 705 !! ! 15-02 (M. Martin) Combined routine for surface types 706 !! ! 17-03 (M. Martin) Added horizontal averaging options 1098 707 !!----------------------------------------------------------------------- 1099 708 … … 1117 726 & psurfmask ! Land-sea mask 1118 727 LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 728 REAL(KIND=wp), INTENT(IN) :: & 729 & plamscl, & ! Diameter in metres of obs footprint in E/W, N/S directions 730 & pphiscl ! This is the full width (rather than half-width) 731 LOGICAL, INTENT(IN) :: & 732 & lindegrees ! T=> plamscl and pphiscl are specified in degrees, F=> in metres 1119 733 1120 734 !! * Local declarations … … 1125 739 INTEGER :: isurf 1126 740 INTEGER :: iobs 741 INTEGER :: imaxifp, imaxjfp 742 INTEGER :: imodi, imodj 1127 743 INTEGER :: idayend 1128 744 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 1129 & igrdi, & 1130 & igrdj 745 & igrdi, & 746 & igrdj, & 747 & igrdip1, & 748 & igrdjp1 1131 749 INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 1132 750 & icount_night, & … … 1136 754 REAL(wp), DIMENSION(1) :: zext, zobsmask 1137 755 REAL(wp) :: zdaystp 1138 REAL(wp), DIMENSION(2,2,1) :: &1139 & zweig1140 756 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 757 & zweig, & 1141 758 & zmask, & 1142 759 & zsurf, & 1143 760 & zsurfm, & 761 & zsurftmp, & 1144 762 & zglam, & 1145 & zgphi 763 & zgphi, & 764 & zglamf, & 765 & zgphif 766 1146 767 REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 1147 768 & zintmp, & … … 1155 776 inrc = kt - kit000 + 2 1156 777 isurf = surfdataqc%nsstp(inrc) 778 779 ! Work out the maximum footprint size for the 780 ! interpolation/averaging in model grid-points - has to be even. 781 782 CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 783 1157 784 1158 785 IF ( ldnightav ) THEN … … 1221 848 1222 849 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) & 850 & zweig(imaxifp,imaxjfp,1), & 851 & igrdi(imaxifp,imaxjfp,isurf), & 852 & igrdj(imaxifp,imaxjfp,isurf), & 853 & zglam(imaxifp,imaxjfp,isurf), & 854 & zgphi(imaxifp,imaxjfp,isurf), & 855 & zmask(imaxifp,imaxjfp,isurf), & 856 & zsurf(imaxifp,imaxjfp,isurf), & 857 & zsurftmp(imaxifp,imaxjfp,isurf), & 858 & zglamf(imaxifp+1,imaxjfp+1,isurf), & 859 & zgphif(imaxifp+1,imaxjfp+1,isurf), & 860 & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 861 & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 1229 862 & ) 1230 863 1231 864 DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 1232 865 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) 866 DO ji = 0, imaxifp 867 imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1 868 869 !Deal with wrap around in longitude 870 IF ( imodi < 1 ) imodi = imodi + jpiglo 871 IF ( imodi > jpiglo ) imodi = imodi - jpiglo 872 873 DO jj = 0, imaxjfp 874 imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1 875 !If model values are out of the domain to the north/south then 876 !set them to be the edge of the domain 877 IF ( imodj < 1 ) imodj = 1 878 IF ( imodj > jpjglo ) imodj = jpjglo 879 880 igrdip1(ji+1,jj+1,iobs) = imodi 881 igrdjp1(ji+1,jj+1,iobs) = imodj 882 883 IF ( ji >= 1 .AND. jj >= 1 ) THEN 884 igrdi(ji,jj,iobs) = imodi 885 igrdj(ji,jj,iobs) = imodj 886 ENDIF 887 888 END DO 889 END DO 1241 890 END DO 1242 891 1243 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &892 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 1244 893 & igrdi, igrdj, glamt, zglam ) 1245 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &894 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 1246 895 & igrdi, igrdj, gphit, zgphi ) 1247 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &896 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 1248 897 & igrdi, igrdj, psurfmask, zmask ) 1249 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &898 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 1250 899 & igrdi, igrdj, psurf, zsurf ) 900 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 901 & igrdip1, igrdjp1, glamf, zglamf ) 902 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 903 & igrdip1, igrdjp1, gphif, zgphif ) 1251 904 1252 905 ! 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 906 IF ( idayend == 0 .AND. ldnightav ) THEN 907 908 ALLOCATE( & 909 & zsurfm(imaxifp,imaxjfp,isurf) & 910 & ) 911 912 CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, & 913 & surfdataqc%vdmean(:,:), zsurfm ) 914 1264 915 ENDIF 1265 916 … … 1290 941 zphi = surfdataqc%rphi(jobs) 1291 942 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 943 IF ( ldnightav .AND. idayend == 0 ) THEN 1299 944 ! Night-time averaged data 1300 CALL obs_int_h2d( 1, 1, zweig, zsurfm(:,:,iobs), zext)945 zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 1301 946 ELSE 1302 CALL obs_int_h2d( 1, 1, zweig, zsurf(:,:,iobs), zext ) 947 zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 948 ENDIF 949 950 IF ( k2dint <= 4 ) THEN 951 952 ! Get weights to interpolate the model value to the observation point 953 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 954 & zglam(:,:,iobs), zgphi(:,:,iobs), & 955 & zmask(:,:,iobs), zweig, zobsmask ) 956 957 ! Interpolate the model value to the observation point 958 CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 959 960 ELSE 961 962 ! Get weights to average the model SLA to the observation footprint 963 CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam, zphi, & 964 & zglam(:,:,iobs), zgphi(:,:,iobs), & 965 & zglamf(:,:,iobs), zgphif(:,:,iobs), & 966 & zmask(:,:,iobs), plamscl, pphiscl, & 967 & lindegrees, zweig, zobsmask ) 968 969 ! Average the model SST to the observation footprint 970 CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 971 & zweig, zsurftmp(:,:,iobs), zext ) 972 1303 973 ENDIF 1304 974 … … 1310 980 surfdataqc%rmod(jobs,1) = zext(1) 1311 981 ENDIF 982 983 IF ( zext(1) == obfillflt ) THEN 984 ! If the observation value is a fill value, set QC flag to bad 985 surfdataqc%nqc(jobs) = 4 986 ENDIF 1312 987 1313 988 END DO … … 1315 990 ! Deallocate the data for interpolation 1316 991 DEALLOCATE( & 992 & zweig, & 1317 993 & igrdi, & 1318 994 & igrdj, & … … 1320 996 & zgphi, & 1321 997 & zmask, & 1322 & zsurf & 998 & zsurf, & 999 & zsurftmp, & 1000 & zglamf, & 1001 & zgphif, & 1002 & igrdip1,& 1003 & igrdjp1 & 1323 1004 & ) 1324 1005 1325 1006 ! 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 1007 IF ( idayend == 0 .AND. ldnightav ) THEN 1008 DEALLOCATE( & 1009 & zsurfm & 1010 & ) 1332 1011 ENDIF 1333 1012
Note: See TracChangeset
for help on using the changeset viewer.
