Changeset 11223 for NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/fldread.F90
- Timestamp:
- 2019-07-05T20:53:14+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/fldread.F90
r10425 r11223 80 80 INTEGER :: nreclast ! last record to be read in the current file 81 81 CHARACTER(len = 256) :: lsmname ! current name of the NetCDF mask file acting as a key 82 ! ! variables related to BDY 82 83 INTEGER :: igrd ! grid type for bdy data 83 84 INTEGER :: ibdy ! bdy set id number 85 INTEGER, POINTER, DIMENSION(:) :: imap ! Array of integer pointers to 1D arrays 86 LOGICAL :: ltotvel ! total velocity or not (T/F) 84 87 END TYPE FLD 85 86 TYPE, PUBLIC :: MAP_POINTER !: Map from input data file to local domain87 INTEGER, POINTER, DIMENSION(:) :: ptr ! Array of integer pointers to 1D arrays88 LOGICAL :: ll_unstruc ! Unstructured (T) or structured (F) boundary data file89 END TYPE MAP_POINTER90 88 91 89 !$AGRIF_DO_NOT_TREAT … … 129 127 CONTAINS 130 128 131 SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset, jpk_bdy, fvl)129 SUBROUTINE fld_read( kt, kn_fsbc, sd, kit, kt_offset ) 132 130 !!--------------------------------------------------------------------- 133 131 !! *** ROUTINE fld_read *** … … 144 142 INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step) 145 143 TYPE(FLD), INTENT(inout), DIMENSION(:) :: sd ! input field related variables 146 TYPE(MAP_POINTER),INTENT(in), OPTIONAL, DIMENSION(:) :: map ! global-to-local mapping indices147 144 INTEGER , INTENT(in ), OPTIONAL :: kit ! subcycle timestep for timesplitting option 148 145 INTEGER , INTENT(in ), OPTIONAL :: kt_offset ! provide fields at time other than "now" … … 150 147 ! ! kt_offset = +1 => fields at "after" time level 151 148 ! ! etc. 152 INTEGER , INTENT(in ), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data153 LOGICAL , INTENT(in ), OPTIONAL :: fvl ! number of vertical levels in the BDY data154 149 !! 155 150 INTEGER :: itmp ! local variable … … 166 161 REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation 167 162 CHARACTER(LEN=1000) :: clfmt ! write format 168 TYPE(MAP_POINTER) :: imap ! global-to-local mapping indices169 163 !!--------------------------------------------------------------------- 170 164 ll_firstcall = kt == nit000 … … 175 169 ENDIF 176 170 IF( PRESENT(kt_offset) ) it_offset = kt_offset 177 178 imap%ptr => NULL()179 171 180 172 ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar … … 188 180 IF( ll_firstcall ) THEN ! initialization 189 181 DO jf = 1, imf 190 IF( PRESENT(map) ) imap = map(jf) 191 IF( PRESENT(jpk_bdy) ) THEN 192 CALL fld_init( kn_fsbc, sd(jf), imap, jpk_bdy, fvl ) ! read each before field (put them in after as they will be swapped) 193 ELSE 194 CALL fld_init( kn_fsbc, sd(jf), imap ) ! read each before field (put them in after as they will be swapped) 195 ENDIF 182 IF( TRIM(sd(jf)%clrootname) == 'NOT USED' ) CYCLE 183 CALL fld_init( kn_fsbc, sd(jf) ) ! read each before field (put them in after as they will be swapped) 196 184 END DO 197 185 IF( lwp ) CALL wgt_print() ! control print … … 202 190 ! 203 191 DO jf = 1, imf ! --- loop over field --- ! 204 192 193 IF( TRIM(sd(jf)%clrootname) == 'NOT USED' ) CYCLE 194 205 195 IF( isecsbc > sd(jf)%nrec_a(2) .OR. ll_firstcall ) THEN ! read/update the after data? 206 207 IF( PRESENT(map) ) imap = map(jf) ! temporary definition of map208 196 209 197 sd(jf)%nrec_b(:) = sd(jf)%nrec_a(:) ! swap before record informations … … 222 210 itmp = sd(jf)%nrec_a(1) ! temporary storage 223 211 sd(jf)%nrec_a(1) = sd(jf)%nreclast ! read the last record of the file currently opened 224 CALL fld_get( sd(jf) , imap )! read after data212 CALL fld_get( sd(jf) ) ! read after data 225 213 sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2) ! re-swap before record field 226 214 sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1) ! update before record informations … … 240 228 & .AND. sd(jf)%nrec_b(1) /= sd(jf)%nrec_a(1) - 1 ) THEN 241 229 sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) - 1 ! move back to before record 242 CALL fld_get( sd(jf) , imap )! read after data230 CALL fld_get( sd(jf) ) ! read after data 243 231 sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2) ! re-swap before record field 244 232 sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1) ! update before record informations … … 285 273 286 274 ! read after data 287 IF( PRESENT(jpk_bdy) ) THEN 288 CALL fld_get( sd(jf), imap, jpk_bdy, fvl) 289 ELSE 290 CALL fld_get( sd(jf), imap ) 291 ENDIF 275 CALL fld_get( sd(jf) ) 276 292 277 ENDIF ! read new data? 293 278 END DO ! --- end loop over field --- ! … … 296 281 297 282 DO jf = 1, imf ! --- loop over field --- ! 283 ! 284 IF( TRIM(sd(jf)%clrootname) == 'NOT USED' ) CYCLE 298 285 ! 299 286 IF( sd(jf)%ln_tint ) THEN ! temporal interpolation … … 327 314 328 315 329 SUBROUTINE fld_init( kn_fsbc, sdjf , map , jpk_bdy, fvl)316 SUBROUTINE fld_init( kn_fsbc, sdjf ) 330 317 !!--------------------------------------------------------------------- 331 318 !! *** ROUTINE fld_init *** … … 336 323 INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step) 337 324 TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables 338 TYPE(MAP_POINTER),INTENT(in) :: map ! global-to-local mapping indices339 INTEGER , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data340 LOGICAL , INTENT(in), OPTIONAL :: fvl ! number of vertical levels in the BDY data341 325 !! 342 326 LOGICAL :: llprevyr ! are we reading previous year file? … … 351 335 CHARACTER(LEN=1000) :: clfmt ! write format 352 336 !!--------------------------------------------------------------------- 337 ! 353 338 llprevyr = .FALSE. 354 339 llprevmth = .FALSE. … … 433 418 ! 434 419 ! read before data in after arrays(as we will swap it later) 435 IF( PRESENT(jpk_bdy) ) THEN 436 CALL fld_get( sdjf, map, jpk_bdy, fvl ) 437 ELSE 438 CALL fld_get( sdjf, map ) 439 ENDIF 420 CALL fld_get( sdjf ) 440 421 ! 441 422 clfmt = "(' fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')" … … 613 594 614 595 615 SUBROUTINE fld_get( sdjf , map, jpk_bdy, fvl)596 SUBROUTINE fld_get( sdjf ) 616 597 !!--------------------------------------------------------------------- 617 598 !! *** ROUTINE fld_get *** … … 620 601 !!---------------------------------------------------------------------- 621 602 TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables 622 TYPE(MAP_POINTER), INTENT(in ) :: map ! global-to-local mapping indices623 INTEGER , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the bdy data624 LOGICAL , INTENT(in), OPTIONAL :: fvl ! number of vertical levels in the bdy data625 603 ! 626 604 INTEGER :: ipk ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) … … 634 612 ipk = SIZE( sdjf%fnow, 3 ) 635 613 ! 636 IF( ASSOCIATED(map%ptr) ) THEN 637 IF( PRESENT(jpk_bdy) ) THEN 638 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), & 639 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl ) 640 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), & 641 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl ) 642 ENDIF 643 ELSE 644 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map ) 645 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1), map ) 646 ENDIF 647 ENDIF 614 IF( ASSOCIATED(sdjf%imap) ) THEN 615 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), & 616 & sdjf%nrec_a(1), sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel ) 617 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), & 618 & sdjf%nrec_a(1), sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel ) 619 ENDIF 648 620 ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 649 621 CALL wgt_list( sdjf, iw ) … … 700 672 END SUBROUTINE fld_get 701 673 702 SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy, fvl ) 674 675 SUBROUTINE fld_map( knum, cdvar, pdta, krec, kmap, kgrd, kbdy, ldtotvel ) 703 676 !!--------------------------------------------------------------------- 704 677 !! *** ROUTINE fld_map *** … … 707 680 !! using a general mapping (for open boundaries) 708 681 !!---------------------------------------------------------------------- 709 710 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 711 712 INTEGER , INTENT(in ) :: num ! stream number 713 CHARACTER(LEN=*) , INTENT(in ) :: clvar ! variable name 714 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional) 715 INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice) 716 TYPE(MAP_POINTER) , INTENT(in ) :: map ! global-to-local mapping indices 717 INTEGER , INTENT(in), OPTIONAL :: igrd, ibdy, jpk_bdy ! grid type, set number and number of vertical levels in the bdy data 718 LOGICAL , INTENT(in), OPTIONAL :: fvl ! grid type, set number and number of vertical levels in the bdy data 719 INTEGER :: jpkm1_bdy! number of vertical levels in the bdy data minus 1 720 !! 721 INTEGER :: ipi ! length of boundary data on local process 722 INTEGER :: ipj ! length of dummy dimension ( = 1 ) 723 INTEGER :: ipk ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 724 INTEGER :: ilendta ! length of data in file 725 INTEGER :: idvar ! variable ID 726 INTEGER :: ib, ik, ji, jj ! loop counters 727 INTEGER :: ierr 728 REAL(wp) :: fv ! fillvalue 729 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read ! work space for global data 730 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read_z ! work space for global data 731 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read_dz ! work space for global data 732 !!--------------------------------------------------------------------- 733 ! 734 ipi = SIZE( dta, 1 ) 735 ipj = 1 736 ipk = SIZE( dta, 3 ) 737 ! 738 idvar = iom_varid( num, clvar ) 739 ilendta = iom_file(num)%dimsz(1,idvar) 740 741 IF ( ln_bdy ) THEN 742 ipj = iom_file(num)%dimsz(2,idvar) 743 IF( map%ll_unstruc) THEN ! unstructured open boundary data file 744 dta_read => dta_global 745 IF( PRESENT(jpk_bdy) ) THEN 746 IF( jpk_bdy>0 ) THEN 747 dta_read_z => dta_global_z 748 dta_read_dz => dta_global_dz 749 jpkm1_bdy = jpk_bdy-1 750 ENDIF 751 ENDIF 752 ELSE ! structured open boundary file 753 dta_read => dta_global2 754 IF( PRESENT(jpk_bdy) ) THEN 755 IF( jpk_bdy>0 ) THEN 756 dta_read_z => dta_global2_z 757 dta_read_dz => dta_global2_dz 758 jpkm1_bdy = jpk_bdy-1 759 ENDIF 760 ENDIF 761 ENDIF 762 ENDIF 763 764 IF(lwp) WRITE(numout,*) 'Dim size for ', TRIM(clvar),' is ', ilendta 765 IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk 766 ! 767 SELECT CASE( ipk ) 768 CASE(1) ; 769 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1 ), nrec ) 770 IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 771 DO ib = 1, ipi 772 DO ik = 1, ipk 773 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ik) 774 END DO 775 END DO 776 ELSE ! we assume that this is a structured open boundary file 777 DO ib = 1, ipi 778 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 779 ji=map%ptr(ib)-(jj-1)*ilendta 780 DO ik = 1, ipk 781 dta(ib,1,ik) = dta_read(ji,jj,ik) 782 END DO 783 END DO 784 ENDIF 682 INTEGER , INTENT(in ) :: knum ! stream number 683 CHARACTER(LEN=*) , INTENT(in ) :: cdvar ! variable name 684 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pdta ! bdy output field on model grid 685 INTEGER , INTENT(in ) :: krec ! record number to read (ie time slice) 686 INTEGER , DIMENSION(:) , INTENT(in ) :: kmap ! global-to-local bdy mapping indices 687 ! optional variables used for vertical interpolation: 688 INTEGER, OPTIONAL , INTENT(in ) :: kgrd ! grid type (t, u, v) 689 INTEGER, OPTIONAL , INTENT(in ) :: kbdy ! bdy number 690 LOGICAL, OPTIONAL , INTENT(in ) :: ldtotvel ! true if toal ( = barotrop + barocline) velocity 691 !! 692 INTEGER :: ipi ! length of boundary data on local process 693 INTEGER :: ipj ! length of dummy dimension ( = 1 ) 694 INTEGER :: ipk ! number of vertical levels of pdta ( 2D: ipk=1 ; 3D: ipk=jpk ) 695 INTEGER :: ipkb ! number of vertical levels in boundary data file 696 INTEGER :: idvar ! variable ID 697 INTEGER :: indims ! number of dimensions of the variable 698 INTEGER, DIMENSION(4) :: idimsz ! size of variable dimensions 699 REAL(wp) :: zfv ! fillvalue 700 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zz_read ! work space for global boundary data 701 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zdta_read ! work space local data requiring vertical interpolation 702 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zdta_read_z ! work space local data requiring vertical interpolation 703 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zdta_read_dz ! work space local data requiring vertical interpolation 704 CHARACTER(LEN=1),DIMENSION(3) :: clgrid 705 LOGICAL :: lluld ! is the variable using the unlimited dimension 706 !!--------------------------------------------------------------------- 707 ! 708 clgrid = (/'t','u','v'/) 709 ! 710 ipi = SIZE( pdta, 1 ) 711 ipj = SIZE( pdta, 2 ) ! must be equal to 1 712 ipk = SIZE( pdta, 3 ) 713 ! 714 idvar = iom_varid( knum, cdvar, kndims = indims, kdimsz = idimsz, lduld = lluld ) 715 IF( indims == 4 .OR. ( indims == 3 .AND. .NOT. lluld ) ) THEN ; ipkb = idimsz(3) ! xy(zl)t or xy(zl) 716 ELSE ; ipkb = 1 ! xy or xyt 717 ENDIF 718 ! 719 ALLOCATE( zz_read( idimsz(1), idimsz(2), ipkb ) ) ! ++++++++ !!! this can be very big... 720 ! 721 IF( ipk == 1 ) THEN 722 723 IF( ipkb /= 1 ) CALL ctl_stop( 'fld_map : we must have ipkb = 1 to read surface data' ) 724 CALL iom_get ( knum, jpdom_unknown, cdvar, zz_read(:,:,1), krec ) ! call iom_get with a 2D file 725 CALL fld_map_core( zz_read, kmap, pdta ) 785 726 786 727 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 787 728 ! Do we include something here to adjust barotropic velocities ! 788 729 ! in case of a depth difference between bdy files and ! 789 ! bathymetry in the case ln_ full_vel = .false. and jpk_bdy>0?!730 ! bathymetry in the case ln_totvel = .false. and ipkb>0? ! 790 731 ! [as the enveloping and parital cells could change H] ! 791 732 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 792 733 793 CASE DEFAULT ; 794 795 IF( PRESENT(jpk_bdy) .AND. jpk_bdy>0 ) THEN ! boundary data not on model grid: vertical interpolation 796 CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar ) 797 dta_read(:,:,:) = -ABS(fv) 798 dta_read_z(:,:,:) = 0._wp 799 dta_read_dz(:,:,:) = 0._wp 800 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:jpk_bdy), nrec ) 801 SELECT CASE( igrd ) 802 CASE(1) 803 CALL iom_get ( num, jpdom_unknown, 'gdept', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 804 CALL iom_get ( num, jpdom_unknown, 'e3t', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 805 CASE(2) 806 CALL iom_get ( num, jpdom_unknown, 'gdepu', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 807 CALL iom_get ( num, jpdom_unknown, 'e3u', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 808 CASE(3) 809 CALL iom_get ( num, jpdom_unknown, 'gdepv', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 810 CALL iom_get ( num, jpdom_unknown, 'e3v', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 811 END SELECT 812 813 IF ( ln_bdy ) & 814 CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta) 815 816 ELSE ! boundary data assumed to be on model grid 817 818 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec ) 819 IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 820 DO ib = 1, ipi 821 DO ik = 1, ipk 822 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ik) 823 END DO 734 ELSE 735 ! 736 CALL iom_get ( knum, jpdom_unknown, cdvar, zz_read(:,:,:), krec ) ! call iom_get with a 3D file 737 ! 738 IF( ipkb /= ipk ) THEN ! boundary data not on model vertical grid : vertical interpolation 739 ! 740 IF( ipk == jpk .AND. iom_varid(knum,'gdep'//clgrid(kgrd)) /= -1 .AND. iom_varid(knum,'e3'//clgrid(kgrd)) /= -1 ) THEN 741 742 ALLOCATE( zdta_read(ipi,ipj,ipkb), zdta_read_z(ipi,ipj,ipkb), zdta_read_dz(ipi,ipj,ipkb) ) 743 744 CALL fld_map_core( zz_read, kmap, zdta_read ) 745 CALL iom_get ( knum, jpdom_unknown, 'gdep'//clgrid(kgrd), zz_read ) ! read only once? Potential temporal evolution? 746 CALL fld_map_core( zz_read, kmap, zdta_read_z ) 747 CALL iom_get ( knum, jpdom_unknown, 'e3'//clgrid(kgrd), zz_read ) ! read only once? Potential temporal evolution? 748 CALL fld_map_core( zz_read, kmap, zdta_read_dz ) 749 750 CALL iom_getatt(knum, '_FillValue', zfv, cdvar=cdvar ) 751 CALL fld_bdy_interp(zdta_read, zdta_read_z, zdta_read_dz, pdta, kgrd, kbdy, zfv, ldtotvel) 752 DEALLOCATE( zdta_read, zdta_read_z, zdta_read_dz ) 753 754 ELSE 755 IF( ipk /= jpk ) CALL ctl_stop( 'fld_map : this should be an impossible case...' ) 756 WRITE(ctmp1,*) 'fld_map : vertical interpolation for bdy variable '//TRIM(cdvar)//' requires ' 757 IF( iom_varid(knum, 'gdep'//clgrid(kgrd)) == -1 ) CALL ctl_stop( ctmp1//'gdep'//clgrid(kgrd)//' variable' ) 758 IF( iom_varid(knum, 'e3'//clgrid(kgrd)) == -1 ) CALL ctl_stop( ctmp1// 'e3'//clgrid(kgrd)//' variable' ) 759 760 ENDIF 761 ! 762 ELSE ! bdy data assumed to be the same levels as bdy variables 763 ! 764 CALL fld_map_core( zz_read, kmap, pdta ) 765 ! 766 ENDIF ! ipkb /= ipk 767 ENDIF ! ipk == 1 768 769 DEALLOCATE( zz_read ) 770 771 END SUBROUTINE fld_map 772 773 774 SUBROUTINE fld_map_core( pdta_read, kmap, pdta_bdy ) 775 !!--------------------------------------------------------------------- 776 !! *** ROUTINE fld_map_core *** 777 !! 778 !! ** Purpose : inner core of fld_map 779 !!---------------------------------------------------------------------- 780 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdta_read ! global boundary data 781 INTEGER, DIMENSION(: ), INTENT(in ) :: kmap ! global-to-local bdy mapping indices 782 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pdta_bdy ! bdy output field on model grid 783 !! 784 INTEGER, DIMENSION(3) :: idim_read, idim_bdy ! arrays dimensions 785 INTEGER :: ji, jj, jk, jb ! loop counters 786 INTEGER :: im1 787 !!--------------------------------------------------------------------- 788 ! 789 idim_read = SHAPE( pdta_read ) 790 idim_bdy = SHAPE( pdta_bdy ) 791 ! 792 ! in all cases: idim_bdy(2) == 1 .AND. idim_read(1) * idim_read(2) == idim_bdy(1) 793 ! structured BDY with rimwidth > 1 : idim_read(2) == rimwidth /= 1 794 ! structured BDY with rimwidth == 1 or unstructured BDY: idim_read(2) == 1 795 ! 796 IF( idim_read(2) > 1 ) THEN ! structured BDY with rimwidth > 1 797 DO jk = 1, idim_bdy(3) 798 DO jb = 1, idim_bdy(1) 799 im1 = kmap(jb) - 1 800 jj = im1 / idim_read(1) + 1 801 ji = MOD( im1, idim_read(1) ) + 1 802 pdta_bdy(jb,1,jk) = pdta_read(ji,jj,jk) 824 803 END DO 825 ELSE ! we assume that this is a structured open boundary file 826 DO ib = 1, ipi 827 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 828 ji=map%ptr(ib)-(jj-1)*ilendta 829 DO ik = 1, ipk 830 dta(ib,1,ik) = dta_read(ji,jj,ik) 831 END DO 804 END DO 805 ELSE 806 DO jk = 1, idim_bdy(3) 807 DO jb = 1, idim_bdy(1) ! horizontal remap of bdy data on the local bdy 808 pdta_bdy(jb,1,jk) = pdta_read(kmap(jb),1,jk) 832 809 END DO 833 ENDIF 834 ENDIF ! PRESENT(jpk_bdy) 835 END SELECT 836 837 END SUBROUTINE fld_map 810 END DO 811 ENDIF 812 813 END SUBROUTINE fld_map_core 838 814 839 SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta)840 815 816 SUBROUTINE fld_bdy_interp(pdta_read, pdta_read_z, pdta_read_dz, pdta, kgrd, kbdy, pfv, ldtotvel) 841 817 !!--------------------------------------------------------------------- 842 818 !! *** ROUTINE fld_bdy_interp *** … … 847 823 USE bdy_oce, ONLY: idx_bdy ! indexing for map <-> ij transformation 848 824 849 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read ! work space for global data 850 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read_z ! work space for global data 851 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read_dz ! work space for global data 852 REAL(wp) , INTENT(in) :: fv ! fillvalue and alternative -ABS(fv) 853 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional) 854 TYPE(MAP_POINTER) , INTENT(in ) :: map ! global-to-local mapping indices 855 LOGICAL , INTENT(in), OPTIONAL :: fvl ! grid type, set number and number of vertical levels in the bdy data 856 INTEGER , INTENT(in) :: igrd, ibdy, jpk_bdy ! number of levels in bdy data 857 INTEGER , INTENT(in) :: ilendta ! length of data in file 858 !! 859 INTEGER :: ipi ! length of boundary data on local process 860 INTEGER :: ipj ! length of dummy dimension ( = 1 ) 861 INTEGER :: ipk ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 862 INTEGER :: jpkm1_bdy ! number of levels in bdy data minus 1 863 INTEGER :: ib, ik, ikk ! loop counters 864 INTEGER :: ji, jj, zij, zjj ! temporary indices 865 REAL(wp) :: zl, zi, zh ! tmp variable for current depth and interpolation factor 866 REAL(wp) :: fv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(fv) 867 CHARACTER (LEN=10) :: ibstr 825 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdta_read ! data read in bdy file 826 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdta_read_z ! depth of the data read in bdy file 827 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdta_read_dz ! thickness of the levels in bdy file 828 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pdta ! output field on model grid (2 dimensional) 829 REAL(wp) , INTENT(in ) :: pfv ! fillvalue of the data read in bdy file 830 LOGICAL , INTENT(in ) :: ldtotvel ! true if toal ( = barotrop + barocline) velocity 831 INTEGER , INTENT(in ) :: kgrd ! grid type (t, u, v) 832 INTEGER , INTENT(in ) :: kbdy ! bdy number 833 !! 834 INTEGER :: ipi ! length of boundary data on local process 835 INTEGER :: ipkb ! number of vertical levels in boundary data file 836 INTEGER :: jb, ji, jj, jk, jkb ! loop counters 837 REAL(wp) :: zcoef 838 REAL(wp) :: zl, zi, zh ! tmp variable for current depth and interpolation factor 839 REAL(wp) :: zfv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(pfv) 840 REAL(wp), DIMENSION(jpk) :: zdepth, zdhalf ! level and half-level depth 868 841 !!--------------------------------------------------------------------- 869 842 870 871 ipi = SIZE( dta, 1 ) 872 ipj = SIZE( dta_read, 2 ) 873 ipk = SIZE( dta, 3 ) 874 jpkm1_bdy = jpk_bdy-1 843 ipi = SIZE( pdta, 1 ) 844 ipkb = SIZE( pdta_read, 3 ) 875 845 876 fv_alt = -ABS(fv) ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 877 DO ib = 1, ipi 878 zij = idx_bdy(ibdy)%nbi(ib,igrd) 879 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 880 IF(narea==2) WRITE(*,*) 'MAPI', ib, igrd, map%ptr(ib), narea-1, zij, zjj 881 ENDDO 882 ! 883 IF ( map%ll_unstruc ) THEN ! unstructured open boundary data file 884 885 DO ib = 1, ipi 886 DO ik = 1, jpk_bdy 887 IF( ( dta_read(map%ptr(ib),1,ik) == fv ) ) THEN 888 dta_read_z(map%ptr(ib),1,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data 889 dta_read_dz(map%ptr(ib),1,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct 846 zfv_alt = -ABS(pfv) ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 847 ! 848 WHERE( pdta_read == pfv ) 849 pdta_read_z = zfv_alt ! safety: put fillvalue into external depth field so consistent with data 850 pdta_read_dz = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct 851 ENDWHERE 852 853 DO jb = 1, ipi 854 ji = idx_bdy(kbdy)%nbi(jb,kgrd) 855 jj = idx_bdy(kbdy)%nbj(jb,kgrd) 856 zh = SUM(pdta_read_dz(jb,1,:) ) 857 ! 858 ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 859 SELECT CASE( kgrd ) 860 CASE(1) 861 IF( ABS( (zh - ht_n(ji,jj)) / ht_n(ji,jj)) * tmask(ji,jj,1) > 0.01_wp ) THEN 862 WRITE(ctmp1,"(I10.10)") jb 863 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%') 864 ! IF(lwp) WRITE(numout,*) 'DEPTHT', zh, sum(e3t_n(ji,jj,:), mask=tmask(ji,jj,:)==1), ht_n(ji,jj), jb, jb, ji, jj 865 ENDIF 866 CASE(2) 867 IF( ABS( (zh - hu_n(ji,jj)) * r1_hu_n(ji,jj)) * umask(ji,jj,1) > 0.01_wp ) THEN 868 WRITE(ctmp1,"(I10.10)") jb 869 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%') 870 ! IF(lwp) WRITE(numout,*) 'DEPTHU', zh, SUM(e3u_n(ji,jj,:), mask=umask(ji,jj,:)==1), SUM(umask(ji,jj,:)), & 871 ! & hu_n(ji,jj), jb, jb, ji, jj, narea-1, pdta_read(jb,1,:) 872 ENDIF 873 CASE(3) 874 IF( ABS( (zh - hv_n(ji,jj)) * r1_hv_n(ji,jj)) * vmask(ji,jj,1) > 0.01_wp ) THEN 875 WRITE(ctmp1,"(I10.10)") jb 876 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%') 877 ENDIF 878 END SELECT 879 ! 880 SELECT CASE( kgrd ) 881 CASE(1) 882 ! depth of T points: 883 zdepth(:) = gdept_n(ji,jj,:) 884 CASE(2) 885 ! depth of U points: we must not use gdept_n as we don't want to do a communication 886 ! --> copy what is done for gdept_n in domvvl... 887 zdhalf(1) = 0.0_wp 888 zdepth(1) = 0.5_wp * e3uw_n(ji,jj,1) 889 DO jk = 2, jpk ! vertical sum 890 ! zcoef = umask - wumask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 891 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 892 ! ! 0.5 where jk = mikt 893 !!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ?? 894 zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) ) 895 zdhalf(jk) = zdhalf(jk-1) + e3u_n(ji,jj,jk-1) 896 zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5 * e3uw_n(ji,jj,jk)) & 897 & + (1-zcoef) * ( zdepth(jk-1) + e3uw_n(ji,jj,jk)) 898 END DO 899 CASE(3) 900 ! depth of V points: we must not use gdept_n as we don't want to do a communication 901 ! --> copy what is done for gdept_n in domvvl... 902 zdhalf(1) = 0.0_wp 903 zdepth(1) = 0.5_wp * e3vw_n(ji,jj,1) 904 DO jk = 2, jpk ! vertical sum 905 ! zcoef = vmask - wvmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 906 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 907 ! ! 0.5 where jk = mikt 908 !!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ?? 909 zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) ) 910 zdhalf(jk) = zdhalf(jk-1) + e3v_n(ji,jj,jk-1) 911 zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5 * e3vw_n(ji,jj,jk)) & 912 & + (1-zcoef) * ( zdepth(jk-1) + e3vw_n(ji,jj,jk)) 913 END DO 914 END SELECT 915 ! 916 DO jk = 1, jpk 917 IF( zdepth(jk) < pdta_read_z(jb,1, 1) ) THEN ! above the first level of external data 918 pdta(jb,1,jk) = pdta_read(jb,1,1) 919 ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkb) ) THEN ! below the last level of external data 920 pdta(jb,1,jk) = pdta_read(jb,1,MAXLOC(pdta_read_z(jb,1,:),1)) 921 ELSE ! inbetween: vertical interpolation between jkb & jkb+1 922 DO jkb = 1, ipkb-1 ! when gdept_n(jkb) < zdepth(jk) < gdept_n(jkb+1) 923 IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) <= 0._wp ) & 924 & .AND. ( pdta_read_z(jb,1,jkb+1) /= zfv_alt) ) THEN ! linear interpolation between 2 levels 925 zi = ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) / ( pdta_read_z(jb,1,jkb+1) - pdta_read_z(jb,1,jkb) ) 926 pdta(jb,1,jk) = pdta_read(jb,1,jkb) + ( pdta_read (jb,1,jkb+1) - pdta_read (jb,1,jkb) ) * zi 927 ENDIF 928 END DO 929 ENDIF 930 END DO ! jpk 931 ! 932 END DO ! ipi 933 934 IF(kgrd == 2) THEN ! do we need to adjust the transport term? 935 DO jb = 1, ipi 936 ji = idx_bdy(kbdy)%nbi(jb,kgrd) 937 jj = idx_bdy(kbdy)%nbj(jb,kgrd) 938 zh = SUM(pdta_read_dz(jb,1,:) ) 939 ztrans = 0._wp 940 ztrans_new = 0._wp 941 DO jkb = 1, ipkb ! calculate transport on input grid 942 ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb) 943 ENDDO 944 DO jk = 1, jpk ! calculate transport on model grid 945 ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3u_n(ji,jj,jk ) * umask(ji,jj,jk) 946 ENDDO 947 DO jk = 1, jpk ! make transport correction 948 IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 949 pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu_n(ji,jj) ) * umask(ji,jj,jk) 950 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 951 pdta(jb,1,jk) = pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hu_n(ji,jj) * umask(ji,jj,jk) 890 952 ENDIF 891 953 ENDDO 892 ENDDO 893 894 DO ib = 1, ipi 895 zij = idx_bdy(ibdy)%nbi(ib,igrd) 896 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 897 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 898 ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 899 SELECT CASE( igrd ) 900 CASE(1) 901 IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 902 WRITE(ibstr,"(I10.10)") map%ptr(ib) 903 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 904 ! 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 905 ENDIF 906 CASE(2) 907 IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 908 WRITE(ibstr,"(I10.10)") map%ptr(ib) 909 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 910 IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u_n(zij,zjj,:), mask=umask(zij,zjj,:)==1), sum(umask(zij,zjj,:)), & 911 & hu_n(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1 , & 912 & dta_read(map%ptr(ib),1,:) 913 ENDIF 914 CASE(3) 915 IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 916 WRITE(ibstr,"(I10.10)") map%ptr(ib) 917 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 918 ENDIF 919 END SELECT 920 DO ik = 1, ipk 921 SELECT CASE( igrd ) 922 CASE(1) 923 zl = gdept_n(zij,zjj,ik) ! if using in step could use fsdept instead of gdept_n? 924 CASE(2) 925 IF(ln_sco) THEN 926 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? 927 ELSE 928 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 929 ENDIF 930 CASE(3) 931 IF(ln_sco) THEN 932 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? 933 ELSE 934 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) ) 935 ENDIF 936 END SELECT 937 IF( zl < dta_read_z(map%ptr(ib),1,1) ) THEN ! above the first level of external data 938 dta(ib,1,ik) = dta_read(map%ptr(ib),1,1) 939 ELSEIF( zl > MAXVAL(dta_read_z(map%ptr(ib),1,:),1) ) THEN ! below the last level of external data 940 dta(ib,1,ik) = dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1)) 941 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 942 DO ikk = 1, jpkm1_bdy ! when gdept_n(ikk) < zl < gdept_n(ikk+1) 943 IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) & 944 & .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN 945 zi = ( zl - dta_read_z(map%ptr(ib),1,ikk) ) / & 946 & ( dta_read_z(map%ptr(ib),1,ikk+1) - dta_read_z(map%ptr(ib),1,ikk) ) 947 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ikk) + & 948 & ( dta_read(map%ptr(ib),1,ikk+1) - dta_read(map%ptr(ib),1,ikk) ) * zi 949 ENDIF 950 END DO 951 ENDIF 952 END DO 953 END DO 954 955 IF(igrd == 2) THEN ! do we need to adjust the transport term? 956 DO ib = 1, ipi 957 zij = idx_bdy(ibdy)%nbi(ib,igrd) 958 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 959 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 960 ztrans = 0._wp 961 ztrans_new = 0._wp 962 DO ik = 1, jpk_bdy ! calculate transport on input grid 963 ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 964 ENDDO 965 DO ik = 1, ipk ! calculate transport on model grid 966 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik) 967 ENDDO 968 DO ik = 1, ipk ! make transport correction 969 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 970 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 971 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 972 IF( ABS(ztrans * r1_hu_n(zij,zjj)) > 0.01_wp ) & 973 & CALL ctl_warn('fld_bdy_interp: barotropic component of > 0.01 ms-1 found in baroclinic velocities at') 974 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) * umask(zij,zjj,ik) 975 ENDIF 976 ENDDO 954 ENDDO 955 ENDIF 956 957 IF(kgrd == 3) THEN ! do we need to adjust the transport term? 958 DO jb = 1, ipi 959 ji = idx_bdy(kbdy)%nbi(jb,kgrd) 960 jj = idx_bdy(kbdy)%nbj(jb,kgrd) 961 zh = SUM(pdta_read_dz(jb,1,:) ) 962 ztrans = 0._wp 963 ztrans_new = 0._wp 964 DO jkb = 1, ipkb ! calculate transport on input grid 965 ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb) 977 966 ENDDO 978 ENDIF 979 980 IF(igrd == 3) THEN ! do we need to adjust the transport term? 981 DO ib = 1, ipi 982 zij = idx_bdy(ibdy)%nbi(ib,igrd) 983 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 984 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 985 ztrans = 0._wp 986 ztrans_new = 0._wp 987 DO ik = 1, jpk_bdy ! calculate transport on input grid 988 ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 989 ENDDO 990 DO ik = 1, ipk ! calculate transport on model grid 991 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik) 992 ENDDO 993 DO ik = 1, ipk ! make transport correction 994 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 995 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 996 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 997 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) * vmask(zij,zjj,ik) 998 ENDIF 999 ENDDO 967 DO jk = 1, jpk ! calculate transport on model grid 968 ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3v_n(ji,jj,jk ) * vmask(ji,jj,jk) 1000 969 ENDDO 1001 ENDIF 1002 1003 ELSE ! structured open boundary file 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 DO ik = 1, jpk_bdy 1009 IF( ( dta_read(ji,jj,ik) == fv ) ) THEN 1010 dta_read_z(ji,jj,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data 1011 dta_read_dz(ji,jj,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct 970 DO jk = 1, jpk ! make transport correction 971 IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 972 pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv_n(ji,jj) ) * vmask(ji,jj,jk) 973 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 974 pdta(jb,1,jk) = pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hv_n(ji,jj) * vmask(ji,jj,jk) 1012 975 ENDIF 1013 976 ENDDO 1014 ENDDO 1015 1016 1017 DO ib = 1, ipi 1018 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1019 ji=map%ptr(ib)-(jj-1)*ilendta 1020 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1021 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1022 zh = SUM(dta_read_dz(ji,jj,:) ) 1023 ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 1024 SELECT CASE( igrd ) 1025 CASE(1) 1026 IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 1027 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1028 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1029 ! 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 1030 ENDIF 1031 CASE(2) 1032 IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 1033 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1034 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1035 ENDIF 1036 CASE(3) 1037 IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 1038 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1039 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1040 ENDIF 1041 END SELECT 1042 DO ik = 1, ipk 1043 SELECT CASE( igrd ) ! coded for sco - need zco and zps option using min 1044 CASE(1) 1045 zl = gdept_n(zij,zjj,ik) ! if using in step could use fsdept instead of gdept_n? 1046 CASE(2) 1047 IF(ln_sco) THEN 1048 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? 1049 ELSE 1050 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 1051 ENDIF 1052 CASE(3) 1053 IF(ln_sco) THEN 1054 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? 1055 ELSE 1056 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) ) 1057 ENDIF 1058 END SELECT 1059 IF( zl < dta_read_z(ji,jj,1) ) THEN ! above the first level of external data 1060 dta(ib,1,ik) = dta_read(ji,jj,1) 1061 ELSEIF( zl > MAXVAL(dta_read_z(ji,jj,:),1) ) THEN ! below the last level of external data 1062 dta(ib,1,ik) = dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1)) 1063 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 1064 DO ikk = 1, jpkm1_bdy ! when gdept_n(ikk) < zl < gdept_n(ikk+1) 1065 IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) & 1066 & .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN 1067 zi = ( zl - dta_read_z(ji,jj,ikk) ) / & 1068 & ( dta_read_z(ji,jj,ikk+1) - dta_read_z(ji,jj,ikk) ) 1069 dta(ib,1,ik) = dta_read(ji,jj,ikk) + & 1070 & ( dta_read(ji,jj,ikk+1) - dta_read(ji,jj,ikk) ) * zi 1071 ENDIF 1072 END DO 1073 ENDIF 1074 END DO 1075 END DO 1076 1077 IF(igrd == 2) THEN ! do we need to adjust the transport term? 1078 DO ib = 1, ipi 1079 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1080 ji=map%ptr(ib)-(jj-1)*ilendta 1081 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1082 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1083 zh = SUM(dta_read_dz(ji,jj,:) ) 1084 ztrans = 0._wp 1085 ztrans_new = 0._wp 1086 DO ik = 1, jpk_bdy ! calculate transport on input grid 1087 ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 1088 ENDDO 1089 DO ik = 1, ipk ! calculate transport on model grid 1090 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik) 1091 ENDDO 1092 DO ik = 1, ipk ! make transport correction 1093 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 1094 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 1095 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 1096 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 1097 ENDIF 1098 ENDDO 1099 ENDDO 1100 ENDIF 1101 1102 IF(igrd == 3) THEN ! do we need to adjust the transport term? 1103 DO ib = 1, ipi 1104 jj = 1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1105 ji = map%ptr(ib)-(jj-1)*ilendta 1106 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1107 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1108 zh = SUM(dta_read_dz(ji,jj,:) ) 1109 ztrans = 0._wp 1110 ztrans_new = 0._wp 1111 DO ik = 1, jpk_bdy ! calculate transport on input grid 1112 ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 1113 ENDDO 1114 DO ik = 1, ipk ! calculate transport on model grid 1115 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik) 1116 ENDDO 1117 DO ik = 1, ipk ! make transport correction 1118 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 1119 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 1120 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 1121 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 1122 ENDIF 1123 ENDDO 1124 ENDDO 1125 ENDIF 1126 1127 ENDIF ! endif unstructured or structured 1128 977 ENDDO 978 ENDIF 979 1129 980 END SUBROUTINE fld_bdy_interp 1130 981 … … 1151 1002 imf = SIZE( sd ) 1152 1003 DO ju = 1, imf 1004 IF( TRIM(sd(ju)%clrootname) == 'NOT USED' ) CYCLE 1153 1005 ill = LEN_TRIM( sd(ju)%vcomp ) 1154 1006 DO jn = 2-COUNT((/sd(ju)%ln_tint/)), 2 … … 1159 1011 iv = -1 1160 1012 DO jv = 1, imf 1013 IF( TRIM(sd(jv)%clrootname) == 'NOT USED' ) CYCLE 1161 1014 IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) ) iv = jv 1162 1015 END DO … … 1299 1152 ! 1300 1153 DO jf = 1, SIZE(sdf) 1301 sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname ) 1154 sdf(jf)%clrootname = sdf_n(jf)%clname 1155 IF( TRIM(sdf_n(jf)%clname) /= 'NOT USED' ) sdf(jf)%clrootname = TRIM( cdir )//sdf(jf)%clrootname 1302 1156 sdf(jf)%clname = "not yet defined" 1303 1157 sdf(jf)%nfreqh = sdf_n(jf)%nfreqh … … 1308 1162 sdf(jf)%num = -1 1309 1163 sdf(jf)%wgtname = " " 1310 IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) sdf(jf)%wgtname = TRIM( cdir )// TRIM( sdf_n(jf)%wname )1164 IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) sdf(jf)%wgtname = TRIM( cdir )//sdf_n(jf)%wname 1311 1165 sdf(jf)%lsmname = " " 1312 IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 ) sdf(jf)%lsmname = TRIM( cdir )// TRIM( sdf_n(jf)%lname )1166 IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 ) sdf(jf)%lsmname = TRIM( cdir )//sdf_n(jf)%lname 1313 1167 sdf(jf)%vcomp = sdf_n(jf)%vcomp 1314 1168 sdf(jf)%rotn(:) = .TRUE. ! pretend to be rotated -> won't try to rotate data before the first call to fld_get … … 1317 1171 IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim ) & 1318 1172 & CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.') 1319 sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn 1173 sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn 1174 sdf(jf)%igrd = 0 1175 sdf(jf)%ibdy = 0 1176 sdf(jf)%imap => NULL() 1177 sdf(jf)%ltotvel = .FALSE. 1320 1178 END DO 1321 1179 !
Note: See TracChangeset
for help on using the changeset viewer.
