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 11822 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/fldread.F90 – NEMO

Ignore:
Timestamp:
2019-10-29T11:41:36+01:00 (4 years ago)
Author:
acc
Message:

Branch 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. Sette tested updates to branch to align with trunk changes between 10721 and 11740. Sette tests are passing but results differ from branch before these changes (except for GYRE_PISCES and VORTEX) and branch results already differed from trunk because of algorithmic fixes. Will need more checks to confirm correctness.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/fldread.F90

    r11480 r11822  
    4646   PUBLIC   fld_clopn 
    4747 
    48    INTEGER :: nfld_Nnn = 1 
    4948   TYPE, PUBLIC ::   FLD_N      !: Namelist field informations 
    5049      CHARACTER(len = 256) ::   clname      ! generic name of the NetCDF flux file 
    51       REAL(wp)             ::   nfreqh      ! frequency of each flux file 
     50      REAL(wp)             ::   freqh       ! frequency of each flux file 
    5251      CHARACTER(len = 34)  ::   clvar       ! generic name of the variable in the NetCDF flux file 
    5352      LOGICAL              ::   ln_tint     ! time interpolation or not (T/F) 
     
    6564      CHARACTER(len = 256)            ::   clrootname   ! generic name of the NetCDF file 
    6665      CHARACTER(len = 256)            ::   clname       ! current name of the NetCDF file 
    67       REAL(wp)                        ::   nfreqh       ! frequency of each flux file 
     66      REAL(wp)                        ::   freqh        ! frequency of each flux file 
    6867      CHARACTER(len = 34)             ::   clvar        ! generic name of the variable in the NetCDF flux file 
    6968      LOGICAL                         ::   ln_tint      ! time interpolation or not (T/F) 
     
    8180      INTEGER                         ::   nreclast     ! last record to be read in the current file 
    8281      CHARACTER(len = 256)            ::   lsmname      ! current name of the NetCDF mask file acting as a key 
    83       INTEGER                         ::   igrd         ! grid type for bdy data 
    84       INTEGER                         ::   ibdy         ! bdy set id number 
     82      !                                                 !  
     83      !                                                 ! Variables related to BDY 
     84      INTEGER                         ::   igrd         !   grid type for bdy data 
     85      INTEGER                         ::   ibdy         !   bdy set id number 
     86      INTEGER, POINTER, DIMENSION(:)  ::   imap         !   Array of integer pointers to 1D arrays 
     87      LOGICAL                         ::   ltotvel      !   total velocity or not (T/F) 
     88      LOGICAL                         ::   lzint        !   T if it requires a vertical interpolation 
    8589   END TYPE FLD 
    86  
    87    TYPE, PUBLIC ::   MAP_POINTER      !: Map from input data file to local domain 
    88       INTEGER, POINTER, DIMENSION(:)  ::  ptr           ! Array of integer pointers to 1D arrays 
    89       LOGICAL                         ::  ll_unstruc    ! Unstructured (T) or structured (F) boundary data file 
    90    END TYPE MAP_POINTER 
    9190 
    9291!$AGRIF_DO_NOT_TREAT 
     
    130129CONTAINS 
    131130 
    132    SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset, jpk_bdy, fvl, Kmm ) 
     131   SUBROUTINE fld_read( kt, kn_fsbc, sd, kit, kt_offset, Kmm ) 
    133132      !!--------------------------------------------------------------------- 
    134133      !!                    ***  ROUTINE fld_read  *** 
     
    145144      INTEGER  , INTENT(in   )               ::   kn_fsbc   ! sbc computation period (in time step)  
    146145      TYPE(FLD), INTENT(inout), DIMENSION(:) ::   sd        ! input field related variables 
    147       TYPE(MAP_POINTER),INTENT(in), OPTIONAL, DIMENSION(:) ::   map   ! global-to-local mapping indices 
    148146      INTEGER  , INTENT(in   ), OPTIONAL     ::   kit       ! subcycle timestep for timesplitting option 
    149147      INTEGER  , INTENT(in   ), OPTIONAL     ::   kt_offset ! provide fields at time other than "now" 
     
    151149      !                                                     !   kt_offset = +1 => fields at "after"  time level 
    152150      !                                                     !   etc. 
    153       INTEGER  , INTENT(in   ), OPTIONAL     ::   jpk_bdy   ! number of vertical levels in the BDY data 
    154       LOGICAL  , INTENT(in   ), OPTIONAL     ::   fvl   ! number of vertical levels in the BDY data 
    155151      INTEGER  , INTENT(in   ), OPTIONAL     ::   Kmm   ! ocean time level index 
    156152      !! 
     
    168164      REAL(wp) ::   ztintb       ! ratio applied to before records when doing time interpolation 
    169165      CHARACTER(LEN=1000) ::   clfmt  ! write format 
    170       TYPE(MAP_POINTER)   ::   imap   ! global-to-local mapping indices 
    171166      !!--------------------------------------------------------------------- 
    172167      ll_firstcall = kt == nit000 
     
    177172      ENDIF 
    178173      IF( PRESENT(kt_offset) )   it_offset = kt_offset 
    179  
    180       imap%ptr => NULL() 
    181174 
    182175      ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar  
     
    190183      IF( ll_firstcall ) THEN                      ! initialization 
    191184         DO jf = 1, imf  
    192             IF( PRESENT(map) ) imap = map(jf) 
    193                IF( PRESENT(jpk_bdy) ) THEN 
    194                   CALL fld_init( kn_fsbc, sd(jf), imap, jpk_bdy, fvl )  ! read each before field (put them in after as they will be swapped) 
    195                ELSE 
    196                   CALL fld_init( kn_fsbc, sd(jf), imap )  ! read each before field (put them in after as they will be swapped) 
    197                ENDIF 
     185            IF( TRIM(sd(jf)%clrootname) == 'NOT USED' )   CYCLE 
     186            CALL fld_init( kn_fsbc, sd(jf) )       ! read each before field (put them in after as they will be swapped) 
    198187         END DO 
    199188         IF( lwp ) CALL wgt_print()                ! control print 
     
    204193         ! 
    205194         DO jf = 1, imf                            ! ---   loop over field   --- ! 
    206              
     195 
     196            IF( TRIM(sd(jf)%clrootname) == 'NOT USED' )   CYCLE 
     197                       
    207198            IF( isecsbc > sd(jf)%nrec_a(2) .OR. ll_firstcall ) THEN    ! read/update the after data? 
    208  
    209                IF( PRESENT(map) )   imap = map(jf)   ! temporary definition of map 
    210199 
    211200               sd(jf)%nrec_b(:) = sd(jf)%nrec_a(:)                                  ! swap before record informations 
     
    215204               CALL fld_rec( kn_fsbc, sd(jf), kt_offset = it_offset, kit = kit )    ! update after record informations 
    216205 
    217                ! if kn_fsbc*rdt is larger than nfreqh (which is kind of odd), 
     206               ! if kn_fsbc*rdt is larger than freqh (which is kind of odd), 
    218207               ! it is possible that the before value is no more the good one... we have to re-read it 
    219208               ! if before is not the last record of the file currently opened and after is the first record to be read 
     
    224213                  itmp = sd(jf)%nrec_a(1)                       ! temporary storage 
    225214                  sd(jf)%nrec_a(1) = sd(jf)%nreclast            ! read the last record of the file currently opened 
    226                   CALL fld_get( sd(jf), imap )                  ! read after data 
     215                  CALL fld_get( sd(jf) )                        ! read after data 
    227216                  sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)   ! re-swap before record field 
    228217                  sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1)           ! update before record informations 
    229                   sd(jf)%nrec_b(2) = sd(jf)%nrec_a(2) - NINT( sd(jf)%nfreqh * 3600 )  ! assume freq to be in hours in this case 
     218                  sd(jf)%nrec_b(2) = sd(jf)%nrec_a(2) - NINT( sd(jf)%freqh * 3600. )  ! assume freq to be in hours in this case 
    230219                  sd(jf)%rotn(1)   = sd(jf)%rotn(2)             ! update before rotate informations 
    231220                  sd(jf)%nrec_a(1) = itmp                       ! move back to after record  
     
    236225               IF( sd(jf)%ln_tint ) THEN 
    237226                   
    238                   ! if kn_fsbc*rdt is larger than nfreqh (which is kind of odd), 
     227                  ! if kn_fsbc*rdt is larger than freqh (which is kind of odd), 
    239228                  ! it is possible that the before value is no more the good one... we have to re-read it 
    240229                  ! if before record is not just just before the after record... 
     
    242231                     &                   .AND. sd(jf)%nrec_b(1) /= sd(jf)%nrec_a(1) - 1 ) THEN    
    243232                     sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) - 1       ! move back to before record 
    244                      CALL fld_get( sd(jf), imap )                  ! read after data 
     233                     CALL fld_get( sd(jf) )                        ! read after data 
    245234                     sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)   ! re-swap before record field 
    246235                     sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1)           ! update before record informations 
    247                      sd(jf)%nrec_b(2) = sd(jf)%nrec_a(2) - NINT( sd(jf)%nfreqh * 3600 )  ! assume freq to be in hours in this case 
     236                     sd(jf)%nrec_b(2) = sd(jf)%nrec_a(2) - NINT( sd(jf)%freqh * 3600. )  ! assume freq to be in hours in this case 
    248237                     sd(jf)%rotn(1)   = sd(jf)%rotn(2)             ! update before rotate informations 
    249238                     sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) + 1       ! move back to after record 
     
    270259                     ! year/month/week/day, next year/month/week/day file must exist 
    271260                     isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdt)   ! second at the end of the run 
    272                      llstop = isecend > sd(jf)%nrec_a(2)                                   ! read more than 1 record of next year 
     261                     llstop = isecend > sd(jf)%nrec_a(2)                             ! read more than 1 record of next year 
    273262                     ! we suppose that the date of next file is next day (should be ok even for weekly files...) 
    274263                     CALL fld_clopn( sd(jf), nyear  + COUNT((/llnxtyr /))                                           ,         & 
     
    279268                        CALL ctl_warn('next year/month/week/day file: '//TRIM(sd(jf)%clname)//     & 
    280269                           &     ' not present -> back to current year/month/day') 
    281                         CALL fld_clopn( sd(jf) )       ! back to the current year/month/day 
     270                        CALL fld_clopn( sd(jf) )               ! back to the current year/month/day 
    282271                        sd(jf)%nrec_a(1) = sd(jf)%nreclast     ! force to read the last record in the current year file 
    283272                     ENDIF 
     
    287276                   
    288277               ! read after data 
    289                IF( PRESENT(jpk_bdy) ) THEN 
    290                   CALL fld_get( sd(jf), imap, jpk_bdy, fvl, Kmm ) 
    291                ELSE 
    292                   CALL fld_get( sd(jf), imap ) 
    293                ENDIF 
     278 
     279               CALL fld_get( sd(jf), Kmm ) 
     280                
    294281            ENDIF   ! read new data? 
    295282         END DO                                    ! --- end loop over field --- ! 
     
    298285 
    299286         DO jf = 1, imf                            ! ---   loop over field   --- ! 
     287            ! 
     288            IF( TRIM(sd(jf)%clrootname) == 'NOT USED' )   CYCLE 
    300289            ! 
    301290            IF( sd(jf)%ln_tint ) THEN              ! temporal interpolation 
     
    329318 
    330319 
    331    SUBROUTINE fld_init( kn_fsbc, sdjf, map , jpk_bdy, fvl) 
     320   SUBROUTINE fld_init( kn_fsbc, sdjf ) 
    332321      !!--------------------------------------------------------------------- 
    333322      !!                    ***  ROUTINE fld_init  *** 
     
    338327      INTEGER  , INTENT(in   ) ::   kn_fsbc      ! sbc computation period (in time step)  
    339328      TYPE(FLD), INTENT(inout) ::   sdjf         ! input field related variables 
    340       TYPE(MAP_POINTER),INTENT(in) ::   map      ! global-to-local mapping indices 
    341       INTEGER  , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data 
    342       LOGICAL  , INTENT(in), OPTIONAL :: fvl     ! number of vertical levels in the BDY data 
    343329      !! 
    344330      LOGICAL :: llprevyr              ! are we reading previous year  file? 
     
    353339      CHARACTER(LEN=1000) ::   clfmt   ! write format 
    354340      !!--------------------------------------------------------------------- 
     341      ! 
    355342      llprevyr   = .FALSE. 
    356343      llprevmth  = .FALSE. 
     
    367354         ! 
    368355         IF( sdjf%nrec_a(1) == 0  ) THEN   ! we redefine record sdjf%nrec_a(1) with the last record of previous year file 
    369             IF    ( sdjf%nfreqh == -12 ) THEN   ! yearly mean 
     356            IF    ( NINT(sdjf%freqh) == -12 ) THEN   ! yearly mean 
    370357               IF( sdjf%cltype == 'yearly' ) THEN             ! yearly file 
    371358                  sdjf%nrec_a(1) = 1                                                       ! force to read the unique record 
     
    374361                  CALL ctl_stop( "fld_init: yearly mean file must be in a yearly type of file: "//TRIM(sdjf%clrootname) ) 
    375362               ENDIF 
    376             ELSEIF( sdjf%nfreqh ==  -1 ) THEN   ! monthly mean 
     363            ELSEIF( NINT(sdjf%freqh) ==  -1 ) THEN   ! monthly mean 
    377364               IF( sdjf%cltype == 'monthly' ) THEN            ! monthly file 
    378365                  sdjf%nrec_a(1) = 1                                                       ! force to read the unique record 
     
    383370                  llprevyr = .NOT. sdjf%ln_clim                                            ! use previous year  file? 
    384371               ENDIF 
    385             ELSE                                ! higher frequency mean (in hours)  
     372            ELSE                                     ! higher frequency mean (in hours)  
    386373               IF    ( sdjf%cltype      == 'monthly' ) THEN   ! monthly file 
    387                   sdjf%nrec_a(1) = NINT( 24 * nmonth_len(nmonth-1) / sdjf%nfreqh )         ! last record of previous month 
     374                  sdjf%nrec_a(1) = NINT( 24. * REAL(nmonth_len(nmonth-1),wp) / sdjf%freqh )! last record of previous month 
    388375                  llprevmth = .TRUE.                                                       ! use previous month file? 
    389376                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file? 
    390377               ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ! weekly file 
    391378                  llprevweek = .TRUE.                                                      ! use previous week  file? 
    392                   sdjf%nrec_a(1) = NINT( 24 * 7 / sdjf%nfreqh )                            ! last record of previous week 
     379                  sdjf%nrec_a(1) = NINT( 24. * 7. / sdjf%freqh )                           ! last record of previous week 
    393380                  isec_week = NINT(rday) * 7                                               ! add a shift toward previous week 
    394381               ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ! daily file 
    395                   sdjf%nrec_a(1) = NINT( 24 / sdjf%nfreqh )                                ! last record of previous day 
     382                  sdjf%nrec_a(1) = NINT( 24. / sdjf%freqh )                                ! last record of previous day 
    396383                  llprevday = .TRUE.                                                       ! use previous day   file? 
    397384                  llprevmth = llprevday .AND. nday   == 1                                  ! use previous month file? 
    398385                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file? 
    399386               ELSE                                           ! yearly file 
    400                   sdjf%nrec_a(1) = NINT( 24 * nyear_len(0) / sdjf%nfreqh )                 ! last record of previous year  
     387                  sdjf%nrec_a(1) = NINT( 24. * REAL(nyear_len(0),wp) / sdjf%freqh )        ! last record of previous year  
    401388                  llprevyr = .NOT. sdjf%ln_clim                                            ! use previous year  file? 
    402389               ENDIF 
     
    435422         ! 
    436423         ! read before data in after arrays(as we will swap it later) 
    437          IF( PRESENT(jpk_bdy) ) THEN 
    438             CALL fld_get( sdjf, map, jpk_bdy, fvl ) 
    439          ELSE 
    440             CALL fld_get( sdjf, map ) 
    441          ENDIF 
     424         CALL fld_get( sdjf ) 
    442425         ! 
    443426         clfmt = "('   fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')" 
     
    458441      !!              if sdjf%ln_tint = .FALSE. 
    459442      !!                  nrec_a(1): record number 
    460       !!                  nrec_b(2) and nrec_a(2): time of the beginning and end of the record (for print only) 
     443      !!                  nrec_b(2) and nrec_a(2): time of the beginning and end of the record 
    461444      !!---------------------------------------------------------------------- 
    462445      INTEGER  , INTENT(in   )           ::   kn_fsbc   ! sbc computation period (in time step)  
     
    486469      ELSE                                      ;   it_offset = 0 
    487470      ENDIF 
    488       IF( PRESENT(kt_offset) )   it_offset = kt_offset 
     471      IF( PRESENT(kt_offset) )      it_offset = kt_offset 
    489472      IF( PRESENT(kit) ) THEN   ;   it_offset = ( kit + it_offset ) * NINT( rdt/REAL(nn_baro,wp) ) 
    490473      ELSE                      ;   it_offset =         it_offset   * NINT(       rdt            ) 
    491474      ENDIF 
    492475      ! 
    493       !                                      ! =========== ! 
    494       IF    ( sdjf%nfreqh == -12 ) THEN      ! yearly mean 
    495          !                                   ! =========== ! 
    496          ! 
    497          IF( sdjf%ln_tint ) THEN                 ! time interpolation, shift by 1/2 record 
     476      !                                           ! =========== ! 
     477      IF    ( NINT(sdjf%freqh) == -12 ) THEN      ! yearly mean 
     478         !                                        ! =========== ! 
     479         ! 
     480         IF( sdjf%ln_tint ) THEN                  ! time interpolation, shift by 1/2 record 
    498481            ! 
    499482            !                  INT( ztmp ) 
     
    507490            !       forcing record :    1  
    508491            !                             
    509             ztmp = REAL( nsec_year, wp ) / ( REAL( nyear_len(1), wp ) * rday ) + 0.5 & 
    510            &       + REAL( it_offset, wp ) / ( REAL( nyear_len(1), wp ) * rday ) 
     492            ztmp =  REAL( nsec_year, wp ) / ( REAL( nyear_len(1), wp ) * rday ) + 0.5 & 
     493               &  + REAL( it_offset, wp ) / ( REAL( nyear_len(1), wp ) * rday ) 
    511494            sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/)) 
    512495            ! swap at the middle of the year 
     
    516499                                    & INT(ztmp) * INT(rday) * nyear_len(1) + INT(ztmp) * NINT( 0.5 * rday) * nyear_len(2)  
    517500            ENDIF 
    518          ELSE                                    ! no time interpolation 
     501         ELSE                                     ! no time interpolation 
    519502            sdjf%nrec_a(1) = 1 
    520503            sdjf%nrec_a(2) = NINT(rday) * nyear_len(1) + nsec1jan000   ! swap at the end    of the year 
     
    522505         ENDIF 
    523506         ! 
    524          !                                   ! ============ ! 
    525       ELSEIF( sdjf%nfreqh ==  -1 ) THEN      ! monthly mean ! 
    526          !                                   ! ============ ! 
    527          ! 
    528          IF( sdjf%ln_tint ) THEN                 ! time interpolation, shift by 1/2 record 
     507         !                                        ! ============ ! 
     508      ELSEIF( NINT(sdjf%freqh) ==  -1 ) THEN      ! monthly mean ! 
     509         !                                        ! ============ ! 
     510         ! 
     511         IF( sdjf%ln_tint ) THEN                  ! time interpolation, shift by 1/2 record 
    529512            ! 
    530513            !                  INT( ztmp ) 
     
    538521            !       forcing record :  nmonth  
    539522            !                             
    540             ztmp = REAL( nsec_month, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) + 0.5 & 
    541            &       + REAL( it_offset, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) 
     523            ztmp =  REAL( nsec_month, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) + 0.5 & 
     524           &      + REAL( it_offset, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) 
    542525            imth = nmonth + INT( ztmp ) - COUNT((/llbefore/)) 
    543526            IF( sdjf%cltype == 'monthly' ) THEN   ;   sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/)) 
     
    553536         ENDIF 
    554537         ! 
    555          !                                   ! ================================ ! 
    556       ELSE                                   ! higher frequency mean (in hours) 
    557          !                                   ! ================================ ! 
    558          ! 
    559          ifreq_sec = NINT( sdjf%nfreqh * 3600 )                                         ! frequency mean (in seconds) 
     538         !                                        ! ================================ ! 
     539      ELSE                                        ! higher frequency mean (in hours) 
     540         !                                        ! ================================ ! 
     541         ! 
     542         ifreq_sec = NINT( sdjf%freqh * 3600. )                                         ! frequency mean (in seconds) 
    560543         IF( sdjf%cltype(1:4) == 'week' )   isec_week = ksec_week( sdjf%cltype(6:8) )   ! since the first day of the current week 
    561544         ! number of second since the beginning of the file 
     
    567550         ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdt + REAL( it_offset, wp )        ! centrered in the middle of sbc time step 
    568551         ztmp = ztmp + 0.01 * rdt                                                       ! avoid truncation error  
    569          IF( sdjf%ln_tint ) THEN                ! time interpolation, shift by 1/2 record 
     552         IF( sdjf%ln_tint ) THEN                 ! time interpolation, shift by 1/2 record 
    570553            ! 
    571554            !          INT( ztmp/ifreq_sec + 0.5 ) 
     
    581564            !                    
    582565            ztmp= ztmp / REAL(ifreq_sec, wp) + 0.5 
    583          ELSE                                   ! no time interpolation 
     566         ELSE                                    ! no time interpolation 
    584567            ! 
    585568            !           INT( ztmp/ifreq_sec ) 
     
    612595      ENDIF 
    613596      ! 
     597      IF( .NOT. sdjf%ln_tint ) sdjf%nrec_a(2) = sdjf%nrec_a(2) - 1   ! last second belongs to bext record : *----( 
     598      ! 
    614599   END SUBROUTINE fld_rec 
    615600 
    616601 
    617    SUBROUTINE fld_get( sdjf, map, jpk_bdy, fvl, Kmm ) 
     602   SUBROUTINE fld_get( sdjf, Kmm ) 
    618603      !!--------------------------------------------------------------------- 
    619604      !!                    ***  ROUTINE fld_get  *** 
     
    622607      !!---------------------------------------------------------------------- 
    623608      TYPE(FLD)        , INTENT(inout) ::   sdjf   ! input field related variables 
    624       TYPE(MAP_POINTER), INTENT(in   ) ::   map    ! global-to-local mapping indices 
    625       INTEGER  , INTENT(in), OPTIONAL  ::   jpk_bdy ! number of vertical levels in the bdy data 
    626       LOGICAL  , INTENT(in), OPTIONAL  ::   fvl     ! number of vertical levels in the bdy data 
    627609      INTEGER  , INTENT(in), OPTIONAL  ::   Kmm     ! ocean time level index 
    628610      ! 
     
    637619      ipk = SIZE( sdjf%fnow, 3 ) 
    638620      ! 
    639       IF( ASSOCIATED(map%ptr) ) THEN 
    640          IF( PRESENT(jpk_bdy) ) THEN 
    641             IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2),                & 
    642                                                         sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl, Kmm ) 
    643             ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ),                & 
    644                                                         sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl, Kmm ) 
    645             ENDIF 
    646          ELSE 
    647             IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map ) 
    648             ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1), map ) 
    649             ENDIF 
    650          ENDIF         
     621      IF( ASSOCIATED(sdjf%imap) ) THEN 
     622         IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1),   & 
     623            &                                        sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel, sdjf%lzint ) 
     624         ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1),   & 
     625            &                                        sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel, sdjf%lzint ) 
     626         ENDIF 
    651627      ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 
    652628         CALL wgt_list( sdjf, iw ) 
     
    703679   END SUBROUTINE fld_get 
    704680 
    705    SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy, fvl, Kmm ) 
     681   SUBROUTINE fld_map( knum, cdvar, pdta, krec, kmap, kgrd, kbdy, ldtotvel, ldzint, Kmm ) 
    706682      !!--------------------------------------------------------------------- 
    707683      !!                    ***  ROUTINE fld_map  *** 
     
    710686      !!                using a general mapping (for open boundaries) 
    711687      !!---------------------------------------------------------------------- 
    712  
    713       USE bdy_oce, ONLY: ln_bdy, idx_bdy, dta_global, dta_global_z, dta_global_dz, dta_global2, dta_global2_z, dta_global2_dz                 ! workspace to read in global data arrays 
    714  
    715       INTEGER                   , INTENT(in ) ::   num     ! stream number 
    716       CHARACTER(LEN=*)          , INTENT(in ) ::   clvar   ! variable name 
    717       REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta     ! output field on model grid (2 dimensional) 
    718       INTEGER                   , INTENT(in ) ::   nrec    ! record number to read (ie time slice) 
    719       TYPE(MAP_POINTER)         , INTENT(in ) ::   map     ! global-to-local mapping indices 
    720       INTEGER  , INTENT(in), OPTIONAL         ::   igrd, ibdy, jpk_bdy  ! grid type, set number and number of vertical levels in the bdy data 
    721       LOGICAL  , INTENT(in), OPTIONAL         ::   fvl     ! grid type, set number and number of vertical levels in the bdy data 
    722       INTEGER  , INTENT(in), OPTIONAL         ::   Kmm     ! ocean time level index  
    723       INTEGER                                 ::   jpkm1_bdy! number of vertical levels in the bdy data minus 1 
    724       !! 
    725       INTEGER                                 ::   ipi      ! length of boundary data on local process 
    726       INTEGER                                 ::   ipj      ! length of dummy dimension ( = 1 ) 
    727       INTEGER                                 ::   ipk      ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
    728       INTEGER                                 ::   ilendta  ! length of data in file 
    729       INTEGER                                 ::   idvar    ! variable ID 
    730       INTEGER                                 ::   ib, ik, ji, jj   ! loop counters 
    731       INTEGER                                 ::   ierr 
    732       REAL(wp)                                ::   fv          ! fillvalue  
    733       REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read    ! work space for global data 
    734       REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read_z  ! work space for global data 
    735       REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read_dz ! work space for global data 
    736       !!--------------------------------------------------------------------- 
    737       ! 
    738       ipi = SIZE( dta, 1 ) 
    739       ipj = 1 
    740       ipk = SIZE( dta, 3 ) 
    741       ! 
    742       idvar   = iom_varid( num, clvar ) 
    743       ilendta = iom_file(num)%dimsz(1,idvar) 
    744  
    745       IF ( ln_bdy ) THEN 
    746          ipj = iom_file(num)%dimsz(2,idvar) 
    747          IF( map%ll_unstruc) THEN   ! unstructured open boundary data file 
    748             dta_read => dta_global 
    749             IF( PRESENT(jpk_bdy) ) THEN 
    750                IF( jpk_bdy>0 ) THEN 
    751                   dta_read_z => dta_global_z 
    752                   dta_read_dz => dta_global_dz 
    753                   jpkm1_bdy = jpk_bdy-1 
    754                ENDIF 
    755             ENDIF 
    756          ELSE                       ! structured open boundary file 
    757             dta_read => dta_global2 
    758             IF( PRESENT(jpk_bdy) ) THEN 
    759                IF( jpk_bdy>0 ) THEN 
    760                   dta_read_z => dta_global2_z 
    761                   dta_read_dz => dta_global2_dz 
    762                   jpkm1_bdy = jpk_bdy-1 
    763                ENDIF 
    764             ENDIF 
    765          ENDIF 
    766       ENDIF 
    767  
    768       IF(lwp) WRITE(numout,*) 'Dim size for ',        TRIM(clvar),' is ', ilendta 
    769       IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk 
    770       ! 
    771       SELECT CASE( ipk ) 
    772       CASE(1)        ;    
    773       CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1    ), nrec ) 
    774          IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 
    775             DO ib = 1, ipi 
    776               DO ik = 1, ipk 
    777                 dta(ib,1,ik) =  dta_read(map%ptr(ib),1,ik) 
    778               END DO 
    779             END DO 
    780          ELSE ! we assume that this is a structured open boundary file 
    781             DO ib = 1, ipi 
    782                jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    783                ji=map%ptr(ib)-(jj-1)*ilendta 
    784                DO ik = 1, ipk 
    785                   dta(ib,1,ik) =  dta_read(ji,jj,ik) 
    786                END DO 
    787             END DO 
    788          ENDIF 
     688      INTEGER                   , INTENT(in   ) ::   knum         ! stream number 
     689      CHARACTER(LEN=*)          , INTENT(in   ) ::   cdvar        ! variable name 
     690      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta         ! bdy output field on model grid 
     691      INTEGER                   , INTENT(in   ) ::   krec         ! record number to read (ie time slice) 
     692      INTEGER , DIMENSION(:)    , INTENT(in   ) ::   kmap         ! global-to-local bdy mapping indices 
     693      ! optional variables used for vertical interpolation: 
     694      INTEGER, OPTIONAL         , INTENT(in   ) ::   kgrd         ! grid type (t, u, v) 
     695      INTEGER, OPTIONAL         , INTENT(in   ) ::   kbdy         ! bdy number 
     696      LOGICAL, OPTIONAL         , INTENT(in   ) ::   ldtotvel     ! true if total ( = barotrop + barocline) velocity 
     697      LOGICAL, OPTIONAL         , INTENT(in   ) ::   ldzint       ! true if 3D variable requires a vertical interpolation 
     698      INTEGER, OPTIONAL         , INTENT(in   ) ::   Kmm          ! ocean time level index  
     699      !! 
     700      INTEGER                                   ::   ipi          ! length of boundary data on local process 
     701      INTEGER                                   ::   ipj          ! length of dummy dimension ( = 1 ) 
     702      INTEGER                                   ::   ipk          ! number of vertical levels of pdta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
     703      INTEGER                                   ::   ipkb         ! number of vertical levels in boundary data file 
     704      INTEGER                                   ::   idvar        ! variable ID 
     705      INTEGER                                   ::   indims       ! number of dimensions of the variable 
     706      INTEGER, DIMENSION(4)                     ::   idimsz       ! size of variable dimensions  
     707      REAL(wp)                                  ::   zfv          ! fillvalue  
     708      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zz_read      ! work space for global boundary data 
     709      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zdta_read    ! work space local data requiring vertical interpolation 
     710      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zdta_read_z  ! work space local data requiring vertical interpolation 
     711      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zdta_read_dz ! work space local data requiring vertical interpolation 
     712      CHARACTER(LEN=1),DIMENSION(3)             ::   clgrid 
     713      LOGICAL                                   ::   lluld        ! is the variable using the unlimited dimension 
     714      LOGICAL                                   ::   llzint       ! local value of ldzint 
     715      !!--------------------------------------------------------------------- 
     716      ! 
     717      clgrid = (/'t','u','v'/) 
     718      ! 
     719      ipi = SIZE( pdta, 1 ) 
     720      ipj = SIZE( pdta, 2 )   ! must be equal to 1 
     721      ipk = SIZE( pdta, 3 ) 
     722      ! 
     723      llzint = .FALSE. 
     724      IF( PRESENT(ldzint) )   llzint = ldzint 
     725      ! 
     726      idvar = iom_varid( knum, cdvar, kndims = indims, kdimsz = idimsz, lduld = lluld  ) 
     727      IF( indims == 4 .OR. ( indims == 3 .AND. .NOT. lluld ) ) THEN   ;   ipkb = idimsz(3)   ! xy(zl)t or xy(zl) 
     728      ELSE                                                            ;   ipkb = 1           ! xy or xyt 
     729      ENDIF 
     730      ! 
     731      ALLOCATE( zz_read( idimsz(1), idimsz(2), ipkb ) )  ! ++++++++ !!! this can be very big...          
     732      ! 
     733      IF( ipk == 1 ) THEN 
     734 
     735         IF( ipkb /= 1 ) CALL ctl_stop( 'fld_map : we must have ipkb = 1 to read surface data' ) 
     736         CALL iom_get ( knum, jpdom_unknown, cdvar, zz_read(:,:,1), krec )   ! call iom_get with a 2D file 
     737         CALL fld_map_core( zz_read, kmap, pdta ) 
    789738 
    790739      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    791740      ! Do we include something here to adjust barotropic velocities ! 
    792741      ! in case of a depth difference between bdy files and          ! 
    793       ! bathymetry in the case ln_full_vel = .false. and jpk_bdy>0?  ! 
     742      ! bathymetry in the case ln_totvel = .false. and ipkb>0?       ! 
    794743      ! [as the enveloping and parital cells could change H]         ! 
    795744      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    796745 
    797       CASE DEFAULT   ; 
    798  
    799       IF( PRESENT(jpk_bdy) .AND. jpk_bdy>0 ) THEN       ! boundary data not on model grid: vertical interpolation 
    800          CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar ) 
    801          dta_read(:,:,:) = -ABS(fv) 
    802          dta_read_z(:,:,:) = 0._wp 
    803          dta_read_dz(:,:,:) = 0._wp 
    804          CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:jpk_bdy), nrec ) 
    805          SELECT CASE( igrd )                   
    806             CASE(1) 
    807                CALL iom_get ( num, jpdom_unknown, 'gdept', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 
    808                CALL iom_get ( num, jpdom_unknown, 'e3t',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 
    809             CASE(2)   
    810                CALL iom_get ( num, jpdom_unknown, 'gdepu', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 
    811                CALL iom_get ( num, jpdom_unknown, 'e3u',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 
    812             CASE(3) 
    813                CALL iom_get ( num, jpdom_unknown, 'gdepv', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 
    814                CALL iom_get ( num, jpdom_unknown, 'e3v',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 
    815          END SELECT 
    816  
    817       IF ( ln_bdy ) &  
    818          CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta, Kmm) 
    819  
    820       ELSE ! boundary data assumed to be on model grid 
    821           
    822          CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec )                     
    823          IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 
    824             DO ib = 1, ipi 
    825               DO ik = 1, ipk 
    826                 dta(ib,1,ik) =  dta_read(map%ptr(ib),1,ik) 
    827               END DO 
     746      ELSE 
     747         ! 
     748         CALL iom_get ( knum, jpdom_unknown, cdvar, zz_read(:,:,:), krec )   ! call iom_get with a 3D file 
     749         ! 
     750         IF( ipkb /= ipk .OR. llzint ) THEN   ! boundary data not on model vertical grid : vertical interpolation 
     751            ! 
     752            IF( ipk == jpk .AND. iom_varid(knum,'gdep'//clgrid(kgrd)) /= -1 .AND. iom_varid(knum,'e3'//clgrid(kgrd)) /= -1 ) THEN 
     753                
     754               ALLOCATE( zdta_read(ipi,ipj,ipkb), zdta_read_z(ipi,ipj,ipkb), zdta_read_dz(ipi,ipj,ipkb) ) 
     755                 
     756               CALL fld_map_core( zz_read, kmap, zdta_read ) 
     757               CALL iom_get ( knum, jpdom_unknown, 'gdep'//clgrid(kgrd), zz_read )   ! read only once? Potential temporal evolution? 
     758               CALL fld_map_core( zz_read, kmap, zdta_read_z ) 
     759               CALL iom_get ( knum, jpdom_unknown,   'e3'//clgrid(kgrd), zz_read )   ! read only once? Potential temporal evolution? 
     760               CALL fld_map_core( zz_read, kmap, zdta_read_dz ) 
     761                
     762               CALL iom_getatt(knum, '_FillValue', zfv, cdvar=cdvar ) 
     763               CALL fld_bdy_interp(zdta_read, zdta_read_z, zdta_read_dz, pdta, kgrd, kbdy, zfv, ldtotvel) 
     764               DEALLOCATE( zdta_read, zdta_read_z, zdta_read_dz ) 
     765                
     766            ELSE 
     767               IF( ipk /= jpk ) CALL ctl_stop( 'fld_map : this should be an impossible case...' ) 
     768               WRITE(ctmp1,*) 'fld_map : vertical interpolation for bdy variable '//TRIM(cdvar)//' requires '  
     769               IF( iom_varid(knum, 'gdep'//clgrid(kgrd)) == -1 ) CALL ctl_stop( ctmp1//'gdep'//clgrid(kgrd)//' variable' ) 
     770               IF( iom_varid(knum,   'e3'//clgrid(kgrd)) == -1 ) CALL ctl_stop( ctmp1//  'e3'//clgrid(kgrd)//' variable' ) 
     771 
     772            ENDIF 
     773            ! 
     774         ELSE                            ! bdy data assumed to be the same levels as bdy variables 
     775            ! 
     776            CALL fld_map_core( zz_read, kmap, pdta ) 
     777            ! 
     778         ENDIF   ! ipkb /= ipk 
     779      ENDIF   ! ipk == 1 
     780       
     781      DEALLOCATE( zz_read ) 
     782 
     783   END SUBROUTINE fld_map 
     784 
     785      
     786   SUBROUTINE fld_map_core( pdta_read, kmap, pdta_bdy ) 
     787      !!--------------------------------------------------------------------- 
     788      !!                    ***  ROUTINE fld_map_core  *** 
     789      !! 
     790      !! ** Purpose :  inner core of fld_map 
     791      !!---------------------------------------------------------------------- 
     792      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read    ! global boundary data 
     793      INTEGER,  DIMENSION(:    ), INTENT(in   ) ::   kmap         ! global-to-local bdy mapping indices 
     794      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta_bdy     ! bdy output field on model grid 
     795      !! 
     796      INTEGER,  DIMENSION(3) ::   idim_read,  idim_bdy            ! arrays dimensions 
     797      INTEGER                ::   ji, jj, jk, jb                  ! loop counters 
     798      INTEGER                ::   im1 
     799      !!--------------------------------------------------------------------- 
     800      ! 
     801      idim_read = SHAPE( pdta_read ) 
     802      idim_bdy  = SHAPE( pdta_bdy  ) 
     803      ! 
     804      ! in all cases: idim_bdy(2) == 1 .AND. idim_read(1) * idim_read(2) == idim_bdy(1) 
     805      ! structured BDY with rimwidth > 1                     : idim_read(2) == rimwidth /= 1 
     806      ! structured BDY with rimwidth == 1 or unstructured BDY: idim_read(2) == 1 
     807      ! 
     808      IF( idim_read(2) > 1 ) THEN    ! structured BDY with rimwidth > 1   
     809         DO jk = 1, idim_bdy(3) 
     810            DO jb = 1, idim_bdy(1) 
     811               im1 = kmap(jb) - 1 
     812               jj = im1 / idim_read(1) + 1 
     813               ji = MOD( im1, idim_read(1) ) + 1 
     814               pdta_bdy(jb,1,jk) =  pdta_read(ji,jj,jk) 
    828815            END DO 
    829          ELSE ! we assume that this is a structured open boundary file 
    830             DO ib = 1, ipi 
    831                jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    832                ji=map%ptr(ib)-(jj-1)*ilendta 
    833                DO ik = 1, ipk 
    834                   dta(ib,1,ik) =  dta_read(ji,jj,ik) 
    835                END DO 
     816         END DO 
     817      ELSE 
     818         DO jk = 1, idim_bdy(3) 
     819            DO jb = 1, idim_bdy(1)   ! horizontal remap of bdy data on the local bdy  
     820               pdta_bdy(jb,1,jk) = pdta_read(kmap(jb),1,jk) 
    836821            END DO 
    837          ENDIF 
    838       ENDIF ! PRESENT(jpk_bdy) 
    839       END SELECT 
    840  
    841    END SUBROUTINE fld_map 
     822         END DO 
     823      ENDIF 
     824       
     825   END SUBROUTINE fld_map_core 
    842826    
    843    SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta, Kmm) 
    844  
     827   SUBROUTINE fld_bdy_interp(pdta_read, pdta_read_z, pdta_read_dz, pdta, kgrd, kbdy, pfv, ldtotvel, Kmm ) 
    845828      !!--------------------------------------------------------------------- 
    846829      !!                    ***  ROUTINE fld_bdy_interp  *** 
     
    851834      USE bdy_oce, ONLY:  idx_bdy         ! indexing for map <-> ij transformation 
    852835 
    853       REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read      ! work space for global data 
    854       REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_z    ! work space for global data 
    855       REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_dz   ! work space for global data 
    856       REAL(wp) , INTENT(in)                                ::   fv            ! fillvalue and alternative -ABS(fv) 
    857       REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta                        ! output field on model grid (2 dimensional) 
    858       TYPE(MAP_POINTER)         , INTENT(in ) ::   map                        ! global-to-local mapping indices 
    859       LOGICAL  , INTENT(in), OPTIONAL         ::   fvl                        ! grid type, set number and number of vertical levels in the bdy data 
    860       INTEGER  , INTENT(in)                   ::   igrd, ibdy, jpk_bdy        ! number of levels in bdy data 
    861       INTEGER  , INTENT(in)                   ::   ilendta                    ! length of data in file 
    862       INTEGER  , INTENT(in), OPTIONAL         ::   Kmm                        ! ocean time level index 
    863       !! 
    864       INTEGER                                 ::   ipi                        ! length of boundary data on local process 
    865       INTEGER                                 ::   ipj                        ! length of dummy dimension ( = 1 ) 
    866       INTEGER                                 ::   ipk                        ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
    867       INTEGER                                 ::   jpkm1_bdy                  ! number of levels in bdy data minus 1 
    868       INTEGER                                 ::   ib, ik, ikk                ! loop counters 
    869       INTEGER                                 ::   ji, jj, zij, zjj           ! temporary indices 
    870       REAL(wp)                                ::   zl, zi, zh                 ! tmp variable for current depth and interpolation factor 
    871       REAL(wp)                                ::   fv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(fv) 
    872       CHARACTER (LEN=10)                      ::   ibstr 
     836      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read       ! data read in bdy file 
     837      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdta_read_z     ! depth of the data read in bdy file 
     838      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdta_read_dz    ! thickness of the levels in bdy file 
     839      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta            ! output field on model grid (2 dimensional) 
     840      REAL(wp)                  , INTENT(in   ) ::   pfv             ! fillvalue of the data read in bdy file 
     841      LOGICAL                   , INTENT(in   ) ::   ldtotvel        ! true if toal ( = barotrop + barocline) velocity 
     842      INTEGER                   , INTENT(in   ) ::   kgrd            ! grid type (t, u, v) 
     843      INTEGER                   , INTENT(in   ) ::   kbdy            ! bdy number 
     844      INTEGER, OPTIONAL         , INTENT(in   ) ::   Kmm             ! ocean time level index 
     845      !! 
     846      INTEGER                                   ::   ipi             ! length of boundary data on local process 
     847      INTEGER                                   ::   ipkb            ! number of vertical levels in boundary data file 
     848      INTEGER                                   ::   jb, ji, jj, jk, jkb   ! loop counters 
     849      REAL(wp)                                  ::   zcoef 
     850      REAL(wp)                                  ::   zl, zi, zh      ! tmp variable for current depth and interpolation factor 
     851      REAL(wp)                                  ::   zfv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(pfv) 
     852      REAL(wp), DIMENSION(jpk)                  ::   zdepth, zdhalf  ! level and half-level depth 
    873853      !!--------------------------------------------------------------------- 
    874854      
    875  
    876       ipi       = SIZE( dta, 1 ) 
    877       ipj       = SIZE( dta_read, 2 ) 
    878       ipk       = SIZE( dta, 3 ) 
    879       jpkm1_bdy = jpk_bdy-1 
     855      ipi  = SIZE( pdta, 1 ) 
     856      ipkb = SIZE( pdta_read, 3 ) 
    880857       
    881       fv_alt = -ABS(fv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 
    882       DO ib = 1, ipi 
    883             zij = idx_bdy(ibdy)%nbi(ib,igrd) 
    884             zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    885          IF(narea==2) WRITE(*,*) 'MAPI', ib, igrd, map%ptr(ib), narea-1, zij, zjj 
    886       ENDDO 
    887       ! 
    888       IF ( map%ll_unstruc ) THEN                            ! unstructured open boundary data file 
    889  
    890          DO ib = 1, ipi 
    891             DO ik = 1, jpk_bdy 
    892                IF( ( dta_read(map%ptr(ib),1,ik) == fv ) ) THEN 
    893                   dta_read_z(map%ptr(ib),1,ik)  = fv_alt ! safety: put fillvalue into external depth field so consistent with data 
    894                   dta_read_dz(map%ptr(ib),1,ik) = 0._wp  ! safety: put 0._wp into external thickness factors to ensure transport is correct 
     858      zfv_alt = -ABS(pfv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 
     859      ! 
     860      WHERE( pdta_read == pfv ) 
     861         pdta_read_z  = zfv_alt ! safety: put fillvalue into external depth field so consistent with data 
     862         pdta_read_dz = 0._wp   ! safety: put 0._wp into external thickness factors to ensure transport is correct 
     863      ENDWHERE 
     864       
     865      DO jb = 1, ipi 
     866         ji = idx_bdy(kbdy)%nbi(jb,kgrd) 
     867         jj = idx_bdy(kbdy)%nbj(jb,kgrd) 
     868         zh  = SUM(pdta_read_dz(jb,1,:) ) 
     869         ! 
     870         ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 
     871         SELECT CASE( kgrd )                          
     872         CASE(1) 
     873            IF( ABS( (zh - ht(ji,jj)) / ht(ji,jj)) * tmask(ji,jj,1) > 0.01_wp ) THEN 
     874               WRITE(ctmp1,"(I10.10)") jb  
     875               CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%') 
     876               !   IF(lwp) WRITE(numout,*) 'DEPTHT', zh, sum(e3t(ji,jj,:,Kmm), mask=tmask(ji,jj,:)==1),  ht(ji,jj), jb, jb, ji, jj 
     877            ENDIF 
     878         CASE(2) 
     879            IF( ABS( (zh - hu(ji,jj,Kmm)) * r1_hu(ji,jj,Kmm)) * umask(ji,jj,1) > 0.01_wp ) THEN 
     880               WRITE(ctmp1,"(I10.10)") jb  
     881               CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%') 
     882               !   IF(lwp) WRITE(numout,*) 'DEPTHU', zh, SUM(e3u(ji,jj,:,Kmm), mask=umask(ji,jj,:)==1),  SUM(umask(ji,jj,:)), & 
     883               !      &                hu(ji,jj,Kmm), jb, jb, ji, jj, narea-1, pdta_read(jb,1,:) 
     884            ENDIF 
     885         CASE(3) 
     886            IF( ABS( (zh - hv(ji,jj,Kmm)) * r1_hv(ji,jj,Kmm)) * vmask(ji,jj,1) > 0.01_wp ) THEN 
     887               WRITE(ctmp1,"(I10.10)") jb 
     888               CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%') 
     889            ENDIF 
     890         END SELECT 
     891         ! 
     892         SELECT CASE( kgrd )                          
     893         CASE(1) 
     894            ! depth of T points: 
     895            zdepth(:) = gdept(ji,jj,:,Kmm) 
     896         CASE(2) 
     897            ! depth of U points: we must not use gdept_n as we don't want to do a communication 
     898            !   --> copy what is done for gdept_n in domvvl... 
     899            zdhalf(1) = 0.0_wp 
     900            zdepth(1) = 0.5_wp * e3uw(ji,jj,1,Kmm) 
     901            DO jk = 2, jpk                               ! vertical sum 
     902               !    zcoef = umask - wumask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     903               !                              ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
     904               !                              ! 0.5 where jk = mikt      
     905               !!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ?? 
     906               zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) ) 
     907               zdhalf(jk) = zdhalf(jk-1) + e3u(ji,jj,jk-1,Kmm) 
     908               zdepth(jk) =      zcoef  * ( zdhalf(jk  ) + 0.5 * e3uw(ji,jj,jk,Kmm))  & 
     909                  &         + (1-zcoef) * ( zdepth(jk-1) + e3uw(ji,jj,jk,Kmm)) 
     910            END DO 
     911         CASE(3) 
     912            ! depth of V points: we must not use gdept_n as we don't want to do a communication 
     913            !   --> copy what is done for gdept_n in domvvl... 
     914            zdhalf(1) = 0.0_wp 
     915            zdepth(1) = 0.5_wp * e3vw(ji,jj,1,Kmm) 
     916            DO jk = 2, jpk                               ! vertical sum 
     917               !    zcoef = vmask - wvmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     918               !                              ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
     919               !                              ! 0.5 where jk = mikt      
     920               !!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ?? 
     921               zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) ) 
     922               zdhalf(jk) = zdhalf(jk-1) + e3v(ji,jj,jk-1,Kmm) 
     923               zdepth(jk) =      zcoef  * ( zdhalf(jk  ) + 0.5 * e3vw(ji,jj,jk,Kmm))  & 
     924                  &         + (1-zcoef) * ( zdepth(jk-1) + e3vw(ji,jj,jk,Kmm)) 
     925            END DO 
     926         END SELECT 
     927         ! 
     928         DO jk = 1, jpk                       
     929            IF(     zdepth(jk) < pdta_read_z(jb,1,          1) ) THEN                ! above the first level of external data 
     930               pdta(jb,1,jk) =  pdta_read(jb,1,1) 
     931            ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkb) ) THEN                       ! below the last level of external data  
     932               pdta(jb,1,jk) =  pdta_read(jb,1,MAXLOC(pdta_read_z(jb,1,:),1)) 
     933            ELSE                                                             ! inbetween: vertical interpolation between jkb & jkb+1 
     934               DO jkb = 1, ipkb-1                                            ! when  gdept_n(jkb) < zdepth(jk) < gdept_n(jkb+1) 
     935                  IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) <= 0._wp ) & 
     936                     &    .AND. ( pdta_read_z(jb,1,jkb+1) /= zfv_alt) ) THEN   ! linear interpolation between 2 levels 
     937                     zi = ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) / ( pdta_read_z(jb,1,jkb+1) - pdta_read_z(jb,1,jkb) ) 
     938                     pdta(jb,1,jk) = pdta_read(jb,1,jkb) + ( pdta_read  (jb,1,jkb+1) - pdta_read  (jb,1,jkb) ) * zi 
     939                  ENDIF 
     940               END DO 
     941            ENDIF 
     942         END DO   ! jpk 
     943         ! 
     944      END DO   ! ipi 
     945       
     946      IF(kgrd == 2) THEN                                  ! do we need to adjust the transport term? 
     947         DO jb = 1, ipi 
     948            ji = idx_bdy(kbdy)%nbi(jb,kgrd) 
     949            jj = idx_bdy(kbdy)%nbj(jb,kgrd) 
     950            zh  = SUM(pdta_read_dz(jb,1,:) ) 
     951            ztrans = 0._wp 
     952            ztrans_new = 0._wp 
     953            DO jkb = 1, ipkb                              ! calculate transport on input grid 
     954               ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb) 
     955            ENDDO 
     956            DO jk = 1, jpk                                ! calculate transport on model grid 
     957               ztrans_new = ztrans_new +      pdta(jb,1,jk ) *        e3u(ji,jj,jk,Kmm ) * umask(ji,jj,jk) 
     958            ENDDO 
     959            DO jk = 1, jpk                                ! make transport correction 
     960               IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
     961                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu(ji,jj,Kmm) ) * umask(ji,jj,jk) 
     962               ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
     963                  pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hu(ji,jj,Kmm)   * umask(ji,jj,jk) 
    895964               ENDIF 
    896965            ENDDO 
    897          ENDDO  
    898  
    899          DO ib = 1, ipi 
    900             zij = idx_bdy(ibdy)%nbi(ib,igrd) 
    901             zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    902             zh  = SUM(dta_read_dz(map%ptr(ib),1,:) ) 
    903             ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 
    904             SELECT CASE( igrd )                          
    905                CASE(1) 
    906                   IF( ABS( (zh - ht(zij,zjj)) / ht(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 
    907                      WRITE(ibstr,"(I10.10)") map%ptr(ib)  
    908                      CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
    909                  !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t(zij,zjj,:,Kmm), mask=tmask(zij,zjj,:)==1),  ht(zij,zjj), map%ptr(ib), ib, zij, zjj 
    910                   ENDIF 
    911                CASE(2) 
    912                   IF( ABS( (zh - hu(zij,zjj,Kmm)) * r1_hu(zij,zjj,Kmm)) * umask(zij,zjj,1) > 0.01_wp ) THEN 
    913                      WRITE(ibstr,"(I10.10)") map%ptr(ib)  
    914                      CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
    915                      IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u(zij,zjj,:,Kmm), mask=umask(zij,zjj,:)==1),  sum(umask(zij,zjj,:)), & 
    916                        &                hu(zij,zjj,Kmm), map%ptr(ib), ib, zij, zjj, narea-1  , & 
    917                         &                dta_read(map%ptr(ib),1,:) 
    918                   ENDIF 
    919                CASE(3) 
    920                   IF( ABS( (zh - hv(zij,zjj,Kmm)) * r1_hv(zij,zjj,Kmm)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 
    921                      WRITE(ibstr,"(I10.10)") map%ptr(ib)  
    922                      CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
    923                   ENDIF 
    924             END SELECT 
    925             DO ik = 1, ipk                       
    926                SELECT CASE( igrd )                        
    927                   CASE(1) 
    928                      zl =  gdept(zij,zjj,ik,Kmm)                                          ! if using in step could use fsdept instead of gdept_n? 
    929                   CASE(2) 
    930                      IF(ln_sco) THEN 
    931                         zl =  ( gdept(zij,zjj,ik,Kmm) + gdept(zij+1,zjj,ik,Kmm) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n? 
    932                      ELSE 
    933                         zl =  MIN( gdept(zij,zjj,ik,Kmm), gdept(zij+1,zjj,ik,Kmm) )  
    934                      ENDIF 
    935                   CASE(3) 
    936                      IF(ln_sco) THEN 
    937                         zl =  ( gdept(zij,zjj,ik,Kmm) + gdept(zij,zjj+1,ik,Kmm) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n? 
    938                      ELSE 
    939                         zl =  MIN( gdept(zij,zjj,ik,Kmm), gdept(zij,zjj+1,ik,Kmm) ) 
    940                      ENDIF 
    941                END SELECT 
    942                IF( zl < dta_read_z(map%ptr(ib),1,1) ) THEN                                         ! above the first level of external data 
    943                   dta(ib,1,ik) =  dta_read(map%ptr(ib),1,1) 
    944                ELSEIF( zl > MAXVAL(dta_read_z(map%ptr(ib),1,:),1) ) THEN                           ! below the last level of external data  
    945                   dta(ib,1,ik) =  dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1)) 
    946                ELSE                                                                                ! inbetween : vertical interpolation between ikk & ikk+1 
    947                   DO ikk = 1, jpkm1_bdy                                                            ! when  gdept(ikk,Kmm) < zl < gdept(ikk+1,Kmm) 
    948                      IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) & 
    949                     &    .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN 
    950                         zi           = ( zl - dta_read_z(map%ptr(ib),1,ikk) ) / & 
    951                        &               ( dta_read_z(map%ptr(ib),1,ikk+1) - dta_read_z(map%ptr(ib),1,ikk) ) 
    952                         dta(ib,1,ik) = dta_read(map%ptr(ib),1,ikk) + & 
    953                        &               ( dta_read(map%ptr(ib),1,ikk+1) - dta_read(map%ptr(ib),1,ikk) ) * zi 
    954                      ENDIF 
    955                   END DO 
    956                ENDIF 
    957             END DO 
    958          END DO 
    959  
    960          IF(igrd == 2) THEN                                 ! do we need to adjust the transport term? 
    961             DO ib = 1, ipi 
    962               zij = idx_bdy(ibdy)%nbi(ib,igrd) 
    963               zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    964               zh  = SUM(dta_read_dz(map%ptr(ib),1,:) ) 
    965               ztrans = 0._wp 
    966               ztrans_new = 0._wp 
    967               DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
    968                   ztrans     = ztrans     + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 
    969               ENDDO 
    970               DO ik = 1, ipk                                ! calculate transport on model grid 
    971                   ztrans_new = ztrans_new + dta(ib,1,ik) * e3u(zij,zjj,ik,Kmm) * umask(zij,zjj,ik) 
    972               ENDDO 
    973               DO ik = 1, ipk                                ! make transport correction 
    974                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    975                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu(zij,zjj,Kmm) ) * umask(zij,zjj,ik) 
    976                  ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
    977                     IF( ABS(ztrans * r1_hu(zij,zjj,Kmm)) > 0.01_wp ) & 
    978                    &   CALL ctl_warn('fld_bdy_interp: barotropic component of > 0.01 ms-1 found in baroclinic velocities at') 
    979                     dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu(zij,zjj,Kmm) * umask(zij,zjj,ik) 
    980                  ENDIF 
    981               ENDDO 
     966         ENDDO 
     967      ENDIF 
     968       
     969      IF(kgrd == 3) THEN                                  ! do we need to adjust the transport term? 
     970         DO jb = 1, ipi 
     971            ji = idx_bdy(kbdy)%nbi(jb,kgrd) 
     972            jj = idx_bdy(kbdy)%nbj(jb,kgrd) 
     973            zh  = SUM(pdta_read_dz(jb,1,:) ) 
     974            ztrans = 0._wp 
     975            ztrans_new = 0._wp 
     976            DO jkb = 1, ipkb                              ! calculate transport on input grid 
     977               ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb) 
    982978            ENDDO 
    983          ENDIF 
    984  
    985          IF(igrd == 3) THEN                                 ! do we need to adjust the transport term? 
    986             DO ib = 1, ipi 
    987               zij = idx_bdy(ibdy)%nbi(ib,igrd) 
    988               zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    989               zh  = SUM(dta_read_dz(map%ptr(ib),1,:) ) 
    990               ztrans = 0._wp 
    991               ztrans_new = 0._wp 
    992               DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
    993                   ztrans     = ztrans     + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 
    994               ENDDO 
    995               DO ik = 1, ipk                                ! calculate transport on model grid 
    996                   ztrans_new = ztrans_new + dta(ib,1,ik) * e3v(zij,zjj,ik,Kmm) * vmask(zij,zjj,ik) 
    997               ENDDO 
    998               DO ik = 1, ipk                                ! make transport correction 
    999                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    1000                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv(zij,zjj,Kmm) ) * vmask(zij,zjj,ik) 
    1001                  ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
    1002                     dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv(zij,zjj,Kmm) * vmask(zij,zjj,ik) 
    1003                  ENDIF 
    1004               ENDDO 
     979            DO jk = 1, jpk                                ! calculate transport on model grid 
     980               ztrans_new = ztrans_new +      pdta(jb,1,jk ) *        e3v(ji,jj,jk,Kmm ) * vmask(ji,jj,jk) 
    1005981            ENDDO 
    1006          ENDIF 
    1007    
    1008       ELSE ! structured open boundary file 
    1009  
    1010          DO ib = 1, ipi 
    1011             jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    1012             ji=map%ptr(ib)-(jj-1)*ilendta 
    1013             DO ik = 1, jpk_bdy                       
    1014                IF( ( dta_read(ji,jj,ik) == fv ) ) THEN 
    1015                   dta_read_z(ji,jj,ik)  = fv_alt ! safety: put fillvalue into external depth field so consistent with data 
    1016                   dta_read_dz(ji,jj,ik) = 0._wp  ! safety: put 0._wp into external thickness factors to ensure transport is correct 
     982            DO jk = 1, jpk                                ! make transport correction 
     983               IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
     984                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv(ji,jj,Kmm) ) * vmask(ji,jj,jk) 
     985               ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
     986                  pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hv(ji,jj,Kmm)   * vmask(ji,jj,jk) 
    1017987               ENDIF 
    1018988            ENDDO 
    1019          ENDDO  
    1020         
    1021  
    1022          DO ib = 1, ipi 
    1023             jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    1024             ji=map%ptr(ib)-(jj-1)*ilendta 
    1025             zij = idx_bdy(ibdy)%nbi(ib,igrd) 
    1026             zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    1027             zh  = SUM(dta_read_dz(ji,jj,:) ) 
    1028             ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 
    1029             SELECT CASE( igrd )                          
    1030                CASE(1) 
    1031                   IF( ABS( (zh - ht(zij,zjj)) / ht(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 
    1032                      WRITE(ibstr,"(I10.10)") map%ptr(ib)  
    1033                      CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
    1034                  !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t(zij,zjj,:,Kmm), mask=tmask(zij,zjj,:)==1),  ht(zij,zjj), map%ptr(ib), ib, zij, zjj 
    1035                   ENDIF 
    1036                CASE(2) 
    1037                   IF( ABS( (zh - hu(zij,zjj,Kmm)) * r1_hu(zij,zjj,Kmm)) * umask(zij,zjj,1) > 0.01_wp ) THEN 
    1038                      WRITE(ibstr,"(I10.10)") map%ptr(ib)  
    1039                      CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
    1040                   ENDIF 
    1041                CASE(3) 
    1042                   IF( ABS( (zh - hv(zij,zjj,Kmm)) * r1_hv(zij,zjj,Kmm)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 
    1043                      WRITE(ibstr,"(I10.10)") map%ptr(ib)  
    1044                      CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
    1045                   ENDIF 
    1046             END SELECT 
    1047             DO ik = 1, ipk                       
    1048                SELECT CASE( igrd )                          ! coded for sco - need zco and zps option using min 
    1049                   CASE(1) 
    1050                      zl =  gdept(zij,zjj,ik,Kmm)              ! if using in step could use fsdept instead of gdept_n? 
    1051                   CASE(2) 
    1052                      IF(ln_sco) THEN 
    1053                         zl =  ( gdept(zij,zjj,ik,Kmm) + gdept(zij+1,zjj,ik,Kmm) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n? 
    1054                      ELSE 
    1055                         zl =  MIN( gdept(zij,zjj,ik,Kmm), gdept(zij+1,zjj,ik,Kmm) ) 
    1056                      ENDIF 
    1057                   CASE(3) 
    1058                      IF(ln_sco) THEN 
    1059                         zl =  ( gdept(zij,zjj,ik,Kmm) + gdept(zij,zjj+1,ik,Kmm) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n? 
    1060                      ELSE 
    1061                         zl =  MIN( gdept(zij,zjj,ik,Kmm), gdept(zij,zjj+1,ik,Kmm) ) 
    1062                      ENDIF 
    1063                END SELECT 
    1064                IF( zl < dta_read_z(ji,jj,1) ) THEN                                      ! above the first level of external data 
    1065                   dta(ib,1,ik) =  dta_read(ji,jj,1) 
    1066                ELSEIF( zl > MAXVAL(dta_read_z(ji,jj,:),1) ) THEN                        ! below the last level of external data  
    1067                   dta(ib,1,ik) =  dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1)) 
    1068                ELSE                                                                     ! inbetween : vertical interpolation between ikk & ikk+1 
    1069                   DO ikk = 1, jpkm1_bdy                                                 ! when  gdept(ikk,Kmm) < zl < gdept(ikk+1,Kmm) 
    1070                      IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) & 
    1071                     &    .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN 
    1072                         zi           = ( zl - dta_read_z(ji,jj,ikk) ) / & 
    1073                        &               ( dta_read_z(ji,jj,ikk+1) - dta_read_z(ji,jj,ikk) ) 
    1074                         dta(ib,1,ik) = dta_read(ji,jj,ikk) + & 
    1075                        &               ( dta_read(ji,jj,ikk+1) - dta_read(ji,jj,ikk) ) * zi 
    1076                      ENDIF 
    1077                   END DO 
    1078                ENDIF 
    1079             END DO 
    1080          END DO 
    1081  
    1082          IF(igrd == 2) THEN                                 ! do we need to adjust the transport term? 
    1083             DO ib = 1, ipi 
    1084                jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    1085                ji=map%ptr(ib)-(jj-1)*ilendta 
    1086                zij = idx_bdy(ibdy)%nbi(ib,igrd) 
    1087                zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    1088                zh = SUM(dta_read_dz(ji,jj,:) ) 
    1089                ztrans = 0._wp 
    1090                ztrans_new = 0._wp 
    1091                DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
    1092                   ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 
    1093                ENDDO 
    1094                DO ik = 1, ipk                                ! calculate transport on model grid 
    1095                   ztrans_new = ztrans_new + dta(ib,1,ik) * e3u(zij,zjj,ik,Kmm) * umask(zij,zjj,ik) 
    1096                ENDDO 
    1097                DO ik = 1, ipk                                ! make transport correction 
    1098                   IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    1099                      dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu(zij,zjj,Kmm) ) * umask(zij,zjj,ik) 
    1100                   ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
    1101                      dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * r1_hu(zij,zjj,Kmm) ) * umask(zij,zjj,ik) 
    1102                   ENDIF 
    1103                ENDDO 
    1104             ENDDO 
    1105          ENDIF 
    1106  
    1107          IF(igrd == 3) THEN                                 ! do we need to adjust the transport term? 
    1108             DO ib = 1, ipi 
    1109                jj  = 1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    1110                ji  = map%ptr(ib)-(jj-1)*ilendta 
    1111                zij = idx_bdy(ibdy)%nbi(ib,igrd) 
    1112                zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    1113                zh  = SUM(dta_read_dz(ji,jj,:) ) 
    1114                ztrans = 0._wp 
    1115                ztrans_new = 0._wp 
    1116                DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
    1117                   ztrans     = ztrans     + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 
    1118                ENDDO 
    1119                DO ik = 1, ipk                                ! calculate transport on model grid 
    1120                   ztrans_new = ztrans_new + dta(ib,1,ik) * e3v(zij,zjj,ik,Kmm) * vmask(zij,zjj,ik) 
    1121                ENDDO 
    1122                DO ik = 1, ipk                                ! make transport correction 
    1123                   IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    1124                      dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv(zij,zjj,Kmm) ) * vmask(zij,zjj,ik) 
    1125                   ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
    1126                      dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * r1_hv(zij,zjj,Kmm) ) * vmask(zij,zjj,ik) 
    1127                   ENDIF 
    1128                ENDDO 
    1129             ENDDO 
    1130          ENDIF 
    1131  
    1132       ENDIF ! endif unstructured or structured 
    1133  
     989         ENDDO 
     990      ENDIF 
     991       
    1134992   END SUBROUTINE fld_bdy_interp 
    1135993 
     
    11561014      imf = SIZE( sd ) 
    11571015      DO ju = 1, imf 
     1016         IF( TRIM(sd(ju)%clrootname) == 'NOT USED' )   CYCLE 
    11581017         ill = LEN_TRIM( sd(ju)%vcomp ) 
    11591018         DO jn = 2-COUNT((/sd(ju)%ln_tint/)), 2 
     
    11641023                  iv = -1 
    11651024                  DO jv = 1, imf 
     1025                     IF( TRIM(sd(jv)%clrootname) == 'NOT USED' )   CYCLE 
    11661026                     IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) )   iv = jv 
    11671027                  END DO 
     
    12021062      LOGICAL, OPTIONAL, INTENT(in   ) ::   ldstop   ! stop if open to read a non-existing file (default = .TRUE.) 
    12031063      ! 
    1204       LOGICAL :: llprevyr              ! are we reading previous year  file? 
    1205       LOGICAL :: llprevmth             ! are we reading previous month file? 
    1206       INTEGER :: iyear, imonth, iday   ! first day of the current file in yyyy mm dd 
    1207       INTEGER :: isec_week             ! number of seconds since start of the weekly file 
    1208       INTEGER :: indexyr               ! year undex (O/1/2: previous/current/next) 
    1209       INTEGER :: iyear_len, imonth_len ! length (days) of iyear and imonth             !  
    1210       CHARACTER(len = 256)::   clname  ! temporary file name 
     1064      LOGICAL  :: llprevyr              ! are we reading previous year  file? 
     1065      LOGICAL  :: llprevmth             ! are we reading previous month file? 
     1066      INTEGER  :: iyear, imonth, iday   ! first day of the current file in yyyy mm dd 
     1067      INTEGER  :: isec_week             ! number of seconds since start of the weekly file 
     1068      INTEGER  :: indexyr               ! year undex (O/1/2: previous/current/next) 
     1069      REAL(wp) :: zyear_len, zmonth_len ! length (days) of iyear and imonth             !  
     1070      CHARACTER(len = 256) ::   clname  ! temporary file name 
    12111071      !!---------------------------------------------------------------------- 
    12121072      IF( PRESENT(kyear) ) THEN                             ! use given values  
     
    12591119         ! find the last record to be read -> update sdjf%nreclast 
    12601120         indexyr = iyear - nyear + 1 
    1261          iyear_len = nyear_len( indexyr ) 
     1121         zyear_len = REAL(nyear_len( indexyr ), wp) 
    12621122         SELECT CASE ( indexyr ) 
    1263          CASE ( 0 )   ;   imonth_len = 31   ! previous year -> imonth = 12 
    1264          CASE ( 1 )   ;   imonth_len = nmonth_len(imonth)  
    1265          CASE ( 2 )   ;   imonth_len = 31   ! next     year -> imonth = 1 
     1123         CASE ( 0 )   ;   zmonth_len = 31.   ! previous year -> imonth = 12 
     1124         CASE ( 1 )   ;   zmonth_len = REAL(nmonth_len(imonth), wp) 
     1125         CASE ( 2 )   ;   zmonth_len = 31.   ! next     year -> imonth = 1 
    12661126         END SELECT 
    12671127         ! 
    12681128         ! last record to be read in the current file 
    1269          IF    ( sdjf%nfreqh == -12 ) THEN                 ;   sdjf%nreclast = 1    !  yearly mean 
    1270          ELSEIF( sdjf%nfreqh ==  -1 ) THEN                                          ! monthly mean 
     1129         IF    ( sdjf%freqh == -12. ) THEN                 ;   sdjf%nreclast = 1    !  yearly mean 
     1130         ELSEIF( sdjf%freqh ==  -1. ) THEN                                          ! monthly mean 
    12711131            IF(     sdjf%cltype      == 'monthly' ) THEN   ;   sdjf%nreclast = 1 
    12721132            ELSE                                           ;   sdjf%nreclast = 12 
    12731133            ENDIF 
    12741134         ELSE                                                                       ! higher frequency mean (in hours) 
    1275             IF(     sdjf%cltype      == 'monthly' ) THEN   ;   sdjf%nreclast = NINT( 24 * imonth_len / sdjf%nfreqh ) 
    1276             ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   sdjf%nreclast = NINT( 24 * 7          / sdjf%nfreqh ) 
    1277             ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   sdjf%nreclast = NINT( 24              / sdjf%nfreqh ) 
    1278             ELSE                                           ;   sdjf%nreclast = NINT( 24 * iyear_len  / sdjf%nfreqh ) 
     1135            IF(     sdjf%cltype      == 'monthly' ) THEN   ;   sdjf%nreclast = NINT( 24. * zmonth_len / sdjf%freqh ) 
     1136            ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   sdjf%nreclast = NINT( 24. * 7.         / sdjf%freqh ) 
     1137            ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   sdjf%nreclast = NINT( 24.              / sdjf%freqh ) 
     1138            ELSE                                           ;   sdjf%nreclast = NINT( 24. * zyear_len  / sdjf%freqh ) 
    12791139            ENDIF 
    12801140         ENDIF 
     
    13041164      ! 
    13051165      DO jf = 1, SIZE(sdf) 
    1306          sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname ) 
     1166         sdf(jf)%clrootname = sdf_n(jf)%clname 
     1167         IF( TRIM(sdf_n(jf)%clname) /= 'NOT USED' )   sdf(jf)%clrootname = TRIM( cdir )//sdf(jf)%clrootname 
    13071168         sdf(jf)%clname     = "not yet defined" 
    1308          sdf(jf)%nfreqh     = sdf_n(jf)%nfreqh 
     1169         sdf(jf)%freqh      = sdf_n(jf)%freqh 
    13091170         sdf(jf)%clvar      = sdf_n(jf)%clvar 
    13101171         sdf(jf)%ln_tint    = sdf_n(jf)%ln_tint 
     
    13131174         sdf(jf)%num        = -1 
    13141175         sdf(jf)%wgtname    = " " 
    1315          IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname ) 
     1176         IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//sdf_n(jf)%wname 
    13161177         sdf(jf)%lsmname = " " 
    1317          IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 )   sdf(jf)%lsmname = TRIM( cdir )//TRIM( sdf_n(jf)%lname ) 
     1178         IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 )   sdf(jf)%lsmname = TRIM( cdir )//sdf_n(jf)%lname 
    13181179         sdf(jf)%vcomp      = sdf_n(jf)%vcomp 
    13191180         sdf(jf)%rotn(:)    = .TRUE.   ! pretend to be rotated -> won't try to rotate data before the first call to fld_get 
     
    13221183         IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim )   & 
    13231184            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.') 
    1324          sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn 
     1185         sdf(jf)%nreclast   = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn 
     1186         sdf(jf)%igrd       = 0 
     1187         sdf(jf)%ibdy       = 0 
     1188         sdf(jf)%imap       => NULL() 
     1189         sdf(jf)%ltotvel    = .FALSE. 
     1190         sdf(jf)%lzint      = .FALSE. 
    13251191      END DO 
    13261192      ! 
     
    13361202         DO jf = 1, SIZE(sdf) 
    13371203            WRITE(numout,*) '      root filename: '  , TRIM( sdf(jf)%clrootname ), '   variable name: ', TRIM( sdf(jf)%clvar ) 
    1338             WRITE(numout,*) '         frequency: '      ,       sdf(jf)%nfreqh      ,   & 
     1204            WRITE(numout,*) '         frequency: '      ,       sdf(jf)%freqh       ,   & 
    13391205               &                  '   time interp: '    ,       sdf(jf)%ln_tint     ,   & 
    13401206               &                  '   climatology: '    ,       sdf(jf)%ln_clim 
Note: See TracChangeset for help on using the changeset viewer.