- Timestamp:
- 2016-10-14T11:10:43+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
r6140 r7029 7 7 !! ! 05-2008 (S. Alderson) Modified for Interpolation in memory from input grid to model grid 8 8 !! ! 10-2013 (D. Delrosso, P. Oddo) suppression of land point prior to interpolation 9 !! ! 12-2015 (J. Harle) Adding BDY on-the-fly interpolation 9 10 !!---------------------------------------------------------------------- 10 11 … … 67 68 INTEGER :: nreclast ! last record to be read in the current file 68 69 CHARACTER(len = 256) :: lsmname ! current name of the NetCDF mask file acting as a key 70 INTEGER :: igrd ! grid type for bdy data 71 INTEGER :: ibdy ! bdy set id number 69 72 END TYPE FLD 70 73 … … 114 117 CONTAINS 115 118 116 SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset )119 SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset, jpk_bdy, fvl ) 117 120 !!--------------------------------------------------------------------- 118 121 !! *** ROUTINE fld_read *** … … 135 138 ! ! kt_offset = +1 => fields at "after" time level 136 139 ! ! etc. 140 INTEGER , INTENT(in ), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data 141 LOGICAL , INTENT(in ), OPTIONAL :: fvl ! number of vertical levels in the BDY data 142 !! 137 143 INTEGER :: itmp ! local variable 138 144 INTEGER :: imf ! size of the structure sd … … 171 177 DO jf = 1, imf 172 178 IF( PRESENT(map) ) imap = map(jf) 173 CALL fld_init( kn_fsbc, sd(jf), imap ) ! read each before field (put them in after as they will be swapped) 179 IF( PRESENT(jpk_bdy) ) THEN 180 CALL fld_init( kn_fsbc, sd(jf), imap, jpk_bdy, fvl ) ! read each before field (put them in after as they will be swapped) 181 ELSE 182 CALL fld_init( kn_fsbc, sd(jf), imap ) ! read each before field (put them in after as they will be swapped) 183 ENDIF 174 184 END DO 175 185 IF( lwp ) CALL wgt_print() ! control print … … 263 273 264 274 ! read after data 265 CALL fld_get( sd(jf), imap ) 266 275 IF( PRESENT(jpk_bdy) ) THEN 276 CALL fld_get( sd(jf), imap, jpk_bdy, fvl) 277 ELSE 278 CALL fld_get( sd(jf), imap ) 279 ENDIF 267 280 ENDIF ! read new data? 268 281 END DO ! --- end loop over field --- ! … … 302 315 303 316 304 SUBROUTINE fld_init( kn_fsbc, sdjf, map )317 SUBROUTINE fld_init( kn_fsbc, sdjf, map , jpk_bdy, fvl) 305 318 !!--------------------------------------------------------------------- 306 319 !! *** ROUTINE fld_init *** … … 309 322 !! - if time interpolation, read before data 310 323 !!---------------------------------------------------------------------- 311 INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step) 312 TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables 313 TYPE(MAP_POINTER),INTENT(in) :: map ! global-to-local mapping indices 324 INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step) 325 TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables 326 TYPE(MAP_POINTER),INTENT(in) :: map ! global-to-local mapping indices 327 INTEGER , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data 328 LOGICAL , INTENT(in), OPTIONAL :: fvl ! number of vertical levels in the BDY data 314 329 !! 315 330 LOGICAL :: llprevyr ! are we reading previous year file? … … 405 420 ENDIF 406 421 ! 407 CALL fld_get( sdjf, map ) ! read before data in after arrays(as we will swap it later) 422 ! read before data in after arrays(as we will swap it later) 423 IF( PRESENT(jpk_bdy) ) THEN 424 CALL fld_get( sdjf, map, jpk_bdy, fvl ) 425 ELSE 426 CALL fld_get( sdjf, map ) 427 ENDIF 408 428 ! 409 429 clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')" … … 581 601 582 602 583 SUBROUTINE fld_get( sdjf, map )603 SUBROUTINE fld_get( sdjf, map, jpk_bdy, fvl ) 584 604 !!--------------------------------------------------------------------- 585 605 !! *** ROUTINE fld_get *** … … 589 609 TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables 590 610 TYPE(MAP_POINTER), INTENT(in ) :: map ! global-to-local mapping indices 611 INTEGER , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the bdy data 612 LOGICAL , INTENT(in), OPTIONAL :: fvl ! number of vertical levels in the bdy data 591 613 ! 592 614 INTEGER :: ipk ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) … … 601 623 ! 602 624 IF( ASSOCIATED(map%ptr) ) THEN 603 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map ) 604 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1), map ) 605 ENDIF 625 IF( PRESENT(jpk_bdy) ) THEN 626 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), & 627 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl ) 628 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), & 629 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl ) 630 ENDIF 631 ELSE 632 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map ) 633 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1), map ) 634 ENDIF 635 ENDIF 606 636 ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 607 637 CALL wgt_list( sdjf, iw ) … … 658 688 END SUBROUTINE fld_get 659 689 660 661 SUBROUTINE fld_map( num, clvar, dta, nrec, map ) 690 SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy, fvl ) 662 691 !!--------------------------------------------------------------------- 663 692 !! *** ROUTINE fld_map *** … … 667 696 !!---------------------------------------------------------------------- 668 697 #if defined key_bdy 669 USE bdy_oce, ONLY: dta_global, dta_global2! workspace to read in global data arrays698 USE bdy_oce, ONLY: idx_bdy, dta_global, dta_global_z, dta_global_dz, dta_global2, dta_global2_z, dta_global2_dz ! workspace to read in global data arrays 670 699 #endif 671 700 INTEGER , INTENT(in ) :: num ! stream number 672 701 CHARACTER(LEN=*) , INTENT(in ) :: clvar ! variable name 673 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional)702 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional) 674 703 INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice) 675 704 TYPE(MAP_POINTER) , INTENT(in ) :: map ! global-to-local mapping indices 705 INTEGER , INTENT(in), OPTIONAL :: igrd, ibdy, jpk_bdy ! grid type, set number and number of vertical levels in the bdy data 706 LOGICAL , INTENT(in), OPTIONAL :: fvl ! grid type, set number and number of vertical levels in the bdy data 707 INTEGER :: jpkm1_bdy! number of vertical levels in the bdy data minus 1 676 708 !! 677 709 INTEGER :: ipi ! length of boundary data on local process … … 682 714 INTEGER :: ib, ik, ji, jj ! loop counters 683 715 INTEGER :: ierr 684 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read ! work space for global data 716 REAL(wp) :: fv ! fillvalue 717 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read ! work space for global data 718 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read_z ! work space for global data 719 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read_dz ! work space for global data 685 720 !!--------------------------------------------------------------------- 686 721 ! … … 696 731 IF( map%ll_unstruc) THEN ! unstructured open boundary data file 697 732 dta_read => dta_global 698 ELSE ! structured open boundary data file 733 IF( PRESENT(jpk_bdy) ) THEN 734 IF( jpk_bdy>0 ) THEN 735 dta_read_z => dta_global_z 736 dta_read_dz => dta_global_dz 737 jpkm1_bdy = jpk_bdy-1 738 ENDIF 739 ENDIF 740 ELSE ! structured open boundary file 699 741 dta_read => dta_global2 742 IF( PRESENT(jpk_bdy) ) THEN 743 IF( jpk_bdy>0 ) THEN 744 dta_read_z => dta_global2_z 745 dta_read_dz => dta_global2_dz 746 jpkm1_bdy = jpk_bdy-1 747 ENDIF 748 ENDIF 700 749 ENDIF 701 750 #endif … … 705 754 ! 706 755 SELECT CASE( ipk ) 707 CASE(1) ; CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1 ), nrec ) 708 CASE DEFAULT ; CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec ) 756 CASE(1) ; 757 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1 ), nrec ) 758 IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 759 DO ib = 1, ipi 760 DO ik = 1, ipk 761 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ik) 762 END DO 763 END DO 764 ELSE ! we assume that this is a structured open boundary file 765 DO ib = 1, ipi 766 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 767 ji=map%ptr(ib)-(jj-1)*ilendta 768 DO ik = 1, ipk 769 dta(ib,1,ik) = dta_read(ji,jj,ik) 770 END DO 771 END DO 772 ENDIF 773 774 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 775 ! Do we include something here to adjust barotropic velocities ! 776 ! in case of a depth difference between bdy files and ! 777 ! bathymetry in the case ln_full_vel = .false. and jpk_bdy>0? ! 778 ! [as the enveloping and parital cells could change H] ! 779 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 780 781 CASE DEFAULT ; 782 783 IF( PRESENT(jpk_bdy) .AND. jpk_bdy>0 ) THEN ! boundary data not on model grid: vertical interpolation 784 CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar ) 785 dta_read(:,:,:) = -ABS(fv) 786 dta_read_z(:,:,:) = 0._wp 787 dta_read_dz(:,:,:) = 0._wp 788 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:jpk_bdy), nrec ) 789 SELECT CASE( igrd ) 790 CASE(1) 791 CALL iom_get ( num, jpdom_unknown, 'gdept', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 792 CALL iom_get ( num, jpdom_unknown, 'e3t', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 793 CASE(2) 794 CALL iom_get ( num, jpdom_unknown, 'gdepu', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 795 CALL iom_get ( num, jpdom_unknown, 'e3u', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 796 CASE(3) 797 CALL iom_get ( num, jpdom_unknown, 'gdepv', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 798 CALL iom_get ( num, jpdom_unknown, 'e3v', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 799 END SELECT 800 801 #if defined key_bdy 802 CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta) 803 #endif 804 ELSE ! boundary data assumed to be on model grid 805 806 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec ) 807 IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 808 DO ib = 1, ipi 809 DO ik = 1, ipk 810 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ik) 811 END DO 812 END DO 813 ELSE ! we assume that this is a structured open boundary file 814 DO ib = 1, ipi 815 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 816 ji=map%ptr(ib)-(jj-1)*ilendta 817 DO ik = 1, ipk 818 dta(ib,1,ik) = dta_read(ji,jj,ik) 819 END DO 820 END DO 821 ENDIF 822 ENDIF ! PRESENT(jpk_bdy) 709 823 END SELECT 710 ! 711 IF( map%ll_unstruc ) THEN ! unstructured open boundary data file 824 825 END SUBROUTINE fld_map 826 827 #if defined key_bdy 828 SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta) 829 830 !!--------------------------------------------------------------------- 831 !! *** ROUTINE fld_bdy_interp *** 832 !! 833 !! ** Purpose : on the fly vertical interpolation to allow the use of 834 !! boundary data from non-native vertical grid 835 !!---------------------------------------------------------------------- 836 USE bdy_oce, ONLY: idx_bdy ! indexing for map <-> ij transformation 837 838 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read ! work space for global data 839 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read_z ! work space for global data 840 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read_dz ! work space for global data 841 REAL(wp) , INTENT(in) :: fv ! fillvalue and alternative -ABS(fv) 842 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional) 843 TYPE(MAP_POINTER) , INTENT(in ) :: map ! global-to-local mapping indices 844 LOGICAL , INTENT(in), OPTIONAL :: fvl ! grid type, set number and number of vertical levels in the bdy data 845 INTEGER , INTENT(in) :: igrd, ibdy, jpk_bdy ! number of levels in bdy data 846 INTEGER , INTENT(in) :: ilendta ! length of data in file 847 !! 848 INTEGER :: ipi ! length of boundary data on local process 849 INTEGER :: ipj ! length of dummy dimension ( = 1 ) 850 INTEGER :: ipk ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 851 INTEGER :: jpkm1_bdy ! number of levels in bdy data minus 1 852 INTEGER :: ib, ik, ikk ! loop counters 853 INTEGER :: ji, jj, zij, zjj ! temporary indices 854 REAL(wp) :: zl, zi, zh ! tmp variable for current depth and interpolation factor 855 REAL(wp) :: fv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(fv) 856 CHARACTER (LEN=10) :: ibstr 857 !!--------------------------------------------------------------------- 858 859 860 ipi = SIZE( dta, 1 ) 861 ipj = SIZE( dta_read, 2 ) 862 ipk = SIZE( dta, 3 ) 863 jpkm1_bdy = jpk_bdy-1 864 865 fv_alt = -ABS(fv) ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 866 DO ib = 1, ipi 867 zij = idx_bdy(ibdy)%nbi(ib,igrd) 868 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 869 IF(narea==2) WRITE(*,*) 'MAPI', ib, igrd, map%ptr(ib), narea-1, zij, zjj 870 ENDDO 871 ! 872 IF ( map%ll_unstruc ) THEN ! unstructured open boundary data file 873 712 874 DO ib = 1, ipi 713 DO ik = 1, ipk 714 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ik) 875 DO ik = 1, jpk_bdy 876 IF( ( dta_read(map%ptr(ib),1,ik) == fv ) ) THEN 877 dta_read_z(map%ptr(ib),1,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data 878 dta_read_dz(map%ptr(ib),1,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct 879 ENDIF 880 ENDDO 881 ENDDO 882 883 DO ib = 1, ipi 884 zij = idx_bdy(ibdy)%nbi(ib,igrd) 885 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 886 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 887 ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 888 SELECT CASE( igrd ) 889 CASE(1) 890 IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 891 WRITE(ibstr,"(I10.10)") map%ptr(ib) 892 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 893 ! IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1), ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj 894 ENDIF 895 CASE(2) 896 IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 897 WRITE(ibstr,"(I10.10)") map%ptr(ib) 898 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 899 IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u_n(zij,zjj,:), mask=umask(zij,zjj,:)==1), sum(umask(zij,zjj,:)), & 900 & hu_n(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1 , & 901 & dta_read(map%ptr(ib),1,:) 902 ENDIF 903 CASE(3) 904 IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 905 WRITE(ibstr,"(I10.10)") map%ptr(ib) 906 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 907 ENDIF 908 END SELECT 909 DO ik = 1, ipk 910 SELECT CASE( igrd ) 911 CASE(1) 912 zl = gdept_n(zij,zjj,ik) ! if using in step could use fsdept instead of gdept_n? 913 CASE(2) 914 IF(ln_sco) THEN 915 zl = ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 916 ELSE 917 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 918 ENDIF 919 CASE(3) 920 IF(ln_sco) THEN 921 zl = ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 922 ELSE 923 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) ) 924 ENDIF 925 END SELECT 926 IF( zl < dta_read_z(map%ptr(ib),1,1) ) THEN ! above the first level of external data 927 dta(ib,1,ik) = dta_read(map%ptr(ib),1,1) 928 ELSEIF( zl > MAXVAL(dta_read_z(map%ptr(ib),1,:),1) ) THEN ! below the last level of external data 929 dta(ib,1,ik) = dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1)) 930 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 931 DO ikk = 1, jpkm1_bdy ! when gdept_n(ikk) < zl < gdept_n(ikk+1) 932 IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) & 933 & .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN 934 zi = ( zl - dta_read_z(map%ptr(ib),1,ikk) ) / & 935 & ( dta_read_z(map%ptr(ib),1,ikk+1) - dta_read_z(map%ptr(ib),1,ikk) ) 936 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ikk) + & 937 & ( dta_read(map%ptr(ib),1,ikk+1) - dta_read(map%ptr(ib),1,ikk) ) * zi 938 ENDIF 939 END DO 940 ENDIF 715 941 END DO 716 942 END DO 717 ELSE ! structured open boundary data file 943 944 IF(igrd == 2) THEN ! do we need to adjust the transport term? 945 DO ib = 1, ipi 946 zij = idx_bdy(ibdy)%nbi(ib,igrd) 947 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 948 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 949 ztrans = 0._wp 950 ztrans_new = 0._wp 951 DO ik = 1, jpk_bdy ! calculate transport on input grid 952 ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 953 ENDDO 954 DO ik = 1, ipk ! calculate transport on model grid 955 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik) 956 ENDDO 957 DO ik = 1, ipk ! make transport correction 958 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 959 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 960 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 961 IF( ABS(ztrans * r1_hu_n(zij,zjj)) > 0.01_wp ) & 962 & CALL ctl_warn('fld_bdy_interp: barotropic component of > 0.01 ms-1 found in baroclinic velocities at') 963 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) * umask(zij,zjj,ik) 964 ENDIF 965 ENDDO 966 ENDDO 967 ENDIF 968 969 IF(igrd == 3) THEN ! do we need to adjust the transport term? 970 DO ib = 1, ipi 971 zij = idx_bdy(ibdy)%nbi(ib,igrd) 972 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 973 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 974 ztrans = 0._wp 975 ztrans_new = 0._wp 976 DO ik = 1, jpk_bdy ! calculate transport on input grid 977 ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 978 ENDDO 979 DO ik = 1, ipk ! calculate transport on model grid 980 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik) 981 ENDDO 982 DO ik = 1, ipk ! make transport correction 983 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 984 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 985 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 986 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) * vmask(zij,zjj,ik) 987 ENDIF 988 ENDDO 989 ENDDO 990 ENDIF 991 992 ELSE ! structured open boundary file 993 718 994 DO ib = 1, ipi 719 995 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 720 996 ji=map%ptr(ib)-(jj-1)*ilendta 721 DO ik = 1, ipk 722 dta(ib,1,ik) = dta_read(ji,jj,ik) 997 DO ik = 1, jpk_bdy 998 IF( ( dta_read(ji,jj,ik) == fv ) ) THEN 999 dta_read_z(ji,jj,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data 1000 dta_read_dz(ji,jj,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct 1001 ENDIF 1002 ENDDO 1003 ENDDO 1004 1005 1006 DO ib = 1, ipi 1007 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1008 ji=map%ptr(ib)-(jj-1)*ilendta 1009 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1010 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1011 zh = SUM(dta_read_dz(ji,jj,:) ) 1012 ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 1013 SELECT CASE( igrd ) 1014 CASE(1) 1015 IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 1016 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1017 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1018 ! IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1), ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj 1019 ENDIF 1020 CASE(2) 1021 IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 1022 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1023 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1024 ENDIF 1025 CASE(3) 1026 IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 1027 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1028 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1029 ENDIF 1030 END SELECT 1031 DO ik = 1, ipk 1032 SELECT CASE( igrd ) ! coded for sco - need zco and zps option using min 1033 CASE(1) 1034 zl = gdept_n(zij,zjj,ik) ! if using in step could use fsdept instead of gdept_n? 1035 CASE(2) 1036 IF(ln_sco) THEN 1037 zl = ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 1038 ELSE 1039 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 1040 ENDIF 1041 CASE(3) 1042 IF(ln_sco) THEN 1043 zl = ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 1044 ELSE 1045 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) ) 1046 ENDIF 1047 END SELECT 1048 IF( zl < dta_read_z(ji,jj,1) ) THEN ! above the first level of external data 1049 dta(ib,1,ik) = dta_read(ji,jj,1) 1050 ELSEIF( zl > MAXVAL(dta_read_z(ji,jj,:),1) ) THEN ! below the last level of external data 1051 dta(ib,1,ik) = dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1)) 1052 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 1053 DO ikk = 1, jpkm1_bdy ! when gdept_n(ikk) < zl < gdept_n(ikk+1) 1054 IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) & 1055 & .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN 1056 zi = ( zl - dta_read_z(ji,jj,ikk) ) / & 1057 & ( dta_read_z(ji,jj,ikk+1) - dta_read_z(ji,jj,ikk) ) 1058 dta(ib,1,ik) = dta_read(ji,jj,ikk) + & 1059 & ( dta_read(ji,jj,ikk+1) - dta_read(ji,jj,ikk) ) * zi 1060 ENDIF 1061 END DO 1062 ENDIF 723 1063 END DO 724 1064 END DO 725 ENDIF 726 ! 727 END SUBROUTINE fld_map 728 1065 1066 IF(igrd == 2) THEN ! do we need to adjust the transport term? 1067 DO ib = 1, ipi 1068 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1069 ji=map%ptr(ib)-(jj-1)*ilendta 1070 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1071 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1072 zh = SUM(dta_read_dz(ji,jj,:) ) 1073 ztrans = 0._wp 1074 ztrans_new = 0._wp 1075 DO ik = 1, jpk_bdy ! calculate transport on input grid 1076 ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 1077 ENDDO 1078 DO ik = 1, ipk ! calculate transport on model grid 1079 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik) 1080 ENDDO 1081 DO ik = 1, ipk ! make transport correction 1082 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 1083 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 1084 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 1085 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 1086 ENDIF 1087 ENDDO 1088 ENDDO 1089 ENDIF 1090 1091 IF(igrd == 3) THEN ! do we need to adjust the transport term? 1092 DO ib = 1, ipi 1093 jj = 1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1094 ji = map%ptr(ib)-(jj-1)*ilendta 1095 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1096 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1097 zh = SUM(dta_read_dz(ji,jj,:) ) 1098 ztrans = 0._wp 1099 ztrans_new = 0._wp 1100 DO ik = 1, jpk_bdy ! calculate transport on input grid 1101 ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 1102 ENDDO 1103 DO ik = 1, ipk ! calculate transport on model grid 1104 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik) 1105 ENDDO 1106 DO ik = 1, ipk ! make transport correction 1107 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 1108 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 1109 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 1110 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 1111 ENDIF 1112 ENDDO 1113 ENDDO 1114 ENDIF 1115 1116 ENDIF ! endif unstructured or structured 1117 1118 END SUBROUTINE fld_bdy_interp 1119 1120 ! SUBROUTINE fld_bdy_conserve(dta_read, dta_read_z, map, jpk_bdy, igrd, ibdy, fv, dta) 1121 1122 ! END SUBROUTINE fld_bdy_conserve 1123 1124 #endif 729 1125 730 1126 SUBROUTINE fld_rot( kt, sd )
Note: See TracChangeset
for help on using the changeset viewer.
