Changeset 4746 for branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP
- Timestamp:
- 2014-08-15T11:26:21+02:00 (10 years ago)
- Location:
- branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS
- Files:
-
- 2 edited
-
diaobs.F90 (modified) (5 diffs)
-
obs_oper.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
r4624 r4746 118 118 !!---------------------------------------------------------------------- 119 119 120 !! * Substitutions 121 # include "domzgr_substitute.h90" 120 122 CONTAINS 121 123 … … 207 209 ! Read namelist parameters 208 210 !----------------------------------------------------------------------- 209 211 212 !Initalise all values in namelist arrays 213 enactfiles(:) = '' 214 coriofiles(:) = '' 215 profbfiles(:) = '' 216 slafilesact(:) = '' 217 slafilespas(:) = '' 218 slafbfiles(:) = '' 219 sstfiles(:) = '' 220 sstfbfiles(:) = '' 221 seaicefiles(:) = '' 210 222 velcurfiles(:) = '' 211 223 veladcpfiles(:) = '' 224 velavcurfiles(:) = '' 225 velhrcurfiles(:) = '' 226 velavadcpfiles(:) = '' 227 velhradcpfiles(:) = '' 228 velfbfiles(:) = '' 229 velcurfiles(:) = '' 230 veladcpfiles(:) = '' 231 endailyavtypes(:) = -1 232 endailyavtypes(1) = 820 233 ln_profb_ena(:) = .FALSE. 234 ln_profb_enatim(:) = .TRUE. 235 ln_velfb_av(:) = .FALSE. 236 ln_ignmis = .FALSE. 212 237 CALL ini_date( dobsini ) 213 238 CALL fin_date( dobsend ) … … 972 997 !! ! 07-04 (G. Smith) Generalized surface operators 973 998 !! ! 08-10 (M. Valdivieso) obs operator for velocity profiles 999 !! ! 14-08 (J. While) observation operator for profiles in 1000 !! generalised vertical coordinates 974 1001 !!---------------------------------------------------------------------- 975 1002 !! * Modules used … … 977 1004 & rdt, & 978 1005 & gdept_1d, & 1006 #if defined key_vvl 1007 & gdept_n, & 1008 #else 1009 & gdept_1d, & 1010 #endif 979 1011 & tmask, umask, vmask 980 1012 USE phycst, ONLY : & ! Physical constants … … 1034 1066 IF ( ln_t3d .OR. ln_s3d ) THEN 1035 1067 DO jprofset = 1, nprofsets 1036 IF ( ld_enact(jprofset) ) THEN 1037 CALL obs_pro_opt( prodatqc(jprofset), & 1038 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1039 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1040 & gdept_1d, tmask, n1dint, n2dint, & 1041 & kdailyavtypes = endailyavtypes ) 1068 IF( ln_zco .OR. ln_zps ) THEN 1069 IF ( ld_enact(jprofset) ) THEN 1070 CALL obs_pro_opt( prodatqc(jprofset), & 1071 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1072 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1073 & gdept_1d, tmask, n1dint, n2dint, & 1074 & kdailyavtypes = endailyavtypes ) 1075 ELSE 1076 CALL obs_pro_opt( prodatqc(jprofset), & 1077 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1078 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1079 & gdept_1d, tmask, n1dint, n2dint ) 1080 ENDIF 1042 1081 ELSE 1043 CALL obs_pro_opt( prodatqc(jprofset), & 1044 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1045 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1046 & gdept_1d, tmask, n1dint, n2dint ) 1082 IF ( ld_enact(jprofset) ) THEN 1083 CALL obs_pro_sco_opt( prodatqc(jprofset), & 1084 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1085 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1086 & fsdept(:,:,:), tmask, n1dint, n2dint, & 1087 & kdailyavtypes = endailyavtypes ) 1088 ELSE 1089 CALL obs_pro_sco_opt( prodatqc(jprofset), & 1090 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1091 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1092 & fsdept(:,:,:), tmask, n1dint, n2dint ) 1093 ENDIF 1047 1094 ENDIF 1048 1095 END DO -
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
r4245 r4746 9 9 !! obs_pro_opt : Compute the model counterpart of temperature and 10 10 !! salinity observations from profiles 11 !! obs_pro_sco_opt: Compute the model counterpart of temperature and 12 !! salinity observations from profiles in generalised 13 !! vertical coordinates 11 14 !! obs_sla_opt : Compute the model counterpart of sea level anomaly 12 15 !! observations … … 37 40 USE dom_oce, ONLY : & 38 41 & glamt, glamu, glamv, & 39 & gphit, gphiu, gphiv 42 & gphit, gphiu, gphiv, & 43 #if defined key_vvl 44 & gdept_n 45 #else 46 & gdept_0 47 #endif 40 48 USE lib_mpp, ONLY : & 41 49 & ctl_warn, ctl_stop 42 50 USE obs_grid, ONLY : & 51 & obs_level_search 52 43 53 IMPLICIT NONE 44 54 … … 47 57 48 58 PUBLIC obs_pro_opt, & ! Compute the model counterpart of profile observations 59 & obs_pro_sco_opt, & ! Compute the model counterpart of profile observations 60 ! in generalised vertical coordinates 49 61 & obs_sla_opt, & ! Compute the model counterpart of SLA observations 50 62 & obs_sst_opt, & ! Compute the model counterpart of SST observations … … 61 73 !!---------------------------------------------------------------------- 62 74 75 !! * Substitutions 76 # include "domzgr_substitute.h90" 63 77 CONTAINS 64 78 … … 449 463 END SUBROUTINE obs_pro_opt 450 464 465 SUBROUTINE obs_pro_sco_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 466 & ptn, psn, pgdept, ptmask, k1dint, k2dint, & 467 & kdailyavtypes ) 468 !!----------------------------------------------------------------------- 469 !! 470 !! *** ROUTINE obs_pro_opt *** 471 !! 472 !! ** Purpose : Compute the model counterpart of profiles 473 !! data by interpolating from the model grid to the 474 !! observation point. Generalised vertical coordinate version 475 !! 476 !! ** Method : Linearly interpolate to each observation point using 477 !! the model values at the corners of the surrounding grid box. 478 !! 479 !! First, model values on the model grid are interpolated vertically to the 480 !! Depths of the profile observations. Two vertical interpolation schemes are 481 !! available: 482 !! - linear (k1dint = 0) 483 !! - Cubic spline (k1dint = 1) 484 !! 485 !! 486 !! Secondly the interpolated values are interpolated horizontally to the 487 !! obs (lon, lat) point. 488 !! Several horizontal interpolation schemes are available: 489 !! - distance-weighted (great circle) (k2dint = 0) 490 !! - distance-weighted (small angle) (k2dint = 1) 491 !! - bilinear (geographical grid) (k2dint = 2) 492 !! - bilinear (quadrilateral grid) (k2dint = 3) 493 !! - polynomial (quadrilateral grid) (k2dint = 4) 494 !! 495 !! For the cubic spline the 2nd derivative of the interpolating 496 !! polynomial is computed before entering the vertical interpolation 497 !! routine. 498 !! 499 !! For ENACT moored buoy data (e.g., TAO), the model equivalent is 500 !! a daily mean model temperature field. So, we first compute 501 !! the mean, then interpolate only at the end of the day. 502 !! 503 !! This is the procedure to be used with generalised vertical model 504 !! coordinates (ie s-coordinates. It is ~4x slower than the equivalent 505 !! horizontal then vertical interpolation algorithm, but can deal with situations 506 !! where the model levels are not flat. 507 !! ONLY PERFORMED if ln_sco=.TRUE. 508 !! 509 !! Note: the in situ temperature observations must be converted 510 !! to potential temperature (the model variable) prior to 511 !! assimilation. 512 !!?????????????????????????????????????????????????????????????? 513 !! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR??? 514 !!?????????????????????????????????????????????????????????????? 515 !! 516 !! ** Action : 517 !! 518 !! History : 519 !! ! 2014-08 (J. While) Adapted from obs_pro_opt to handel generalised 520 !! vertical coordinates 521 !!----------------------------------------------------------------------- 522 523 !! * Modules used 524 USE obs_profiles_def ! Definition of storage space for profile obs. 525 USE dom_oce, ONLY : & 526 #if defined key_vvl 527 gdepw_n 528 #else 529 gdepw_0 530 #endif 531 532 IMPLICIT NONE 533 534 !! * Arguments 535 TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening 536 INTEGER, INTENT(IN) :: kt ! Time step 537 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 538 INTEGER, INTENT(IN) :: kpj 539 INTEGER, INTENT(IN) :: kpk 540 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 541 ! (kit000-1 = restart time) 542 INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header) 543 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 544 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 545 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 546 & ptn, & ! Model temperature field 547 & psn, & ! Model salinity field 548 & ptmask ! Land-sea mask 549 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,jpj,kpk) :: & 550 & pgdept ! Model array of depth levels 551 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 552 & kdailyavtypes ! Types for daily averages 553 !! * Local declarations 554 INTEGER :: ji 555 INTEGER :: jj 556 INTEGER :: jk 557 INTEGER :: iico, ijco 558 INTEGER :: jobs 559 INTEGER :: inrc 560 INTEGER :: ipro 561 INTEGER :: idayend 562 INTEGER :: ista 563 INTEGER :: iend 564 INTEGER :: iobs 565 INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes 566 INTEGER, DIMENSION(imaxavtypes) :: & 567 & idailyavtypes 568 REAL(KIND=wp) :: zlam 569 REAL(KIND=wp) :: zphi 570 REAL(KIND=wp) :: zdaystp 571 REAL(KIND=wp), DIMENSION(kpk) :: & 572 & zobsmask, & 573 & zobsk, & 574 & zobs2k 575 REAL(KIND=wp), DIMENSION(2,2,1) :: & 576 & zweig, & 577 & l_zweig 578 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 579 & zmask, & 580 & zintt, & 581 & zints, & 582 & zinmt, & 583 & zgdept,& 584 & zgdepw,& 585 & zinms 586 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 587 & zglam, & 588 & zgphi 589 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 590 & igrdi, & 591 & igrdj 592 INTEGER :: & 593 & inum_obs 594 REAL(KIND=wp), DIMENSION(1) :: zmsk_1 595 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 596 INTEGER, ALLOCATABLE, DIMENSION(:) :: v_indic 597 598 !------------------------------------------------------------------------ 599 ! Local initialization 600 !------------------------------------------------------------------------ 601 ! ... Record and data counters 602 inrc = kt - kit000 + 2 603 ipro = prodatqc%npstp(inrc) 604 605 ! Daily average types 606 IF ( PRESENT(kdailyavtypes) ) THEN 607 idailyavtypes(:) = kdailyavtypes(:) 608 ELSE 609 idailyavtypes(:) = -1 610 ENDIF 611 612 ! Initialize daily mean for first time-step 613 idayend = MOD( kt - kit000 + 1, kdaystp ) 614 615 ! Added kt == 0 test to catch restart case 616 IF ( idayend == 1 .OR. kt == 0) THEN 617 618 IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 619 DO jk = 1, jpk 620 DO jj = 1, jpj 621 DO ji = 1, jpi 622 prodatqc%vdmean(ji,jj,jk,1) = 0.0 623 prodatqc%vdmean(ji,jj,jk,2) = 0.0 624 END DO 625 END DO 626 END DO 627 628 ENDIF 629 630 DO jk = 1, jpk 631 DO jj = 1, jpj 632 DO ji = 1, jpi 633 ! Increment the temperature field for computing daily mean 634 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 635 & + ptn(ji,jj,jk) 636 ! Increment the salinity field for computing daily mean 637 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 638 & + psn(ji,jj,jk) 639 END DO 640 END DO 641 END DO 642 643 ! Compute the daily mean at the end of day 644 zdaystp = 1.0 / REAL( kdaystp ) 645 IF ( idayend == 0 ) THEN 646 DO jk = 1, jpk 647 DO jj = 1, jpj 648 DO ji = 1, jpi 649 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 650 & * zdaystp 651 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 652 & * zdaystp 653 END DO 654 END DO 655 END DO 656 ENDIF 657 658 ! Get the data for interpolation 659 ALLOCATE( & 660 & igrdi(2,2,ipro), & 661 & igrdj(2,2,ipro), & 662 & zglam(2,2,ipro), & 663 & zgphi(2,2,ipro), & 664 & zmask(2,2,kpk,ipro), & 665 & zintt(2,2,kpk,ipro), & 666 & zints(2,2,kpk,ipro), & 667 & zgdept(2,2,kpk,ipro), & 668 & zgdepw(2,2,kpk,ipro) & 669 & ) 670 671 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 672 iobs = jobs - prodatqc%nprofup 673 igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1 674 igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1 675 igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1 676 igrdj(1,2,iobs) = prodatqc%mj(jobs,1) 677 igrdi(2,1,iobs) = prodatqc%mi(jobs,1) 678 igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1 679 igrdi(2,2,iobs) = prodatqc%mi(jobs,1) 680 igrdj(2,2,iobs) = prodatqc%mj(jobs,1) 681 END DO 682 683 CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam ) 684 CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi ) 685 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask ) 686 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn, zintt ) 687 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn, zints ) 688 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, fsdept(:,:,:), & 689 & zgdept ) 690 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, fsdepw(:,:,:), & 691 & zgdepw ) 692 693 ! At the end of the day also get interpolated means 694 IF ( idayend == 0 ) THEN 695 696 ALLOCATE( & 697 & zinmt(2,2,kpk,ipro), & 698 & zinms(2,2,kpk,ipro) & 699 & ) 700 701 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, & 702 & prodatqc%vdmean(:,:,:,1), zinmt ) 703 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, & 704 & prodatqc%vdmean(:,:,:,2), zinms ) 705 706 ENDIF 707 708 ! Return if no observations to process 709 ! Has to be done after comm commands to ensure processors 710 ! stay in sync 711 IF ( ipro == 0 ) RETURN 712 713 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 714 715 iobs = jobs - prodatqc%nprofup 716 717 IF ( kt /= prodatqc%mstp(jobs) ) THEN 718 719 IF(lwp) THEN 720 WRITE(numout,*) 721 WRITE(numout,*) ' E R R O R : Observation', & 722 & ' time step is not consistent with the', & 723 & ' model time step' 724 WRITE(numout,*) ' =========' 725 WRITE(numout,*) 726 WRITE(numout,*) ' Record = ', jobs, & 727 & ' kt = ', kt, & 728 & ' mstp = ', prodatqc%mstp(jobs), & 729 & ' ntyp = ', prodatqc%ntyp(jobs) 730 ENDIF 731 CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 732 ENDIF 733 734 zlam = prodatqc%rlam(jobs) 735 zphi = prodatqc%rphi(jobs) 736 737 ! Horizontal weights 738 ! Only calculated once, for both T and S. 739 ! Masked values are calculated later. 740 741 IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 742 & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 743 744 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 745 & zglam(:,:,iobs), zgphi(:,:,iobs), & 746 & zmask(:,:,1,iobs), zweig, zmsk_1 ) 747 748 ENDIF 749 750 ! IF zmsk_1 = 0; then ob is on land 751 IF (zmsk_1(1) < 0.1) THEN 752 WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask' 753 754 ELSE 755 756 ! Temperature 757 758 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 759 760 zobsk(:) = obfillflt 761 762 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 763 764 IF ( idayend == 0 ) THEN 765 766 ! Daily averaged moored buoy (MRB) data 767 768 ! vertically interpolate all 4 corners 769 ista = prodatqc%npvsta(jobs,1) 770 iend = prodatqc%npvend(jobs,1) 771 inum_obs = iend - ista + 1 772 ALLOCATE(interp_corner(2,2,inum_obs),v_indic(inum_obs)) 773 774 DO iin=1,2 775 DO ijn=1,2 776 777 778 779 IF ( k1dint == 1 ) THEN 780 CALL obs_int_z1d_spl( kpk, & 781 & zinmt(iin,ijn,:,jobs), & 782 & zobs2k, zgdept(iin,ijn,:,jobs), & 783 & zmask(iin,ijn,:,jobs)) 784 ENDIF 785 786 CALL obs_level_search(kpk, & 787 & zgdept(iin,ijn,:,jobs), & 788 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 789 & v_indic) 790 CALL obs_int_z1d(kpk, v_indic, k1dint, inum_obs, & 791 & prodatqc%var(1)%vdep(ista:iend), & 792 & zinmt(iin,ijn,:,jobs), & 793 & zobs2k, interp_corner(iin,ijn,:), & 794 & zgdept(iin,ijn,:,jobs), & 795 & zmask(iin,ijn,:,jobs)) 796 797 ENDDO 798 ENDDO 799 800 801 ELSE 802 803 CALL ctl_stop( ' A nonzero' // & 804 & ' number of profile T BUOY data should' // & 805 & ' only occur at the end of a given day' ) 806 807 ENDIF 808 809 ELSE 810 811 ! Point data 812 813 ! vertically interpolate all 4 corners 814 ista = prodatqc%npvsta(jobs,1) 815 iend = prodatqc%npvend(jobs,1) 816 inum_obs = iend - ista + 1 817 ALLOCATE(interp_corner(2,2,inum_obs), v_indic(inum_obs)) 818 DO iin=1,2 819 DO ijn=1,2 820 821 822 IF ( k1dint == 1 ) THEN 823 CALL obs_int_z1d_spl( kpk, & 824 & zintt(iin,ijn,:,jobs),& 825 & zobs2k, zgdept(iin,ijn,:,jobs), & 826 & zmask(iin,ijn,:,jobs)) 827 828 ENDIF 829 830 CALL obs_level_search(kpk, & 831 & zgdept(iin,ijn,:,jobs),& 832 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 833 & v_indic) 834 CALL obs_int_z1d(kpk, v_indic, k1dint, inum_obs, & 835 & prodatqc%var(1)%vdep(ista:iend), & 836 & zintt(iin,ijn,:,jobs), & 837 & zobs2k,interp_corner(iin,ijn,:), & 838 & zgdept(iin,ijn,:,jobs), & 839 & zmask(iin,ijn,:,jobs) ) 840 841 ENDDO 842 ENDDO 843 844 ENDIF 845 846 !------------------------------------------------------------- 847 ! Compute the horizontal interpolation for every profile level 848 !------------------------------------------------------------- 849 850 DO ikn=1,inum_obs 851 iend=ista+ikn-1 852 853 ! This code forces the horizontal weights to be 854 ! zero IF the observation is below the bottom of the 855 ! corners of the interpolation nodes, Or if it is in 856 ! the mask. This is important for observations are near 857 ! steep bathymetry 858 DO iin=1,2 859 DO ijn=1,2 860 861 depth_loop1: DO ik=kpk,2,-1 862 IF(zmask(iin,ijn,ik-1,jobs ) > 0.9 )THEN 863 864 l_zweig(iin,ijn,1) = & 865 & zweig(iin,ijn,1) * & 866 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,jobs) ) & 867 & - prodatqc%var(1)%vdep(iend)),0._wp) 868 869 EXIT depth_loop1 870 ENDIF 871 ENDDO depth_loop1 872 873 ENDDO 874 ENDDO 875 876 CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 877 & prodatqc%var(1)%vmod(iend:iend) ) 878 879 ENDDO 880 881 882 DEALLOCATE(interp_corner,v_indic) 883 884 ENDIF 885 886 887 ! Salinity 888 889 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 890 891 zobsk(:) = obfillflt 892 893 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 894 895 IF ( idayend == 0 ) THEN 896 897 ! Daily averaged moored buoy (MRB) data 898 899 ! vertically interpolate all 4 corners 900 ista = prodatqc%npvsta(jobs,2) 901 iend = prodatqc%npvend(jobs,2) 902 inum_obs = iend - ista + 1 903 ALLOCATE(interp_corner(2,2,inum_obs),v_indic(inum_obs)) 904 905 DO iin=1,2 906 DO ijn=1,2 907 908 909 910 IF ( k1dint == 1 ) THEN 911 CALL obs_int_z1d_spl( kpk, & 912 & zinms(iin,ijn,:,jobs), & 913 & zobs2k, zgdept(iin,ijn,:,jobs), & 914 & zmask(iin,ijn,:,jobs)) 915 ENDIF 916 917 CALL obs_level_search(kpk, & 918 & zgdept(iin,ijn,:,jobs), & 919 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 920 & v_indic) 921 CALL obs_int_z1d(kpk, v_indic, k1dint, inum_obs, & 922 & prodatqc%var(2)%vdep(ista:iend), & 923 & zinms(iin,ijn,:,jobs), & 924 & zobs2k, interp_corner(iin,ijn,:), & 925 & zgdept(iin,ijn,:,jobs), & 926 & zmask(iin,ijn,:,jobs)) 927 928 ENDDO 929 ENDDO 930 931 932 ELSE 933 934 CALL ctl_stop( ' A nonzero' // & 935 & ' number of profile T BUOY data should' // & 936 & ' only occur at the end of a given day' ) 937 938 ENDIF 939 940 ELSE 941 942 ! Point data 943 944 ! vertically interpolate all 4 corners 945 ista = prodatqc%npvsta(jobs,2) 946 iend = prodatqc%npvend(jobs,2) 947 inum_obs = iend - ista + 1 948 ALLOCATE(interp_corner(2,2,inum_obs), v_indic(inum_obs)) 949 950 DO iin=1,2 951 DO ijn=1,2 952 953 954 IF ( k1dint == 1 ) THEN 955 CALL obs_int_z1d_spl( kpk, & 956 & zints(iin,ijn,:,jobs),& 957 & zobs2k, zgdept(iin,ijn,:,jobs), & 958 & zmask(iin,ijn,:,jobs)) 959 960 ENDIF 961 962 CALL obs_level_search(kpk, & 963 & zgdept(iin,ijn,:,jobs),& 964 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 965 & v_indic) 966 CALL obs_int_z1d(kpk, v_indic, k1dint, inum_obs, & 967 & prodatqc%var(2)%vdep(ista:iend), & 968 & zints(iin,ijn,:,jobs), & 969 & zobs2k,interp_corner(iin,ijn,:), & 970 & zgdept(iin,ijn,:,jobs), & 971 & zmask(iin,ijn,:,jobs) ) 972 973 ENDDO 974 ENDDO 975 976 ENDIF 977 978 !------------------------------------------------------------- 979 ! Compute the horizontal interpolation for every profile level 980 !------------------------------------------------------------- 981 982 DO ikn=1,inum_obs 983 iend=ista+ikn-1 984 985 ! This code forces the horizontal weights to be 986 ! zero IF the observation is below the bottom of the 987 ! corners of the interpolation nodes, Or if it is in 988 ! the mask. This is important for observations are near 989 ! steep bathymetry 990 DO iin=1,2 991 DO ijn=1,2 992 993 depth_loop2: DO ik=kpk,2,-1 994 IF(zmask(iin,ijn,ik-1,jobs ) > 0.9 )THEN 995 996 l_zweig(iin,ijn,1) = & 997 & zweig(iin,ijn,1) * & 998 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,jobs) ) & 999 & - prodatqc%var(2)%vdep(iend)),0._wp) 1000 1001 EXIT depth_loop2 1002 ENDIF 1003 ENDDO depth_loop2 1004 1005 ENDDO 1006 ENDDO 1007 1008 CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 1009 & prodatqc%var(2)%vmod(iend:iend) ) 1010 1011 ENDDO 1012 1013 1014 DEALLOCATE(interp_corner,v_indic) 1015 1016 ENDIF 1017 1018 ENDIF 1019 1020 END DO 1021 1022 ! Deallocate the data for interpolation 1023 DEALLOCATE( & 1024 & igrdi, & 1025 & igrdj, & 1026 & zglam, & 1027 & zgphi, & 1028 & zmask, & 1029 & zintt, & 1030 & zints & 1031 & ) 1032 ! At the end of the day also get interpolated means 1033 IF ( idayend == 0 ) THEN 1034 DEALLOCATE( & 1035 & zinmt, & 1036 & zinms & 1037 & ) 1038 ENDIF 1039 1040 prodatqc%nprofup = prodatqc%nprofup + ipro 1041 1042 END SUBROUTINE obs_pro_sco_opt 1043 451 1044 SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, & 452 1045 & psshn, psshmask, k2dint ) … … 1293 1886 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 1294 1887 & zglamv(:,:,iobs), zgphiv(:,:,iobs), & 1295 & zvmask(:,:,:,iobs), zweigv, zobsmask u)1888 & zvmask(:,:,:,iobs), zweigv, zobsmaskv ) 1296 1889 1297 1890 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
