New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 7029 for branches/NERC – NEMO

Changeset 7029 for branches/NERC


Ignore:
Timestamp:
2016-10-14T11:10:43+02:00 (7 years ago)
Author:
jamesharle
Message:

Adding ORCHESTRA configuration
Merging with branches/2016/dev_r5549_BDY_ZEROGRAD
Merging with branches/2016/dev_r5840_BDY_MSK
Merging with branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP

Location:
branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM
Files:
7 added
14 edited

Legend:

Unmodified
Added
Removed
  • branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM/CONFIG/SHARED/namelist_ref

    r6497 r7029  
    651651    ln_vol        = .false.               !  total volume correction (see nn_volctl parameter) 
    652652    nn_volctl     = 1                     !  = 0, the total water flux across open boundaries is zero 
     653    nb_jpk_bdy    = -1                    ! number of levels in the bdy data (set < 0 if consistent with planned run) 
    653654/ 
    654655!----------------------------------------------------------------------- 
  • branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM/CONFIG/cfg.txt

    r6403 r7029  
    1111GYRE OPA_SRC 
    1212ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
     13ORCHESTRA OPA_SRC LIM_SRC_3 TOP_SRC 
  • branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r6140 r7029  
    4949      LOGICAL                           ::  ll_tem 
    5050      LOGICAL                           ::  ll_sal 
     51      LOGICAL                           ::  ll_fvl 
    5152      REAL(wp), POINTER, DIMENSION(:)   ::  ssh 
    5253      REAL(wp), POINTER, DIMENSION(:)   ::  u2d 
     
    9192   ! 
    9293   INTEGER                    ::   nb_bdy                   !: number of open boundary sets 
     94   INTEGER                    ::   nb_jpk_bdy               !: number of levels in the bdy data (set < 0 if consistent with planned run) 
    9395   INTEGER, DIMENSION(jp_bdy) ::   nn_rimwidth              !: boundary rim width for Flow Relaxation Scheme 
    9496   INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P  
     
    134136                                                                          !: =1 => some data to be read in from data files 
    135137   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays (unstr.  bdy) 
     138   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global_z      !: workspace for reading in global depth arrays (unstr.  bdy) 
     139   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global_dz     !: workspace for reading in global depth arrays (unstr.  bdy) 
    136140   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2       !: workspace for reading in global data arrays (struct. bdy) 
     141   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2_z     !: workspace for reading in global depth arrays (struct. bdy) 
     142   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2_dz    !: workspace for reading in global depth arrays (struct. bdy) 
    137143!$AGRIF_DO_NOT_TREAT 
    138144   TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process) 
  • branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r6140 r7029  
    267267 
    268268                        jend = jstart + dta%nread(2) - 1 
    269                         CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  & 
    270                                      & kit=jit, kt_offset=time_offset ) 
     269                        IF( ln_full_vel_array(ib_bdy) ) THEN 
     270                           CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  & 
     271                                     & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(ib_bdy)  ) 
     272                        ELSE 
     273                           CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  & 
     274                                     & kit=jit, kt_offset=time_offset  ) 
     275                        ENDIF 
    271276 
    272277                        ! If full velocities in boundary data then extract barotropic velocities from 3D fields 
     
    333338                     jend = jstart + dta%nread(1) - 1 
    334339                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), & 
    335                                   & map=nbmap_ptr(jstart:jend), kt_offset=time_offset ) 
     340                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(ib_bdy) ) 
    336341                  ENDIF 
    337342                  ! If full velocities in boundary data then split into barotropic and baroclinic data 
     
    459464      NAMELIST/nambdy_dta/ bn_a_i, bn_ht_i, bn_ht_s 
    460465#endif 
    461       NAMELIST/nambdy_dta/ ln_full_vel 
     466      NAMELIST/nambdy_dta/ ln_full_vel, nb_jpk_bdy 
    462467      !!--------------------------------------------------------------------------- 
    463468      ! 
  • branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90

    r6140 r7029  
    5757         CASE('orlanski' )   ;   CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. ) 
    5858         CASE('orlanski_npo');   CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. ) 
     59         CASE('zerograd')    ;   CALL bdy_dyn3d_zgrad( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
     60         CASE('neumann')     ;   CALL bdy_dyn3d_nmn( idx_bdy(ib_bdy), ib_bdy ) 
    5961         CASE DEFAULT        ;   CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' ) 
    6062         END SELECT 
     
    110112   END SUBROUTINE bdy_dyn3d_spe 
    111113 
     114   SUBROUTINE bdy_dyn3d_zgrad( idx, dta, kt , ib_bdy ) 
     115      !!---------------------------------------------------------------------- 
     116      !!                  ***  SUBROUTINE bdy_dyn3d_zgrad  *** 
     117      !! 
     118      !! ** Purpose : - Enforce a zero gradient of normal velocity 
     119      !! 
     120      !!---------------------------------------------------------------------- 
     121      INTEGER                     ::   kt 
     122      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
     123      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
     124      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index 
     125      !! 
     126      INTEGER  ::   jb, jk         ! dummy loop indices 
     127      INTEGER  ::   ii, ij, igrd   ! local integers 
     128      REAL(wp) ::   zwgt           ! boundary weight 
     129      INTEGER  ::   fu, fv 
     130      !!---------------------------------------------------------------------- 
     131      ! 
     132      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zgrad') 
     133      ! 
     134      igrd = 2                      ! Copying tangential velocity into bdy points 
     135      DO jb = 1, idx%nblenrim(igrd) 
     136         DO jk = 1, jpkm1 
     137            ii   = idx%nbi(jb,igrd) 
     138            ij   = idx%nbj(jb,igrd) 
     139            fu   = ABS( ABS (NINT( idx%flagu(jb,igrd) ) ) - 1 ) 
     140            ua(ii,ij,jk) = ua(ii,ij,jk) * REAL( 1 - fu) + ( ua(ii,ij+fu,jk) * umask(ii,ij+fu,jk) & 
     141                        &+ ua(ii,ij-fu,jk) * umask(ii,ij-fu,jk) ) * umask(ii,ij,jk) * REAL( fu ) 
     142         END DO 
     143      END DO 
     144      ! 
     145      igrd = 3                      ! Copying tangential velocity into bdy points 
     146      DO jb = 1, idx%nblenrim(igrd) 
     147         DO jk = 1, jpkm1 
     148            ii   = idx%nbi(jb,igrd) 
     149            ij   = idx%nbj(jb,igrd) 
     150            fv   = ABS( ABS (NINT( idx%flagv(jb,igrd) ) ) - 1 ) 
     151            va(ii,ij,jk) = va(ii,ij,jk) * REAL( 1 - fv ) + ( va(ii+fv,ij,jk) * vmask(ii+fv,ij,jk) & 
     152                        &+ va(ii-fv,ij,jk) * vmask(ii-fv,ij,jk) ) * vmask(ii,ij,jk) * REAL( fv ) 
     153         END DO 
     154      END DO 
     155      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ! Boundary points should be updated   
     156      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )    
     157      ! 
     158      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
     159 
     160      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zgrad') 
     161 
     162   END SUBROUTINE bdy_dyn3d_zgrad 
    112163 
    113164   SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy ) 
     
    296347   END SUBROUTINE bdy_dyn3d_dmp 
    297348 
     349   SUBROUTINE bdy_dyn3d_nmn( idx, ib_bdy ) 
     350      !!---------------------------------------------------------------------- 
     351      !!                 ***  SUBROUTINE bdy_dyn3d_nmn  *** 
     352      !!              
     353      !!              - Apply Neumann condition to baroclinic velocities.  
     354      !!              - Wrapper routine for bdy_nmn 
     355      !!  
     356      !! 
     357      !!---------------------------------------------------------------------- 
     358      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices 
     359      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index 
     360 
     361      INTEGER  ::   jb, igrd                               ! dummy loop indices 
     362      !!---------------------------------------------------------------------- 
     363 
     364      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_nmn') 
     365      ! 
     366      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities.  
     367      ! 
     368      igrd = 2      ! Neumann bc on u-velocity;  
     369      !             
     370      CALL bdy_nmn( idx, igrd, ua ) 
     371 
     372      igrd = 3      ! Neumann bc on v-velocity 
     373      !   
     374      CALL bdy_nmn( idx, igrd, va ) 
     375      ! 
     376      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )    ! Boundary points should be updated 
     377      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy ) 
     378      ! 
     379      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_nmn') 
     380      ! 
     381   END SUBROUTINE bdy_dyn3d_nmn 
     382 
    298383#else 
    299384   !!---------------------------------------------------------------------- 
  • branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r6140 r7029  
    100100         &             cn_ice_lim, nn_ice_lim_dta,                           & 
    101101         &             rn_ice_tem, rn_ice_sal, rn_ice_age,                 & 
    102          &             ln_vol, nn_volctl, nn_rimwidth 
     102         &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy 
    103103         ! 
    104104      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend 
     
    209209             dta_bdy(ib_bdy)%ll_u3d = .true. 
    210210             dta_bdy(ib_bdy)%ll_v3d = .true. 
     211          CASE('neumann') 
     212             IF(lwp) WRITE(numout,*) '      Neumann conditions' 
     213             dta_bdy(ib_bdy)%ll_u3d = .false. 
     214             dta_bdy(ib_bdy)%ll_v3d = .false. 
     215          CASE('zerograd') 
     216             IF(lwp) WRITE(numout,*) '      Zero gradient for baroclinic velocities' 
     217             dta_bdy(ib_bdy)%ll_u3d = .false. 
     218             dta_bdy(ib_bdy)%ll_v3d = .false. 
    211219          CASE('zero') 
    212220             IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
     
    377385          IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 
    378386          IF(lwp) WRITE(numout,*) 
     387        ENDIF 
     388        IF( nb_jpk_bdy > 0 ) THEN 
     389           IF(lwp) WRITE(numout,*) '*** open boundary will be interpolate in the vertical onto the native grid ***' 
     390        ELSE 
     391           IF(lwp) WRITE(numout,*) '*** open boundary will be read straight onto the native grid without vertical interpolation ***' 
    379392        ENDIF 
    380393     ENDIF 
     
    499512            &      nbrdta(jpbdta, jpbgrd, nb_bdy) ) 
    500513 
    501          ALLOCATE( dta_global(jpbdtau, 1, jpk) ) 
    502          IF ( icount>0 ) ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) ) 
     514         IF( nb_jpk_bdy>0 ) THEN 
     515            ALLOCATE( dta_global(jpbdtau, 1, nb_jpk_bdy) ) 
     516            ALLOCATE( dta_global_z(jpbdtau, 1, nb_jpk_bdy) ) 
     517            ALLOCATE( dta_global_dz(jpbdtau, 1, nb_jpk_bdy) ) 
     518         ELSE 
     519            ALLOCATE( dta_global(jpbdtau, 1, jpk) ) 
     520            ALLOCATE( dta_global_z(jpbdtau, 1, jpk) ) ! needed ?? TODO 
     521            ALLOCATE( dta_global_dz(jpbdtau, 1, jpk) )! needed ?? TODO 
     522         ENDIF 
     523 
     524         IF ( icount>0 ) THEN 
     525            IF( nb_jpk_bdy>0 ) THEN 
     526               ALLOCATE( dta_global2(jpbdtas, nrimmax, nb_jpk_bdy) ) 
     527               ALLOCATE( dta_global2_z(jpbdtas, nrimmax, nb_jpk_bdy) ) 
     528               ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, nb_jpk_bdy) ) 
     529            ELSE 
     530               ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) ) 
     531               ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk) ) ! needed ?? TODO 
     532               ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk) )! needed ?? TODO   
     533            ENDIF 
     534         ENDIF 
    503535         !  
    504536      ENDIF 
     
    10921124      !          = 0  elsewhere    
    10931125  
     1126      bdytmask(:,:) = ssmask(:,:) 
     1127 
    10941128      IF( ln_mask_file ) THEN 
    10951129         CALL iom_open( cn_mask_file, inum ) 
     
    11081142         CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond. 
    11091143 
    1110  
    1111          ! Mask corrections 
    1112          ! ---------------- 
    1113          DO ik = 1, jpkm1 
    1114             DO ij = 1, jpj 
    1115                DO ii = 1, jpi 
    1116                   tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij) 
    1117                   umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij) 
    1118                   vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij) 
    1119                END DO       
    1120             END DO 
    1121             DO ij = 2, jpjm1 
    1122                DO ii = 2, jpim1 
    1123                   fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   & 
    1124                      &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1) 
    1125                END DO       
    1126             END DO 
    1127          END DO 
    1128          tmask_i (:,:) = ssmask(:,:) * tmask_i(:,:) 
    1129          ! 
    11301144      ENDIF ! ln_mask_file=.TRUE. 
    11311145       
    1132       bdytmask(:,:) = ssmask(:,:) 
    11331146      IF( .NOT.ln_mask_file ) THEN 
    11341147         ! If .not. ln_mask_file then we need to derive mask on U and V grid from mask on T grid here. 
  • branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90

    r6140 r7029  
    2727   PUBLIC   bdy_orlanski_2d     ! routine called where? 
    2828   PUBLIC   bdy_orlanski_3d     ! routine called where? 
     29   PUBLIC   bdy_nmn             ! routine called where? 
    2930 
    3031   !!---------------------------------------------------------------------- 
     
    355356   END SUBROUTINE bdy_orlanski_3d 
    356357 
     358   SUBROUTINE bdy_nmn( idx, igrd, phia ) 
     359      !!---------------------------------------------------------------------- 
     360      !!                 ***  SUBROUTINE bdy_nmn  *** 
     361      !!                     
     362      !! ** Purpose : Duplicate the value at open boundaries, zero gradient. 
     363      !!  
     364      !!---------------------------------------------------------------------- 
     365      INTEGER,                    INTENT(in)     ::   igrd     ! grid index 
     366      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   phia     ! model after 3D field (to be updated) 
     367      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
     368      !!  
     369      REAL(wp) ::   zcoef, zcoef1, zcoef2 
     370      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask      ! land/sea mask for field 
     371      REAL(wp), POINTER, DIMENSION(:,:)        :: bdypmask      ! land/sea mask for field 
     372      INTEGER  ::   ib, ik   ! dummy loop indices 
     373      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses 
     374      !!---------------------------------------------------------------------- 
     375      ! 
     376      IF( nn_timing == 1 ) CALL timing_start('bdy_nmn') 
     377      ! 
     378      SELECT CASE(igrd) 
     379         CASE(1) 
     380            pmask => tmask(:,:,:) 
     381            bdypmask => bdytmask(:,:) 
     382         CASE(2) 
     383            pmask => umask(:,:,:) 
     384            bdypmask => bdyumask(:,:) 
     385         CASE(3) 
     386            pmask => vmask(:,:,:) 
     387            bdypmask => bdyvmask(:,:) 
     388         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_nmn' ) 
     389      END SELECT 
     390      DO ib = 1, idx%nblenrim(igrd) 
     391         ii = idx%nbi(ib,igrd) 
     392         ij = idx%nbj(ib,igrd) 
     393         DO ik = 1, jpkm1 
     394            ! search the sense of the gradient 
     395            zcoef1 = bdypmask(ii-1,ij  )*pmask(ii-1,ij,ik) +  bdypmask(ii+1,ij  )*pmask(ii+1,ij,ik) 
     396            zcoef2 = bdypmask(ii  ,ij-1)*pmask(ii,ij-1,ik) +  bdypmask(ii  ,ij+1)*pmask(ii,ij+1,ik) 
     397            IF ( nint(zcoef1+zcoef2) == 0) THEN 
     398               ! corner **** we probably only want to set the tangentail component for the dynamics here 
     399               zcoef = pmask(ii-1,ij,ik) + pmask(ii+1,ij,ik) +  pmask(ii,ij-1,ik) +  pmask(ii,ij+1,ik) 
     400               IF (zcoef > .5_wp) THEN ! Only set none isolated points. 
     401                 phia(ii,ij,ik) = phia(ii-1,ij  ,ik) * pmask(ii-1,ij  ,ik) + & 
     402                   &              phia(ii+1,ij  ,ik) * pmask(ii+1,ij  ,ik) + & 
     403                   &              phia(ii  ,ij-1,ik) * pmask(ii  ,ij-1,ik) + & 
     404                   &              phia(ii  ,ij+1,ik) * pmask(ii  ,ij+1,ik) 
     405                 phia(ii,ij,ik) = ( phia(ii,ij,ik) / zcoef ) * pmask(ii,ij,ik) 
     406               ELSE 
     407                 phia(ii,ij,ik) = phia(ii,ij  ,ik) * pmask(ii,ij  ,ik) 
     408               ENDIF 
     409            ELSEIF ( nint(zcoef1+zcoef2) == 2) THEN 
     410               ! oblique corner **** we probably only want to set the normal component for the dynamics here 
     411               zcoef = pmask(ii-1,ij,ik)*bdypmask(ii-1,ij  ) + pmask(ii+1,ij,ik)*bdypmask(ii+1,ij  ) + & 
     412                   &   pmask(ii,ij-1,ik)*bdypmask(ii,ij -1 ) +  pmask(ii,ij+1,ik)*bdypmask(ii,ij+1  ) 
     413               phia(ii,ij,ik) = phia(ii-1,ij  ,ik) * pmask(ii-1,ij  ,ik)*bdypmask(ii-1,ij  ) + & 
     414                   &            phia(ii+1,ij  ,ik) * pmask(ii+1,ij  ,ik)*bdypmask(ii+1,ij  )  + & 
     415                   &            phia(ii  ,ij-1,ik) * pmask(ii  ,ij-1,ik)*bdypmask(ii,ij -1 ) + & 
     416                   &            phia(ii  ,ij+1,ik) * pmask(ii  ,ij+1,ik)*bdypmask(ii,ij+1  ) 
     417  
     418               phia(ii,ij,ik) = ( phia(ii,ij,ik) / MAX(1._wp, zcoef) ) * pmask(ii,ij,ik) 
     419            ELSE 
     420               ip = nint(bdypmask(ii+1,ij  )*pmask(ii+1,ij,ik) - bdypmask(ii-1,ij  )*pmask(ii-1,ij,ik)) 
     421               jp = nint(bdypmask(ii  ,ij+1)*pmask(ii,ij+1,ik) - bdypmask(ii  ,ij-1)*pmask(ii,ij-1,ik)) 
     422               phia(ii,ij,ik) = phia(ii+ip,ij+jp,ik) * pmask(ii+ip,ij+jp,ik) 
     423            ENDIF 
     424         END DO 
     425      END DO 
     426      ! 
     427      IF( nn_timing == 1 ) CALL timing_stop('bdy_nmn') 
     428      ! 
     429   END SUBROUTINE bdy_nmn 
    357430 
    358431#else 
     
    367440      WRITE(*,*) 'bdy_orlanski_3d: You should not have seen this print! error?', kt 
    368441   END SUBROUTINE bdy_orlanski_3d 
     442   SUBROUTINE bdy_nmn( idx, igrd, phia )      ! Empty routine 
     443      WRITE(*,*) 'bdy_nmn: You should not have seen this print! error?', kt 
     444   END SUBROUTINE bdy_nmn 
    369445#endif 
    370446 
  • branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90

    r6140 r7029  
    155155      ! 
    156156      REAL(wp) ::   zwgt           ! boundary weight 
    157       INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
    158       INTEGER  ::   ii, ij,zcoef, zcoef1,zcoef2, ip, jp   ! 2D addresses 
     157      REAL(wp) ::   zcoef, zcoef1,zcoef2 
     158      INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
     159      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses 
    159160      !!---------------------------------------------------------------------- 
    160161      ! 
     
    169170            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  ) 
    170171            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1) 
    171             IF ( zcoef1+zcoef2 == 0) THEN 
     172            IF ( NINT(zcoef1+zcoef2) == 0) THEN 
    172173               ! corner 
    173174               zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) +  tmask(ii,ij-1,ik) +  tmask(ii,ij+1,ik) 
     
    176177                 &                    tsa(ii  ,ij-1,ik,jp_tem) * tmask(ii  ,ij-1,ik) + & 
    177178                 &                    tsa(ii  ,ij+1,ik,jp_tem) * tmask(ii  ,ij+1,ik) 
    178                tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1, zcoef) ) * tmask(ii,ij,ik) 
     179               tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1._wp, zcoef) ) * tmask(ii,ij,ik) 
    179180               tsa(ii,ij,ik,jp_sal) = tsa(ii-1,ij  ,ik,jp_sal) * tmask(ii-1,ij  ,ik) + & 
    180181                 &                    tsa(ii+1,ij  ,ik,jp_sal) * tmask(ii+1,ij  ,ik) + & 
    181182                 &                    tsa(ii  ,ij-1,ik,jp_sal) * tmask(ii  ,ij-1,ik) + & 
    182183                 &                    tsa(ii  ,ij+1,ik,jp_sal) * tmask(ii  ,ij+1,ik) 
    183                tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1, zcoef) ) * tmask(ii,ij,ik) 
     184               tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1._wp, zcoef) ) * tmask(ii,ij,ik) 
    184185            ELSE 
    185                ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  ) 
    186                jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1) 
     186               ip = NINT(bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )) 
     187               jp = NINT(bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)) 
    187188               tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii+ip,ij+jp,ik) 
    188189               tsa(ii,ij,ik,jp_sal) = tsa(ii+ip,ij+jp,ik,jp_sal) * tmask(ii+ip,ij+jp,ik) 
  • branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r6140 r7029  
    2424   USE oce             ! ocean dynamics and tracers 
    2525   USE dom_oce         ! ocean space and time domain 
    26    ! 
     26   USE bdy_oce       
    2727   USE in_out_manager  ! I/O manager 
    2828   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2929   USE lib_mpp         ! 
     30   USE iom 
    3031   USE wrk_nemo        ! Memory allocation 
    3132   USE timing          ! Timing 
     
    102103      INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers 
    103104      INTEGER  ::   ijf, ijl, ij0, ij1       !   -       - 
    104       INTEGER  ::   ios 
     105      INTEGER  ::   ios, inum 
    105106      INTEGER  ::   isrow                    ! index for ORCA1 starting row 
    106107      INTEGER , POINTER, DIMENSION(:,:) ::  imsk 
     
    108109      !! 
    109110      NAMELIST/namlbc/ rn_shlat, ln_vorlat 
     111#if defined key_bdy  
     112      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,                 & 
     113         &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     & 
     114         &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             & 
     115         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, & 
     116         &             cn_ice_lim, nn_ice_lim_dta,                           & 
     117         &             rn_ice_tem, rn_ice_sal, rn_ice_age,                 & 
     118         &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy 
     119#endif 
    110120      !!--------------------------------------------------------------------- 
    111121      ! 
     
    155165      END DO   
    156166       
     167#if defined key_bdy  
     168      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist :Unstructured open boundaries   
     169      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903) 
     170903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp ) 
     171 
     172      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries 
     173      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 ) 
     174904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp ) 
     175      IF(lwm) WRITE ( numond, nambdy ) 
     176 
     177     IF( ln_mask_file ) THEN ! correct for bdy mask 
     178         CALL iom_open( cn_mask_file, inum ) 
     179         CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) ) 
     180         CALL iom_close( inum ) 
     181 
     182         ! Mask corrections 
     183         ! ---------------- 
     184         DO jk = 1, jpkm1 
     185            DO jj = 1, jpj 
     186               DO ji = 1, jpi 
     187                  tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj) 
     188               END DO 
     189            END DO 
     190         END DO 
     191      ENDIF 
     192#endif 
    157193      ! (ISF) define barotropic mask and mask the ice shelf point 
    158194      ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked 
  • branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90

    r5328 r7029  
    2424   USE wrk_nemo        ! Memory Allocation 
    2525   USE timing          ! Timing 
     26#if defined key_bdy  
     27   USE bdy_oce         ! ocean open boundary conditions 
     28#endif 
    2629 
    2730   IMPLICIT NONE 
     
    7881      REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke 
    7982      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv  
     83#if defined key_bdy 
     84      INTEGER  ::   jb                 ! dummy loop indices 
     85      INTEGER  ::   ii, ij, igrd, ib_bdy   ! local integers 
     86      INTEGER  ::   fu, fv 
     87#endif 
    8088      !!---------------------------------------------------------------------- 
    8189      ! 
     
    98106      zhke(:,:,jpk) = 0._wp 
    99107       
     108#if defined key_bdy 
     109      ! Maria Luneva & Fred Wobus: July-2016 
     110      ! compensate for lack of turbulent kinetic energy on liquid bdy points 
     111      DO ib_bdy = 1, nb_bdy 
     112         IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 
     113            igrd = 2           ! Copying normal velocity into points outside bdy 
     114            DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     115               DO jk = 1, jpkm1 
     116                  ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
     117                  ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
     118                  fu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) ) 
     119                  un(ii-fu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk) 
     120               END DO 
     121            END DO 
     122            ! 
     123            igrd = 3           ! Copying normal velocity into points outside bdy 
     124            DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     125               DO jk = 1, jpkm1 
     126                  ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
     127                  ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
     128                  fv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) ) 
     129                  vn(ii,ij-fv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk) 
     130               END DO 
     131            END DO 
     132         ENDIF 
     133      ENDDO   
     134#endif  
     135 
    100136      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==! 
    101137      ! 
     
    133169         ! 
    134170      END SELECT 
     171 
     172#if defined key_bdy 
     173      ! restore velocity masks at points outside boundary 
     174      un(:,:,:) = un(:,:,:) * umask(:,:,:) 
     175      vn(:,:,:) = vn(:,:,:) * vmask(:,:,:) 
     176#endif       
     177 
     178 
    135179      ! 
    136180      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==! 
  • branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90

    r6519 r7029  
    99   !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  add C1D case   
    1010   !!            3.6  ! 2014-15  DIMG format removed 
     11   !!            3.6  ! 2015-15  (J. Harle) Added procedure to read REAL attributes 
    1112   !!-------------------------------------------------------------------- 
    1213 
     
    6768   END INTERFACE 
    6869   INTERFACE iom_getatt 
    69       MODULE PROCEDURE iom_g0d_intatt 
     70      MODULE PROCEDURE iom_g0d_intatt, iom_g0d_ratt 
    7071   END INTERFACE 
    7172   INTERFACE iom_rstput 
     
    979980         IF( iom_file(kiomid)%nfid > 0 ) THEN 
    980981            SELECT CASE (iom_file(kiomid)%iolib) 
    981             CASE (jpnf90   )   ;   CALL iom_nf90_getatt( kiomid, cdatt, pvar ) 
     982            CASE (jpnf90   )   ;   CALL iom_nf90_getatt( kiomid, cdatt, pv_i0d=pvar ) 
    982983            CASE DEFAULT 
    983984               CALL ctl_stop( 'iom_g0d_att: accepted IO library is only jpnf90' ) 
     
    987988   END SUBROUTINE iom_g0d_intatt 
    988989 
     990   SUBROUTINE iom_g0d_ratt( kiomid, cdatt, pvar, cdvar ) 
     991      INTEGER         , INTENT(in   )                 ::   kiomid    ! Identifier of the file 
     992      CHARACTER(len=*), INTENT(in   )                 ::   cdatt     ! Name of the attribute 
     993      REAL(wp)        , INTENT(  out)                 ::   pvar      ! written field 
     994      CHARACTER(len=*), INTENT(in   ), OPTIONAL       ::   cdvar     ! Name of the variable 
     995      ! 
     996      IF( kiomid > 0 ) THEN 
     997         IF( iom_file(kiomid)%nfid > 0 ) THEN 
     998            SELECT CASE (iom_file(kiomid)%iolib) 
     999            CASE (jpnf90   )   ;   IF( PRESENT(cdvar) ) THEN 
     1000                                      CALL iom_nf90_getatt( kiomid, cdatt, pv_r0d=pvar, cdvar=cdvar ) 
     1001                                   ELSE 
     1002                                      CALL iom_nf90_getatt( kiomid, cdatt, pv_r0d=pvar ) 
     1003                                   ENDIF 
     1004            CASE DEFAULT     
     1005               CALL ctl_stop( 'iom_g0d_att: accepted IO library is only jpnf90' ) 
     1006            END SELECT 
     1007         ENDIF 
     1008      ENDIF 
     1009   END SUBROUTINE iom_g0d_ratt 
    9891010 
    9901011   !!---------------------------------------------------------------------- 
  • branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90

    r6140 r7029  
    77   !!            9.0  ! 06 02  (S. Masson) Adaptation to NEMO 
    88   !!             "   ! 07 07  (D. Storkey) Changes to iom_nf90_gettime 
     9   !!            3.6  ! 2015-15  (J. Harle) Added procedure to read REAL attributes 
    910   !!-------------------------------------------------------------------- 
    1011   !!gm  caution add !DIR nec: improved performance to be checked as well as no result changes 
     
    3536   END INTERFACE 
    3637   INTERFACE iom_nf90_getatt 
    37       MODULE PROCEDURE iom_nf90_intatt 
     38      MODULE PROCEDURE iom_nf90_att 
    3839   END INTERFACE 
    3940   INTERFACE iom_nf90_rstput 
     
    313314 
    314315 
    315    SUBROUTINE iom_nf90_intatt( kiomid, cdatt, pvar ) 
    316       !!----------------------------------------------------------------------- 
    317       !!                  ***  ROUTINE  iom_nf90_intatt  *** 
     316   SUBROUTINE iom_nf90_att( kiomid, cdatt, pv_i0d, pv_r0d, cdvar) 
     317      !!----------------------------------------------------------------------- 
     318      !!                  ***  ROUTINE  iom_nf90_att  *** 
    318319      !! 
    319320      !! ** Purpose : read an integer attribute with NF90 
     
    321322      INTEGER         , INTENT(in   ) ::   kiomid   ! Identifier of the file 
    322323      CHARACTER(len=*), INTENT(in   ) ::   cdatt    ! attribute name 
    323       INTEGER         , INTENT(  out) ::   pvar     ! read field 
     324      INTEGER         , INTENT(  out), OPTIONAL       ::   pv_i0d    ! read field 
     325      REAL(wp),         INTENT(  out), OPTIONAL       ::   pv_r0d    ! read field  
     326      CHARACTER(len=*), INTENT(in   ), OPTIONAL       ::   cdvar     ! name of the variable 
    324327      ! 
    325328      INTEGER                         ::   if90id   ! temporary integer 
     329      INTEGER                         ::   ivarid           ! NetCDF  variable Id 
    326330      LOGICAL                         ::   llok     ! temporary logical 
    327331      CHARACTER(LEN=100)              ::   clinfo   ! info character 
     
    329333      !  
    330334      if90id = iom_file(kiomid)%nfid 
    331       llok = NF90_Inquire_attribute(if90id, NF90_GLOBAL, cdatt) == nf90_noerr 
     335      IF( PRESENT(cdvar) ) THEN 
     336         llok = NF90_INQ_VARID( if90id, TRIM(cdvar), ivarid ) == nf90_noerr   ! does the variable exist in the file 
     337         IF( llok ) THEN 
     338            llok = NF90_Inquire_attribute(if90id, ivarid, cdatt) == nf90_noerr 
     339         ELSE 
     340            CALL ctl_warn('iom_nf90_getatt: no variable '//cdvar//' found') 
     341         ENDIF 
     342      ELSE 
     343         llok = NF90_Inquire_attribute(if90id, NF90_GLOBAL, cdatt) == nf90_noerr 
     344      ENDIF  
     345! 
    332346      IF( llok) THEN 
    333347         clinfo = 'iom_nf90_getatt, file: '//TRIM(iom_file(kiomid)%name)//', att: '//TRIM(cdatt) 
    334          CALL iom_nf90_check(NF90_GET_ATT(if90id, NF90_GLOBAL, cdatt, values=pvar), clinfo) 
     348         IF(     PRESENT(pv_r0d) ) THEN 
     349            IF( PRESENT(cdvar) ) THEN 
     350               CALL iom_nf90_check(NF90_GET_ATT(if90id, ivarid, cdatt, values=pv_r0d), clinfo) 
     351            ELSE 
     352               CALL iom_nf90_check(NF90_GET_ATT(if90id, NF90_GLOBAL, cdatt, values=pv_r0d), clinfo) 
     353            ENDIF 
     354         ELSE 
     355            IF( PRESENT(cdvar) ) THEN 
     356               CALL iom_nf90_check(NF90_GET_ATT(if90id, ivarid, cdatt, values=pv_i0d), clinfo) 
     357            ELSE 
     358               CALL iom_nf90_check(NF90_GET_ATT(if90id, NF90_GLOBAL, cdatt, values=pv_i0d), clinfo) 
     359            ENDIF 
     360         ENDIF 
    335361      ELSE 
    336362         CALL ctl_warn('iom_nf90_getatt: no attribute '//cdatt//' found') 
    337          pvar = -999 
     363         IF(     PRESENT(pv_r0d) ) THEN 
     364            pv_r0d = -999._wp 
     365         ELSE 
     366            pv_i0d = -999 
     367         ENDIF 
    338368      ENDIF 
    339369      !  
    340    END SUBROUTINE iom_nf90_intatt 
     370   END SUBROUTINE iom_nf90_att 
    341371 
    342372 
  • branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90

    r6412 r7029  
    4141      USE in_out_manager  ! I/O Manager 
    4242      USE iom 
     43#if defined key_bdy 
     44      USE bdy_oce 
     45#endif 
    4346      !!  
    4447      INTEGER :: ji, jj, jn, jproc, jarea     ! dummy loop indices 
     
    7376      ! read namelist for ln_zco 
    7477      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh 
    75  
     78#if defined key_bdy 
     79      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,                 & 
     80         &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     & 
     81         &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &   
     82         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, & 
     83         &             cn_ice_lim, nn_ice_lim_dta,                           & 
     84         &             rn_ice_tem, rn_ice_sal, rn_ice_age,                 & 
     85         &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy 
     86#endif 
    7687      !!---------------------------------------------------------------------- 
    7788      !!  OPA 9.0 , LOCEAN-IPSL (2005)  
     
    137148      imask(:,:)=1 
    138149      WHERE ( zdta(:,:) - zdtaisf(:,:) <= rn_isfhmin ) imask = 0 
     150 
     151#if defined key_bdy 
     152      ! Adjust imask with bdy_msk if exists 
     153 
     154      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist : BDY 
     155      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903) 
     156903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist (mppini_2)', lwp ) 
     157 
     158      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist : BDY 
     159      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 ) 
     160904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist (mppini_2)', lwp ) 
     161      IF(lwm) WRITE ( numond, namzgr ) 
     162 
     163      IF( ln_mask_file ) THEN 
     164         CALL iom_open( cn_mask_file, inum ) 
     165         CALL iom_get ( inum, jpdom_unknown, 'bdy_msk', zdta(:,:), kstart=(/jpizoom,jpjzoom/), kcount=(/jpiglo,jpjglo/) ) 
     166         CALL iom_close( inum ) 
     167         WHERE ( zdta(:,:) <= 0. ) imask = 0 
     168      ENDIF 
     169#endif 
    139170 
    140171      !  1. Dimension arrays for subdomains 
  • branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90

    r6140 r7029  
    77   !!                 !  05-2008  (S. Alderson)  Modified for Interpolation in memory from input grid to model grid 
    88   !!                 !  10-2013  (D. Delrosso, P. Oddo)  suppression of land point prior to interpolation 
     9   !!                 !  12-2015  (J. Harle) Adding BDY on-the-fly interpolation 
    910   !!---------------------------------------------------------------------- 
    1011 
     
    6768      INTEGER                         ::   nreclast     ! last record to be read in the current file 
    6869      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 
    6972   END TYPE FLD 
    7073 
     
    114117CONTAINS 
    115118 
    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 ) 
    117120      !!--------------------------------------------------------------------- 
    118121      !!                    ***  ROUTINE fld_read  *** 
     
    135138      !                                                     !   kt_offset = +1 => fields at "after"  time level 
    136139      !                                                     !   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      !! 
    137143      INTEGER  ::   itmp         ! local variable 
    138144      INTEGER  ::   imf          ! size of the structure sd 
     
    171177         DO jf = 1, imf  
    172178            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 
    174184         END DO 
    175185         IF( lwp ) CALL wgt_print()                ! control print 
     
    263273 
    264274               ! 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 
    267280            ENDIF   ! read new data? 
    268281         END DO                                    ! --- end loop over field --- ! 
     
    302315 
    303316 
    304    SUBROUTINE fld_init( kn_fsbc, sdjf, map ) 
     317   SUBROUTINE fld_init( kn_fsbc, sdjf, map , jpk_bdy, fvl) 
    305318      !!--------------------------------------------------------------------- 
    306319      !!                    ***  ROUTINE fld_init  *** 
     
    309322      !!               - if time interpolation, read before data  
    310323      !!---------------------------------------------------------------------- 
    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 
    314329      !! 
    315330      LOGICAL :: llprevyr              ! are we reading previous year  file? 
     
    405420         ENDIF 
    406421         ! 
    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 
    408428         ! 
    409429         clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')" 
     
    581601 
    582602 
    583    SUBROUTINE fld_get( sdjf, map ) 
     603   SUBROUTINE fld_get( sdjf, map, jpk_bdy, fvl ) 
    584604      !!--------------------------------------------------------------------- 
    585605      !!                    ***  ROUTINE fld_get  *** 
     
    589609      TYPE(FLD)        , INTENT(inout) ::   sdjf   ! input field related variables 
    590610      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 
    591613      ! 
    592614      INTEGER ::   ipk      ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
     
    601623      ! 
    602624      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         
    606636      ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 
    607637         CALL wgt_list( sdjf, iw ) 
     
    658688   END SUBROUTINE fld_get 
    659689 
    660  
    661    SUBROUTINE fld_map( num, clvar, dta, nrec, map ) 
     690   SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy, fvl ) 
    662691      !!--------------------------------------------------------------------- 
    663692      !!                    ***  ROUTINE fld_map  *** 
     
    667696      !!---------------------------------------------------------------------- 
    668697#if defined key_bdy 
    669       USE bdy_oce, ONLY:  dta_global, dta_global2         ! workspace to read in global data arrays 
     698      USE bdy_oce, ONLY:  idx_bdy, dta_global, dta_global_z, dta_global_dz, dta_global2, dta_global2_z, dta_global2_dz                 ! workspace to read in global data arrays 
    670699#endif  
    671700      INTEGER                   , INTENT(in ) ::   num     ! stream number 
    672701      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) 
    674703      INTEGER                   , INTENT(in ) ::   nrec    ! record number to read (ie time slice) 
    675704      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 
    676708      !! 
    677709      INTEGER                                 ::   ipi      ! length of boundary data on local process 
     
    682714      INTEGER                                 ::   ib, ik, ji, jj   ! loop counters 
    683715      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 
    685720      !!--------------------------------------------------------------------- 
    686721      ! 
     
    696731      IF( map%ll_unstruc) THEN   ! unstructured open boundary data file 
    697732         dta_read => dta_global 
    698       ELSE                       ! structured open boundary data file 
     733         IF( PRESENT(jpk_bdy) ) THEN 
     734            IF( jpk_bdy>0 ) THEN 
     735               dta_read_z => dta_global_z 
     736               dta_read_dz => dta_global_dz 
     737               jpkm1_bdy = jpk_bdy-1 
     738            ENDIF 
     739         ENDIF 
     740      ELSE                       ! structured open boundary file 
    699741         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 
    700749      ENDIF 
    701750#endif 
     
    705754      ! 
    706755      SELECT CASE( ipk ) 
    707       CASE(1)        ;   CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1    ), nrec ) 
    708       CASE DEFAULT   ;   CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec ) 
     756      CASE(1)        ;    
     757      CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1    ), nrec ) 
     758         IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 
     759            DO ib = 1, ipi 
     760              DO ik = 1, ipk 
     761                dta(ib,1,ik) =  dta_read(map%ptr(ib),1,ik) 
     762              END DO 
     763            END DO 
     764         ELSE ! we assume that this is a structured open boundary file 
     765            DO ib = 1, ipi 
     766               jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
     767               ji=map%ptr(ib)-(jj-1)*ilendta 
     768               DO ik = 1, ipk 
     769                  dta(ib,1,ik) =  dta_read(ji,jj,ik) 
     770               END DO 
     771            END DO 
     772         ENDIF 
     773 
     774      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     775      ! Do we include something here to adjust barotropic velocities ! 
     776      ! in case of a depth difference between bdy files and          ! 
     777      ! bathymetry in the case ln_full_vel = .false. and jpk_bdy>0?  ! 
     778      ! [as the enveloping and parital cells could change H]         ! 
     779      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     780 
     781      CASE DEFAULT   ; 
     782 
     783      IF( PRESENT(jpk_bdy) .AND. jpk_bdy>0 ) THEN       ! boundary data not on model grid: vertical interpolation 
     784         CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar ) 
     785         dta_read(:,:,:) = -ABS(fv) 
     786         dta_read_z(:,:,:) = 0._wp 
     787         dta_read_dz(:,:,:) = 0._wp 
     788         CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:jpk_bdy), nrec ) 
     789         SELECT CASE( igrd )                   
     790            CASE(1) 
     791               CALL iom_get ( num, jpdom_unknown, 'gdept', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 
     792               CALL iom_get ( num, jpdom_unknown, 'e3t',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 
     793            CASE(2)   
     794               CALL iom_get ( num, jpdom_unknown, 'gdepu', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 
     795               CALL iom_get ( num, jpdom_unknown, 'e3u',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 
     796            CASE(3) 
     797               CALL iom_get ( num, jpdom_unknown, 'gdepv', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 
     798               CALL iom_get ( num, jpdom_unknown, 'e3v',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 
     799         END SELECT 
     800 
     801#if defined key_bdy 
     802         CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta) 
     803#endif 
     804      ELSE ! boundary data assumed to be on model grid 
     805          
     806         CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec )                     
     807         IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 
     808            DO ib = 1, ipi 
     809              DO ik = 1, ipk 
     810                dta(ib,1,ik) =  dta_read(map%ptr(ib),1,ik) 
     811              END DO 
     812            END DO 
     813         ELSE ! we assume that this is a structured open boundary file 
     814            DO ib = 1, ipi 
     815               jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
     816               ji=map%ptr(ib)-(jj-1)*ilendta 
     817               DO ik = 1, ipk 
     818                  dta(ib,1,ik) =  dta_read(ji,jj,ik) 
     819               END DO 
     820            END DO 
     821         ENDIF 
     822      ENDIF ! PRESENT(jpk_bdy) 
    709823      END SELECT 
    710       ! 
    711       IF( map%ll_unstruc ) THEN  ! unstructured open boundary data file 
     824 
     825   END SUBROUTINE fld_map 
     826    
     827#if defined key_bdy 
     828   SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta) 
     829 
     830      !!--------------------------------------------------------------------- 
     831      !!                    ***  ROUTINE fld_bdy_interp  *** 
     832      !! 
     833      !! ** Purpose :   on the fly vertical interpolation to allow the use of 
     834      !!                boundary data from non-native vertical grid 
     835      !!---------------------------------------------------------------------- 
     836      USE bdy_oce, ONLY:  idx_bdy         ! indexing for map <-> ij transformation 
     837 
     838      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read      ! work space for global data 
     839      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_z    ! work space for global data 
     840      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_dz   ! work space for global data 
     841      REAL(wp) , INTENT(in)                                ::   fv            ! fillvalue and alternative -ABS(fv) 
     842      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta                        ! output field on model grid (2 dimensional) 
     843      TYPE(MAP_POINTER)         , INTENT(in ) ::   map                        ! global-to-local mapping indices 
     844      LOGICAL  , INTENT(in), OPTIONAL         ::   fvl                        ! grid type, set number and number of vertical levels in the bdy data 
     845      INTEGER  , INTENT(in)                   ::   igrd, ibdy, jpk_bdy        ! number of levels in bdy data 
     846      INTEGER  , INTENT(in)                   ::   ilendta                    ! length of data in file 
     847      !! 
     848      INTEGER                                 ::   ipi                        ! length of boundary data on local process 
     849      INTEGER                                 ::   ipj                        ! length of dummy dimension ( = 1 ) 
     850      INTEGER                                 ::   ipk                        ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
     851      INTEGER                                 ::   jpkm1_bdy                  ! number of levels in bdy data minus 1 
     852      INTEGER                                 ::   ib, ik, ikk                ! loop counters 
     853      INTEGER                                 ::   ji, jj, zij, zjj           ! temporary indices 
     854      REAL(wp)                                ::   zl, zi, zh                 ! tmp variable for current depth and interpolation factor 
     855      REAL(wp)                                ::   fv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(fv) 
     856      CHARACTER (LEN=10)                      ::   ibstr 
     857      !!--------------------------------------------------------------------- 
     858      
     859 
     860      ipi       = SIZE( dta, 1 ) 
     861      ipj       = SIZE( dta_read, 2 ) 
     862      ipk       = SIZE( dta, 3 ) 
     863      jpkm1_bdy = jpk_bdy-1 
     864       
     865      fv_alt = -ABS(fv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 
     866      DO ib = 1, ipi 
     867            zij = idx_bdy(ibdy)%nbi(ib,igrd) 
     868            zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
     869         IF(narea==2) WRITE(*,*) 'MAPI', ib, igrd, map%ptr(ib), narea-1, zij, zjj 
     870      ENDDO 
     871      ! 
     872      IF ( map%ll_unstruc ) THEN                            ! unstructured open boundary data file 
     873 
    712874         DO ib = 1, ipi 
    713             DO ik = 1, ipk 
    714                dta(ib,1,ik) =  dta_read(map%ptr(ib),1,ik) 
     875            DO ik = 1, jpk_bdy 
     876               IF( ( dta_read(map%ptr(ib),1,ik) == fv ) ) THEN 
     877                  dta_read_z(map%ptr(ib),1,ik)  = fv_alt ! safety: put fillvalue into external depth field so consistent with data 
     878                  dta_read_dz(map%ptr(ib),1,ik) = 0._wp  ! safety: put 0._wp into external thickness factors to ensure transport is correct 
     879               ENDIF 
     880            ENDDO 
     881         ENDDO  
     882 
     883         DO ib = 1, ipi 
     884            zij = idx_bdy(ibdy)%nbi(ib,igrd) 
     885            zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
     886            zh  = SUM(dta_read_dz(map%ptr(ib),1,:) ) 
     887            ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 
     888            SELECT CASE( igrd )                          
     889               CASE(1) 
     890                  IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 
     891                     WRITE(ibstr,"(I10.10)") map%ptr(ib)  
     892                     CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
     893                 !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1),  ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj 
     894                  ENDIF 
     895               CASE(2) 
     896                  IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 
     897                     WRITE(ibstr,"(I10.10)") map%ptr(ib)  
     898                     CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
     899                     IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u_n(zij,zjj,:), mask=umask(zij,zjj,:)==1),  sum(umask(zij,zjj,:)), & 
     900                       &                hu_n(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1  , & 
     901                        &                dta_read(map%ptr(ib),1,:) 
     902                  ENDIF 
     903               CASE(3) 
     904                  IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 
     905                     WRITE(ibstr,"(I10.10)") map%ptr(ib)  
     906                     CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
     907                  ENDIF 
     908            END SELECT 
     909            DO ik = 1, ipk                       
     910               SELECT CASE( igrd )                        
     911                  CASE(1) 
     912                     zl =  gdept_n(zij,zjj,ik)                                          ! if using in step could use fsdept instead of gdept_n? 
     913                  CASE(2) 
     914                     IF(ln_sco) THEN 
     915                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n? 
     916                     ELSE 
     917                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) )  
     918                     ENDIF 
     919                  CASE(3) 
     920                     IF(ln_sco) THEN 
     921                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n? 
     922                     ELSE 
     923                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) ) 
     924                     ENDIF 
     925               END SELECT 
     926               IF( zl < dta_read_z(map%ptr(ib),1,1) ) THEN                                         ! above the first level of external data 
     927                  dta(ib,1,ik) =  dta_read(map%ptr(ib),1,1) 
     928               ELSEIF( zl > MAXVAL(dta_read_z(map%ptr(ib),1,:),1) ) THEN                           ! below the last level of external data  
     929                  dta(ib,1,ik) =  dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1)) 
     930               ELSE                                                                                ! inbetween : vertical interpolation between ikk & ikk+1 
     931                  DO ikk = 1, jpkm1_bdy                                                            ! when  gdept_n(ikk) < zl < gdept_n(ikk+1) 
     932                     IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) & 
     933                    &    .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN 
     934                        zi           = ( zl - dta_read_z(map%ptr(ib),1,ikk) ) / & 
     935                       &               ( dta_read_z(map%ptr(ib),1,ikk+1) - dta_read_z(map%ptr(ib),1,ikk) ) 
     936                        dta(ib,1,ik) = dta_read(map%ptr(ib),1,ikk) + & 
     937                       &               ( dta_read(map%ptr(ib),1,ikk+1) - dta_read(map%ptr(ib),1,ikk) ) * zi 
     938                     ENDIF 
     939                  END DO 
     940               ENDIF 
    715941            END DO 
    716942         END DO 
    717       ELSE                       ! structured open boundary data file 
     943 
     944         IF(igrd == 2) THEN                                 ! do we need to adjust the transport term? 
     945            DO ib = 1, ipi 
     946              zij = idx_bdy(ibdy)%nbi(ib,igrd) 
     947              zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
     948              zh  = SUM(dta_read_dz(map%ptr(ib),1,:) ) 
     949              ztrans = 0._wp 
     950              ztrans_new = 0._wp 
     951              DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
     952                  ztrans     = ztrans     + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 
     953              ENDDO 
     954              DO ik = 1, ipk                                ! calculate transport on model grid 
     955                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik) 
     956              ENDDO 
     957              DO ik = 1, ipk                                ! make transport correction 
     958                 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
     959                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 
     960                 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
     961                    IF( ABS(ztrans * r1_hu_n(zij,zjj)) > 0.01_wp ) & 
     962                   &   CALL ctl_warn('fld_bdy_interp: barotropic component of > 0.01 ms-1 found in baroclinic velocities at') 
     963                    dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) * umask(zij,zjj,ik) 
     964                 ENDIF 
     965              ENDDO 
     966            ENDDO 
     967         ENDIF 
     968 
     969         IF(igrd == 3) THEN                                 ! do we need to adjust the transport term? 
     970            DO ib = 1, ipi 
     971              zij = idx_bdy(ibdy)%nbi(ib,igrd) 
     972              zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
     973              zh  = SUM(dta_read_dz(map%ptr(ib),1,:) ) 
     974              ztrans = 0._wp 
     975              ztrans_new = 0._wp 
     976              DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
     977                  ztrans     = ztrans     + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 
     978              ENDDO 
     979              DO ik = 1, ipk                                ! calculate transport on model grid 
     980                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik) 
     981              ENDDO 
     982              DO ik = 1, ipk                                ! make transport correction 
     983                 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
     984                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 
     985                 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
     986                    dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) * vmask(zij,zjj,ik) 
     987                 ENDIF 
     988              ENDDO 
     989            ENDDO 
     990         ENDIF 
     991   
     992      ELSE ! structured open boundary file 
     993 
    718994         DO ib = 1, ipi 
    719995            jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    720996            ji=map%ptr(ib)-(jj-1)*ilendta 
    721             DO ik = 1, ipk 
    722                dta(ib,1,ik) =  dta_read(ji,jj,ik) 
     997            DO ik = 1, jpk_bdy                       
     998               IF( ( dta_read(ji,jj,ik) == fv ) ) THEN 
     999                  dta_read_z(ji,jj,ik)  = fv_alt ! safety: put fillvalue into external depth field so consistent with data 
     1000                  dta_read_dz(ji,jj,ik) = 0._wp  ! safety: put 0._wp into external thickness factors to ensure transport is correct 
     1001               ENDIF 
     1002            ENDDO 
     1003         ENDDO  
     1004        
     1005 
     1006         DO ib = 1, ipi 
     1007            jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
     1008            ji=map%ptr(ib)-(jj-1)*ilendta 
     1009            zij = idx_bdy(ibdy)%nbi(ib,igrd) 
     1010            zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
     1011            zh  = SUM(dta_read_dz(ji,jj,:) ) 
     1012            ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 
     1013            SELECT CASE( igrd )                          
     1014               CASE(1) 
     1015                  IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 
     1016                     WRITE(ibstr,"(I10.10)") map%ptr(ib)  
     1017                     CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
     1018                 !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1),  ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj 
     1019                  ENDIF 
     1020               CASE(2) 
     1021                  IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 
     1022                     WRITE(ibstr,"(I10.10)") map%ptr(ib)  
     1023                     CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
     1024                  ENDIF 
     1025               CASE(3) 
     1026                  IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 
     1027                     WRITE(ibstr,"(I10.10)") map%ptr(ib)  
     1028                     CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
     1029                  ENDIF 
     1030            END SELECT 
     1031            DO ik = 1, ipk                       
     1032               SELECT CASE( igrd )                          ! coded for sco - need zco and zps option using min 
     1033                  CASE(1) 
     1034                     zl =  gdept_n(zij,zjj,ik)              ! if using in step could use fsdept instead of gdept_n? 
     1035                  CASE(2) 
     1036                     IF(ln_sco) THEN 
     1037                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n? 
     1038                     ELSE 
     1039                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 
     1040                     ENDIF 
     1041                  CASE(3) 
     1042                     IF(ln_sco) THEN 
     1043                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n? 
     1044                     ELSE 
     1045                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) ) 
     1046                     ENDIF 
     1047               END SELECT 
     1048               IF( zl < dta_read_z(ji,jj,1) ) THEN                                      ! above the first level of external data 
     1049                  dta(ib,1,ik) =  dta_read(ji,jj,1) 
     1050               ELSEIF( zl > MAXVAL(dta_read_z(ji,jj,:),1) ) THEN                        ! below the last level of external data  
     1051                  dta(ib,1,ik) =  dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1)) 
     1052               ELSE                                                                     ! inbetween : vertical interpolation between ikk & ikk+1 
     1053                  DO ikk = 1, jpkm1_bdy                                                 ! when  gdept_n(ikk) < zl < gdept_n(ikk+1) 
     1054                     IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) & 
     1055                    &    .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN 
     1056                        zi           = ( zl - dta_read_z(ji,jj,ikk) ) / & 
     1057                       &               ( dta_read_z(ji,jj,ikk+1) - dta_read_z(ji,jj,ikk) ) 
     1058                        dta(ib,1,ik) = dta_read(ji,jj,ikk) + & 
     1059                       &               ( dta_read(ji,jj,ikk+1) - dta_read(ji,jj,ikk) ) * zi 
     1060                     ENDIF 
     1061                  END DO 
     1062               ENDIF 
    7231063            END DO 
    7241064         END DO 
    725       ENDIF 
    726       ! 
    727    END SUBROUTINE fld_map 
    728  
     1065 
     1066         IF(igrd == 2) THEN                                 ! do we need to adjust the transport term? 
     1067            DO ib = 1, ipi 
     1068               jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
     1069               ji=map%ptr(ib)-(jj-1)*ilendta 
     1070               zij = idx_bdy(ibdy)%nbi(ib,igrd) 
     1071               zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
     1072               zh = SUM(dta_read_dz(ji,jj,:) ) 
     1073               ztrans = 0._wp 
     1074               ztrans_new = 0._wp 
     1075               DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
     1076                  ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 
     1077               ENDDO 
     1078               DO ik = 1, ipk                                ! calculate transport on model grid 
     1079                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik) 
     1080               ENDDO 
     1081               DO ik = 1, ipk                                ! make transport correction 
     1082                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
     1083                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 
     1084                  ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
     1085                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 
     1086                  ENDIF 
     1087               ENDDO 
     1088            ENDDO 
     1089         ENDIF 
     1090 
     1091         IF(igrd == 3) THEN                                 ! do we need to adjust the transport term? 
     1092            DO ib = 1, ipi 
     1093               jj  = 1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
     1094               ji  = map%ptr(ib)-(jj-1)*ilendta 
     1095               zij = idx_bdy(ibdy)%nbi(ib,igrd) 
     1096               zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
     1097               zh  = SUM(dta_read_dz(ji,jj,:) ) 
     1098               ztrans = 0._wp 
     1099               ztrans_new = 0._wp 
     1100               DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
     1101                  ztrans     = ztrans     + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 
     1102               ENDDO 
     1103               DO ik = 1, ipk                                ! calculate transport on model grid 
     1104                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik) 
     1105               ENDDO 
     1106               DO ik = 1, ipk                                ! make transport correction 
     1107                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
     1108                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 
     1109                  ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
     1110                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 
     1111                  ENDIF 
     1112               ENDDO 
     1113            ENDDO 
     1114         ENDIF 
     1115 
     1116      ENDIF ! endif unstructured or structured 
     1117 
     1118   END SUBROUTINE fld_bdy_interp 
     1119 
     1120!  SUBROUTINE fld_bdy_conserve(dta_read, dta_read_z, map, jpk_bdy, igrd, ibdy, fv, dta) 
     1121 
     1122!  END SUBROUTINE fld_bdy_conserve 
     1123 
     1124#endif 
    7291125 
    7301126   SUBROUTINE fld_rot( kt, sd ) 
Note: See TracChangeset for help on using the changeset viewer.