Changeset 7412 for branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC
- Timestamp:
- 2016-12-01T11:30:29+01:00 (8 years ago)
- Location:
- branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC
- Files:
-
- 8 edited
-
fldread.F90 (modified) (17 diffs)
-
sbcice_lim.F90 (modified) (3 diffs)
-
sbcice_lim_2.F90 (modified) (2 diffs)
-
sbcmod.F90 (modified) (1 diff)
-
sbctide.F90 (modified) (3 diffs)
-
sbcwave.F90 (modified) (1 diff)
-
tideini.F90 (modified) (4 diffs)
-
updtide.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
r6140 r7412 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 *** … … 666 695 !! using a general mapping (for open boundaries) 667 696 !!---------------------------------------------------------------------- 668 #if defined key_bdy 669 USE bdy_oce, ONLY: dta_global, dta_global2! workspace to read in global data arrays670 #endif 697 698 USE bdy_oce, ONLY: ln_bdy, 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 699 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 ! … … 692 727 ilendta = iom_file(num)%dimsz(1,idvar) 693 728 694 #if defined key_bdy 695 ipj = iom_file(num)%dimsz(2,idvar) 696 IF( map%ll_unstruc) THEN ! unstructured open boundary data file 697 dta_read => dta_global 698 ELSE ! structured open boundary data file 699 dta_read => dta_global2 729 IF ( ln_bdy ) THEN 730 ipj = iom_file(num)%dimsz(2,idvar) 731 IF( map%ll_unstruc) THEN ! unstructured open boundary data file 732 dta_read => dta_global 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 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 749 ENDIF 700 750 ENDIF 701 #endif702 751 703 752 IF(lwp) WRITE(numout,*) 'Dim size for ', TRIM(clvar),' is ', ilendta … … 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 ( ln_bdy ) & 802 CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta) 803 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 SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta) 828 829 !!--------------------------------------------------------------------- 830 !! *** ROUTINE fld_bdy_interp *** 831 !! 832 !! ** Purpose : on the fly vertical interpolation to allow the use of 833 !! boundary data from non-native vertical grid 834 !!---------------------------------------------------------------------- 835 USE bdy_oce, ONLY: idx_bdy ! indexing for map <-> ij transformation 836 837 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read ! work space for global data 838 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read_z ! work space for global data 839 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read_dz ! work space for global data 840 REAL(wp) , INTENT(in) :: fv ! fillvalue and alternative -ABS(fv) 841 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional) 842 TYPE(MAP_POINTER) , INTENT(in ) :: map ! global-to-local mapping indices 843 LOGICAL , INTENT(in), OPTIONAL :: fvl ! grid type, set number and number of vertical levels in the bdy data 844 INTEGER , INTENT(in) :: igrd, ibdy, jpk_bdy ! number of levels in bdy data 845 INTEGER , INTENT(in) :: ilendta ! length of data in file 846 !! 847 INTEGER :: ipi ! length of boundary data on local process 848 INTEGER :: ipj ! length of dummy dimension ( = 1 ) 849 INTEGER :: ipk ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 850 INTEGER :: jpkm1_bdy ! number of levels in bdy data minus 1 851 INTEGER :: ib, ik, ikk ! loop counters 852 INTEGER :: ji, jj, zij, zjj ! temporary indices 853 REAL(wp) :: zl, zi, zh ! tmp variable for current depth and interpolation factor 854 REAL(wp) :: fv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(fv) 855 CHARACTER (LEN=10) :: ibstr 856 !!--------------------------------------------------------------------- 857 858 859 ipi = SIZE( dta, 1 ) 860 ipj = SIZE( dta_read, 2 ) 861 ipk = SIZE( dta, 3 ) 862 jpkm1_bdy = jpk_bdy-1 863 864 fv_alt = -ABS(fv) ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 865 DO ib = 1, ipi 866 zij = idx_bdy(ibdy)%nbi(ib,igrd) 867 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 868 IF(narea==2) WRITE(*,*) 'MAPI', ib, igrd, map%ptr(ib), narea-1, zij, zjj 869 ENDDO 870 ! 871 IF ( map%ll_unstruc ) THEN ! unstructured open boundary data file 872 712 873 DO ib = 1, ipi 713 DO ik = 1, ipk 714 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ik) 874 DO ik = 1, jpk_bdy 875 IF( ( dta_read(map%ptr(ib),1,ik) == fv ) ) THEN 876 dta_read_z(map%ptr(ib),1,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data 877 dta_read_dz(map%ptr(ib),1,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct 878 ENDIF 879 ENDDO 880 ENDDO 881 882 DO ib = 1, ipi 883 zij = idx_bdy(ibdy)%nbi(ib,igrd) 884 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 885 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 886 ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 887 SELECT CASE( igrd ) 888 CASE(1) 889 IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 890 WRITE(ibstr,"(I10.10)") map%ptr(ib) 891 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 892 ! 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 893 ENDIF 894 CASE(2) 895 IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 896 WRITE(ibstr,"(I10.10)") map%ptr(ib) 897 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 898 IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u_n(zij,zjj,:), mask=umask(zij,zjj,:)==1), sum(umask(zij,zjj,:)), & 899 & hu_n(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1 , & 900 & dta_read(map%ptr(ib),1,:) 901 ENDIF 902 CASE(3) 903 IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 904 WRITE(ibstr,"(I10.10)") map%ptr(ib) 905 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 906 ENDIF 907 END SELECT 908 DO ik = 1, ipk 909 SELECT CASE( igrd ) 910 CASE(1) 911 zl = gdept_n(zij,zjj,ik) ! if using in step could use fsdept instead of gdept_n? 912 CASE(2) 913 IF(ln_sco) THEN 914 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? 915 ELSE 916 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 917 ENDIF 918 CASE(3) 919 IF(ln_sco) THEN 920 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? 921 ELSE 922 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) ) 923 ENDIF 924 END SELECT 925 IF( zl < dta_read_z(map%ptr(ib),1,1) ) THEN ! above the first level of external data 926 dta(ib,1,ik) = dta_read(map%ptr(ib),1,1) 927 ELSEIF( zl > MAXVAL(dta_read_z(map%ptr(ib),1,:),1) ) THEN ! below the last level of external data 928 dta(ib,1,ik) = dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1)) 929 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 930 DO ikk = 1, jpkm1_bdy ! when gdept_n(ikk) < zl < gdept_n(ikk+1) 931 IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) & 932 & .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN 933 zi = ( zl - dta_read_z(map%ptr(ib),1,ikk) ) / & 934 & ( dta_read_z(map%ptr(ib),1,ikk+1) - dta_read_z(map%ptr(ib),1,ikk) ) 935 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ikk) + & 936 & ( dta_read(map%ptr(ib),1,ikk+1) - dta_read(map%ptr(ib),1,ikk) ) * zi 937 ENDIF 938 END DO 939 ENDIF 715 940 END DO 716 941 END DO 717 ELSE ! structured open boundary data file 942 943 IF(igrd == 2) THEN ! do we need to adjust the transport term? 944 DO ib = 1, ipi 945 zij = idx_bdy(ibdy)%nbi(ib,igrd) 946 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 947 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 948 ztrans = 0._wp 949 ztrans_new = 0._wp 950 DO ik = 1, jpk_bdy ! calculate transport on input grid 951 ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 952 ENDDO 953 DO ik = 1, ipk ! calculate transport on model grid 954 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik) 955 ENDDO 956 DO ik = 1, ipk ! make transport correction 957 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 958 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 959 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 960 IF( ABS(ztrans * r1_hu_n(zij,zjj)) > 0.01_wp ) & 961 & CALL ctl_warn('fld_bdy_interp: barotropic component of > 0.01 ms-1 found in baroclinic velocities at') 962 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) * umask(zij,zjj,ik) 963 ENDIF 964 ENDDO 965 ENDDO 966 ENDIF 967 968 IF(igrd == 3) THEN ! do we need to adjust the transport term? 969 DO ib = 1, ipi 970 zij = idx_bdy(ibdy)%nbi(ib,igrd) 971 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 972 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 973 ztrans = 0._wp 974 ztrans_new = 0._wp 975 DO ik = 1, jpk_bdy ! calculate transport on input grid 976 ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 977 ENDDO 978 DO ik = 1, ipk ! calculate transport on model grid 979 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik) 980 ENDDO 981 DO ik = 1, ipk ! make transport correction 982 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 983 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 984 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 985 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) * vmask(zij,zjj,ik) 986 ENDIF 987 ENDDO 988 ENDDO 989 ENDIF 990 991 ELSE ! structured open boundary file 992 718 993 DO ib = 1, ipi 719 994 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 720 995 ji=map%ptr(ib)-(jj-1)*ilendta 721 DO ik = 1, ipk 722 dta(ib,1,ik) = dta_read(ji,jj,ik) 996 DO ik = 1, jpk_bdy 997 IF( ( dta_read(ji,jj,ik) == fv ) ) THEN 998 dta_read_z(ji,jj,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data 999 dta_read_dz(ji,jj,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct 1000 ENDIF 1001 ENDDO 1002 ENDDO 1003 1004 1005 DO ib = 1, ipi 1006 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1007 ji=map%ptr(ib)-(jj-1)*ilendta 1008 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1009 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1010 zh = SUM(dta_read_dz(ji,jj,:) ) 1011 ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 1012 SELECT CASE( igrd ) 1013 CASE(1) 1014 IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 1015 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1016 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1017 ! 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 1018 ENDIF 1019 CASE(2) 1020 IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 1021 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1022 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1023 ENDIF 1024 CASE(3) 1025 IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 1026 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1027 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1028 ENDIF 1029 END SELECT 1030 DO ik = 1, ipk 1031 SELECT CASE( igrd ) ! coded for sco - need zco and zps option using min 1032 CASE(1) 1033 zl = gdept_n(zij,zjj,ik) ! if using in step could use fsdept instead of gdept_n? 1034 CASE(2) 1035 IF(ln_sco) THEN 1036 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? 1037 ELSE 1038 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 1039 ENDIF 1040 CASE(3) 1041 IF(ln_sco) THEN 1042 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? 1043 ELSE 1044 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) ) 1045 ENDIF 1046 END SELECT 1047 IF( zl < dta_read_z(ji,jj,1) ) THEN ! above the first level of external data 1048 dta(ib,1,ik) = dta_read(ji,jj,1) 1049 ELSEIF( zl > MAXVAL(dta_read_z(ji,jj,:),1) ) THEN ! below the last level of external data 1050 dta(ib,1,ik) = dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1)) 1051 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 1052 DO ikk = 1, jpkm1_bdy ! when gdept_n(ikk) < zl < gdept_n(ikk+1) 1053 IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) & 1054 & .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN 1055 zi = ( zl - dta_read_z(ji,jj,ikk) ) / & 1056 & ( dta_read_z(ji,jj,ikk+1) - dta_read_z(ji,jj,ikk) ) 1057 dta(ib,1,ik) = dta_read(ji,jj,ikk) + & 1058 & ( dta_read(ji,jj,ikk+1) - dta_read(ji,jj,ikk) ) * zi 1059 ENDIF 1060 END DO 1061 ENDIF 723 1062 END DO 724 1063 END DO 725 ENDIF 726 ! 727 END SUBROUTINE fld_map 1064 1065 IF(igrd == 2) THEN ! do we need to adjust the transport term? 1066 DO ib = 1, ipi 1067 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1068 ji=map%ptr(ib)-(jj-1)*ilendta 1069 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1070 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1071 zh = SUM(dta_read_dz(ji,jj,:) ) 1072 ztrans = 0._wp 1073 ztrans_new = 0._wp 1074 DO ik = 1, jpk_bdy ! calculate transport on input grid 1075 ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 1076 ENDDO 1077 DO ik = 1, ipk ! calculate transport on model grid 1078 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik) 1079 ENDDO 1080 DO ik = 1, ipk ! make transport correction 1081 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 1082 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 1083 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 1084 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 1085 ENDIF 1086 ENDDO 1087 ENDDO 1088 ENDIF 1089 1090 IF(igrd == 3) THEN ! do we need to adjust the transport term? 1091 DO ib = 1, ipi 1092 jj = 1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1093 ji = map%ptr(ib)-(jj-1)*ilendta 1094 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1095 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1096 zh = SUM(dta_read_dz(ji,jj,:) ) 1097 ztrans = 0._wp 1098 ztrans_new = 0._wp 1099 DO ik = 1, jpk_bdy ! calculate transport on input grid 1100 ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 1101 ENDDO 1102 DO ik = 1, ipk ! calculate transport on model grid 1103 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik) 1104 ENDDO 1105 DO ik = 1, ipk ! make transport correction 1106 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 1107 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 1108 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 1109 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 1110 ENDIF 1111 ENDDO 1112 ENDDO 1113 ENDIF 1114 1115 ENDIF ! endif unstructured or structured 1116 1117 END SUBROUTINE fld_bdy_interp 728 1118 729 1119 -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r6416 r7412 62 62 USE timing ! Timing 63 63 64 #if defined key_bdy 65 USE bdyice_lim ! unstructured open boundary data (bdy_ice_lim routine) 66 #endif 64 USE bdy_oce , ONLY: ln_bdy 65 USE bdyice_lim ! unstructured open boundary data (bdy_ice_lim routine) 67 66 68 67 IMPLICIT NONE … … 166 165 IF( nn_monocat /= 2 ) CALL lim_itd_me ! Mechanical redistribution ! (ridging/rafting) 167 166 ! 168 #if defined key_bdy 169 CALL bdy_ice_lim( kt ) ! bdy ice thermo 170 IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) 171 #endif 167 IF( ln_bdy ) CALL bdy_ice_lim( kt ) ! bdy ice thermo 168 IF( ln_bdy .AND. ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) 172 169 ! 173 170 CALL lim_update1( kt ) ! Corrections … … 380 377 r1_nlay_s = 1._wp / REAL( nlay_s, wp ) 381 378 ! 382 #if defined key_bdy 383 IF( lwp .AND. ln_limdiahsb ) CALL ctl_warn('online conservation check activated but it does not work with BDY') 384 #endif 379 IF( lwp .AND. ln_bdy .AND. ln_limdiahsb ) & 380 & CALL ctl_warn('online conservation check activated but it does not work with BDY') 385 381 ! 386 382 END SUBROUTINE ice_run -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90
r6140 r7412 54 54 # endif 55 55 56 #if defined key_bdy 56 USE bdy_oce , ONLY: ln_bdy 57 57 USE bdyice_lim ! unstructured open boundary data (bdy_ice_lim routine) 58 #endif59 58 60 59 IMPLICIT NONE … … 230 229 CALL lim_trp_2 ( kt ) ! Ice transport ( Advection/diffusion ) 231 230 IF( ln_limdmp ) CALL lim_dmp_2 ( kt ) ! Ice damping 232 #if defined key_bdy 233 CALL bdy_ice_lim( kt ) ! bdy ice thermo 234 #endif 231 IF( ln_bdy ) CALL bdy_ice_lim( kt ) ! bdy ice thermo 235 232 END IF 236 233 ! ! Ice surface fluxes in coupled mode -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
r7403 r7412 47 47 USE traqsr ! active tracers: light penetration 48 48 USE sbcwave ! Wave module 49 USE bdy_ par ! Require lk_bdy49 USE bdy_oce , ONLY: ln_bdy 50 50 ! 51 51 USE prtctl ! Print control (prt_ctl routine) -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90
r6140 r7412 22 22 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: pot_astro ! 23 23 24 #if defined key_tide25 24 !!---------------------------------------------------------------------- 26 !! 'key_tide' :tidal potential25 !! tidal potential 27 26 !!---------------------------------------------------------------------- 28 27 !! sbc_tide : … … 30 29 !!---------------------------------------------------------------------- 31 30 32 LOGICAL, PUBLIC, PARAMETER :: lk_tide = .TRUE.33 31 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: amp_pot, phi_pot 34 32 … … 125 123 END SUBROUTINE tide_init_potential 126 124 127 #else128 !!----------------------------------------------------------------------129 !! Default case : Empty module130 !!----------------------------------------------------------------------131 LOGICAL, PUBLIC, PARAMETER :: lk_tide = .FALSE.132 CONTAINS133 SUBROUTINE sbc_tide( kt ) ! Empty routine134 INTEGER , INTENT(in) :: kt ! ocean time-step135 WRITE(*,*) 'sbc_tide: You should not have seen this print! error?', kt136 END SUBROUTINE sbc_tide137 #endif138 139 125 !!====================================================================== 140 126 END MODULE sbctide -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
r7403 r7412 141 141 wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - e3t_n(:,:,jk) * ze3hdiv(:,:,jk) 142 142 END DO 143 #if defined key_bdy 144 IF( l k_bdy ) THEN143 ! 144 IF( ln_bdy ) THEN 145 145 DO jk = 1, jpkm1 146 146 wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 147 147 END DO 148 148 ENDIF 149 #endif 149 150 150 CALL wrk_dealloc( jpi,jpj,jpk, ze3hdiv ) 151 151 ! -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90
r6140 r7412 25 25 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: ftide !: 26 26 27 LOGICAL , PUBLIC :: ln_tide !: 27 28 LOGICAL , PUBLIC :: ln_tide_pot !: 28 29 LOGICAL , PUBLIC :: ln_tide_ramp !: … … 48 49 INTEGER :: ios ! Local integer output status for namelist read 49 50 ! 50 NAMELIST/nam_tide/ln_tide _pot, ln_tide_ramp, rdttideramp, clname51 NAMELIST/nam_tide/ln_tide, ln_tide_pot, ln_tide_ramp, rdttideramp, clname 51 52 !!---------------------------------------------------------------------- 52 !53 IF(lwp) THEN54 WRITE(numout,*)55 WRITE(numout,*) 'tide_init : Initialization of the tidal components'56 WRITE(numout,*) '~~~~~~~~~ '57 ENDIF58 !59 CALL tide_init_Wave60 53 ! 61 54 ! Read Namelist nam_tide … … 69 62 IF(lwm) WRITE ( numond, nam_tide ) 70 63 ! 64 IF (ln_tide) THEN 65 IF (lwp) THEN 66 WRITE(numout,*) 67 WRITE(numout,*) 'tide_init : Initialization of the tidal components' 68 WRITE(numout,*) '~~~~~~~~~ ' 69 WRITE(numout,*) ' Namelist nam_tide' 70 WRITE(numout,*) ' Use tidal components : ln_tide = ', ln_tide 71 WRITE(numout,*) ' Apply astronomical potential : ln_tide_pot = ', ln_tide_pot 72 WRITE(numout,*) ' nb_harmo = ', nb_harmo 73 WRITE(numout,*) ' ln_tide_ramp = ', ln_tide_ramp 74 WRITE(numout,*) ' rdttideramp = ', rdttideramp 75 ENDIF 76 ELSE 77 IF(lwp) WRITE(numout,*) 78 IF(lwp) WRITE(numout,*) 'tide_init : tidal components not used (ln_tide = F)' 79 IF(lwp) WRITE(numout,*) '~~~~~~~~~ ' 80 RETURN 81 ENDIF 82 ! 83 CALL tide_init_Wave 84 ! 71 85 nb_harmo=0 72 86 DO jk = 1, jpmax_harmo … … 79 93 IF( nb_harmo == 0 ) CALL ctl_stop( 'tide_init : No tidal components set in nam_tide' ) 80 94 ! 81 IF(lwp) THEN82 WRITE(numout,*) ' Namelist nam_tide'83 WRITE(numout,*) ' Apply astronomical potential : ln_tide_pot =', ln_tide_pot84 WRITE(numout,*) ' nb_harmo = ', nb_harmo85 WRITE(numout,*) ' ln_tide_ramp = ', ln_tide_ramp86 WRITE(numout,*) ' rdttideramp = ', rdttideramp87 ENDIF88 95 IF( ln_tide_ramp.AND.((nitend-nit000+1)*rdt/rday < rdttideramp) ) & 89 96 & CALL ctl_stop('rdttideramp must be lower than run duration') -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90
r5913 r7412 5 5 !!====================================================================== 6 6 !! History : 9.0 ! 07 (O. Le Galloudec) Original code 7 !!----------------------------------------------------------------------8 #if defined key_tide9 !!----------------------------------------------------------------------10 !! 'key_tide' : tidal potential11 7 !!---------------------------------------------------------------------- 12 8 !! upd_tide : update tidal potential … … 81 77 END SUBROUTINE upd_tide 82 78 83 #else84 !!----------------------------------------------------------------------85 !! Dummy module : NO TIDE86 !!----------------------------------------------------------------------87 CONTAINS88 SUBROUTINE upd_tide( kt, kit, time_offset ) ! Empty routine89 INTEGER, INTENT(in) :: kt ! integer arg, dummy routine90 INTEGER, INTENT(in), OPTIONAL :: kit ! optional arg, dummy routine91 INTEGER, INTENT(in), OPTIONAL :: time_offset ! optional arg, dummy routine92 WRITE(*,*) 'upd_tide: You should not have seen this print! error?', kt93 END SUBROUTINE upd_tide94 95 #endif96 97 79 !!====================================================================== 98 80
Note: See TracChangeset
for help on using the changeset viewer.
