Changeset 5963
- Timestamp:
- 2015-12-01T16:06:40+01:00 (8 years ago)
- Location:
- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS
- Files:
-
- 3 edited
-
diaobs.F90 (modified) (6 diffs)
-
obs_oper.F90 (modified) (6 diffs)
-
obs_prep.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
r4990 r5963 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 210 213 enactfiles(:) = '' 211 214 coriofiles(:) = '' … … 232 235 ln_velfb_av(:) = .FALSE. 233 236 ln_ignmis = .FALSE. 234 237 235 238 CALL ini_date( dobsini ) 236 239 CALL fin_date( dobsend ) … … 991 994 !! ! 07-04 (G. Smith) Generalized surface operators 992 995 !! ! 08-10 (M. Valdivieso) obs operator for velocity profiles 996 !! ! 14-08 (J. While) observation operator for profiles in 997 !! generalised vertical coordinates 993 998 !!---------------------------------------------------------------------- 994 999 !! * Modules used … … 996 1001 & rdt, & 997 1002 & gdept_1d, & 1003 #if defined key_vvl 1004 & gdept_n, & 1005 #else 1006 & gdept_1d, & 1007 #endif 998 1008 & tmask, umask, vmask 999 1009 USE phycst, ONLY : & ! Physical constants … … 1053 1063 IF ( ln_t3d .OR. ln_s3d ) THEN 1054 1064 DO jprofset = 1, nprofsets 1055 IF ( ld_enact(jprofset) ) THEN 1056 CALL obs_pro_opt( prodatqc(jprofset), & 1057 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1058 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1059 & gdept_1d, tmask, n1dint, n2dint, & 1060 & kdailyavtypes = endailyavtypes ) 1065 IF( ln_zco .OR. ln_zps ) THEN 1066 IF ( ld_enact(jprofset) ) THEN 1067 CALL obs_pro_opt( prodatqc(jprofset), & 1068 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1069 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1070 & gdept_1d, tmask, n1dint, n2dint, & 1071 & kdailyavtypes = endailyavtypes ) 1072 ELSE 1073 CALL obs_pro_opt( prodatqc(jprofset), & 1074 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1075 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1076 & gdept_1d, tmask, n1dint, n2dint ) 1077 ENDIF 1061 1078 ELSE 1062 CALL obs_pro_opt( prodatqc(jprofset), & 1063 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1064 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1065 & gdept_1d, tmask, n1dint, n2dint ) 1079 IF ( ld_enact(jprofset) ) THEN 1080 CALL obs_pro_sco_opt( prodatqc(jprofset), & 1081 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1082 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1083 & fsdept(:,:,:), tmask, n1dint, n2dint, & 1084 & kdailyavtypes = endailyavtypes ) 1085 ELSE 1086 CALL obs_pro_sco_opt( prodatqc(jprofset), & 1087 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1088 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1089 & fsdept(:,:,:), tmask, n1dint, n2dint ) 1090 ENDIF 1066 1091 ENDIF 1067 1092 END DO -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
r4245 r5963 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 554 !! * Local declarations 555 INTEGER :: ji 556 INTEGER :: jj 557 INTEGER :: jk 558 INTEGER :: iico, ijco 559 INTEGER :: jobs 560 INTEGER :: inrc 561 INTEGER :: ipro 562 INTEGER :: idayend 563 INTEGER :: ista 564 INTEGER :: iend 565 INTEGER :: iobs 566 INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes 567 INTEGER, DIMENSION(imaxavtypes) :: & 568 & idailyavtypes 569 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 570 & igrdi, & 571 & igrdj 572 INTEGER :: & 573 & inum_obs 574 INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 575 REAL(KIND=wp) :: zlam 576 REAL(KIND=wp) :: zphi 577 REAL(KIND=wp) :: zdaystp 578 REAL(KIND=wp), DIMENSION(kpk) :: & 579 & zobsmask, & 580 & zobsk, & 581 & zobs2k 582 REAL(KIND=wp), DIMENSION(2,2,1) :: & 583 & zweig, & 584 & l_zweig 585 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 586 & zmask, & 587 & zintt, & 588 & zints, & 589 & zinmt, & 590 & zgdept,& 591 & zgdepw,& 592 & zinms 593 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 594 & zglam, & 595 & zgphi 596 REAL(KIND=wp), DIMENSION(1) :: zmsk_1 597 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 598 599 !------------------------------------------------------------------------ 600 ! Local initialization 601 !------------------------------------------------------------------------ 602 ! ... Record and data counters 603 inrc = kt - kit000 + 2 604 ipro = prodatqc%npstp(inrc) 605 606 ! Daily average types 607 IF ( PRESENT(kdailyavtypes) ) THEN 608 idailyavtypes(:) = kdailyavtypes(:) 609 ELSE 610 idailyavtypes(:) = -1 611 ENDIF 612 613 ! Initialize daily mean for first time-step 614 idayend = MOD( kt - kit000 + 1, kdaystp ) 615 616 ! Added kt == 0 test to catch restart case 617 IF ( idayend == 1 .OR. kt == 0) THEN 618 619 IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 620 DO jk = 1, jpk 621 DO jj = 1, jpj 622 DO ji = 1, jpi 623 prodatqc%vdmean(ji,jj,jk,1) = 0.0 624 prodatqc%vdmean(ji,jj,jk,2) = 0.0 625 END DO 626 END DO 627 END DO 628 629 ENDIF 630 631 DO jk = 1, jpk 632 DO jj = 1, jpj 633 DO ji = 1, jpi 634 ! Increment the temperature field for computing daily mean 635 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 636 & + ptn(ji,jj,jk) 637 ! Increment the salinity field for computing daily mean 638 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 639 & + psn(ji,jj,jk) 640 END DO 641 END DO 642 END DO 643 644 ! Compute the daily mean at the end of day 645 zdaystp = 1.0 / REAL( kdaystp ) 646 IF ( idayend == 0 ) THEN 647 DO jk = 1, jpk 648 DO jj = 1, jpj 649 DO ji = 1, jpi 650 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 651 & * zdaystp 652 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 653 & * zdaystp 654 END DO 655 END DO 656 END DO 657 ENDIF 658 659 ! Get the data for interpolation 660 ALLOCATE( & 661 & igrdi(2,2,ipro), & 662 & igrdj(2,2,ipro), & 663 & zglam(2,2,ipro), & 664 & zgphi(2,2,ipro), & 665 & zmask(2,2,kpk,ipro), & 666 & zintt(2,2,kpk,ipro), & 667 & zints(2,2,kpk,ipro), & 668 & zgdept(2,2,kpk,ipro), & 669 & zgdepw(2,2,kpk,ipro) & 670 & ) 671 672 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 673 iobs = jobs - prodatqc%nprofup 674 igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1 675 igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1 676 igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1 677 igrdj(1,2,iobs) = prodatqc%mj(jobs,1) 678 igrdi(2,1,iobs) = prodatqc%mi(jobs,1) 679 igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1 680 igrdi(2,2,iobs) = prodatqc%mi(jobs,1) 681 igrdj(2,2,iobs) = prodatqc%mj(jobs,1) 682 END DO 683 684 CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam ) 685 CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi ) 686 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask ) 687 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn, zintt ) 688 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn, zints ) 689 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, fsdept(:,:,:), & 690 & zgdept ) 691 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, fsdepw(:,:,:), & 692 & zgdepw ) 693 694 ! At the end of the day also get interpolated means 695 IF ( idayend == 0 ) THEN 696 697 ALLOCATE( & 698 & zinmt(2,2,kpk,ipro), & 699 & zinms(2,2,kpk,ipro) & 700 & ) 701 702 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, & 703 & prodatqc%vdmean(:,:,:,1), zinmt ) 704 CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, & 705 & prodatqc%vdmean(:,:,:,2), zinms ) 706 707 ENDIF 708 709 ! Return if no observations to process 710 ! Has to be done after comm commands to ensure processors 711 ! stay in sync 712 IF ( ipro == 0 ) RETURN 713 714 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 715 716 iobs = jobs - prodatqc%nprofup 717 718 IF ( kt /= prodatqc%mstp(jobs) ) THEN 719 720 IF(lwp) THEN 721 WRITE(numout,*) 722 WRITE(numout,*) ' E R R O R : Observation', & 723 & ' time step is not consistent with the', & 724 & ' model time step' 725 WRITE(numout,*) ' =========' 726 WRITE(numout,*) 727 WRITE(numout,*) ' Record = ', jobs, & 728 & ' kt = ', kt, & 729 & ' mstp = ', prodatqc%mstp(jobs), & 730 & ' ntyp = ', prodatqc%ntyp(jobs) 731 ENDIF 732 CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 733 ENDIF 734 735 zlam = prodatqc%rlam(jobs) 736 zphi = prodatqc%rphi(jobs) 737 738 ! Horizontal weights 739 ! Only calculated once, for both T and S. 740 ! Masked values are calculated later. 741 742 IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 743 & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 744 745 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 746 & zglam(:,:,iobs), zgphi(:,:,iobs), & 747 & zmask(:,:,1,iobs), zweig, zmsk_1 ) 748 749 ENDIF 750 751 ! IF zmsk_1 = 0; then ob is on land 752 IF (zmsk_1(1) < 0.1) THEN 753 WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask' 754 755 ELSE 756 757 ! Temperature 758 759 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 760 761 zobsk(:) = obfillflt 762 763 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 764 765 IF ( idayend == 0 ) THEN 766 767 ! Daily averaged moored buoy (MRB) data 768 769 ! vertically interpolate all 4 corners 770 ista = prodatqc%npvsta(jobs,1) 771 iend = prodatqc%npvend(jobs,1) 772 inum_obs = iend - ista + 1 773 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 774 775 DO iin=1,2 776 DO ijn=1,2 777 778 779 780 IF ( k1dint == 1 ) THEN 781 CALL obs_int_z1d_spl( kpk, & 782 & zinmt(iin,ijn,:,jobs), & 783 & zobs2k, zgdept(iin,ijn,:,jobs), & 784 & zmask(iin,ijn,:,jobs)) 785 ENDIF 786 787 CALL obs_level_search(kpk, & 788 & zgdept(iin,ijn,:,jobs), & 789 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 790 & iv_indic) 791 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 792 & prodatqc%var(1)%vdep(ista:iend), & 793 & zinmt(iin,ijn,:,jobs), & 794 & zobs2k, interp_corner(iin,ijn,:), & 795 & zgdept(iin,ijn,:,jobs), & 796 & zmask(iin,ijn,:,jobs)) 797 798 ENDDO 799 ENDDO 800 801 802 ELSE 803 804 CALL ctl_stop( ' A nonzero' // & 805 & ' number of profile T BUOY data should' // & 806 & ' only occur at the end of a given day' ) 807 808 ENDIF 809 810 ELSE 811 812 ! Point data 813 814 ! vertically interpolate all 4 corners 815 ista = prodatqc%npvsta(jobs,1) 816 iend = prodatqc%npvend(jobs,1) 817 inum_obs = iend - ista + 1 818 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 819 DO iin=1,2 820 DO ijn=1,2 821 822 823 IF ( k1dint == 1 ) THEN 824 CALL obs_int_z1d_spl( kpk, & 825 & zintt(iin,ijn,:,jobs),& 826 & zobs2k, zgdept(iin,ijn,:,jobs), & 827 & zmask(iin,ijn,:,jobs)) 828 829 ENDIF 830 831 CALL obs_level_search(kpk, & 832 & zgdept(iin,ijn,:,jobs),& 833 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 834 & iv_indic) 835 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 836 & prodatqc%var(1)%vdep(ista:iend), & 837 & zintt(iin,ijn,:,jobs), & 838 & zobs2k,interp_corner(iin,ijn,:), & 839 & zgdept(iin,ijn,:,jobs), & 840 & zmask(iin,ijn,:,jobs) ) 841 842 ENDDO 843 ENDDO 844 845 ENDIF 846 847 !------------------------------------------------------------- 848 ! Compute the horizontal interpolation for every profile level 849 !------------------------------------------------------------- 850 851 DO ikn=1,inum_obs 852 iend=ista+ikn-1 853 854 ! This code forces the horizontal weights to be 855 ! zero IF the observation is below the bottom of the 856 ! corners of the interpolation nodes, Or if it is in 857 ! the mask. This is important for observations are near 858 ! steep bathymetry 859 DO iin=1,2 860 DO ijn=1,2 861 862 depth_loop1: DO ik=kpk,2,-1 863 IF(zmask(iin,ijn,ik-1,jobs ) > 0.9 )THEN 864 865 l_zweig(iin,ijn,1) = & 866 & zweig(iin,ijn,1) * & 867 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,jobs) ) & 868 & - prodatqc%var(1)%vdep(iend)),0._wp) 869 870 EXIT depth_loop1 871 ENDIF 872 ENDDO depth_loop1 873 874 ENDDO 875 ENDDO 876 877 CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 878 & prodatqc%var(1)%vmod(iend:iend) ) 879 880 ENDDO 881 882 883 DEALLOCATE(interp_corner,iv_indic) 884 885 ENDIF 886 887 888 ! Salinity 889 890 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 891 892 zobsk(:) = obfillflt 893 894 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 895 896 IF ( idayend == 0 ) THEN 897 898 ! Daily averaged moored buoy (MRB) data 899 900 ! vertically interpolate all 4 corners 901 ista = prodatqc%npvsta(jobs,2) 902 iend = prodatqc%npvend(jobs,2) 903 inum_obs = iend - ista + 1 904 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 905 906 DO iin=1,2 907 DO ijn=1,2 908 909 910 911 IF ( k1dint == 1 ) THEN 912 CALL obs_int_z1d_spl( kpk, & 913 & zinms(iin,ijn,:,jobs), & 914 & zobs2k, zgdept(iin,ijn,:,jobs), & 915 & zmask(iin,ijn,:,jobs)) 916 ENDIF 917 918 CALL obs_level_search(kpk, & 919 & zgdept(iin,ijn,:,jobs), & 920 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 921 & iv_indic) 922 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 923 & prodatqc%var(2)%vdep(ista:iend), & 924 & zinms(iin,ijn,:,jobs), & 925 & zobs2k, interp_corner(iin,ijn,:), & 926 & zgdept(iin,ijn,:,jobs), & 927 & zmask(iin,ijn,:,jobs)) 928 929 ENDDO 930 ENDDO 931 932 933 ELSE 934 935 CALL ctl_stop( ' A nonzero' // & 936 & ' number of profile T BUOY data should' // & 937 & ' only occur at the end of a given day' ) 938 939 ENDIF 940 941 ELSE 942 943 ! Point data 944 945 ! vertically interpolate all 4 corners 946 ista = prodatqc%npvsta(jobs,2) 947 iend = prodatqc%npvend(jobs,2) 948 inum_obs = iend - ista + 1 949 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 950 951 DO iin=1,2 952 DO ijn=1,2 953 954 955 IF ( k1dint == 1 ) THEN 956 CALL obs_int_z1d_spl( kpk, & 957 & zints(iin,ijn,:,jobs),& 958 & zobs2k, zgdept(iin,ijn,:,jobs), & 959 & zmask(iin,ijn,:,jobs)) 960 961 ENDIF 962 963 CALL obs_level_search(kpk, & 964 & zgdept(iin,ijn,:,jobs),& 965 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 966 & iv_indic) 967 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 968 & prodatqc%var(2)%vdep(ista:iend), & 969 & zints(iin,ijn,:,jobs), & 970 & zobs2k,interp_corner(iin,ijn,:), & 971 & zgdept(iin,ijn,:,jobs), & 972 & zmask(iin,ijn,:,jobs) ) 973 974 ENDDO 975 ENDDO 976 977 ENDIF 978 979 !------------------------------------------------------------- 980 ! Compute the horizontal interpolation for every profile level 981 !------------------------------------------------------------- 982 983 DO ikn=1,inum_obs 984 iend=ista+ikn-1 985 986 ! This code forces the horizontal weights to be 987 ! zero IF the observation is below the bottom of the 988 ! corners of the interpolation nodes, Or if it is in 989 ! the mask. This is important for observations are near 990 ! steep bathymetry 991 DO iin=1,2 992 DO ijn=1,2 993 994 depth_loop2: DO ik=kpk,2,-1 995 IF(zmask(iin,ijn,ik-1,jobs ) > 0.9 )THEN 996 997 l_zweig(iin,ijn,1) = & 998 & zweig(iin,ijn,1) * & 999 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,jobs) ) & 1000 & - prodatqc%var(2)%vdep(iend)),0._wp) 1001 1002 EXIT depth_loop2 1003 ENDIF 1004 ENDDO depth_loop2 1005 1006 ENDDO 1007 ENDDO 1008 1009 CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 1010 & prodatqc%var(2)%vmod(iend:iend) ) 1011 1012 ENDDO 1013 1014 1015 DEALLOCATE(interp_corner,iv_indic) 1016 1017 ENDIF 1018 1019 ENDIF 1020 1021 END DO 1022 1023 ! Deallocate the data for interpolation 1024 DEALLOCATE( & 1025 & igrdi, & 1026 & igrdj, & 1027 & zglam, & 1028 & zgphi, & 1029 & zmask, & 1030 & zintt, & 1031 & zints & 1032 & ) 1033 ! At the end of the day also get interpolated means 1034 IF ( idayend == 0 ) THEN 1035 DEALLOCATE( & 1036 & zinmt, & 1037 & zinms & 1038 & ) 1039 ENDIF 1040 1041 prodatqc%nprofup = prodatqc%nprofup + ipro 1042 1043 END SUBROUTINE obs_pro_sco_opt 1044 451 1045 SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, & 452 1046 & psshn, psshmask, k2dint ) … … 1293 1887 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 1294 1888 & zglamv(:,:,iobs), zgphiv(:,:,iobs), & 1295 & zvmask(:,:,:,iobs), zweigv, zobsmask u)1889 & zvmask(:,:,:,iobs), zweigv, zobsmaskv ) 1296 1890 1297 1891 ENDIF -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
r4292 r5963 48 48 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 49 49 !!---------------------------------------------------------------------- 50 51 !! * Substitutions 52 # include "domzgr_substitute.h90" 50 53 51 54 CONTAINS … … 1709 1712 !! * Modules used 1710 1713 USE dom_oce, ONLY : & ! Geographical information 1711 & gdepw_1d 1712 1714 & gdepw_1d, & 1715 & gdepw_0, & 1716 #if defined key_vvl 1717 & gdepw_n, & 1718 & gdept_n, & 1719 #endif 1720 & ln_zco, & 1721 & ln_zps 1722 1713 1723 !! * Arguments 1714 1724 INTEGER, INTENT(IN) :: kprofno ! Number of profiles … … 1754 1764 & igrdj 1755 1765 LOGICAL :: lgridobs ! Is observation on a model grid point. 1766 LOGICAL :: ll_next_to_land ! Is a profile next to land 1756 1767 INTEGER :: iig, ijg ! i,j of observation on model grid point. 1757 1768 INTEGER :: jobs, jobsp, jk, ji, jj … … 1816 1827 END DO 1817 1828 1829 ! Check if next to land 1830 IF ( ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN 1831 ll_next_to_land=.TRUE. 1832 ELSE 1833 ll_next_to_land=.FALSE. 1834 ENDIF 1835 1818 1836 ! Reject observations 1819 1837 … … 1832 1850 ENDIF 1833 1851 1834 ! Flag if the observation falls with a model land cell 1835 IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1836 & == 0.0_wp ) THEN 1837 kobsqc(jobsp) = kobsqc(jobsp) + 12 1838 klanobs = klanobs + 1 1839 CYCLE 1852 ! To check if an observations falls within land there are two cases: 1853 ! 1: z-coordibnates, where the check uses the mask 1854 ! 2: terrain following (eg s-coordinates), 1855 ! where we use the depth of the bottom cell to mask observations 1856 1857 IF( ln_zps .OR. ln_zco ) THEN !(CASE 1) 1858 1859 ! Flag if the observation falls with a model land cell 1860 IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1861 & == 0.0_wp ) THEN 1862 kobsqc(jobsp) = kobsqc(jobsp) + 12 1863 klanobs = klanobs + 1 1864 CYCLE 1865 ENDIF 1866 1867 ! Flag if the observation is close to land 1868 IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 1869 & 0.0_wp) THEN 1870 knlaobs = knlaobs + 1 1871 IF (ld_nea) THEN 1872 kobsqc(jobsp) = kobsqc(jobsp) + 14 1873 ENDIF 1874 ENDIF 1875 1876 ELSE ! Case 2 1877 1878 ! Flag if the observation is deeper than the bathymetry 1879 ! Or if it is within the mask 1880 IF ( ALL( fsdepw(iig-1:iig+1,ijg-1:ijg+1,kpk) < pobsdep(jobsp) ) & 1881 & .OR. & 1882 & ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1883 & == 0.0_wp) ) THEN 1884 kobsqc(jobsp) = kobsqc(jobsp) + 12 1885 klanobs = klanobs + 1 1886 CYCLE 1887 ENDIF 1888 1889 ! Flag if the observation is close to land 1890 IF ( ll_next_to_land ) THEN 1891 knlaobs = knlaobs + 1 1892 IF (ld_nea) THEN 1893 kobsqc(jobsp) = kobsqc(jobsp) + 14 1894 ENDIF 1895 ENDIF 1896 1840 1897 ENDIF 1841 1898 … … 1850 1907 ENDIF 1851 1908 ENDIF 1852 1853 ! Flag if the observation falls is close to land 1854 IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 1855 & 0.0_wp) THEN 1856 IF (ld_nea) kobsqc(jobsp) = kobsqc(jobsp) + 14 1857 knlaobs = knlaobs + 1 1858 ENDIF 1909 1859 1910 1860 1911 ! Set observation depth equal to that of the first model depth
Note: See TracChangeset
for help on using the changeset viewer.
