Changeset 5705


Ignore:
Timestamp:
2015-08-25T14:17:45+02:00 (5 years ago)
Author:
jamesharle
Message:

Update to fld_bdy_interp - structuring, typos and minor bugfixes

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90

    r5700 r5705  
    106106   INTEGER                            ::   nxt_wgt = 1  ! point to next available space in ref_wgts array 
    107107   REAL(wp), PARAMETER                ::   undeff_lsm = -999.00_wp 
     108#  include "domzgr_substitute.h90" 
    108109 
    109110!$AGRIF_END_DO_NOT_TREAT 
     
    839840      USE bdy_oce, ONLY:  idx_bdy         ! indexing for map <-> ij transformation 
    840841 
    841       REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read    ! work space for global data 
    842       REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_z  ! work space for global data 
    843       REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_dz  ! work space for global data 
    844       REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta        ! output field on model grid (2 dimensional) 
    845       TYPE(MAP_POINTER)         , INTENT(in ) ::   map     ! global-to-local mapping indices 
    846       LOGICAL  , INTENT(in), OPTIONAL         ::   fvl  ! grid type, set number and number of vertical levels in the bdy data 
    847       INTEGER  , INTENT(in)                   ::   igrd, ibdy, jpk_bdy      ! number of levels in bdy data 
    848       INTEGER                                 ::   jpkm1_bdy    ! number of levels in bdy data minus 1 
    849       REAL(wp) , INTENT(in)                                ::   fv ! fillvalue and alternative -ABS(fv) 
    850       !! 
    851       INTEGER                                 ::   ipi        ! length of boundary data on local process 
    852       INTEGER                                 ::   ipj        ! length of dummy dimension ( = 1 ) 
    853       INTEGER                                 ::   ipk ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
    854       INTEGER, INTENT(in)                                 ::   ilendta    ! length of data in file 
    855       INTEGER                                 ::   ib, ik, ikk! loop counters 
    856       INTEGER                                 ::   ji, jj ! loop counters 
    857       REAL(wp)                                ::   zl, zi, zh, zz, zdz    ! tmp variable for current depth and interpolation factor 
     842      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read      ! work space for global data 
     843      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_z    ! work space for global data 
     844      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_dz   ! work space for global data 
     845      REAL(wp) , INTENT(in)                                ::   fv            ! fillvalue and alternative -ABS(fv) 
     846      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta                        ! output field on model grid (2 dimensional) 
     847      TYPE(MAP_POINTER)         , INTENT(in ) ::   map                        ! global-to-local mapping indices 
     848      LOGICAL  , INTENT(in), OPTIONAL         ::   fvl                        ! grid type, set number and number of vertical levels in the bdy data 
     849      INTEGER  , INTENT(in)                   ::   igrd, ibdy, jpk_bdy        ! number of levels in bdy data 
     850      INTEGER  , INTENT(in)                   ::   ilendta                    ! length of data in file 
     851      !! 
     852      INTEGER                                 ::   ipi                        ! length of boundary data on local process 
     853      INTEGER                                 ::   ipj                        ! length of dummy dimension ( = 1 ) 
     854      INTEGER                                 ::   ipk                        ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
     855      INTEGER                                 ::   jpkm1_bdy                  ! number of levels in bdy data minus 1 
     856      INTEGER                                 ::   ib, ik, ikk                ! loop counters 
     857      INTEGER                                 ::   ji, jj, zij, zjj           ! temporary indices 
     858      REAL(wp)                                ::   zl, zi, zh                 ! tmp variable for current depth and interpolation factor 
    858859      REAL(wp)                                ::   fv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(fv) 
     860      CHARACTER (LEN=10)                      ::   ibstr 
    859861      !!--------------------------------------------------------------------- 
    860862      
     
    868870       
    869871      ! 
    870       IF ( map%ll_unstruc ) THEN ! unstructured open boundary data file 
     872      IF ( map%ll_unstruc ) THEN                            ! unstructured open boundary data file 
    871873 
    872874         DO ib = 1, ipi 
    873875            DO ik = 1, jpk_bdy 
    874876               IF( ( dta_read(map%ptr(ib),1,ik) == fv ) ) THEN 
    875                   dta_read_z(map%ptr(ib),1,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data 
    876                   dta_read_dz(map%ptr(ib),1,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct 
     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 
    877879               ENDIF 
    878880            ENDDO 
     
    880882 
    881883         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,:) ) 
    882887            DO ik = 1, ipk                       
    883                zl =  gdept_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)   ! if using in step could use fsdept instead of gdept_0? 
     888               SELECT CASE( igrd )                        
     889                  CASE(1) 
     890                     zl =  gdept_0(zij,zjj,ik)                                          ! if using in step could use fsdept instead of gdept_0? 
     891                  CASE(2) 
     892                     IF(ln_sco) THEN 
     893                        zl =  ( gdept_0(zij,zjj,ik) + gdept_0(zij+1,zjj,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_0? 
     894                     ELSE 
     895                        zl =  MIN( gdept_0(zij,zjj,ik), gdept_0(zij+1,zjj,ik) )  
     896                     ENDIF 
     897                  CASE(3) 
     898                     IF(ln_sco) THEN 
     899                        zl =  ( gdept_0(zij,zjj,ik) + gdept_0(zij,zjj+1,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_0? 
     900                     ELSE 
     901                        zl =  MIN( gdept_0(zij,zjj,ik), gdept_0(zij,zjj+1,ik) ) 
     902                     ENDIF 
     903               END SELECT 
    884904               IF( zl < dta_read_z(map%ptr(ib),1,1) ) THEN                                         ! above the first level of external data 
    885905                  dta(ib,1,ik) =  dta_read(map%ptr(ib),1,1) 
    886906               ELSEIF( zl > MAXVAL(dta_read_z(map%ptr(ib),1,:),1) ) THEN                           ! below the last level of external data  
    887907                  dta(ib,1,ik) =  dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1)) 
    888                ELSE                                                                          ! inbetween : vertical interpolation between ikk & ikk+1 
    889                   DO ikk = 1, jpkm1_bdy                                                          ! when  gdept_0(ikk) < zl < gdept_0(ikk+1) 
    890                      IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp)   & 
     908               ELSE                                                                                ! inbetween : vertical interpolation between ikk & ikk+1 
     909                  DO ikk = 1, jpkm1_bdy                                                            ! when  gdept_0(ikk) < zl < gdept_0(ikk+1) 
     910                     IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) & 
    891911                    &    .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN 
    892                         zi = ( zl - dta_read_z(map%ptr(ib),1,ikk) ) / (dta_read_z(map%ptr(ib),1,ikk+1)-dta_read_z(map%ptr(ib),1,ikk)) 
     912                        zi           = ( zl - dta_read_z(map%ptr(ib),1,ikk) ) / & 
     913                       &               ( dta_read_z(map%ptr(ib),1,ikk+1) - dta_read_z(map%ptr(ib),1,ikk) ) 
    893914                        dta(ib,1,ik) = dta_read(map%ptr(ib),1,ikk) + & 
    894                       &                ( dta_read(map%ptr(ib),1,ikk+1) - dta_read(map%ptr(ib),1,ikk) ) * zi 
     915                       &               ( dta_read(map%ptr(ib),1,ikk+1) - dta_read(map%ptr(ib),1,ikk) ) * zi 
    895916                     ENDIF 
    896917                  END DO 
     
    899920         END DO 
    900921 
    901          IF(igrd == 2) THEN ! do we need to adjust the transport term? 
     922         IF(igrd == 2) THEN                                 ! do we need to adjust the transport term? 
    902923            DO ib = 1, ipi 
    903               zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 
     924              zij = idx_bdy(ibdy)%nbi(ib,igrd) 
     925              zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
     926              zh  = SUM(dta_read_dz(map%ptr(ib),1,:) ) 
    904927              ztrans = 0._wp 
    905928              ztrans_new = 0._wp 
    906               DO ik = 1, jpk_bdy 
    907                   ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 
     929              DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
     930                  ztrans     = ztrans    + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 
    908931              ENDDO 
    909               DO ik = 1, ipk 
    910                   zdz =  e3u_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)  
    911                   ztrans_new = ztrans_new + dta(ib,1,ik) * zdz 
     932              DO ik = 1, ipk                                ! calculate transport on model grid 
     933                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_0(zij,zjj,ik) * umask(zij,zjj,ik) 
    912934              ENDDO 
    913               DO ik = 1, ipk 
    914                  zdz =  e3u_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)    
    915                  zz  =  hur(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd))       ! Is array HUR the correct one to use???? 
     935              DO ik = 1, ipk                                ! make transport correction 
    916936                 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    917                     dta(ib,1,ik) = dta(ib,1,ik) + ( ztrans - ztrans_new ) * ( zdz * zz ) 
     937                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) / hu_0(zij,zjj) ) * umask(zij,zjj,ik) 
    918938                 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
    919                     dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * ( zdz * zz ) 
     939                    IF( ABS(ztrans / hu_0(zij,zjj)) > 0.01_wp ) & 
     940                   &   CALL ctl_warn('fld_bdy_interp: barotropic component of > 0.01 ms-1 found in baroclinic velocities at') 
     941                    dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) / hu_0(zij,zjj) * umask(zij,zjj,ik) 
    920942                 ENDIF 
    921943              ENDDO 
     
    923945         ENDIF 
    924946 
    925          IF(igrd == 3) THEN ! do we need to adjust the transport term? 
     947         IF(igrd == 3) THEN                                 ! do we need to adjust the transport term? 
    926948            DO ib = 1, ipi 
    927               zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 
     949              zij = idx_bdy(ibdy)%nbi(ib,igrd) 
     950              zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
     951              zh  = SUM(dta_read_dz(map%ptr(ib),1,:) ) 
    928952              ztrans = 0._wp 
    929953              ztrans_new = 0._wp 
    930               DO ik = 1, jpk_bdy 
    931                   ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 
     954              DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
     955                  ztrans     = ztrans    + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 
    932956              ENDDO 
    933               DO ik = 1, ipk 
    934                   zdz =  e3v_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)  
    935                   ztrans_new = ztrans_new + dta(ib,1,ik) * zdz 
     957              DO ik = 1, ipk                                ! calculate transport on model grid 
     958                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_0(zij,zjj,ik) * vmask(zij,zjj,ik) 
    936959              ENDDO 
    937               DO ik = 1, ipk 
    938                  zdz =  e3v_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)    
    939                  zz  =  hvr(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd))    
     960              DO ik = 1, ipk                                ! make transport correction 
    940961                 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    941                     dta(ib,1,ik) = dta(ib,1,ik) + ( ztrans - ztrans_new ) * ( zdz * zz ) 
     962                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) / hv_0(zij,zjj) ) * vmask(zij,zjj,ik) 
    942963                 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
    943                     dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * ( zdz * zz ) 
     964                    dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) / hv_0(zij,zjj) * vmask(zij,zjj,ik) 
    944965                 ENDIF 
    945966              ENDDO 
     
    947968         ENDIF 
    948969   
    949          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    950          ! At this point write out a single velocity profile/dz/H for model and input data             ! 
    951          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    952  
    953970      ELSE ! structured open boundary file 
    954971 
     
    958975            DO ik = 1, jpk_bdy                       
    959976               IF( ( dta_read(ji,jj,ik) == fv ) ) THEN 
    960                   dta_read_z(ji,jj,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data 
    961                   dta_read_dz(ji,jj,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct 
     977                  dta_read_z(ji,jj,ik)  = fv_alt ! safety: put fillvalue into external depth field so consistent with data 
     978                  dta_read_dz(ji,jj,ik) = 0._wp  ! safety: put 0._wp into external thickness factors to ensure transport is correct 
    962979               ENDIF 
    963980            ENDDO 
     
    968985            jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    969986            ji=map%ptr(ib)-(jj-1)*ilendta 
     987            zij = idx_bdy(ibdy)%nbi(ib,igrd) 
     988            zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
     989            zh  = SUM(dta_read_dz(ji,jj,:) ) 
    970990            DO ik = 1, ipk                       
    971                zl =  gdept_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)   ! if using in step could use fsdept instead of gdept_0? 
    972                IF( zl < dta_read_z(ji,jj,1) ) THEN                                         ! above the first level of external data 
     991               SELECT CASE( igrd )                          ! coded for sco - need zco and zps option using min 
     992                  CASE(1) 
     993                     zl =  gdept_0(zij,zjj,ik)              ! if using in step could use fsdept instead of gdept_0? 
     994                  CASE(2) 
     995                     IF(ln_sco) THEN 
     996                        zl =  ( gdept_0(zij,zjj,ik) + gdept_0(zij+1,zjj,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_0? 
     997                     ELSE 
     998                        zl =  MIN( gdept_0(zij,zjj,ik), gdept_0(zij+1,zjj,ik) ) 
     999                     ENDIF 
     1000                  CASE(3) 
     1001                     IF(ln_sco) THEN 
     1002                        zl =  ( gdept_0(zij,zjj,ik) + gdept_0(zij,zjj+1,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_0? 
     1003                     ELSE 
     1004                        zl =  MIN( gdept_0(zij,zjj,ik), gdept_0(zij,zjj+1,ik) ) 
     1005                     ENDIF 
     1006               END SELECT 
     1007               IF( zl < dta_read_z(ji,jj,1) ) THEN                                      ! above the first level of external data 
    9731008                  dta(ib,1,ik) =  dta_read(ji,jj,1) 
    974                ELSEIF( zl > MAXVAL(dta_read_z(ji,jj,:),1) ) THEN                           ! below the last level of external data  
     1009               ELSEIF( zl > MAXVAL(dta_read_z(ji,jj,:),1) ) THEN                        ! below the last level of external data  
    9751010                  dta(ib,1,ik) =  dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1)) 
    976                ELSE                                                                          ! inbetween : vertical interpolation between ikk & ikk+1 
    977                   DO ikk = 1, jpkm1_bdy                                                          ! when  gdept_0(ikk) < zl < gdept_0(ikk+1) 
    978                      IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp)   & 
     1011               ELSE                                                                     ! inbetween : vertical interpolation between ikk & ikk+1 
     1012                  DO ikk = 1, jpkm1_bdy                                                 ! when  gdept_0(ikk) < zl < gdept_0(ikk+1) 
     1013                     IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) & 
    9791014                    &    .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN 
    980                         zi = ( zl - dta_read_z(ji,jj,ikk) ) / (dta_read_z(ji,jj,ikk+1)-dta_read_z(ji,jj,ikk)) 
     1015                        zi           = ( zl - dta_read_z(ji,jj,ikk) ) / & 
     1016                       &               ( dta_read_z(ji,jj,ikk+1) - dta_read_z(ji,jj,ikk) ) 
    9811017                        dta(ib,1,ik) = dta_read(ji,jj,ikk) + & 
    982                       &                ( dta_read(ji,jj,ikk+1) - dta_read(ji,jj,ikk) ) * zi 
     1018                       &               ( dta_read(ji,jj,ikk+1) - dta_read(ji,jj,ikk) ) * zi 
    9831019                     ENDIF 
    9841020                  END DO 
     
    9871023         END DO 
    9881024 
    989          IF(igrd == 2) THEN ! do we need to adjust the transport term? 
     1025         IF(igrd == 2) THEN                                 ! do we need to adjust the transport term? 
    9901026            jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    9911027            ji=map%ptr(ib)-(jj-1)*ilendta 
     1028            zij = idx_bdy(ibdy)%nbi(ib,igrd) 
     1029            zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    9921030            DO ib = 1, ipi 
    9931031              zh = SUM(dta_read_dz(ji,jj,:) ) 
    9941032              ztrans = 0._wp 
    9951033              ztrans_new = 0._wp 
    996               DO ik = 1, jpk_bdy 
     1034              DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
    9971035                  ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 
    9981036              ENDDO 
    999               DO ik = 1, ipk 
    1000                   zdz =  e3u_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik) 
    1001                   ztrans_new = ztrans_new + dta(ib,1,ik) * zdz 
     1037              DO ik = 1, ipk                                ! calculate transport on model grid 
     1038                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_0(zij,zjj,ik) * umask(zij,zjj,ik) 
    10021039              ENDDO 
    1003               DO ik = 1, ipk 
    1004                  zdz =  e3u_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik) 
    1005                  zz  =  hur(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd)) 
     1040              DO ik = 1, ipk                                ! make transport correction 
    10061041                 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    1007                     dta(ib,1,ik) = dta(ib,1,ik) + ( ztrans - ztrans_new ) * ( zdz * zz ) 
     1042                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * hur(zij,zjj) ) * umask(zij,zjj,ik) 
    10081043                 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
    1009                     dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * ( zdz * zz ) 
     1044                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * hur(zij,zjj) ) * umask(zij,zjj,ik) 
    10101045                 ENDIF 
    10111046              ENDDO 
     
    10131048         ENDIF 
    10141049 
    1015          IF(igrd == 3) THEN ! do we need to adjust the transport term? 
    1016             jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    1017             ji=map%ptr(ib)-(jj-1)*ilendta 
     1050         IF(igrd == 3) THEN                                 ! do we need to adjust the transport term? 
    10181051           DO ib = 1, ipi 
    1019               zh = SUM(dta_read_dz(ji,jj,:) ) 
     1052              jj  = 1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
     1053              ji  = map%ptr(ib)-(jj-1)*ilendta 
     1054              zij = idx_bdy(ibdy)%nbi(ib,igrd) 
     1055              zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
     1056              zh  = SUM(dta_read_dz(ji,jj,:) ) 
    10201057              ztrans = 0._wp 
    10211058              ztrans_new = 0._wp 
    1022               DO ik = 1, jpk_bdy 
    1023                   ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 
     1059              DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
     1060                  ztrans     = ztrans    + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 
    10241061              ENDDO 
    1025               DO ik = 1, ipk 
    1026                   zdz =  e3v_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik) 
    1027                   ztrans_new = ztrans_new + dta(ib,1,ik) * zdz 
     1062              DO ik = 1, ipk                                ! calculate transport on model grid 
     1063                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_0(zij,zjj,ik) * vmask(zij,zjj,ik) 
    10281064              ENDDO 
    1029               DO ik = 1, ipk 
    1030                  zdz =  e3v_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik) 
    1031                  zz  =  hvr(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd)) 
     1065              DO ik = 1, ipk                                ! make transport correction 
    10321066                 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    1033                     dta(ib,1,ik) = dta(ib,1,ik) + ( ztrans - ztrans_new ) * ( zdz * zz ) 
     1067                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * hvr(zij,zjj) ) * vmask(zij,zjj,ik) 
    10341068                 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
    1035                     dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * ( zdz * zz ) 
     1069                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * hvr(zij,zjj) ) * vmask(zij,zjj,ik) 
    10361070                 ENDIF 
    10371071              ENDDO 
Note: See TracChangeset for help on using the changeset viewer.