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 11223 – NEMO

Changeset 11223


Ignore:
Timestamp:
2019-07-05T20:53:14+02:00 (5 years ago)
Author:
smasson
Message:

dev_r10984_HPC-13 : cleaning of rewriting of bdydta, see #2285

Location:
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SHARED/namelist_ref

    r10808 r11223  
    609609   ln_vol        = .false.    !  total volume correction (see nn_volctl parameter) 
    610610   nn_volctl     =  1         !  = 0, the total water flux across open boundaries is zero 
    611    nb_jpk_bdy    = -1         ! number of levels in the bdy data (set < 0 if consistent with planned run) 
    612611/ 
    613612!----------------------------------------------------------------------- 

  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SPITZ12/EXPREF/namelist_cfg

    r10911 r11223  
    177177   nn_rimwidth   = 1          !  width of the relaxation zone 
    178178   ln_vol        = .false.    !  total volume correction (see nn_volctl parameter) 
    179    nb_jpk_bdy    = -1         ! number of levels in the bdy data (set < 0 if consistent with planned run) 
    180179/ 
    181180!----------------------------------------------------------------------- 

  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icevar.F90

    r10994 r11223  
    7373   PUBLIC   ice_var_zapneg 
    7474   PUBLIC   ice_var_roundoff 
    75    PUBLIC   ice_var_itd 
    76    PUBLIC   ice_var_itd2 
    7775   PUBLIC   ice_var_bv            
    7876   PUBLIC   ice_var_enthalpy            
    7977   PUBLIC   ice_var_sshdyn 
     78   PUBLIC   ice_var_itd 
     79 
     80   INTERFACE ice_var_itd 
     81      MODULE PROCEDURE ice_var_itd_1c1c, ice_var_itd_Nc1c, ice_var_itd_1cMc, ice_var_itd_NcMc 
     82   END INTERFACE 
    8083 
    8184   !!---------------------------------------------------------------------- 
    … …  
    656659   END SUBROUTINE ice_var_roundoff 
    657660    
     661 
     662   SUBROUTINE ice_var_bv 
     663      !!------------------------------------------------------------------- 
     664      !!                ***  ROUTINE ice_var_bv *** 
     665      !! 
     666      !! ** Purpose :   computes mean brine volume (%) in sea ice 
     667      !! 
     668      !! ** Method  : e = - 0.054 * S (ppt) / T (C) 
     669      !! 
     670      !! References : Vancoppenolle et al., JGR, 2007 
     671      !!------------------------------------------------------------------- 
     672      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
     673      !!------------------------------------------------------------------- 
     674      ! 
     675!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done 
     676!!   instead of setting everything to zero as just below 
     677      bv_i (:,:,:) = 0._wp 
     678      DO jl = 1, jpl 
     679         DO jk = 1, nlay_i 
     680            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )    
     681               bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 ) 
     682            END WHERE 
     683         END DO 
     684      END DO 
     685      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:) 
     686      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp 
     687      END WHERE 
     688      ! 
     689   END SUBROUTINE ice_var_bv 
     690 
     691 
     692   SUBROUTINE ice_var_enthalpy 
     693      !!------------------------------------------------------------------- 
     694      !!                   ***  ROUTINE ice_var_enthalpy ***  
     695      !!                  
     696      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature 
     697      !! 
     698      !! ** Method  :   Formula (Bitz and Lipscomb, 1999) 
     699      !!------------------------------------------------------------------- 
     700      INTEGER  ::   ji, jk   ! dummy loop indices 
     701      REAL(wp) ::   ztmelts  ! local scalar  
     702      !!------------------------------------------------------------------- 
     703      ! 
     704      DO jk = 1, nlay_i             ! Sea ice energy of melting 
     705         DO ji = 1, npti 
     706            ztmelts      = - rTmlt  * sz_i_1d(ji,jk) 
     707            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point => likely conservation issue 
     708                                                                !   (sometimes zdf scheme produces abnormally high temperatures)    
     709            e_i_1d(ji,jk) = rhoi * ( rcpi  * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           & 
     710               &                   + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   & 
     711               &                   - rcp   * ztmelts ) 
     712         END DO 
     713      END DO 
     714      DO jk = 1, nlay_s             ! Snow energy of melting 
     715         DO ji = 1, npti 
     716            e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) 
     717         END DO 
     718      END DO 
     719      ! 
     720   END SUBROUTINE ice_var_enthalpy 
     721 
    658722    
    659    SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i ) 
    660       !!------------------------------------------------------------------- 
    661       !!                ***  ROUTINE ice_var_itd   *** 
    662       !! 
    663       !! ** Purpose :  converting 1-cat ice to multiple ice categories 
     723   FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b) 
     724      !!--------------------------------------------------------------------- 
     725      !!                   ***  ROUTINE ice_var_sshdyn  *** 
     726      !!                      
     727      !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded 
     728      !! 
     729      !! ** Method  :  ssh_lead = ssh + (Mice + Msnow) / rau0 
     730      !! 
     731      !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira, 
     732      !!                Sea ice-ocean coupling using a rescaled vertical coordinate z*,  
     733      !!                Ocean Modelling, Volume 24, Issues 1-2, 2008 
     734      !!---------------------------------------------------------------------- 
     735      ! 
     736      ! input 
     737      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh            !: ssh [m] 
     738      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass    !: mass of snow and ice at current  ice time step [Kg/m2] 
     739      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b  !: mass of snow and ice at previous ice time step [Kg/m2] 
     740      ! 
     741      ! output 
     742      REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn  ! equivalent ssh in lead [m] 
     743      ! 
     744      ! temporary 
     745      REAL(wp) :: zintn, zintb                     ! time interpolation weights [] 
     746      REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload  ! snow and ice load [m] 
     747      ! 
     748      ! compute ice load used to define the equivalent ssh in lead 
     749      IF( ln_ice_embd ) THEN 
     750         !                                             
     751         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 
     752         !                                               = (1/nn_fsbc)^2 * {SUM[n]        , n=0,nn_fsbc-1} 
     753         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp 
     754         ! 
     755         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   *    {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 
     756         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 
     757         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 
     758         ! 
     759         zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0 
     760         ! 
     761      ELSE 
     762         zsnwiceload(:,:) = 0.0_wp 
     763      ENDIF 
     764      ! compute equivalent ssh in lead 
     765      ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:) 
     766      ! 
     767   END FUNCTION ice_var_sshdyn 
     768 
     769    
     770   !!------------------------------------------------------------------- 
     771   !!                ***  INTERFACE ice_var_itd   *** 
     772   !! 
     773   !! ** Purpose :  converting N-cat ice to jpl ice categories 
     774   !!------------------------------------------------------------------- 
     775   SUBROUTINE ice_var_itd_1c1c( zhti, zhts, zati, zh_i, zh_s, za_i ) 
     776      !!------------------------------------------------------------------- 
     777      !! ** Purpose :  converting 1-cat ice to 1 ice category 
     778      !!------------------------------------------------------------------- 
     779      REAL(wp), DIMENSION(:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables 
     780      REAL(wp), DIMENSION(:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables 
     781      !!------------------------------------------------------------------- 
     782      zh_i(:) = zhti(:) 
     783      zh_s(:) = zhts(:) 
     784      za_i(:) = zati(:) 
     785   END SUBROUTINE ice_var_itd_1c1c 
     786 
     787   SUBROUTINE ice_var_itd_Nc1c( zhti, zhts, zati, zh_i, zh_s, za_i ) 
     788      !!------------------------------------------------------------------- 
     789      !! ** Purpose :  converting N-cat ice to 1 ice category 
     790      !!------------------------------------------------------------------- 
     791      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables 
     792      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables 
     793      !!------------------------------------------------------------------- 
     794      ! 
     795      za_i(:) = SUM( zati(:,:), dim=2 ) 
     796      ! 
     797      WHERE( za_i(:) /= 0._wp ) 
     798         zh_i(:) = SUM( zhti(:,:) * zati(:,:), dim=2 ) / za_i(:) 
     799         zh_s(:) = SUM( zhts(:,:) * zati(:,:), dim=2 ) / za_i(:) 
     800      ELSEWHERE 
     801         zh_i(:) = 0._wp 
     802         zh_s(:) = 0._wp 
     803      END WHERE 
     804      ! 
     805   END SUBROUTINE ice_var_itd_Nc1c 
     806    
     807   SUBROUTINE ice_var_itd_1cMc( zhti, zhts, zati, zh_i, zh_s, za_i ) 
     808      !!------------------------------------------------------------------- 
     809      !! 
     810      !! ** Purpose :  converting 1-cat ice to jpl ice categories 
    664811      !! 
    665812      !!                  ice thickness distribution follows a gaussian law 
    … …  
    705852      zh_s(1:idim,1:jpl) = 0._wp 
    706853      za_i(1:idim,1:jpl) = 0._wp 
     854      ! 
     855      IF( jpl == 1 ) THEN 
     856         CALL ice_var_itd_1c1c( zhti, zhts, zati, zh_i(:,1), zh_s(:,1), za_i(:,1) ) 
     857         RETURN 
     858      ENDIF 
    707859      ! 
    708860      DO ji = 1, idim 
    … …  
    801953      END DO 
    802954      ! 
    803    END SUBROUTINE ice_var_itd 
    804  
    805  
    806    SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i ) 
    807       !!------------------------------------------------------------------- 
    808       !!                ***  ROUTINE ice_var_itd2   *** 
     955   END SUBROUTINE ice_var_itd_1cMc 
     956 
     957   SUBROUTINE ice_var_itd_NcMc( zhti, zhts, zati, zh_i, zh_s, za_i ) 
     958      !!------------------------------------------------------------------- 
    809959      !! 
    810960      !! ** Purpose :  converting N-cat ice to jpl ice categories 
    … …  
    845995      idim = SIZE( zhti, 1 ) 
    846996      icat = SIZE( zhti, 2 ) 
    847       ! 
    848       ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays 
    849       ALLOCATE( jlmin(idim), jlmax(idim) ) 
    850  
    851       ! --- initialize output fields to 0 --- ! 
    852       zh_i(1:idim,1:jpl) = 0._wp 
    853       zh_s(1:idim,1:jpl) = 0._wp 
    854       za_i(1:idim,1:jpl) = 0._wp 
    855       ! 
    856       ! --- fill the categories --- ! 
    857       !     find where cat-input = cat-output and fill cat-output fields   
    858       jlmax(:) = 0 
    859       jlmin(:) = 999 
    860       jlfil(:,:) = 0 
    861       DO jl1 = 1, jpl 
    862          DO jl2 = 1, icat 
     997      !                                 ! ---------------------- ! 
     998      IF( icat == jpl ) THEN            ! input cat = output cat ! 
     999         !                              ! ---------------------- ! 
     1000         zh_i(:,:) = zhti(:,:) 
     1001         zh_s(:,:) = zhts(:,:) 
     1002         za_i(:,:) = zati(:,:) 
     1003         !                              ! ---------------------- ! 
     1004      ELSEIF( icat == 1 ) THEN          ! specific case if N = 1 ! 
     1005         !                              ! ---------------------- ! 
     1006         ! 
     1007         CALL ice_var_itd_1cMc( zhti(:,1), zhts(:,1), zati(:,1), zh_i, zh_s, za_i ) 
     1008         ! 
     1009         !                              ! ---------------------- ! 
     1010      ELSEIF( jpl == 1 ) THEN           ! specific case if M = 1 ! 
     1011         !                              ! ---------------------- ! 
     1012         ! 
     1013         CALL ice_var_itd_Nc1c( zhti, zhts, zati, zh_i(:,1), zh_s(:,1), za_i(:,1) ) 
     1014         ! 
     1015         !                              ! ----------------------- ! 
     1016      ELSE                              ! input cat /= output cat ! 
     1017         !                              ! ----------------------- ! 
     1018          
     1019         ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays 
     1020         ALLOCATE( jlmin(idim), jlmax(idim) ) 
     1021 
     1022         ! --- initialize output fields to 0 --- ! 
     1023         zh_i(1:idim,1:jpl) = 0._wp 
     1024         zh_s(1:idim,1:jpl) = 0._wp 
     1025         za_i(1:idim,1:jpl) = 0._wp 
     1026         ! 
     1027         ! --- fill the categories --- ! 
     1028         !     find where cat-input = cat-output and fill cat-output fields   
     1029         jlmax(:) = 0 
     1030         jlmin(:) = 999 
     1031         jlfil(:,:) = 0 
     1032         DO jl1 = 1, jpl 
     1033            DO jl2 = 1, icat 
     1034               DO ji = 1, idim 
     1035                  IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN 
     1036                     ! fill the right category 
     1037                     zh_i(ji,jl1) = zhti(ji,jl2) 
     1038                     zh_s(ji,jl1) = zhts(ji,jl2) 
     1039                     za_i(ji,jl1) = zati(ji,jl2) 
     1040                     ! record categories that are filled 
     1041                     jlmax(ji) = MAX( jlmax(ji), jl1 ) 
     1042                     jlmin(ji) = MIN( jlmin(ji), jl1 ) 
     1043                     jlfil(ji,jl1) = jl1 
     1044                  ENDIF 
     1045               END DO 
     1046            END DO 
     1047         END DO 
     1048         ! 
     1049         ! --- fill the gaps between categories --- !   
     1050         !     transfer from categories filled at the previous step to the empty ones in between 
     1051         DO ji = 1, idim 
     1052            jl1 = jlmin(ji) 
     1053            jl2 = jlmax(ji) 
     1054            IF( jl1 > 1 ) THEN 
     1055               ! fill the lower cat (jl1-1) 
     1056               za_i(ji,jl1-1) = ztrans * za_i(ji,jl1) 
     1057               zh_i(ji,jl1-1) = hi_mean(jl1-1) 
     1058               ! remove from cat jl1 
     1059               za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1) 
     1060            ENDIF 
     1061            IF( jl2 < jpl ) THEN 
     1062               ! fill the upper cat (jl2+1) 
     1063               za_i(ji,jl2+1) = ztrans * za_i(ji,jl2) 
     1064               zh_i(ji,jl2+1) = hi_mean(jl2+1) 
     1065               ! remove from cat jl2 
     1066               za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2) 
     1067            ENDIF 
     1068         END DO 
     1069         ! 
     1070         jlfil2(:,:) = jlfil(:,:)  
     1071         ! fill categories from low to high 
     1072         DO jl = 2, jpl-1 
    8631073            DO ji = 1, idim 
    864                IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN 
    865                   ! fill the right category 
    866                   zh_i(ji,jl1) = zhti(ji,jl2) 
    867                   zh_s(ji,jl1) = zhts(ji,jl2) 
    868                   za_i(ji,jl1) = zati(ji,jl2) 
    869                   ! record categories that are filled 
    870                   jlmax(ji) = MAX( jlmax(ji), jl1 ) 
    871                   jlmin(ji) = MIN( jlmin(ji), jl1 ) 
    872                   jlfil(ji,jl1) = jl1 
     1074               IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN 
     1075                  ! fill high 
     1076                  za_i(ji,jl) = ztrans * za_i(ji,jl-1) 
     1077                  zh_i(ji,jl) = hi_mean(jl) 
     1078                  jlfil(ji,jl) = jl 
     1079                  ! remove low 
     1080                  za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1) 
    8731081               ENDIF 
    8741082            END DO 
    8751083         END DO 
    876       END DO 
    877       ! 
    878       ! --- fill the gaps between categories --- !   
    879       !     transfer from categories filled at the previous step to the empty ones in between 
    880       DO ji = 1, idim 
    881          jl1 = jlmin(ji) 
    882          jl2 = jlmax(ji) 
    883          IF( jl1 > 1 ) THEN 
    884             ! fill the lower cat (jl1-1) 
    885             za_i(ji,jl1-1) = ztrans * za_i(ji,jl1) 
    886             zh_i(ji,jl1-1) = hi_mean(jl1-1) 
    887             ! remove from cat jl1 
    888             za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1) 
    889          ENDIF 
    890          IF( jl2 < jpl ) THEN 
    891             ! fill the upper cat (jl2+1) 
    892             za_i(ji,jl2+1) = ztrans * za_i(ji,jl2) 
    893             zh_i(ji,jl2+1) = hi_mean(jl2+1) 
    894             ! remove from cat jl2 
    895             za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2) 
    896          ENDIF 
    897       END DO 
    898       ! 
    899       jlfil2(:,:) = jlfil(:,:)  
    900       ! fill categories from low to high 
    901       DO jl = 2, jpl-1 
    902          DO ji = 1, idim 
    903             IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN 
    904                ! fill high 
    905                za_i(ji,jl) = ztrans * za_i(ji,jl-1) 
    906                zh_i(ji,jl) = hi_mean(jl) 
    907                jlfil(ji,jl) = jl 
    908                ! remove low 
    909                za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1) 
    910             ENDIF 
    911          END DO 
    912       END DO 
    913       ! 
    914       ! fill categories from high to low 
    915       DO jl = jpl-1, 2, -1 
    916          DO ji = 1, idim 
    917             IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN 
    918                ! fill low 
    919                za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1) 
    920                zh_i(ji,jl) = hi_mean(jl)  
    921                jlfil2(ji,jl) = jl 
    922                ! remove high 
    923                za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1) 
    924             ENDIF 
    925          END DO 
    926       END DO 
    927       ! 
    928       DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays 
    929       DEALLOCATE( jlmin, jlmax ) 
    930       ! 
    931    END SUBROUTINE ice_var_itd2 
    932  
    933  
    934    SUBROUTINE ice_var_bv 
    935       !!------------------------------------------------------------------- 
    936       !!                ***  ROUTINE ice_var_bv *** 
    937       !! 
    938       !! ** Purpose :   computes mean brine volume (%) in sea ice 
    939       !! 
    940       !! ** Method  : e = - 0.054 * S (ppt) / T (C) 
    941       !! 
    942       !! References : Vancoppenolle et al., JGR, 2007 
    943       !!------------------------------------------------------------------- 
    944       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    945       !!------------------------------------------------------------------- 
    946       ! 
    947 !!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done 
    948 !!   instead of setting everything to zero as just below 
    949       bv_i (:,:,:) = 0._wp 
    950       DO jl = 1, jpl 
    951          DO jk = 1, nlay_i 
    952             WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )    
    953                bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 ) 
    954             END WHERE 
    955          END DO 
    956       END DO 
    957       WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:) 
    958       ELSEWHERE                     ;   bvm_i(:,:) = 0._wp 
    959       END WHERE 
    960       ! 
    961    END SUBROUTINE ice_var_bv 
    962  
    963  
    964    SUBROUTINE ice_var_enthalpy 
    965       !!------------------------------------------------------------------- 
    966       !!                   ***  ROUTINE ice_var_enthalpy ***  
    967       !!                  
    968       !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature 
    969       !! 
    970       !! ** Method  :   Formula (Bitz and Lipscomb, 1999) 
    971       !!------------------------------------------------------------------- 
    972       INTEGER  ::   ji, jk   ! dummy loop indices 
    973       REAL(wp) ::   ztmelts  ! local scalar  
    974       !!------------------------------------------------------------------- 
    975       ! 
    976       DO jk = 1, nlay_i             ! Sea ice energy of melting 
    977          DO ji = 1, npti 
    978             ztmelts      = - rTmlt  * sz_i_1d(ji,jk) 
    979             t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point => likely conservation issue 
    980                                                                 !   (sometimes zdf scheme produces abnormally high temperatures)    
    981             e_i_1d(ji,jk) = rhoi * ( rcpi  * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           & 
    982                &                   + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   & 
    983                &                   - rcp   * ztmelts ) 
    984          END DO 
    985       END DO 
    986       DO jk = 1, nlay_s             ! Snow energy of melting 
    987          DO ji = 1, npti 
    988             e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) 
    989          END DO 
    990       END DO 
    991       ! 
    992    END SUBROUTINE ice_var_enthalpy 
    993  
    994    FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b) 
    995       !!--------------------------------------------------------------------- 
    996       !!                   ***  ROUTINE ice_var_sshdyn  *** 
    997       !!                      
    998       !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded 
    999       !! 
    1000       !! ** Method  :  ssh_lead = ssh + (Mice + Msnow) / rau0 
    1001       !! 
    1002       !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira, 
    1003       !!                Sea ice-ocean coupling using a rescaled vertical coordinate z*,  
    1004       !!                Ocean Modelling, Volume 24, Issues 1-2, 2008 
    1005       !!---------------------------------------------------------------------- 
    1006       ! 
    1007       ! input 
    1008       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh            !: ssh [m] 
    1009       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass    !: mass of snow and ice at current  ice time step [Kg/m2] 
    1010       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b  !: mass of snow and ice at previous ice time step [Kg/m2] 
    1011       ! 
    1012       ! output 
    1013       REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn  ! equivalent ssh in lead [m] 
    1014       ! 
    1015       ! temporary 
    1016       REAL(wp) :: zintn, zintb                     ! time interpolation weights [] 
    1017       REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload  ! snow and ice load [m] 
    1018       ! 
    1019       ! compute ice load used to define the equivalent ssh in lead 
    1020       IF( ln_ice_embd ) THEN 
    1021          !                                             
    1022          ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 
    1023          !                                               = (1/nn_fsbc)^2 * {SUM[n]        , n=0,nn_fsbc-1} 
    1024          zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp 
    1025          ! 
    1026          ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   *    {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 
    1027          !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 
    1028          zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 
    1029          ! 
    1030          zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0 
    1031          ! 
    1032       ELSE 
    1033          zsnwiceload(:,:) = 0.0_wp 
     1084         ! 
     1085         ! fill categories from high to low 
     1086         DO jl = jpl-1, 2, -1 
     1087            DO ji = 1, idim 
     1088               IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN 
     1089                  ! fill low 
     1090                  za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1) 
     1091                  zh_i(ji,jl) = hi_mean(jl)  
     1092                  jlfil2(ji,jl) = jl 
     1093                  ! remove high 
     1094                  za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1) 
     1095               ENDIF 
     1096            END DO 
     1097         END DO 
     1098         ! 
     1099         DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays 
     1100         DEALLOCATE( jlmin, jlmax ) 
     1101         ! 
    10341102      ENDIF 
    1035       ! compute equivalent ssh in lead 
    1036       ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:) 
    1037       ! 
    1038    END FUNCTION ice_var_sshdyn 
    1039  
     1103      ! 
     1104   END SUBROUTINE ice_var_itd_NcMc 
    10401105 
    10411106#else 

  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdy_oce.F90

    r11191 r11223  
    4242   TYPE, PUBLIC ::   OBC_DATA     !: Storage for external data 
    4343      INTEGER          , DIMENSION(2)   ::  nread 
    44       LOGICAL                           ::  ll_ssh 
    45       LOGICAL                           ::  ll_u2d 
    46       LOGICAL                           ::  ll_v2d 
    47       LOGICAL                           ::  ll_u3d 
    48       LOGICAL                           ::  ll_v3d 
    49       LOGICAL                           ::  ll_tem 
    50       LOGICAL                           ::  ll_sal 
    51       LOGICAL                           ::  ll_fvl 
     44      LOGICAL                           ::  lneed_ssh 
     45      LOGICAL                           ::  lneed_dyn2d 
     46      LOGICAL                           ::  lneed_dyn3d 
     47      LOGICAL                           ::  lneed_tra 
     48      LOGICAL                           ::  lneed_ice 
    5249      REAL(wp), POINTER, DIMENSION(:)   ::  ssh 
    5350      REAL(wp), POINTER, DIMENSION(:)   ::  u2d 
    … …  
    5754      REAL(wp), POINTER, DIMENSION(:,:) ::  tem 
    5855      REAL(wp), POINTER, DIMENSION(:,:) ::  sal 
    59 #if defined key_si3 
    60       LOGICAL                           ::   ll_a_i 
    61       LOGICAL                           ::   ll_h_i 
    62       LOGICAL                           ::   ll_h_s 
    63       REAL(wp), POINTER, DIMENSION(:,:) ::   a_i    !: now ice leads fraction climatology 
    64       REAL(wp), POINTER, DIMENSION(:,:) ::   h_i    !: Now ice  thickness climatology 
    65       REAL(wp), POINTER, DIMENSION(:,:) ::   h_s    !: now snow thickness 
    66 #endif 
     56      REAL(wp), POINTER, DIMENSION(:,:) ::  a_i    !: now ice leads fraction climatology 
     57      REAL(wp), POINTER, DIMENSION(:,:) ::  h_i    !: Now ice  thickness climatology 
     58      REAL(wp), POINTER, DIMENSION(:,:) ::  h_s    !: now snow thickness 
    6759#if defined key_top 
    6860      CHARACTER(LEN=20)                   :: cn_obc  !: type of boundary condition to apply 
    … …  
    8779   ! 
    8880   INTEGER                    ::   nb_bdy                   !: number of open boundary sets 
    89    INTEGER, DIMENSION(jp_bdy) ::   nb_jpk_bdy               !: number of levels in the bdy data (set < 0 if consistent with planned run) 
    9081   INTEGER, DIMENSION(jp_bdy) ::   nn_rimwidth              !: boundary rim width for Flow Relaxation Scheme 
    9182   INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P  
    … …  
    130121   INTEGER,  DIMENSION(jp_bdy)                     ::   nn_dta            !: =0 => *all* data is set to initial conditions 
    131122                                                                          !: =1 => some data to be read in from data files 
    132    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays (unstr.  bdy) 
    133    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global_z      !: workspace for reading in global depth arrays (unstr.  bdy) 
    134    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global_dz     !: workspace for reading in global depth arrays (unstr.  bdy) 
    135    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2       !: workspace for reading in global data arrays (struct. bdy) 
    136    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2_z     !: workspace for reading in global depth arrays (struct. bdy) 
    137    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2_dz    !: workspace for reading in global depth arrays (struct. bdy) 
    138123!$AGRIF_DO_NOT_TREAT 
    139124   TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process) 

  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydta.F90

    r10951 r11223  
    4343   PUBLIC   bdy_dta_init     ! routine called by nemogcm.F90 
    4444 
    45    INTEGER, ALLOCATABLE, DIMENSION(:)   ::   nb_bdy_fld        ! Number of fields to update for each boundary set. 
    46    INTEGER                              ::   nb_bdy_fld_sum    ! Total number of fields to update for all boundary sets. 
    47    LOGICAL,           DIMENSION(jp_bdy) ::   ln_full_vel_array ! =T => full velocities in 3D boundary conditions 
    48                                                                ! =F => baroclinic velocities in 3D boundary conditions 
     45   INTEGER , PARAMETER ::   jpbdyfld  = 10    ! maximum number of files to read  
     46   INTEGER , PARAMETER ::   jp_bdyssh = 1     !  
     47   INTEGER , PARAMETER ::   jp_bdyu2d = 2     !  
     48   INTEGER , PARAMETER ::   jp_bdyv2d = 3     ! 
     49   INTEGER , PARAMETER ::   jp_bdyu3d = 4     ! 
     50   INTEGER , PARAMETER ::   jp_bdyv3d = 5     ! 
     51   INTEGER , PARAMETER ::   jp_bdytem = 6     !  
     52   INTEGER , PARAMETER ::   jp_bdysal = 7     !  
     53   INTEGER , PARAMETER ::   jp_bdya_i = 8     !  
     54   INTEGER , PARAMETER ::   jp_bdyh_i = 9     !  
     55   INTEGER , PARAMETER ::   jp_bdyh_S = 10    !  
     56                                                              ! =F => baroclinic velocities in 3D boundary conditions 
    4957!$AGRIF_DO_NOT_TREAT 
    50    TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:), TARGET ::   bf        ! structure of input fields (file informations, fields read) 
     58   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:,:), TARGET ::   bf   ! structure of input fields (file informations, fields read) 
    5159!$AGRIF_END_DO_NOT_TREAT 
    52    TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:) :: nbmap_ptr   ! array of pointers to nbmap 
    53  
    54 #if defined key_si3 
    55    INTEGER ::   nice_cat                      ! number of categories in the input file 
    56    INTEGER ::   jfld_hti, jfld_hts, jfld_ai   ! indices of ice thickness, snow thickness and concentration in bf structure 
    57    INTEGER, DIMENSION(jp_bdy) :: jfld_htit, jfld_htst, jfld_ait 
    58 #endif 
    5960 
    6061   !!---------------------------------------------------------------------- 
    … …  
    6566CONTAINS 
    6667 
    67       SUBROUTINE bdy_dta( kt, jit, time_offset ) 
     68   SUBROUTINE bdy_dta( kt, kit, kt_offset ) 
    6869      !!---------------------------------------------------------------------- 
    6970      !!                   ***  SUBROUTINE bdy_dta  *** 
    … …  
    7576      !!---------------------------------------------------------------------- 
    7677      INTEGER, INTENT(in)           ::   kt           ! ocean time-step index  
    77       INTEGER, INTENT(in), OPTIONAL ::   jit          ! subcycle time-step index (for timesplitting option) 
    78       INTEGER, INTENT(in), OPTIONAL ::   time_offset  ! time offset in units of timesteps. NB. if jit 
     78      INTEGER, INTENT(in), OPTIONAL ::   kit          ! subcycle time-step index (for timesplitting option) 
     79      INTEGER, INTENT(in), OPTIONAL ::   kt_offset    ! time offset in units of timesteps. NB. if kit 
    7980      !                                               ! is present then units = subcycle timesteps. 
    80       !                                               ! time_offset = 0 => get data at "now" time level 
    81       !                                               ! time_offset = -1 => get data at "before" time level 
    82       !                                               ! time_offset = +1 => get data at "after" time level 
     81      !                                               ! kt_offset = 0 => get data at "now" time level 
     82      !                                               ! kt_offset = -1 => get data at "before" time level 
     83      !                                               ! kt_offset = +1 => get data at "after" time level 
    8384      !                                               ! etc. 
    8485      ! 
    85       INTEGER ::  jbdy, jfld, jstart, jend, ib, jl  ! dummy loop indices 
    86       INTEGER ::  ii, ij, ik, igrd                  ! local integers 
    87       INTEGER,          DIMENSION(jpbgrd) ::   ilen1  
    88       INTEGER, POINTER, DIMENSION(:)      ::   nblen, nblenrim  ! short cuts 
    89       TYPE(OBC_DATA), POINTER             ::   dta              ! short cut 
     86      INTEGER ::  jbdy, jfld, jstart, jend, ib, jl    ! dummy loop indices 
     87      INTEGER ::  ii, ij, ik, igrd, ipl               ! local integers 
     88      INTEGER,   DIMENSION(jpbgrd)     ::   ilen1  
     89      INTEGER,   DIMENSION(:), POINTER ::   nblen, nblenrim  ! short cuts 
     90      TYPE(OBC_DATA)         , POINTER ::   dta_alias        ! short cut 
     91      TYPE(FLD), DIMENSION(:), POINTER ::   bf_alias 
    9092      !!--------------------------------------------------------------------------- 
    9193      ! 
    … …  
    9496      ! Initialise data arrays once for all from initial conditions where required 
    9597      !--------------------------------------------------------------------------- 
    96       IF( kt == nit000 .AND. .NOT.PRESENT(jit) ) THEN 
     98      IF( kt == nit000 .AND. .NOT.PRESENT(kit) ) THEN 
    9799 
    98100         ! Calculate depth-mean currents 
    99101         !----------------------------- 
    100           
     102 
    101103         DO jbdy = 1, nb_bdy 
    102104            ! 
    103105            nblen    => idx_bdy(jbdy)%nblen 
    104106            nblenrim => idx_bdy(jbdy)%nblenrim 
    105             dta      => dta_bdy(jbdy) 
    106107            ! 
    107108            IF( nn_dyn2d_dta(jbdy) == 0 ) THEN  
    108109               ilen1(:) = nblen(:) 
    109                IF( dta%ll_ssh ) THEN  
     110               IF( dta_bdy(jbdy)%lneed_ssh ) THEN  
    110111                  igrd = 1 
    111112                  DO ib = 1, ilen1(igrd) 
    … …  
    113114                     ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    114115                     dta_bdy(jbdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)          
    115                   END DO  
    116                ENDIF 
    117                IF( dta%ll_u2d ) THEN  
     116                  END DO 
     117               ENDIF 
     118               IF( dta_bdy(jbdy)%lneed_dyn2d) THEN  
    118119                  igrd = 2 
    119120                  DO ib = 1, ilen1(igrd) 
    … …  
    121122                     ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    122123                     dta_bdy(jbdy)%u2d(ib) = un_b(ii,ij) * umask(ii,ij,1)          
    123                   END DO  
    124                ENDIF 
    125                IF( dta%ll_v2d ) THEN  
     124                  END DO 
    126125                  igrd = 3 
    127126                  DO ib = 1, ilen1(igrd) 
    … …  
    129128                     ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    130129                     dta_bdy(jbdy)%v2d(ib) = vn_b(ii,ij) * vmask(ii,ij,1)          
    131                   END DO  
     130                  END DO 
    132131               ENDIF 
    133132            ENDIF 
    … …  
    135134            IF( nn_dyn3d_dta(jbdy) == 0 ) THEN  
    136135               ilen1(:) = nblen(:) 
    137                IF( dta%ll_u3d ) THEN  
     136               IF( dta_bdy(jbdy)%lneed_dyn3d ) THEN  
    138137                  igrd = 2  
    139138                  DO ib = 1, ilen1(igrd) 
    … …  
    143142                        dta_bdy(jbdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - un_b(ii,ij) ) * umask(ii,ij,ik)          
    144143                     END DO 
    145                   END DO  
    146                ENDIF 
    147                IF( dta%ll_v3d ) THEN  
     144                  END DO 
    148145                  igrd = 3  
    149146                  DO ib = 1, ilen1(igrd) 
    … …  
    152149                        ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    153150                        dta_bdy(jbdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - vn_b(ii,ij) ) * vmask(ii,ij,ik)          
    154                         END DO 
    155                   END DO  
     151                     END DO 
     152                  END DO 
    156153               ENDIF 
    157154            ENDIF 
    … …  
    159156            IF( nn_tra_dta(jbdy) == 0 ) THEN  
    160157               ilen1(:) = nblen(:) 
    161                IF( dta%ll_tem ) THEN 
     158               IF( dta_bdy(jbdy)%lneed_tra ) THEN 
    162159                  igrd = 1  
    163160                  DO ib = 1, ilen1(igrd) 
    … …  
    165162                        ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    166163                        ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    167                         dta_bdy(jbdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)          
     164                        dta_bdy(jbdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_bdytem) * tmask(ii,ij,ik)          
     165                        dta_bdy(jbdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_bdysal) * tmask(ii,ij,ik)          
    168166                     END DO 
    169                   END DO  
    170                ENDIF 
    171                IF( dta%ll_sal ) THEN 
    172                   igrd = 1  
    173                   DO ib = 1, ilen1(igrd) 
    174                      DO ik = 1, jpkm1 
    175                         ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    176                         ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    177                         dta_bdy(jbdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)          
    178                      END DO 
    179                   END DO  
     167                  END DO 
    180168               ENDIF 
    181169            ENDIF 
    … …  
    184172            IF( nn_ice_dta(jbdy) == 0 ) THEN    ! set ice to initial values 
    185173               ilen1(:) = nblen(:) 
    186                IF( dta%ll_a_i ) THEN 
     174               IF( dta_bdy(jbdy)%lneed_ice ) THEN 
    187175                  igrd = 1    
    188176                  DO jl = 1, jpl 
    … …  
    191179                        ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    192180                        dta_bdy(jbdy)%a_i (ib,jl) =  a_i(ii,ij,jl) * tmask(ii,ij,1)  
    193                      END DO 
    194                   END DO 
    195                ENDIF 
    196                IF( dta%ll_h_i ) THEN 
    197                   igrd = 1    
    198                   DO jl = 1, jpl 
    199                      DO ib = 1, ilen1(igrd) 
    200                         ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    201                         ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    202181                        dta_bdy(jbdy)%h_i (ib,jl) =  h_i(ii,ij,jl) * tmask(ii,ij,1)  
    203                      END DO 
    204                   END DO 
    205                ENDIF 
    206                IF( dta%ll_h_s ) THEN 
    207                   igrd = 1    
    208                   DO jl = 1, jpl 
    209                      DO ib = 1, ilen1(igrd) 
    210                         ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    211                         ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    212182                        dta_bdy(jbdy)%h_s (ib,jl) =  h_s(ii,ij,jl) * tmask(ii,ij,1)  
    213183                     END DO 
    … …  
    222192      ! update external data from files 
    223193      !-------------------------------- 
    224       
    225       jstart = 1 
    226       DO jbdy = 1, nb_bdy    
    227          dta => dta_bdy(jbdy) 
    228          IF( nn_dta(jbdy) == 1 ) THEN ! skip this bit if no external data required 
    229        
    230             IF( PRESENT(jit) ) THEN 
    231                ! Update barotropic boundary conditions only 
    232                ! jit is optional argument for fld_read and bdytide_update 
    233                IF( cn_dyn2d(jbdy) /= 'none' ) THEN 
    234                   IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
    235                      IF( dta%ll_ssh ) dta%ssh(:) = 0._wp 
    236                      IF( dta%ll_u2d ) dta%u2d(:) = 0._wp 
    237                      IF( dta%ll_u3d ) dta%v2d(:) = 0._wp 
    238                   ENDIF 
    239                   IF (cn_tra(jbdy) /= 'runoff') THEN 
    240                      IF( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 ) THEN 
    241  
    242                         jend = jstart + dta%nread(2) - 1 
    243                         IF( ln_full_vel_array(jbdy) ) THEN 
    244                            CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  & 
    245                                      & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy(jbdy),   & 
    246                                      & fvl=ln_full_vel_array(jbdy)  ) 
    247                         ELSE 
    248                            CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  & 
    249                                      & kit=jit, kt_offset=time_offset  ) 
    250                         ENDIF 
    251  
    252                         ! If full velocities in boundary data then extract barotropic velocities from 3D fields 
    253                         IF( ln_full_vel_array(jbdy) .AND.                                             & 
    254                           &    ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR.  & 
    255                           &      nn_dyn3d_dta(jbdy) == 1 ) )THEN 
    256  
    257                            igrd = 2                      ! zonal velocity 
    258                            dta%u2d(:) = 0._wp 
    259                            DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    260                               ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
    261                               ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    262                               DO ik = 1, jpkm1 
    263                                  dta%u2d(ib) = dta%u2d(ib) & 
    264                        &                          + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik) 
    265                               END DO 
    266                               dta%u2d(ib) =  dta%u2d(ib) * r1_hu_n(ii,ij) 
    267                            END DO 
    268                            igrd = 3                      ! meridional velocity 
    269                            dta%v2d(:) = 0._wp 
    270                            DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    271                               ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
    272                               ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    273                               DO ik = 1, jpkm1 
    274                                  dta%v2d(ib) = dta%v2d(ib) & 
    275                        &                       + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik) 
    276                               END DO 
    277                               dta%v2d(ib) =  dta%v2d(ib) * r1_hv_n(ii,ij) 
    278                            END DO 
    279                         ENDIF                     
    280                      ENDIF 
    281                      IF( nn_dyn2d_dta(jbdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
    282                         CALL bdytide_update( kt=kt, idx=idx_bdy(jbdy), dta=dta, td=tides(jbdy),   &  
    283                           &                 jit=jit, time_offset=time_offset ) 
    284                      ENDIF 
    285                   ENDIF 
    286                ENDIF 
    287             ELSE 
    288                IF (cn_tra(jbdy) == 'runoff') then      ! runoff condition 
    289                   jend = nb_bdy_fld(jbdy) 
    290                   CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend),  & 
    291                                & map=nbmap_ptr(jstart:jend), kt_offset=time_offset ) 
    292                   ! 
    293                   igrd = 2                      ! zonal velocity 
    294                   DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    295                      ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
    296                      ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    297                      dta%u2d(ib) = dta%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) ) 
    298                   END DO 
    299                   ! 
    300                   igrd = 3                      ! meridional velocity 
    301                   DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    302                      ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
    303                      ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    304                      dta%v2d(ib) = dta%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) ) 
    305                   END DO 
    306                ELSE 
    307                   IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
    308                      IF( dta%ll_ssh ) dta%ssh(:) = 0._wp 
    309                      IF( dta%ll_u2d ) dta%u2d(:) = 0._wp 
    310                      IF( dta%ll_v2d ) dta%v2d(:) = 0._wp 
    311                   ENDIF 
    312                   IF( dta%nread(1) .gt. 0 ) THEN ! update external data 
    313                      jend = jstart + dta%nread(1) - 1 
    314                      CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), & 
    315                                   & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy(jbdy),   & 
    316                                   & fvl=ln_full_vel_array(jbdy) ) 
    317                   ENDIF 
    318                   ! If full velocities in boundary data then split into barotropic and baroclinic data 
    319                   IF( ln_full_vel_array(jbdy) .and.                                             & 
    320                     & ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR. & 
    321                     &   nn_dyn3d_dta(jbdy) == 1 ) ) THEN 
    322                      igrd = 2                      ! zonal velocity 
    323                      dta%u2d(:) = 0._wp 
    324                      DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    325                         ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
    326                         ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    327                         DO ik = 1, jpkm1 
    328                            dta%u2d(ib) = dta%u2d(ib) & 
    329                  &                       + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik) 
    330                         END DO 
    331                         dta%u2d(ib) =  dta%u2d(ib) * r1_hu_n(ii,ij) 
    332                         DO ik = 1, jpkm1 
    333                            dta%u3d(ib,ik) = dta%u3d(ib,ik) - dta%u2d(ib) 
    334                         END DO 
    335                      END DO 
    336                      igrd = 3                      ! meridional velocity 
    337                      dta%v2d(:) = 0._wp 
    338                      DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    339                         ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
    340                         ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    341                         DO ik = 1, jpkm1 
    342                            dta%v2d(ib) = dta%v2d(ib) & 
    343                  &                       + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik) 
    344                         END DO 
    345                         dta%v2d(ib) =  dta%v2d(ib) * r1_hv_n(ii,ij) 
    346                         DO ik = 1, jpkm1 
    347                            dta%v3d(ib,ik) = dta%v3d(ib,ik) - dta%v2d(ib) 
    348                         END DO 
    349                      END DO 
    350                   ENDIF 
    351  
    352                ENDIF 
     194 
     195      DO jbdy = 1, nb_bdy 
     196 
     197         dta_alias => dta_bdy(jbdy) 
     198         bf_alias  => bf(:,jbdy) 
     199 
     200         ! read/update all bdy data 
     201         ! ------------------------ 
     202         CALL fld_read( kt, 1, bf_alias, kit = kit, kt_offset = kt_offset ) 
     203 
     204         ! apply some corrections in some specific cases... 
     205         ! -------------------------------------------------- 
     206         ! 
     207         ! if runoff condition: change river flow we read (in m3/s) into barotropic velocity (m/s) 
     208         IF( cn_tra(jbdy) == 'runoff' .AND. TRIM(bf_alias(jp_bdyu2d)%clrootname) /= 'NOT USED' ) THEN   ! runoff and we read u/v2d 
     209            ! 
     210            igrd = 2                      ! zonal flow (m3/s) to barotropic zonal velocity (m/s) 
     211            DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
     212               ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     213               ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
     214               dta_alias%u2d(ib) = dta_alias%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) ) 
     215            END DO 
     216            igrd = 3                      ! meridional flow (m3/s) to barotropic meridional velocity (m/s) 
     217            DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
     218               ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     219               ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
     220               dta_alias%v2d(ib) = dta_alias%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) ) 
     221            END DO 
     222         ENDIF 
     223 
     224         ! tidal harmonic forcing ONLY: initialise arrays 
     225         IF( nn_dyn2d_dta(jbdy) == 2 ) THEN   ! we did not read ssh, u/v2d  
     226            IF( dta_alias%lneed_ssh   ) dta_alias%ssh(:) = 0._wp 
     227            IF( dta_alias%lneed_dyn2d ) dta_alias%u2d(:) = 0._wp 
     228            IF( dta_alias%lneed_dyn2d ) dta_alias%v2d(:) = 0._wp 
     229         ENDIF 
     230 
     231         ! If full velocities in boundary data, then split it into barotropic and baroclinic component 
     232         IF( bf_alias(jp_bdyu3d)%ltotvel ) THEN     ! if we read 3D total velocity (can be true only if u3d was read) 
     233            ! 
     234            igrd = 2                       ! zonal velocity 
     235            dta_alias%u2d(:) = 0._wp       ! compute barotrope zonal velocity and put it in u2d 
     236            DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
     237               ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     238               ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
     239               DO ik = 1, jpkm1 
     240                  dta_alias%u2d(ib) = dta_alias%u2d(ib) + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik) 
     241               END DO 
     242               dta_alias%u2d(ib) =  dta_alias%u2d(ib) * r1_hu_n(ii,ij) 
     243               DO ik = 1, jpkm1            ! compute barocline zonal velocity and put it in u3d 
     244                  dta_alias%u3d(ib,ik) = dta_alias%u3d(ib,ik) - dta_alias%u2d(ib) 
     245               END DO 
     246            END DO 
     247            igrd = 3                       ! meridional velocity 
     248            dta_alias%v2d(:) = 0._wp       ! compute barotrope meridional velocity and put it in v2d 
     249            DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
     250               ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     251               ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
     252               DO ik = 1, jpkm1 
     253                  dta_alias%v2d(ib) = dta_alias%v2d(ib) + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik) 
     254               END DO 
     255               dta_alias%v2d(ib) =  dta_alias%v2d(ib) * r1_hv_n(ii,ij) 
     256               DO ik = 1, jpkm1            ! compute barocline meridional velocity and put it in v3d 
     257                  dta_alias%v3d(ib,ik) = dta_alias%v3d(ib,ik) - dta_alias%v2d(ib) 
     258               END DO 
     259            END DO 
     260         ENDIF   ! ltotvel 
     261 
     262         ! update tidal harmonic forcing 
     263         IF( PRESENT(kit) .AND. nn_dyn2d_dta(jbdy) .GE. 2 ) THEN 
     264            CALL bdytide_update( kt = kt, idx = idx_bdy(jbdy), dta = dta_alias, td = tides(jbdy),   &  
     265               &                 kit = kit, kt_offset = kt_offset ) 
     266         ENDIF 
     267 
     268         !  atm surface pressure : add inverted barometer effect to ssh if it was read 
     269         IF ( ln_apr_obc .AND. TRIM(bf_alias(jp_bdyssh)%clrootname) /= 'NOT USED' ) THEN 
     270            igrd = 1 
     271            DO ib = 1, idx_bdy(jbdy)%nblenrim(igrd)   ! ssh is used only on the rim 
     272               ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     273               ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
     274               dta_alias%ssh(ib) = dta_alias%ssh(ib) + ssh_ib(ii,ij) 
     275            END DO 
     276         ENDIF 
     277 
    353278#if defined key_si3 
    354                ! convert N-cat fields (input) into jpl-cat (output) 
    355                IF( cn_ice(jbdy) /= 'none' .AND. nn_ice_dta(jbdy) == 1 ) THEN 
    356                   jfld_hti = jfld_htit(jbdy) 
    357                   jfld_hts = jfld_htst(jbdy) 
    358                   jfld_ai  = jfld_ait(jbdy) 
    359                   IF    ( jpl /= 1 .AND. nice_cat == 1 ) THEN                       ! case input cat = 1 
    360                      CALL ice_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), & 
    361                         &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    ) 
    362                   ELSEIF( jpl /= 1 .AND. nice_cat /= 1 .AND. nice_cat /= jpl ) THEN ! case input cat /=1 and /=jpl 
    363                      CALL ice_var_itd2( bf(jfld_hti)%fnow(:,1,:), bf(jfld_hts)%fnow(:,1,:), bf(jfld_ai)%fnow(:,1,:), & 
    364                         &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    ) 
    365                   ENDIF 
    366                ENDIF 
     279         ! ice: convert N-cat fields (input) into jpl-cat (output) 
     280         IF( dta_alias%lneed_ice ) THEN 
     281            ipl = SIZE(bf_alias(jp_bdya_i)%fnow, 3) 
     282            IF( ipl /= jpl ) THEN      ! ice: convert N-cat fields (input) into jpl-cat (output) 
     283               CALL ice_var_itd(bf_alias(jp_bdyh_i)%fnow(:,1,:), bf_alias(jp_bdyh_s)%fnow(:,1,:), bf_alias(jp_bdya_i)%fnow(:,1,:), & 
     284                  &              dta_alias%h_i               , dta_alias%h_s               , dta_alias%a_i                 ) 
     285            ENDIF 
     286         ENDIF 
    367287#endif 
    368             ENDIF 
    369             jstart = jstart + dta%nread(1) 
    370          ENDIF    ! nn_dta(jbdy) = 1 
    371288      END DO  ! jbdy 
    372  
    373       IF ( ln_apr_obc ) THEN 
    374          DO jbdy = 1, nb_bdy 
    375             IF (cn_tra(jbdy) /= 'runoff')THEN 
    376                igrd = 1                      ! meridional velocity 
    377                DO ib = 1, idx_bdy(jbdy)%nblenrim(igrd) 
    378                   ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
    379                   ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    380                   dta_bdy(jbdy)%ssh(ib) = dta_bdy(jbdy)%ssh(ib) + ssh_ib(ii,ij) 
    381                END DO 
    382             ENDIF 
    383          END DO 
    384       ENDIF 
    385289 
    386290      IF ( ln_tide ) THEN 
    387291         IF (ln_dynspg_ts) THEN      ! Fill temporary arrays with slow-varying bdy data                            
    388             DO jbdy = 1, nb_bdy    ! Tidal component added in ts loop 
    389                IF ( nn_dyn2d_dta(jbdy) .ge. 2 ) THEN 
     292            DO jbdy = 1, nb_bdy      ! Tidal component added in ts loop 
     293               IF ( nn_dyn2d_dta(jbdy) .GE. 2 ) THEN 
    390294                  nblen => idx_bdy(jbdy)%nblen 
    391295                  nblenrim => idx_bdy(jbdy)%nblenrim 
    392                   IF( cn_dyn2d(jbdy) == 'frs' ) THEN; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF  
    393                   IF ( dta_bdy(jbdy)%ll_ssh ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1)) 
    394                   IF ( dta_bdy(jbdy)%ll_u2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2)) 
    395                   IF ( dta_bdy(jbdy)%ll_v2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3)) 
    396                ENDIF 
    397             END DO 
    398          ELSE ! Add tides if not split-explicit free surface else this is done in ts loop 
    399             ! 
    400             CALL bdy_dta_tides( kt=kt, time_offset=time_offset ) 
    401          ENDIF 
    402       ENDIF 
    403  
    404       ! 
    405       IF( ln_timing )   CALL timing_stop('bdy_dta') 
    406       ! 
    407    END SUBROUTINE bdy_dta 
     296                  IF( cn_dyn2d(jbdy) == 'frs' ) THEN ; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF  
     297                     IF ( dta_bdy(jbdy)%lneed_ssh   ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1)) 
     298                     IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2)) 
     299                     IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3)) 
     300                  ENDIF 
     301               END DO 
     302            ELSE ! Add tides if not split-explicit free surface else this is done in ts loop 
     303               ! 
     304               CALL bdy_dta_tides( kt=kt, kt_offset=kt_offset ) 
     305            ENDIF 
     306         ENDIF 
     307         ! 
     308         IF( ln_timing )   CALL timing_stop('bdy_dta') 
     309         ! 
     310      END SUBROUTINE bdy_dta 
    408311 
    409312 
    … …  
    418321      !!                 
    419322      !!---------------------------------------------------------------------- 
    420       INTEGER ::   jbdy, jfld, jstart, jend, ierror, ios     ! Local integers 
    421       ! 
     323      INTEGER ::   jbdy, jfld    ! Local integers 
     324      INTEGER ::   ierror, ios     !  
     325      ! 
     326      CHARACTER(len=3)                       ::   cl3           !  
    422327      CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files 
    423       CHARACTER(len=100), DIMENSION(nb_bdy)  ::   cn_dir_array  ! Root directory for location of data files 
    424       CHARACTER(len = 256)::   clname                           ! temporary file name 
    425328      LOGICAL                                ::   ln_full_vel   ! =T => full velocities in 3D boundary data 
    426329      !                                                         ! =F => baroclinic velocities in 3D boundary data 
    427       INTEGER                                ::   ilen_global   ! Max length required for global bdy dta arrays 
    428       INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ilen1, ilen3  ! size of 1st and 3rd dimensions of local arrays 
    429       INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ibdy           ! bdy set for a particular jfld 
    430       INTEGER, ALLOCATABLE, DIMENSION(:)     ::   igrid         ! index for grid type (1,2,3 = T,U,V) 
    431       INTEGER, POINTER, DIMENSION(:)         ::   nblen, nblenrim  ! short cuts 
    432       TYPE(OBC_DATA), POINTER                ::   dta           ! short cut 
    433 #if defined key_si3 
    434       INTEGER               ::   kndims   ! number of dimensions in an array (i.e. 3 = wo ice cat; 4 = w ice cat) 
    435       INTEGER, DIMENSION(4) ::   kdimsz   ! size   of dimensions 
    436       INTEGER               ::   inum,id1 ! local integer 
    437 #endif 
    438       TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   blf_i         !  array of namelist information structures 
    439       TYPE(FLD_N) ::   bn_tem, bn_sal, bn_u3d, bn_v3d   !  
    440       TYPE(FLD_N) ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read 
    441 #if defined key_si3 
    442       TYPE(FLD_N) ::   bn_a_i, bn_h_i, bn_h_s       
    443 #endif 
     330      INTEGER                                ::   ipk,ipl           ! 
     331      INTEGER                                ::   idvar         ! variable ID 
     332      INTEGER                                ::   indims        ! number of dimensions of the variable 
     333      INTEGER                                ::   iszdim        ! number of dimensions of the variable 
     334      INTEGER, DIMENSION(4)                  ::   i4dimsz       ! size of variable dimensions  
     335      INTEGER                                ::   igrd          ! index for grid type (1,2,3 = T,U,V) 
     336      LOGICAL                                ::   lluld         ! is the variable using the unlimited dimension 
     337      LOGICAL                                ::   llneed        ! 
     338      LOGICAL                                ::   llread        ! 
     339      TYPE(FLD_N), DIMENSION(1), TARGET ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! must be an array to be used with fld_fill 
     340      TYPE(FLD_N), DIMENSION(1), TARGET ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read 
     341      TYPE(FLD_N), DIMENSION(1), TARGET ::   bn_a_i, bn_h_i, bn_h_s       
     342      TYPE(FLD_N), DIMENSION(:), POINTER ::   bn_alias                        ! must be an array to be used with fld_fill 
     343      TYPE(FLD  ), DIMENSION(:), POINTER ::   bf_alias 
     344      ! 
    444345      NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d  
    445 #if defined key_si3 
    446346      NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s 
    447 #endif 
    448347      NAMELIST/nambdy_dta/ ln_full_vel 
    449348      !!--------------------------------------------------------------------------- 
    … …  
    454353      IF(lwp) WRITE(numout,*) '' 
    455354 
    456       ! Set nn_dta 
    457       DO jbdy = 1, nb_bdy 
    458          nn_dta(jbdy) = MAX(   nn_dyn2d_dta  (jbdy)    & 
    459             &                , nn_dyn3d_dta  (jbdy)    & 
    460             &                , nn_tra_dta    (jbdy)    & 
    461 #if defined key_si3 
    462             &                , nn_ice_dta    (jbdy)    & 
    463 #endif 
    464                               ) 
    465          IF(nn_dta(jbdy) > 1)   nn_dta(jbdy) = 1 
    466       END DO 
    467  
    468       ! Work out upper bound of how many fields there are to read in and allocate arrays 
    469       ! --------------------------------------------------------------------------- 
    470       ALLOCATE( nb_bdy_fld(nb_bdy) ) 
    471       nb_bdy_fld(:) = 0 
    472       DO jbdy = 1, nb_bdy          
    473          IF( cn_dyn2d(jbdy) /= 'none' .AND. ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) ) THEN 
    474             nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 3 
    475          ENDIF 
    476          IF( cn_dyn3d(jbdy) /= 'none' .AND. nn_dyn3d_dta(jbdy) == 1 ) THEN 
    477             nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 2 
    478          ENDIF 
    479          IF( cn_tra(jbdy) /= 'none' .AND. nn_tra_dta(jbdy) == 1  ) THEN 
    480             nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 2 
    481          ENDIF 
    482 #if defined key_si3 
    483          IF( cn_ice(jbdy) /= 'none' .AND. nn_ice_dta(jbdy) == 1  ) THEN 
    484             nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 3 
    485          ENDIF 
    486 #endif                
    487          IF(lwp) WRITE(numout,*) 'Maximum number of files to open =', nb_bdy_fld(jbdy) 
    488       END DO             
    489  
    490       nb_bdy_fld_sum = SUM( nb_bdy_fld ) 
    491  
    492       ALLOCATE( bf(nb_bdy_fld_sum), STAT=ierror ) 
     355      ALLOCATE( bf(jpbdyfld,nb_bdy), STAT=ierror ) 
    493356      IF( ierror > 0 ) THEN    
    494357         CALL ctl_stop( 'bdy_dta: unable to allocate bf structure' )   ;   RETURN   
    495358      ENDIF 
    496       ALLOCATE( blf_i(nb_bdy_fld_sum), STAT=ierror ) 
    497       IF( ierror > 0 ) THEN    
    498          CALL ctl_stop( 'bdy_dta: unable to allocate blf_i structure' )   ;   RETURN   
    499       ENDIF 
    500       ALLOCATE( nbmap_ptr(nb_bdy_fld_sum), STAT=ierror ) 
    501       IF( ierror > 0 ) THEN    
    502          CALL ctl_stop( 'bdy_dta: unable to allocate nbmap_ptr structure' )   ;   RETURN   
    503       ENDIF 
    504       ALLOCATE( ilen1(nb_bdy_fld_sum), ilen3(nb_bdy_fld_sum) )  
    505       ALLOCATE( ibdy(nb_bdy_fld_sum) )  
    506       ALLOCATE( igrid(nb_bdy_fld_sum) )  
    507  
     359      bf(:,:)%clrootname = 'NOT USED'   ! default definition used as a flag in fld_read to do nothing. 
     360      bf(:,:)%ltotvel = .FALSE.         ! default definition 
     361  
    508362      ! Read namelists 
    509363      ! -------------- 
    510364      REWIND(numnam_cfg) 
    511       jfld = 0  
    512       DO jbdy = 1, nb_bdy          
    513          IF( nn_dta(jbdy) == 1 ) THEN 
    514             REWIND(numnam_ref) 
    515             READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901) 
    516 901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist', lwp ) 
     365      DO jbdy = 1, nb_bdy 
     366 
     367         WRITE(ctmp1, '(a,i2)') 'BDY number ', jbdy 
     368         WRITE(ctmp2, '(a,i2)') 'block nambdy_dta number ', jbdy 
     369 
     370         ! There is only one nambdy_dta block in namelist_ref -> use it for each bdy so we do a rewind  
     371         REWIND(numnam_ref) 
     372         READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901) 
     373901      IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist', lwp ) 
     374 
     375         !   by-pass nambdy_dta reading if no input data used in this bdy    
     376         IF(       ( dta_bdy(jbdy)%lneed_dyn2d .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1 )   & 
     377            & .OR. ( dta_bdy(jbdy)%lneed_dyn3d .AND.     nn_dyn3d_dta(jbdy)    == 1 )   & 
     378            & .OR. ( dta_bdy(jbdy)%lneed_tra   .AND.       nn_tra_dta(jbdy)    == 1 )   & 
     379            & .OR. ( dta_bdy(jbdy)%lneed_ice   .AND.       nn_ice_dta(jbdy)    == 1 )   )   THEN 
     380            ! WARNING: we don't do a rewind here, each bdy reads its own nambdy_dta block one after another 
    517381            READ  ( numnam_cfg, nambdy_dta, IOSTAT = ios, ERR = 902 ) 
    518382902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_dta in configuration namelist', lwp ) 
    519             IF(lwm) WRITE( numond, nambdy_dta ) 
    520  
    521             cn_dir_array(jbdy) = cn_dir 
    522             ln_full_vel_array(jbdy) = ln_full_vel 
    523  
    524             nblen => idx_bdy(jbdy)%nblen 
    525             nblenrim => idx_bdy(jbdy)%nblenrim 
    526             dta => dta_bdy(jbdy) 
    527             dta%nread(2) = 0 
    528  
    529             ! Only read in necessary fields for this set. 
    530             ! Important that barotropic variables come first. 
    531             IF( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) THEN  
    532  
    533                IF( dta%ll_ssh ) THEN  
    534                   if(lwp) write(numout,*) '++++++ reading in ssh field' 
    535                   jfld = jfld + 1 
    536                   blf_i(jfld) = bn_ssh 
    537                   ibdy(jfld) = jbdy 
    538                   igrid(jfld) = 1 
    539                   ilen1(jfld) = nblen(igrid(jfld)) 
    540                   ilen3(jfld) = 1 
    541                   dta%nread(2) = dta%nread(2) + 1 
    542                ENDIF 
    543  
    544                IF( dta%ll_u2d .and. .not. ln_full_vel_array(jbdy) ) THEN 
    545                   if(lwp) write(numout,*) '++++++ reading in u2d field' 
    546                   jfld = jfld + 1 
    547                   blf_i(jfld) = bn_u2d 
    548                   ibdy(jfld) = jbdy 
    549                   igrid(jfld) = 2 
    550                   ilen1(jfld) = nblen(igrid(jfld)) 
    551                   ilen3(jfld) = 1 
    552                   dta%nread(2) = dta%nread(2) + 1 
    553                ENDIF 
    554  
    555                IF( dta%ll_v2d .and. .not. ln_full_vel_array(jbdy) ) THEN 
    556                   if(lwp) write(numout,*) '++++++ reading in v2d field' 
    557                   jfld = jfld + 1 
    558                   blf_i(jfld) = bn_v2d 
    559                   ibdy(jfld) = jbdy 
    560                   igrid(jfld) = 3 
    561                   ilen1(jfld) = nblen(igrid(jfld)) 
    562                   ilen3(jfld) = 1 
    563                   dta%nread(2) = dta%nread(2) + 1 
    564                ENDIF 
    565  
    566             ENDIF 
    567  
    568             ! read 3D velocities if baroclinic velocities require OR if 
    569             ! barotropic velocities required and ln_full_vel set to .true. 
    570             IF( nn_dyn3d_dta(jbdy) == 1 .OR. & 
    571            &  ( ln_full_vel_array(jbdy) .AND. ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 ) ) ) THEN 
    572  
    573                IF( dta%ll_u3d .OR. ( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) ) THEN 
    574                   if(lwp) write(numout,*) '++++++ reading in u3d field' 
    575                   jfld = jfld + 1 
    576                   blf_i(jfld) = bn_u3d 
    577                   ibdy(jfld) = jbdy 
    578                   igrid(jfld) = 2 
    579                   ilen1(jfld) = nblen(igrid(jfld)) 
    580                   ilen3(jfld) = jpk 
    581                   IF( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) dta%nread(2) = dta%nread(2) + 1 
    582                ENDIF 
    583  
    584                IF( dta%ll_v3d .OR. ( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) ) THEN 
    585                   if(lwp) write(numout,*) '++++++ reading in v3d field' 
    586                   jfld = jfld + 1 
    587                   blf_i(jfld) = bn_v3d 
    588                   ibdy(jfld) = jbdy 
    589                   igrid(jfld) = 3 
    590                   ilen1(jfld) = nblen(igrid(jfld)) 
    591                   ilen3(jfld) = jpk 
    592                   IF( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) dta%nread(2) = dta%nread(2) + 1 
    593                ENDIF 
    594  
    595             ENDIF 
    596  
    597             ! temperature and salinity 
    598             IF( nn_tra_dta(jbdy) == 1 ) THEN 
    599  
    600                IF( dta%ll_tem ) THEN 
    601                   if(lwp) write(numout,*) '++++++ reading in tem field' 
    602                   jfld = jfld + 1 
    603                   blf_i(jfld) = bn_tem 
    604                   ibdy(jfld) = jbdy 
    605                   igrid(jfld) = 1 
    606                   ilen1(jfld) = nblen(igrid(jfld)) 
    607                   ilen3(jfld) = jpk 
    608                ENDIF 
    609  
    610                IF( dta%ll_sal ) THEN 
    611                   if(lwp) write(numout,*) '++++++ reading in sal field' 
    612                   jfld = jfld + 1 
    613                   blf_i(jfld) = bn_sal 
    614                   ibdy(jfld) = jbdy 
    615                   igrid(jfld) = 1 
    616                   ilen1(jfld) = nblen(igrid(jfld)) 
    617                   ilen3(jfld) = jpk 
    618                ENDIF 
    619  
    620             ENDIF 
    621  
    622 #if defined key_si3 
    623             ! sea ice 
    624             IF( nn_ice_dta(jbdy) == 1 ) THEN 
    625                ! Test for types of ice input (1cat or Xcat)  
    626                ! Build file name to find dimensions  
    627                clname=TRIM( cn_dir )//TRIM(bn_a_i%clname) 
    628                IF( .NOT. bn_a_i%ln_clim ) THEN    
    629                                                   WRITE(clname, '(a,"_y",i4.4)' ) TRIM( clname ), nyear    ! add year 
    630                   IF( bn_a_i%cltype /= 'yearly' ) WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname ), nmonth   ! add month 
    631                ELSE 
    632                   IF( bn_a_i%cltype /= 'yearly' ) WRITE(clname, '(a,"_m",i2.2)' ) TRIM( clname ), nmonth   ! add month 
    633                ENDIF 
    634                IF( bn_a_i%cltype == 'daily' .OR. bn_a_i%cltype(1:4) == 'week' ) & 
    635                &                                  WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname ), nday     ! add day 
     383            IF(lwm) WRITE( numond, nambdy_dta )            
     384         ENDIF 
     385 
     386         ! get the number of ice categories in bdy data file (use a_i information to do this) 
     387         ipl = jpl   ! default definition 
     388         IF( dta_bdy(jbdy)%lneed_ice ) THEN    ! if we need ice bdy data 
     389            IF( nn_ice_dta(jbdy) == 1 ) THEN   ! if we get ice bdy data from netcdf file 
     390               CALL fld_fill(  bf(jp_bdya_i,jbdy:jbdy), bn_a_i, cn_dir, 'bdy_dta', 'a_i'//' '//ctmp1, ctmp2 )   ! use namelist info 
     391               CALL fld_clopn( bf(jp_bdya_i,jbdy), nyear, nmonth, nday )   ! not a problem when we call it again after 
     392               idvar = iom_varid( bf(jp_bdya_i,jbdy)%num, bf(jp_bdya_i,jbdy)%clvar, kndims=indims, kdimsz=i4dimsz, lduld=lluld ) 
     393               IF( indims == 4 .OR. ( indims == 3 .AND. .NOT. lluld ) ) THEN   ;   ipl = i4dimsz(3)   ! xylt or xyl 
     394               ELSE                                                            ;   ipl = 1            ! xy or xyt 
     395               ENDIF 
     396            ENDIF 
     397         ENDIF 
     398 
     399         DO jfld = 1, jpbdyfld 
     400 
     401            ! ===================== 
     402            !          ssh  
     403            ! ===================== 
     404            IF( jfld == jp_bdyssh ) THEN 
     405               cl3 = 'ssh' 
     406               igrd = 1                                                    ! T point 
     407               ipk = 1                                                     ! surface data 
     408               llneed = dta_bdy(jbdy)%lneed_ssh                            ! dta_bdy(jbdy)%ssh will be needed 
     409               llread = MOD(nn_dyn2d_dta(jbdy),2) == 1                     ! get data from NetCDF file 
     410               bf_alias => bf(jp_bdyssh,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy 
     411               bn_alias => bn_ssh                                          ! alias for ssh structure of nambdy_dta  
     412            ENDIF 
     413            ! ===================== 
     414            !         dyn2d 
     415            ! ===================== 
     416            IF( jfld == jp_bdyu2d ) THEN 
     417               cl3 = 'u2d' 
     418               igrd = 2                                                    ! U point 
     419               ipk = 1                                                     ! surface data 
     420               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed 
     421               llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1   ! don't get u2d from u3d and read NetCDF file 
     422               bf_alias => bf(jp_bdyu2d,jbdy:jbdy)                         ! alias for u2d structure of bdy number jbdy 
     423               bn_alias => bn_u2d                                          ! alias for u2d structure of nambdy_dta  
     424            ENDIF 
     425            IF( jfld == jp_bdyv2d ) THEN 
     426               cl3 = 'v2d' 
     427               igrd = 3                                                    ! V point 
     428               ipk = 1                                                     ! surface data 
     429               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed 
     430               llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1   ! don't get v2d from v3d and read NetCDF file 
     431               bf_alias => bf(jp_bdyv2d,jbdy:jbdy)                         ! alias for v2d structure of bdy number jbdy 
     432               bn_alias => bn_v2d                                          ! alias for v2d structure of nambdy_dta  
     433            ENDIF 
     434            ! ===================== 
     435            !         dyn3d 
     436            ! ===================== 
     437            IF( jfld == jp_bdyu3d ) THEN 
     438               cl3 = 'u3d' 
     439               igrd = 2                                                    ! U point 
     440               ipk = jpk                                                   ! 3d data 
     441               llneed = dta_bdy(jbdy)%lneed_dyn3d .OR.               &     ! dta_bdy(jbdy)%u3d will be needed 
     442                  &   ( dta_bdy(jbdy)%lneed_dyn2d .AND. ln_full_vel )      !   u3d needed to compute u2d 
     443               llread = nn_dyn3d_dta(jbdy) == 1                            ! get data from NetCDF file 
     444               bf_alias => bf(jp_bdyu3d,jbdy:jbdy)                         ! alias for u3d structure of bdy number jbdy 
     445               bn_alias => bn_u3d                                          ! alias for u3d structure of nambdy_dta  
     446            ENDIF 
     447            IF( jfld == jp_bdyv3d ) THEN 
     448               cl3 = 'v3d' 
     449               igrd = 3                                                    ! V point 
     450               ipk = jpk                                                   ! 3d data 
     451               llneed = dta_bdy(jbdy)%lneed_dyn3d .OR.               &     ! dta_bdy(jbdy)%v3d will be needed 
     452                  &   ( dta_bdy(jbdy)%lneed_dyn2d .AND. ln_full_vel )      !   v3d needed to compute v2d 
     453               llread = nn_dyn3d_dta(jbdy) == 1                            ! get data from NetCDF file 
     454               bf_alias => bf(jp_bdyv3d,jbdy:jbdy)                         ! alias for v3d structure of bdy number jbdy 
     455               bn_alias => bn_v3d                                          ! alias for v3d structure of nambdy_dta  
     456            ENDIF 
     457 
     458            ! ===================== 
     459            !          tra 
     460            ! ===================== 
     461            IF( jfld == jp_bdytem ) THEN 
     462               cl3 = 'tem' 
     463               igrd = 1                                                    ! T point 
     464               ipk = jpk                                                   ! 3d data 
     465               llneed = dta_bdy(jbdy)%lneed_tra                            ! dta_bdy(jbdy)%tem will be needed 
     466               llread = nn_tra_dta(jbdy) == 1                              ! get data from NetCDF file 
     467               bf_alias => bf(jp_bdytem,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy 
     468               bn_alias => bn_tem                                          ! alias for ssh structure of nambdy_dta  
     469            ENDIF 
     470            IF( jfld == jp_bdysal ) THEN 
     471               cl3 = 'sal' 
     472               igrd = 1                                                    ! T point 
     473               ipk = jpk                                                   ! 3d data 
     474               llneed = dta_bdy(jbdy)%lneed_tra                            ! dta_bdy(jbdy)%sal will be needed 
     475               llread = nn_tra_dta(jbdy) == 1                              ! get data from NetCDF file 
     476               bf_alias => bf(jp_bdysal,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy 
     477               bn_alias => bn_sal                                          ! alias for ssh structure of nambdy_dta  
     478            ENDIF 
     479 
     480            ! ===================== 
     481            !          ice 
     482            ! ===================== 
     483            IF( jfld == jp_bdya_i ) THEN 
     484               cl3 = 'a_i' 
     485               igrd = 1                                                    ! T point 
     486               ipk = ipl                                                   !  
     487               llneed = dta_bdy(jbdy)%lneed_ice                            ! dta_bdy(jbdy)%a_i will be needed 
     488               llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file 
     489               bf_alias => bf(jp_bdya_i,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy 
     490               bn_alias => bn_a_i                                          ! alias for ssh structure of nambdy_dta  
     491            ENDIF 
     492            IF( jfld == jp_bdyh_i ) THEN 
     493               cl3 = 'h_i' 
     494               igrd = 1                                                    ! T point 
     495               ipk = ipl                                                   !  
     496               llneed = dta_bdy(jbdy)%lneed_ice                            ! dta_bdy(jbdy)%h_i will be needed 
     497               llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file 
     498               bf_alias => bf(jp_bdyh_i,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy 
     499               bn_alias => bn_h_i                                          ! alias for ssh structure of nambdy_dta  
     500            ENDIF 
     501            IF( jfld == jp_bdyh_s ) THEN 
     502               cl3 = 'h_s' 
     503               igrd = 1                                                    ! T point 
     504               ipk = ipl                                                   !  
     505               llneed = dta_bdy(jbdy)%lneed_ice                            ! dta_bdy(jbdy)%h_s will be needed 
     506               llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file 
     507               bf_alias => bf(jp_bdyh_s,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy 
     508               bn_alias => bn_h_s                                          ! alias for ssh structure of nambdy_dta  
     509            ENDIF 
     510 
     511 
     512            IF( llneed ) THEN                                              ! dta_bdy(jbdy)%xxx will be needed 
     513               !                                                           !   -> must be associated with an allocated target 
     514               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus 
     515               ALLOCATE( bf_alias(1)%fnow( iszdim, 1, ipk ) )              ! allocate the target 
    636516               ! 
    637                CALL iom_open  ( clname, inum ) 
    638                id1 = iom_varid( inum, bn_a_i%clvar, kdimsz=kdimsz, kndims=kndims, ldstop = .FALSE. ) 
    639                CALL iom_close ( inum ) 
    640  
    641                 IF ( kndims == 4 ) THEN 
    642                  nice_cat = kdimsz(4)   ! Xcat input 
    643                ELSE 
    644                  nice_cat = 1           ! 1cat input       
    645                ENDIF 
    646                ! End test 
    647  
    648                IF( dta%ll_a_i ) THEN 
    649                   jfld = jfld + 1 
    650                   blf_i(jfld) = bn_a_i 
    651                   ibdy(jfld)  = jbdy 
    652                   igrid(jfld) = 1 
    653                   ilen1(jfld) = nblen(igrid(jfld)) 
    654                   ilen3(jfld) = nice_cat 
    655                ENDIF 
    656  
    657                IF( dta%ll_h_i ) THEN 
    658                   jfld = jfld + 1 
    659                   blf_i(jfld) = bn_h_i 
    660                   ibdy(jfld)  = jbdy 
    661                   igrid(jfld) = 1 
    662                   ilen1(jfld) = nblen(igrid(jfld)) 
    663                   ilen3(jfld) = nice_cat 
    664                ENDIF 
    665  
    666                IF( dta%ll_h_s ) THEN 
    667                   jfld = jfld + 1 
    668                   blf_i(jfld) = bn_h_s 
    669                   ibdy(jfld)  = jbdy 
    670                   igrid(jfld) = 1 
    671                   ilen1(jfld) = nblen(igrid(jfld)) 
    672                   ilen3(jfld) = nice_cat 
    673                ENDIF 
    674  
    675             ENDIF 
    676 #endif 
    677             ! Recalculate field counts 
    678             !------------------------- 
    679             IF( jbdy == 1 ) THEN  
    680                nb_bdy_fld_sum = 0 
    681                nb_bdy_fld(jbdy) = jfld 
    682                nb_bdy_fld_sum     = jfld               
    683             ELSE 
    684                nb_bdy_fld(jbdy) = jfld - nb_bdy_fld_sum 
    685                nb_bdy_fld_sum = nb_bdy_fld_sum + nb_bdy_fld(jbdy) 
    686             ENDIF 
    687  
    688             dta%nread(1) = nb_bdy_fld(jbdy) 
    689  
    690          ENDIF ! nn_dta == 1 
    691       ENDDO ! jbdy 
    692  
    693       DO jfld = 1, nb_bdy_fld_sum 
    694          ALLOCATE( bf(jfld)%fnow(ilen1(jfld),1,ilen3(jfld)) ) 
    695          IF( blf_i(jfld)%ln_tint ) ALLOCATE( bf(jfld)%fdta(ilen1(jfld),1,ilen3(jfld),2) ) 
    696          nbmap_ptr(jfld)%ptr => idx_bdy(ibdy(jfld))%nbmap(:,igrid(jfld)) 
    697          nbmap_ptr(jfld)%ll_unstruc = ln_coords_file(ibdy(jfld)) 
    698       ENDDO 
    699  
    700       ! fill bf with blf_i and control print 
    701       !------------------------------------- 
    702       jstart = 1 
    703       DO jbdy = 1, nb_bdy 
    704          jend = jstart - 1 + nb_bdy_fld(jbdy)  
    705          CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(jbdy), 'bdy_dta',   & 
    706          &              'open boundary conditions', 'nambdy_dta' ) 
    707          jstart = jend + 1 
    708       ENDDO 
    709  
    710       DO jfld = 1, nb_bdy_fld_sum 
    711          bf(jfld)%igrd = igrid(jfld) 
    712          bf(jfld)%ibdy = ibdy(jfld) 
    713       ENDDO 
    714  
    715       ! Initialise local boundary data arrays 
    716       ! nn_xxx_dta=0 : allocate space - will be filled from initial conditions later 
    717       ! nn_xxx_dta=1 : point to "fnow" arrays 
    718       !------------------------------------- 
    719  
    720       jfld = 0 
    721       DO jbdy=1, nb_bdy 
    722  
    723          nblen => idx_bdy(jbdy)%nblen 
    724          dta => dta_bdy(jbdy) 
    725  
    726          if(lwp) then 
    727             write(numout,*) '++++++ dta%ll_ssh = ',dta%ll_ssh 
    728             write(numout,*) '++++++ dta%ll_u2d = ',dta%ll_u2d 
    729             write(numout,*) '++++++ dta%ll_v2d = ',dta%ll_v2d 
    730             write(numout,*) '++++++ dta%ll_u3d = ',dta%ll_u3d 
    731             write(numout,*) '++++++ dta%ll_v3d = ',dta%ll_v3d 
    732             write(numout,*) '++++++ dta%ll_tem = ',dta%ll_tem 
    733             write(numout,*) '++++++ dta%ll_sal = ',dta%ll_sal 
    734          endif 
    735  
    736          IF ( nn_dyn2d_dta(jbdy) == 0 .or. nn_dyn2d_dta(jbdy) == 2 ) THEN 
    737             if(lwp) write(numout,*) '++++++ dta%ssh/u2d/u3d allocated space' 
    738             IF( dta%ll_ssh ) ALLOCATE( dta%ssh(nblen(1)) ) 
    739             IF( dta%ll_u2d ) ALLOCATE( dta%u2d(nblen(2)) ) 
    740             IF( dta%ll_v2d ) ALLOCATE( dta%v2d(nblen(3)) ) 
    741          ENDIF 
    742          IF ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) THEN 
    743             IF( dta%ll_ssh ) THEN 
    744                if(lwp) write(numout,*) '++++++ dta%ssh pointing to fnow' 
    745                jfld = jfld + 1 
    746                dta%ssh => bf(jfld)%fnow(:,1,1) 
    747             ENDIF 
    748             IF ( dta%ll_u2d ) THEN 
    749                IF ( ln_full_vel_array(jbdy) ) THEN 
    750                   if(lwp) write(numout,*) '++++++ dta%u2d allocated space' 
    751                   ALLOCATE( dta%u2d(nblen(2)) ) 
    752                ELSE 
    753                   if(lwp) write(numout,*) '++++++ dta%u2d pointing to fnow' 
    754                   jfld = jfld + 1 
    755                   dta%u2d => bf(jfld)%fnow(:,1,1) 
    756                ENDIF 
    757             ENDIF 
    758             IF ( dta%ll_v2d ) THEN 
    759                IF ( ln_full_vel_array(jbdy) ) THEN 
    760                   if(lwp) write(numout,*) '++++++ dta%v2d allocated space' 
    761                   ALLOCATE( dta%v2d(nblen(3)) ) 
    762                ELSE 
    763                   if(lwp) write(numout,*) '++++++ dta%v2d pointing to fnow' 
    764                   jfld = jfld + 1 
    765                   dta%v2d => bf(jfld)%fnow(:,1,1) 
    766                ENDIF 
    767             ENDIF 
    768          ENDIF 
    769  
    770          IF ( nn_dyn3d_dta(jbdy) == 0 ) THEN 
    771             if(lwp) write(numout,*) '++++++ dta%u3d/v3d allocated space' 
    772             IF( dta%ll_u3d ) ALLOCATE( dta_bdy(jbdy)%u3d(nblen(2),jpk) ) 
    773             IF( dta%ll_v3d ) ALLOCATE( dta_bdy(jbdy)%v3d(nblen(3),jpk) ) 
    774          ENDIF 
    775          IF ( nn_dyn3d_dta(jbdy) == 1 .or. & 
    776            &  ( ln_full_vel_array(jbdy) .and. ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) ) ) THEN 
    777             IF ( dta%ll_u3d .or. ( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) ) THEN 
    778                if(lwp) write(numout,*) '++++++ dta%u3d pointing to fnow' 
    779                jfld = jfld + 1 
    780                dta_bdy(jbdy)%u3d => bf(jfld)%fnow(:,1,:) 
    781             ENDIF 
    782             IF ( dta%ll_v3d .or. ( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) ) THEN 
    783                if(lwp) write(numout,*) '++++++ dta%v3d pointing to fnow' 
    784                jfld = jfld + 1 
    785                dta_bdy(jbdy)%v3d => bf(jfld)%fnow(:,1,:) 
    786             ENDIF 
    787          ENDIF 
    788  
    789          IF( nn_tra_dta(jbdy) == 0 ) THEN 
    790             if(lwp) write(numout,*) '++++++ dta%tem/sal allocated space' 
    791             IF( dta%ll_tem ) ALLOCATE( dta_bdy(jbdy)%tem(nblen(1),jpk) ) 
    792             IF( dta%ll_sal ) ALLOCATE( dta_bdy(jbdy)%sal(nblen(1),jpk) ) 
    793          ELSE 
    794             IF( dta%ll_tem ) THEN 
    795                if(lwp) write(numout,*) '++++++ dta%tem pointing to fnow' 
    796                jfld = jfld + 1 
    797                dta_bdy(jbdy)%tem => bf(jfld)%fnow(:,1,:) 
    798             ENDIF 
    799             IF( dta%ll_sal ) THEN  
    800                if(lwp) write(numout,*) '++++++ dta%sal pointing to fnow' 
    801                jfld = jfld + 1 
    802                dta_bdy(jbdy)%sal => bf(jfld)%fnow(:,1,:) 
    803             ENDIF 
    804          ENDIF 
    805  
    806 #if defined key_si3 
    807          IF (cn_ice(jbdy) /= 'none') THEN 
    808             IF( nn_ice_dta(jbdy) == 0 ) THEN 
    809                ALLOCATE( dta_bdy(jbdy)%a_i(nblen(1),jpl) ) 
    810                ALLOCATE( dta_bdy(jbdy)%h_i(nblen(1),jpl) ) 
    811                ALLOCATE( dta_bdy(jbdy)%h_s(nblen(1),jpl) ) 
    812             ELSE 
    813                IF ( nice_cat == jpl ) THEN ! case input cat = jpl 
    814                   jfld = jfld + 1 
    815                   dta_bdy(jbdy)%a_i => bf(jfld)%fnow(:,1,:) 
    816                   jfld = jfld + 1 
    817                   dta_bdy(jbdy)%h_i => bf(jfld)%fnow(:,1,:) 
    818                   jfld = jfld + 1 
    819                   dta_bdy(jbdy)%h_s => bf(jfld)%fnow(:,1,:) 
    820                ELSE                        ! case input cat = 1 OR (/=1 and /=jpl) 
    821                   jfld_ait(jbdy)  = jfld + 1 
    822                   jfld_htit(jbdy) = jfld + 2 
    823                   jfld_htst(jbdy) = jfld + 3 
    824                   jfld     = jfld + 3 
    825                   ALLOCATE( dta_bdy(jbdy)%a_i(nblen(1),jpl) ) 
    826                   ALLOCATE( dta_bdy(jbdy)%h_i(nblen(1),jpl) ) 
    827                   ALLOCATE( dta_bdy(jbdy)%h_s(nblen(1),jpl) ) 
    828                   dta_bdy(jbdy)%a_i(:,:) = 0._wp 
    829                   dta_bdy(jbdy)%h_i(:,:) = 0._wp 
    830                   dta_bdy(jbdy)%h_s(:,:) = 0._wp 
    831                ENDIF 
    832  
    833             ENDIF 
    834          ENDIF 
    835 #endif 
     517               IF( llread ) THEN                                           ! get data from NetCDF file 
     518                  CALL fld_fill( bf_alias, bn_alias, cn_dir, 'bdy_dta', cl3//' '//ctmp1, ctmp2 )   ! use namelist info 
     519                  IF( bf_alias(1)%ln_tint ) ALLOCATE( bf_alias(1)%fdta( iszdim, 1, ipk, 2 ) ) 
     520                  bf_alias(1)%imap => idx_bdy(jbdy)%nbmap(1:iszdim,igrd)      ! associate the mapping used for this bdy 
     521                  bf_alias(1)%igrd = igrd                                     ! used only for vertical integration of 3D arrays 
     522                  bf_alias(1)%ltotvel = ln_full_vel                           ! T if u3d is full velocity 
     523               ENDIF 
     524 
     525               ! associate the pointer and get rid of the dimensions with a size equal to 1 
     526               IF( jfld == jp_bdyssh           ) dta_bdy(jbdy)%ssh => bf_alias(1)%fnow(:,1,1) 
     527               IF( jfld == jp_bdyu2d           ) dta_bdy(jbdy)%u2d => bf_alias(1)%fnow(:,1,1) 
     528               IF( jfld == jp_bdyv2d           ) dta_bdy(jbdy)%v2d => bf_alias(1)%fnow(:,1,1) 
     529               IF( jfld == jp_bdyu3d           ) dta_bdy(jbdy)%u3d => bf_alias(1)%fnow(:,1,:) 
     530               IF( jfld == jp_bdyv3d           ) dta_bdy(jbdy)%v3d => bf_alias(1)%fnow(:,1,:) 
     531               IF( jfld == jp_bdytem           ) dta_bdy(jbdy)%tem => bf_alias(1)%fnow(:,1,:) 
     532               IF( jfld == jp_bdysal           ) dta_bdy(jbdy)%sal => bf_alias(1)%fnow(:,1,:) 
     533               IF( jfld == jp_bdya_i ) THEN 
     534                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%a_i => bf_alias(1)%fnow(:,1,:) 
     535                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%a_i(iszdim,jpl) ) 
     536                  ENDIF 
     537               ENDIF 
     538               IF( jfld == jp_bdyh_i ) THEN 
     539                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%h_i => bf_alias(1)%fnow(:,1,:) 
     540                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%h_i(iszdim,jpl) ) 
     541                  ENDIF 
     542               ENDIF 
     543               IF( jfld == jp_bdyh_s ) THEN 
     544                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%h_s => bf_alias(1)%fnow(:,1,:) 
     545                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%h_s(iszdim,jpl) ) 
     546                  ENDIF 
     547               ENDIF 
     548            ENDIF 
     549 
     550         END DO   ! jpbdyfld 
    836551         ! 
    837552      END DO ! jbdy  
    838553      ! 
    839554   END SUBROUTINE bdy_dta_init 
    840  
     555    
    841556   !!============================================================================== 
    842557END MODULE bdydta 

  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90

    r11191 r11223  
    3737 
    3838   INTEGER, PARAMETER ::   jp_nseg = 100   !  
    39    INTEGER, PARAMETER ::   nrimmax =  20   ! maximum rimwidth in structured 
    40                                                ! open boundary data files 
    4139   ! Straight open boundary segment parameters: 
    4240   INTEGER  ::   nbdysege, nbdysegw, nbdysegn, nbdysegs  
    … …  
    7068         &             cn_ice, nn_ice_dta,                                     & 
    7169         &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     & 
    72          &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy 
     70         &             ln_vol, nn_volctl, nn_rimwidth 
    7371         ! 
    7472      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    … …  
    8179      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 901) 
    8280901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp ) 
     81      ! make sur that all elements of the namelist variables have a default definition from namelist_ref 
     82      ln_coords_file (2:jp_bdy) = ln_coords_file (1) 
     83      cn_coords_file (2:jp_bdy) = cn_coords_file (1) 
     84      cn_dyn2d       (2:jp_bdy) = cn_dyn2d       (1) 
     85      nn_dyn2d_dta   (2:jp_bdy) = nn_dyn2d_dta   (1) 
     86      cn_dyn3d       (2:jp_bdy) = cn_dyn3d       (1) 
     87      nn_dyn3d_dta   (2:jp_bdy) = nn_dyn3d_dta   (1) 
     88      cn_tra         (2:jp_bdy) = cn_tra         (1) 
     89      nn_tra_dta     (2:jp_bdy) = nn_tra_dta     (1)     
     90      ln_tra_dmp     (2:jp_bdy) = ln_tra_dmp     (1) 
     91      ln_dyn3d_dmp   (2:jp_bdy) = ln_dyn3d_dmp   (1) 
     92      rn_time_dmp    (2:jp_bdy) = rn_time_dmp    (1) 
     93      rn_time_dmp_out(2:jp_bdy) = rn_time_dmp_out(1) 
     94      cn_ice         (2:jp_bdy) = cn_ice         (1) 
     95      nn_ice_dta     (2:jp_bdy) = nn_ice_dta     (1) 
     96      rn_ice_tem     (2:jp_bdy) = rn_ice_tem     (1) 
     97      rn_ice_sal     (2:jp_bdy) = rn_ice_sal     (1) 
     98      rn_ice_age     (2:jp_bdy) = rn_ice_age     (1) 
    8399      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries 
    84100      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 ) 
    … …  
    87103 
    88104      IF( .NOT. Agrif_Root() ) ln_bdy = .FALSE.   ! forced for Agrif children 
     105 
     106      IF( nb_bdy == 0 ) ln_bdy = .FALSE. 
    89107       
    90108      ! ----------------------------------------- 
    … …  
    97115         ! 
    98116         ! Open boundaries definition (arrays and masks) 
    99          CALL bdy_segs 
     117         CALL bdy_def 
    100118         ! 
    101119         ! Open boundaries initialisation of external data arrays 
    … …  
    115133 
    116134 
    117    SUBROUTINE bdy_segs 
     135   SUBROUTINE bdy_def 
    118136      !!---------------------------------------------------------------------- 
    119137      !!                 ***  ROUTINE bdy_init  *** 
    … …  
    130148      INTEGER  ::   ilen1, ibm1                            !   -       - 
    131149      INTEGER  ::   iwe, ies, iso, ino, inum, id_dummy     !   -       - 
    132       INTEGER  ::   jpbdta, jpbdtau, jpbdtas               !   -       - 
     150      INTEGER  ::   jpbdta                                 !   -       - 
    133151      INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       - 
    134152      INTEGER  ::   ii1, ii2, ii3, ij1, ij2, ij3           !   -       - 
    135153      INTEGER  ::   iibe, ijbe, iibi, ijbi                 !   -       - 
    136154      INTEGER  ::   flagu, flagv                           ! short cuts 
     155      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zz_read      ! work space for 2D global boundary data 
    137156      REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask   ! pointer to 2D mask fields 
    138       REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars 
    139       INTEGER, DIMENSION (2)                  ::   kdimsz 
     157      INTEGER, DIMENSION (4)                  ::   kdimsz 
    140158      INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays  
    141159      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta 
    142160      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points 
    143161      CHARACTER(LEN=1),DIMENSION(jpbgrd)      ::   cgrid 
    144       INTEGER :: jpk_max 
    145162      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                  ! temporary fmask array excluding coastal boundary condition (shlat) 
    146163      REAL(wp), DIMENSION(jpi,jpj) ::   ztmask, zumask, zvmask  ! temporary u/v mask array 
    147164      !! 
    148       CHARACTER(LEN=1)                     ::   ctypebdy   !     -        -  
    149165      INTEGER                              ::   nbdyind, nbdybeg, nbdyend 
    150       !! 
    151       NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend 
    152       INTEGER  ::   ios                 ! Local integer output status for namelist read 
    153166      !!---------------------------------------------------------------------- 
    154167      ! 
    … …  
    161174         &                               ' and general open boundary condition are not compatible' ) 
    162175 
    163       IF( nb_bdy == 0 ) THEN  
    164         IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.' 
    165       ELSE 
    166         IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy 
     176      IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy 
     177 
     178      DO ib_bdy = 1,nb_bdy 
     179 
     180         IF(lwp) THEN 
     181            WRITE(numout,*) ' '  
     182            WRITE(numout,*) '------ Open boundary data set ',ib_bdy,' ------'  
     183            IF( ln_coords_file(ib_bdy) ) THEN 
     184               WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy)) 
     185            ELSE 
     186               WRITE(numout,*) 'Boundary defined in namelist.' 
     187            ENDIF 
     188            WRITE(numout,*) 
     189         ENDIF 
     190 
     191         ! barotropic bdy 
     192         !---------------- 
     193         IF(lwp) THEN 
     194            WRITE(numout,*) 'Boundary conditions for barotropic solution:  ' 
     195            SELECT CASE( cn_dyn2d(ib_bdy) )                   
     196            CASE( 'none' )           ;   WRITE(numout,*) '      no open boundary condition'         
     197            CASE( 'frs' )            ;   WRITE(numout,*) '      Flow Relaxation Scheme' 
     198            CASE( 'flather' )        ;   WRITE(numout,*) '      Flather radiation condition' 
     199            CASE( 'orlanski' )       ;   WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging' 
     200            CASE( 'orlanski_npo' )   ;   WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging' 
     201            CASE DEFAULT             ;   CALL ctl_stop( 'unrecognised value for cn_dyn2d' ) 
     202            END SELECT 
     203         ENDIF 
     204 
     205         dta_bdy(ib_bdy)%lneed_ssh   = cn_dyn2d(ib_bdy) == 'flather' 
     206         dta_bdy(ib_bdy)%lneed_dyn2d = cn_dyn2d(ib_bdy) /= 'none' 
     207 
     208         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn2d ) THEN 
     209            SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !  
     210            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'         
     211            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file' 
     212            CASE( 2 )      ;   WRITE(numout,*) '      tidal harmonic forcing taken from file' 
     213            CASE( 3 )      ;   WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files' 
     214            CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' ) 
     215            END SELECT 
     216         ENDIF 
     217         IF ( dta_bdy(ib_bdy)%lneed_dyn2d .AND. nn_dyn2d_dta(ib_bdy) .GE. 2  .AND. .NOT.ln_tide ) THEN 
     218            CALL ctl_stop( 'You must activate with ln_tide to add tidal forcing at open boundaries' ) 
     219         ENDIF 
     220         IF(lwp) WRITE(numout,*) 
     221 
     222         ! baroclinic bdy 
     223         !---------------- 
     224         IF(lwp) THEN 
     225            WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  ' 
     226            SELECT CASE( cn_dyn3d(ib_bdy) )                   
     227            CASE('none')           ;   WRITE(numout,*) '      no open boundary condition'         
     228            CASE('frs')            ;   WRITE(numout,*) '      Flow Relaxation Scheme' 
     229            CASE('specified')      ;   WRITE(numout,*) '      Specified value' 
     230            CASE('neumann')        ;   WRITE(numout,*) '      Neumann conditions' 
     231            CASE('zerograd')       ;   WRITE(numout,*) '      Zero gradient for baroclinic velocities' 
     232            CASE('zero')           ;   WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
     233            CASE('orlanski')       ;   WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging' 
     234            CASE('orlanski_npo')   ;   WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging' 
     235            CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn3d' ) 
     236            END SELECT 
     237         ENDIF 
     238 
     239         dta_bdy(ib_bdy)%lneed_dyn3d = cn_dyn3d(ib_bdy) == 'frs'      .OR. cn_dyn3d(ib_bdy) == 'specified'   & 
     240            &                     .OR. cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' 
     241 
     242         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn3d ) THEN 
     243            SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   !  
     244            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'         
     245            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file' 
     246            CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' ) 
     247            END SELECT 
     248         END IF 
     249 
     250         IF ( ln_dyn3d_dmp(ib_bdy) ) THEN 
     251            IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN 
     252               IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.' 
     253               ln_dyn3d_dmp(ib_bdy) = .false. 
     254            ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN 
     255               CALL ctl_stop( 'Use FRS OR relaxation' ) 
     256            ELSE 
     257               IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone' 
     258               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
     259               IF(rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 
     260               dta_bdy(ib_bdy)%lneed_dyn3d = .TRUE. 
     261            ENDIF 
     262         ELSE 
     263            IF(lwp) WRITE(numout,*) '      NO relaxation on baroclinic velocities' 
     264         ENDIF 
     265         IF(lwp) WRITE(numout,*) 
     266 
     267         !    tra bdy 
     268         !---------------- 
     269         IF(lwp) THEN 
     270            WRITE(numout,*) 'Boundary conditions for temperature and salinity:  ' 
     271            SELECT CASE( cn_tra(ib_bdy) )                   
     272            CASE('none')           ;   WRITE(numout,*) '      no open boundary condition'         
     273            CASE('frs')            ;   WRITE(numout,*) '      Flow Relaxation Scheme' 
     274            CASE('specified')      ;   WRITE(numout,*) '      Specified value' 
     275            CASE('neumann')        ;   WRITE(numout,*) '      Neumann conditions' 
     276            CASE('runoff')         ;   WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity' 
     277            CASE('orlanski')       ;   WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging' 
     278            CASE('orlanski_npo')   ;   WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging' 
     279            CASE DEFAULT           ;   CALL ctl_stop( 'unrecognised value for cn_tra' ) 
     280            END SELECT 
     281         ENDIF 
     282 
     283         dta_bdy(ib_bdy)%lneed_tra = cn_tra(ib_bdy) == 'frs'       .OR. cn_tra(ib_bdy) == 'specified'   & 
     284            &                   .OR. cn_tra(ib_bdy) == 'orlanski'  .OR. cn_tra(ib_bdy) == 'orlanski_npo'  
     285 
     286         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_tra ) THEN 
     287            SELECT CASE( nn_tra_dta(ib_bdy) )                   !  
     288            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'         
     289            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file' 
     290            CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' ) 
     291            END SELECT 
     292         ENDIF 
     293 
     294         IF ( ln_tra_dmp(ib_bdy) ) THEN 
     295            IF ( cn_tra(ib_bdy) == 'none' ) THEN 
     296               IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.' 
     297               ln_tra_dmp(ib_bdy) = .false. 
     298            ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN 
     299               CALL ctl_stop( 'Use FRS OR relaxation' ) 
     300            ELSE 
     301               IF(lwp) WRITE(numout,*) '      + T/S relaxation zone' 
     302               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
     303               IF(lwp) WRITE(numout,*) '      Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days' 
     304               IF(lwp.AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 
     305               dta_bdy(ib_bdy)%lneed_tra = .TRUE. 
     306            ENDIF 
     307         ELSE 
     308            IF(lwp) WRITE(numout,*) '      NO T/S relaxation' 
     309         ENDIF 
     310         IF(lwp) WRITE(numout,*) 
     311 
     312#if defined key_si3 
     313         IF(lwp) THEN 
     314            WRITE(numout,*) 'Boundary conditions for sea ice:  ' 
     315            SELECT CASE( cn_ice(ib_bdy) )                   
     316            CASE('none')   ;   WRITE(numout,*) '      no open boundary condition'         
     317            CASE('frs')    ;   WRITE(numout,*) '      Flow Relaxation Scheme' 
     318            CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice' ) 
     319            END SELECT 
     320         ENDIF 
     321 
     322         dta_bdy(ib_bdy)%lneed_ice = cn_ice(ib_bdy) /= 'none' 
     323 
     324         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_ice ) THEN  
     325            SELECT CASE( nn_ice_dta(ib_bdy) )                   !  
     326            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'         
     327            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file' 
     328            CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_dta must be 0 or 1' ) 
     329            END SELECT 
     330            WRITE(numout,*) 
     331            WRITE(numout,*) '      tem of bdy sea-ice = ', rn_ice_tem(ib_bdy)          
     332            WRITE(numout,*) '      sal of bdy sea-ice = ', rn_ice_sal(ib_bdy)          
     333            WRITE(numout,*) '      age of bdy sea-ice = ', rn_ice_age(ib_bdy)          
     334         ENDIF 
     335#else 
     336         dta_bdy(ib_bdy)%lneed_ice = .FALSE. 
     337#endif 
     338         ! 
     339         IF(lwp) WRITE(numout,*) 
     340         IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy) 
     341         IF(lwp) WRITE(numout,*) 
     342         ! 
     343      END DO   ! nb_bdy 
     344 
     345      IF( lwp ) THEN 
     346         IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value) 
     347            WRITE(numout,*) 'Volume correction applied at open boundaries' 
     348            WRITE(numout,*) 
     349            SELECT CASE ( nn_volctl ) 
     350            CASE( 1 )      ;   WRITE(numout,*) '      The total volume will be constant' 
     351            CASE( 0 )      ;   WRITE(numout,*) '      The total volume will vary according to the surface E-P flux' 
     352            CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 
     353            END SELECT 
     354            WRITE(numout,*) 
     355            ! 
     356            ! sanity check if used with tides         
     357            IF( ln_tide ) THEN  
     358               WRITE(numout,*) ' The total volume correction is not working with tides. ' 
     359               WRITE(numout,*) ' Set ln_vol to .FALSE. ' 
     360               WRITE(numout,*) ' or ' 
     361               WRITE(numout,*) ' equilibriate your bdy input files ' 
     362               CALL ctl_stop( 'The total volume correction is not working with tides.' ) 
     363            END IF 
     364         ELSE 
     365            WRITE(numout,*) 'No volume correction applied at open boundaries' 
     366            WRITE(numout,*) 
     367         ENDIF 
    167368      ENDIF 
    168  
    169       DO ib_bdy = 1,nb_bdy 
    170         IF(lwp) WRITE(numout,*) ' '  
    171         IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------'  
    172  
    173         IF( ln_coords_file(ib_bdy) ) THEN 
    174            IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy)) 
    175         ELSE 
    176            IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.' 
    177         ENDIF 
    178         IF(lwp) WRITE(numout,*) 
    179  
    180         IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  ' 
    181         SELECT CASE( cn_dyn2d(ib_bdy) )                   
    182           CASE( 'none' )          
    183              IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    184              dta_bdy(ib_bdy)%ll_ssh = .false. 
    185              dta_bdy(ib_bdy)%ll_u2d = .false. 
    186              dta_bdy(ib_bdy)%ll_v2d = .false. 
    187           CASE( 'frs' )           
    188              IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    189              dta_bdy(ib_bdy)%ll_ssh = .false. 
    190              dta_bdy(ib_bdy)%ll_u2d = .true. 
    191              dta_bdy(ib_bdy)%ll_v2d = .true. 
    192           CASE( 'flather' )       
    193              IF(lwp) WRITE(numout,*) '      Flather radiation condition' 
    194              dta_bdy(ib_bdy)%ll_ssh = .true. 
    195              dta_bdy(ib_bdy)%ll_u2d = .true. 
    196              dta_bdy(ib_bdy)%ll_v2d = .true. 
    197           CASE( 'orlanski' )      
    198              IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging' 
    199              dta_bdy(ib_bdy)%ll_ssh = .false. 
    200              dta_bdy(ib_bdy)%ll_u2d = .true. 
    201              dta_bdy(ib_bdy)%ll_v2d = .true. 
    202           CASE( 'orlanski_npo' )  
    203              IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging' 
    204              dta_bdy(ib_bdy)%ll_ssh = .false. 
    205              dta_bdy(ib_bdy)%ll_u2d = .true. 
    206              dta_bdy(ib_bdy)%ll_v2d = .true. 
    207           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn2d' ) 
    208         END SELECT 
    209         IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN 
    210            SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !  
    211               CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
    212               CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
    213               CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      tidal harmonic forcing taken from file' 
    214               CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files' 
    215               CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' ) 
    216            END SELECT 
    217            IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(.NOT.ln_tide)) THEN 
    218              CALL ctl_stop( 'You must activate with ln_tide to add tidal forcing at open boundaries' ) 
    219            ENDIF 
    220         ENDIF 
    221         IF(lwp) WRITE(numout,*) 
    222  
    223         IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  ' 
    224         SELECT CASE( cn_dyn3d(ib_bdy) )                   
    225           CASE('none') 
    226              IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    227              dta_bdy(ib_bdy)%ll_u3d = .false. 
    228              dta_bdy(ib_bdy)%ll_v3d = .false. 
    229           CASE('frs')        
    230              IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    231              dta_bdy(ib_bdy)%ll_u3d = .true. 
    232              dta_bdy(ib_bdy)%ll_v3d = .true. 
    233           CASE('specified') 
    234              IF(lwp) WRITE(numout,*) '      Specified value' 
    235              dta_bdy(ib_bdy)%ll_u3d = .true. 
    236              dta_bdy(ib_bdy)%ll_v3d = .true. 
    237           CASE('neumann') 
    238              IF(lwp) WRITE(numout,*) '      Neumann conditions' 
    239              dta_bdy(ib_bdy)%ll_u3d = .false. 
    240              dta_bdy(ib_bdy)%ll_v3d = .false. 
    241           CASE('zerograd') 
    242              IF(lwp) WRITE(numout,*) '      Zero gradient for baroclinic velocities' 
    243              dta_bdy(ib_bdy)%ll_u3d = .false. 
    244              dta_bdy(ib_bdy)%ll_v3d = .false. 
    245           CASE('zero') 
    246              IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
    247              dta_bdy(ib_bdy)%ll_u3d = .false. 
    248              dta_bdy(ib_bdy)%ll_v3d = .false. 
    249           CASE('orlanski') 
    250              IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging' 
    251              dta_bdy(ib_bdy)%ll_u3d = .true. 
    252              dta_bdy(ib_bdy)%ll_v3d = .true. 
    253           CASE('orlanski_npo') 
    254              IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging' 
    255              dta_bdy(ib_bdy)%ll_u3d = .true. 
    256              dta_bdy(ib_bdy)%ll_v3d = .true. 
    257           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn3d' ) 
    258         END SELECT 
    259         IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 
    260            SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   !  
    261               CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
    262               CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
    263               CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' ) 
    264            END SELECT 
    265         ENDIF 
    266  
    267         IF ( ln_dyn3d_dmp(ib_bdy) ) THEN 
    268            IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN 
    269               IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.' 
    270               ln_dyn3d_dmp(ib_bdy)=.false. 
    271            ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN 
    272               CALL ctl_stop( 'Use FRS OR relaxation' ) 
    273            ELSE 
    274               IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone' 
    275               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
    276               IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 
    277               dta_bdy(ib_bdy)%ll_u3d = .true. 
    278               dta_bdy(ib_bdy)%ll_v3d = .true. 
    279            ENDIF 
    280         ELSE 
    281            IF(lwp) WRITE(numout,*) '      NO relaxation on baroclinic velocities' 
    282         ENDIF 
    283         IF(lwp) WRITE(numout,*) 
    284  
    285         IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  ' 
    286         SELECT CASE( cn_tra(ib_bdy) )                   
    287           CASE('none') 
    288              IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    289              dta_bdy(ib_bdy)%ll_tem = .false. 
    290              dta_bdy(ib_bdy)%ll_sal = .false. 
    291           CASE('frs') 
    292              IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    293              dta_bdy(ib_bdy)%ll_tem = .true. 
    294              dta_bdy(ib_bdy)%ll_sal = .true. 
    295           CASE('specified') 
    296              IF(lwp) WRITE(numout,*) '      Specified value' 
    297              dta_bdy(ib_bdy)%ll_tem = .true. 
    298              dta_bdy(ib_bdy)%ll_sal = .true. 
    299           CASE('neumann') 
    300              IF(lwp) WRITE(numout,*) '      Neumann conditions' 
    301              dta_bdy(ib_bdy)%ll_tem = .false. 
    302              dta_bdy(ib_bdy)%ll_sal = .false. 
    303           CASE('runoff') 
    304              IF(lwp) WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity' 
    305              dta_bdy(ib_bdy)%ll_tem = .false. 
    306              dta_bdy(ib_bdy)%ll_sal = .false. 
    307           CASE('orlanski') 
    308              IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging' 
    309              dta_bdy(ib_bdy)%ll_tem = .true. 
    310              dta_bdy(ib_bdy)%ll_sal = .true. 
    311           CASE('orlanski_npo') 
    312              IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging' 
    313              dta_bdy(ib_bdy)%ll_tem = .true. 
    314              dta_bdy(ib_bdy)%ll_sal = .true. 
    315           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_tra' ) 
    316         END SELECT 
    317         IF( cn_tra(ib_bdy) /= 'none' ) THEN 
    318            SELECT CASE( nn_tra_dta(ib_bdy) )                   !  
    319               CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
    320               CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
    321               CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' ) 
    322            END SELECT 
    323         ENDIF 
    324  
    325         IF ( ln_tra_dmp(ib_bdy) ) THEN 
    326            IF ( cn_tra(ib_bdy) == 'none' ) THEN 
    327               IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.' 
    328               ln_tra_dmp(ib_bdy)=.false. 
    329            ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN 
    330               CALL ctl_stop( 'Use FRS OR relaxation' ) 
    331            ELSE 
    332               IF(lwp) WRITE(numout,*) '      + T/S relaxation zone' 
    333               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
    334               IF(lwp) WRITE(numout,*) '      Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days' 
    335               IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 
    336               dta_bdy(ib_bdy)%ll_tem = .true. 
    337               dta_bdy(ib_bdy)%ll_sal = .true. 
    338            ENDIF 
    339         ELSE 
    340            IF(lwp) WRITE(numout,*) '      NO T/S relaxation' 
    341         ENDIF 
    342         IF(lwp) WRITE(numout,*) 
    343  
    344 #if defined key_si3 
    345          IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  ' 
    346          SELECT CASE( cn_ice(ib_bdy) )                   
    347          CASE('none') 
    348              IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    349              dta_bdy(ib_bdy)%ll_a_i = .false. 
    350              dta_bdy(ib_bdy)%ll_h_i = .false. 
    351              dta_bdy(ib_bdy)%ll_h_s = .false. 
    352          CASE('frs') 
    353              IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    354              dta_bdy(ib_bdy)%ll_a_i = .true. 
    355              dta_bdy(ib_bdy)%ll_h_i = .true. 
    356              dta_bdy(ib_bdy)%ll_h_s = .true. 
    357          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice' ) 
    358          END SELECT 
    359         IF( cn_ice(ib_bdy) /= 'none' ) THEN  
    360            SELECT CASE( nn_ice_dta(ib_bdy) )                   !  
    361               CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
    362               CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
    363               CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_dta must be 0 or 1' ) 
    364            END SELECT 
    365         ENDIF 
    366         IF(lwp) WRITE(numout,*) 
    367         IF(lwp) WRITE(numout,*) '      tem of bdy sea-ice = ', rn_ice_tem(ib_bdy)          
    368         IF(lwp) WRITE(numout,*) '      sal of bdy sea-ice = ', rn_ice_sal(ib_bdy)          
    369         IF(lwp) WRITE(numout,*) '      age of bdy sea-ice = ', rn_ice_age(ib_bdy)          
    370 #endif 
    371  
    372         IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy) 
    373         IF(lwp) WRITE(numout,*) 
    374          ! 
    375       END DO 
    376  
    377      IF( nb_bdy > 0 ) THEN 
    378         IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value) 
    379           IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 
    380           IF(lwp) WRITE(numout,*) 
    381           SELECT CASE ( nn_volctl ) 
    382             CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant' 
    383             CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux' 
    384             CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 
    385           END SELECT 
    386           IF(lwp) WRITE(numout,*) 
    387           ! 
    388           ! sanity check if used with tides         
    389           IF( ln_tide ) THEN  
    390              IF(lwp) WRITE(numout,*) ' The total volume correction is not working with tides. ' 
    391              IF(lwp) WRITE(numout,*) ' Set ln_vol to .FALSE. ' 
    392              IF(lwp) WRITE(numout,*) ' or ' 
    393              IF(lwp) WRITE(numout,*) ' equilibriate your bdy input files ' 
    394              CALL ctl_stop( 'The total volume correction is not working with tides.' ) 
    395           END IF 
    396         ELSE 
    397           IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 
    398           IF(lwp) WRITE(numout,*) 
    399         ENDIF 
    400         IF( nb_jpk_bdy(ib_bdy) > 0 ) THEN 
    401            IF(lwp) WRITE(numout,*) '*** open boundary will be interpolate in the vertical onto the native grid ***' 
    402         ELSE 
    403            IF(lwp) WRITE(numout,*) '*** open boundary will be read straight onto the native grid without vertical interpolation ***' 
    404         ENDIF 
    405      ENDIF 
    406369 
    407370      ! ------------------------------------------------- 
    … …  
    409372      ! ------------------------------------------------- 
    410373 
    411       ! Work out global dimensions of boundary data 
    412       ! --------------------------------------------- 
    413374      REWIND( numnam_cfg )      
    414  
    415375      nblendta(:,:) = 0 
    416376      nbdysege = 0 
    … …  
    418378      nbdysegn = 0 
    419379      nbdysegs = 0 
    420       icount   = 0 ! count user defined segments 
    421       ! Dimensions below are used to allocate arrays to read external data 
    422       jpbdtas = 1 ! Maximum size of boundary data (structured case) 
    423       jpbdtau = 1 ! Maximum size of boundary data (unstructured case) 
    424  
     380 
     381      ! Define all boundaries  
     382      ! --------------------- 
    425383      DO ib_bdy = 1, nb_bdy 
    426  
    427          IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters 
    428   
    429             icount = icount + 1 
    430             ! No REWIND here because may need to read more than one nambdy_index namelist. 
    431             ! Read only namelist_cfg to avoid unseccessfull overwrite  
    432             ! keep full control of the configuration namelist 
    433             READ  ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 ) 
    434 904         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_index in configuration namelist', lwp ) 
    435             IF(lwm) WRITE ( numond, nambdy_index ) 
    436  
    437             SELECT CASE ( TRIM(ctypebdy) ) 
    438               CASE( 'N' ) 
    439                  IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
    440                     nbdyind  = jpjglo - 2  ! set boundary to whole side of model domain. 
    441                     nbdybeg  = 2 
    442                     nbdyend  = jpiglo - 1 
    443                  ENDIF 
    444                  nbdysegn = nbdysegn + 1 
    445                  npckgn(nbdysegn) = ib_bdy ! Save bdy package number 
    446                  jpjnob(nbdysegn) = nbdyind 
    447                  jpindt(nbdysegn) = nbdybeg 
    448                  jpinft(nbdysegn) = nbdyend 
    449                  ! 
    450               CASE( 'S' ) 
    451                  IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
    452                     nbdyind  = 2           ! set boundary to whole side of model domain. 
    453                     nbdybeg  = 2 
    454                     nbdyend  = jpiglo - 1 
    455                  ENDIF 
    456                  nbdysegs = nbdysegs + 1 
    457                  npckgs(nbdysegs) = ib_bdy ! Save bdy package number 
    458                  jpjsob(nbdysegs) = nbdyind 
    459                  jpisdt(nbdysegs) = nbdybeg 
    460                  jpisft(nbdysegs) = nbdyend 
    461                  ! 
    462               CASE( 'E' ) 
    463                  IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
    464                     nbdyind  = jpiglo - 2  ! set boundary to whole side of model domain. 
    465                     nbdybeg  = 2 
    466                     nbdyend  = jpjglo - 1 
    467                  ENDIF 
    468                  nbdysege = nbdysege + 1  
    469                  npckge(nbdysege) = ib_bdy ! Save bdy package number 
    470                  jpieob(nbdysege) = nbdyind 
    471                  jpjedt(nbdysege) = nbdybeg 
    472                  jpjeft(nbdysege) = nbdyend 
    473                  ! 
    474               CASE( 'W' ) 
    475                  IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
    476                     nbdyind  = 2           ! set boundary to whole side of model domain. 
    477                     nbdybeg  = 2 
    478                     nbdyend  = jpjglo - 1 
    479                  ENDIF 
    480                  nbdysegw = nbdysegw + 1 
    481                  npckgw(nbdysegw) = ib_bdy ! Save bdy package number 
    482                  jpiwob(nbdysegw) = nbdyind 
    483                  jpjwdt(nbdysegw) = nbdybeg 
    484                  jpjwft(nbdysegw) = nbdyend 
    485                  ! 
    486               CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' ) 
    487             END SELECT 
    488  
    489             ! For simplicity we assume that in case of straight bdy, arrays have the same length 
    490             ! (even if it is true that last tangential velocity points 
    491             ! are useless). This simplifies a little bit boundary data format (and agrees with format 
    492             ! used so far in obc package) 
    493  
    494             nblendta(1:jpbgrd,ib_bdy) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(ib_bdy) 
    495             jpbdtas = MAX(jpbdtas, (nbdyend - nbdybeg + 1)) 
    496             IF (lwp.and.(nn_rimwidth(ib_bdy)>nrimmax)) & 
    497             & CALL ctl_stop( 'rimwidth must be lower than nrimmax' ) 
    498  
    499          ELSE            ! Read size of arrays in boundary coordinates file. 
     384         ! 
     385         IF( .NOT. ln_coords_file(ib_bdy) ) THEN     ! build bdy coordinates with segments defined in namelist 
     386 
     387            CALL bdy_read_seg( ib_bdy, nblendta(:,ib_bdy) ) 
     388            ! Now look for crossings in user (namelist) defined open boundary segments: 
     389            IF( nbdysege > 0 .OR. nbdysegw > 0 .OR. nbdysegn > 0 .OR. nbdysegs > 0)   CALL bdy_ctl_seg 
     390 
     391         ELSE                                        ! Read size of arrays in boundary coordinates file. 
     392             
    500393            CALL iom_open( cn_coords_file(ib_bdy), inum ) 
    501394            DO igrd = 1, jpbgrd 
    502395               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )   
    503396               nblendta(igrd,ib_bdy) = MAXVAL(kdimsz) 
    504                jpbdtau = MAX(jpbdtau, MAXVAL(kdimsz)) 
    505397            END DO 
    506398            CALL iom_close( inum ) 
    507             ! 
    508          ENDIF  
     399         ENDIF 
    509400         ! 
    510401      END DO ! ib_bdy 
    511  
    512       IF (nb_bdy>0) THEN 
    513          jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy)) 
    514  
    515          ! Allocate arrays 
    516          !--------------- 
    517          ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    & 
    518             &      nbrdta(jpbdta, jpbgrd, nb_bdy) ) 
    519           
    520          jpk_max = MAXVAL(nb_jpk_bdy) 
    521          jpk_max = MAX(jpk_max, jpk) 
    522  
    523          ALLOCATE( dta_global(jpbdtau, 1, jpk_max) ) 
    524          ALLOCATE( dta_global_z(jpbdtau, 1, jpk_max) ) ! needed ?? TODO 
    525          ALLOCATE( dta_global_dz(jpbdtau, 1, jpk_max) )! needed ?? TODO 
    526  
    527          IF ( icount>0 ) THEN 
    528             ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk_max) ) 
    529             ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk_max) ) ! needed ?? TODO 
    530             ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk_max) )! needed ?? TODO   
    531          ENDIF 
    532          !  
    533       ENDIF 
    534  
    535       ! Now look for crossings in user (namelist) defined open boundary segments: 
    536       !-------------------------------------------------------------------------- 
    537       IF( icount>0 )   CALL bdy_ctl_seg 
    538  
     402       
     403      ! Allocate arrays 
     404      !--------------- 
     405      jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy)) 
     406      ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), nbrdta(jpbdta, jpbgrd, nb_bdy) ) 
     407     
    539408      ! Calculate global boundary index arrays or read in from file 
    540409      !------------------------------------------------------------                
    … …  
    544413         IF( ln_coords_file(ib_bdy) ) THEN 
    545414            ! 
     415            ALLOCATE( zz_read( MAXVAL(nblendta), 1 ) )           
    546416            CALL iom_open( cn_coords_file(ib_bdy), inum ) 
     417            ! 
    547418            DO igrd = 1, jpbgrd 
    548                CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
     419               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) ) 
    549420               DO ii = 1,nblendta(igrd,ib_bdy) 
    550                   nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) ) 
     421                  nbidta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) 
    551422               END DO 
    552                CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
     423               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) ) 
    553424               DO ii = 1,nblendta(igrd,ib_bdy) 
    554                   nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) ) 
     425                  nbjdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) 
    555426               END DO 
    556                CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
     427               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) ) 
    557428               DO ii = 1,nblendta(igrd,ib_bdy) 
    558                   nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) ) 
     429                  nbrdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) 
    559430               END DO 
    560431               ! 
    … …  
    564435               IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy) 
    565436               IF (ibr_max < nn_rimwidth(ib_bdy))   & 
    566                      CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) ) 
    567             END DO 
     437                  CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) ) 
     438            END DO 
     439            ! 
    568440            CALL iom_close( inum ) 
     441            DEALLOCATE( zz_read ) 
    569442            ! 
    570          ENDIF  
    571          ! 
    572       END DO       
    573      
     443         ENDIF 
     444         ! 
     445      END DO 
     446 
    574447      ! 2. Now fill indices corresponding to straight open boundary arrays: 
    575       ! East 
    576       !----- 
    577       DO iseg = 1, nbdysege 
    578          ib_bdy = npckge(iseg) 
    579          ! 
    580          ! ------------ T points ------------- 
    581          igrd=1 
    582          icount=0 
    583          DO ir = 1, nn_rimwidth(ib_bdy) 
    584             DO ij = jpjedt(iseg), jpjeft(iseg) 
    585                icount = icount + 1 
    586                nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 
    587                nbjdta(icount, igrd, ib_bdy) = ij 
    588                nbrdta(icount, igrd, ib_bdy) = ir 
    589             ENDDO 
    590          ENDDO 
    591          ! 
    592          ! ------------ U points ------------- 
    593          igrd=2 
    594          icount=0 
    595          DO ir = 1, nn_rimwidth(ib_bdy) 
    596             DO ij = jpjedt(iseg), jpjeft(iseg) 
    597                icount = icount + 1 
    598                nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir 
    599                nbjdta(icount, igrd, ib_bdy) = ij 
    600                nbrdta(icount, igrd, ib_bdy) = ir 
    601             ENDDO 
    602          ENDDO 
    603          ! 
    604          ! ------------ V points ------------- 
    605          igrd=3 
    606          icount=0 
    607          DO ir = 1, nn_rimwidth(ib_bdy) 
    608 !            DO ij = jpjedt(iseg), jpjeft(iseg) - 1 
    609             DO ij = jpjedt(iseg), jpjeft(iseg) 
    610                icount = icount + 1 
    611                nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 
    612                nbjdta(icount, igrd, ib_bdy) = ij 
    613                nbrdta(icount, igrd, ib_bdy) = ir 
    614             ENDDO 
    615             nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    616             nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    617          ENDDO 
    618       ENDDO 
    619       ! 
    620       ! West 
    621       !----- 
    622       DO iseg = 1, nbdysegw 
    623          ib_bdy = npckgw(iseg) 
    624          ! 
    625          ! ------------ T points ------------- 
    626          igrd=1 
    627          icount=0 
    628          DO ir = 1, nn_rimwidth(ib_bdy) 
    629             DO ij = jpjwdt(iseg), jpjwft(iseg) 
    630                icount = icount + 1 
    631                nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    632                nbjdta(icount, igrd, ib_bdy) = ij 
    633                nbrdta(icount, igrd, ib_bdy) = ir 
    634             ENDDO 
    635          ENDDO 
    636          ! 
    637          ! ------------ U points ------------- 
    638          igrd=2 
    639          icount=0 
    640          DO ir = 1, nn_rimwidth(ib_bdy) 
    641             DO ij = jpjwdt(iseg), jpjwft(iseg) 
    642                icount = icount + 1 
    643                nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    644                nbjdta(icount, igrd, ib_bdy) = ij 
    645                nbrdta(icount, igrd, ib_bdy) = ir 
    646             ENDDO 
    647          ENDDO 
    648          ! 
    649          ! ------------ V points ------------- 
    650          igrd=3 
    651          icount=0 
    652          DO ir = 1, nn_rimwidth(ib_bdy) 
    653 !            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 
    654             DO ij = jpjwdt(iseg), jpjwft(iseg) 
    655                icount = icount + 1 
    656                nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    657                nbjdta(icount, igrd, ib_bdy) = ij 
    658                nbrdta(icount, igrd, ib_bdy) = ir 
    659             ENDDO 
    660             nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    661             nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    662          ENDDO 
    663       ENDDO 
    664       ! 
    665       ! North 
    666       !----- 
    667       DO iseg = 1, nbdysegn 
    668          ib_bdy = npckgn(iseg) 
    669          ! 
    670          ! ------------ T points ------------- 
    671          igrd=1 
    672          icount=0 
    673          DO ir = 1, nn_rimwidth(ib_bdy) 
    674             DO ii = jpindt(iseg), jpinft(iseg) 
    675                icount = icount + 1 
    676                nbidta(icount, igrd, ib_bdy) = ii 
    677                nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir  
    678                nbrdta(icount, igrd, ib_bdy) = ir 
    679             ENDDO 
    680          ENDDO 
    681          ! 
    682          ! ------------ U points ------------- 
    683          igrd=2 
    684          icount=0 
    685          DO ir = 1, nn_rimwidth(ib_bdy) 
    686 !            DO ii = jpindt(iseg), jpinft(iseg) - 1 
    687             DO ii = jpindt(iseg), jpinft(iseg) 
    688                icount = icount + 1 
    689                nbidta(icount, igrd, ib_bdy) = ii 
    690                nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 
    691                nbrdta(icount, igrd, ib_bdy) = ir 
    692             ENDDO 
    693             nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    694             nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    695          ENDDO 
    696          ! 
    697          ! ------------ V points ------------- 
    698          igrd=3 
    699          icount=0 
    700          DO ir = 1, nn_rimwidth(ib_bdy) 
    701             DO ii = jpindt(iseg), jpinft(iseg) 
    702                icount = icount + 1 
    703                nbidta(icount, igrd, ib_bdy) = ii 
    704                nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir 
    705                nbrdta(icount, igrd, ib_bdy) = ir 
    706             ENDDO 
    707          ENDDO 
    708       ENDDO 
    709       ! 
    710       ! South 
    711       !----- 
    712       DO iseg = 1, nbdysegs 
    713          ib_bdy = npckgs(iseg) 
    714          ! 
    715          ! ------------ T points ------------- 
    716          igrd=1 
    717          icount=0 
    718          DO ir = 1, nn_rimwidth(ib_bdy) 
    719             DO ii = jpisdt(iseg), jpisft(iseg) 
    720                icount = icount + 1 
    721                nbidta(icount, igrd, ib_bdy) = ii 
    722                nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    723                nbrdta(icount, igrd, ib_bdy) = ir 
    724             ENDDO 
    725          ENDDO 
    726          ! 
    727          ! ------------ U points ------------- 
    728          igrd=2 
    729          icount=0 
    730          DO ir = 1, nn_rimwidth(ib_bdy) 
    731 !            DO ii = jpisdt(iseg), jpisft(iseg) - 1 
    732             DO ii = jpisdt(iseg), jpisft(iseg) 
    733                icount = icount + 1 
    734                nbidta(icount, igrd, ib_bdy) = ii 
    735                nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    736                nbrdta(icount, igrd, ib_bdy) = ir 
    737             ENDDO 
    738             nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    739             nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    740          ENDDO 
    741          ! 
    742          ! ------------ V points ------------- 
    743          igrd=3 
    744          icount=0 
    745          DO ir = 1, nn_rimwidth(ib_bdy) 
    746             DO ii = jpisdt(iseg), jpisft(iseg) 
    747                icount = icount + 1 
    748                nbidta(icount, igrd, ib_bdy) = ii 
    749                nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    750                nbrdta(icount, igrd, ib_bdy) = ir 
    751             ENDDO 
    752          ENDDO 
    753       ENDDO 
     448      CALL bdy_coords_seg( nbidta, nbjdta, nbrdta ) 
    754449 
    755450      !  Deal with duplicated points 
    … …  
    765460                     DO ib2 = 1, nblendta(igrd,ib_bdy2) 
    766461                        IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. & 
    767                         &   (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN 
    768 !                           IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &  
    769 !                                                       &              nbidta(ib1, igrd, ib_bdy1),      &  
    770 !                                                       &              nbjdta(ib2, igrd, ib_bdy2) 
     462                           &   (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN 
     463                           !                           IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &  
     464                           !                                                       &              nbidta(ib1, igrd, ib_bdy1),      &  
     465                           !                                                       &              nbjdta(ib2, igrd, ib_bdy2) 
    771466                           ! keep only points with the lowest distance to boundary: 
    772467                           IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN 
    773                              nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2 
    774                              nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2 
     468                              nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2 
     469                              nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2 
    775470                           ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN 
    776                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 
    777                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 
    778                            ! Arbitrary choice if distances are the same: 
     471                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 
     472                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 
     473                              ! Arbitrary choice if distances are the same: 
    779474                           ELSE 
    780                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 
    781                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 
     475                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 
     476                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 
    782477                           ENDIF 
    783478                        END IF 
    … …  
    811506                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN 
    812507                     CALL ctl_stop('bdy_segs : ERROR : boundary data in file must be defined ', & 
    813                           &        ' in order of distance from edge nbr A utility for re-ordering ', & 
    814                           &        ' boundary coordinates and data files exists in the TOOLS/OBC directory') 
    815                   ENDIF     
     508                        &        ' in order of distance from edge nbr A utility for re-ordering ', & 
     509                        &        ' boundary coordinates and data files exists in the TOOLS/OBC directory') 
     510                  ENDIF 
    816511               ENDIF 
    817512               ! check if point is in local domain 
    … …  
    911606               ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) )   ! both rim 0 and rim 1 have the same weights 
    912607               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( ir - 1 ) *0.5 )      ! tanh formulation 
    913 !               idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.  ! quadratic 
    914 !               idx_bdy(ib_bdy)%nbw(ib,igrd) =  REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy))       ! linear 
    915             END DO 
    916          END DO  
     608               !               idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.  ! quadratic 
     609               !               idx_bdy(ib_bdy)%nbw(ib,igrd) =  REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy))       ! linear 
     610            END DO 
     611         END DO 
    917612 
    918613         ! Compute damping coefficients 
    … …  
    922617               ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) )   ! both rim 0 and rim 1 have the same damping coefficients 
    923618               idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) &  
    924                & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic 
     619                  & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic 
    925620               idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) &  
    926                & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic 
    927             END DO 
    928          END DO  
     621                  & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic 
     622            END DO 
     623         END DO 
    929624 
    930625      END DO ! ib_bdy 
    … …  
    935630 
    936631      ! ------------------------------------------ 
    937       ! handel rim0, do as if rim 1 was free ocean 
     632      ! handle rim0, do as if rim 1 was free ocean 
    938633      ! ------------------------------------------ 
    939634 
    … …  
    944639         DO ii = 2, jpim1 
    945640            zfmask(ii,ij) = ztmask(ii,ij  ) * ztmask(ii+1,ij  )   & 
    946            &              * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 
    947          END DO       
     641               &              * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 
     642         END DO 
    948643      END DO 
    949644      CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. ) 
    … …  
    953648      ! bdytmask = 1  on the computational domain AND on open boundaries 
    954649      !          = 0  elsewhere    
    955   
     650 
    956651      bdytmask(:,:) = ssmask(:,:) 
    957652 
    … …  
    981676 
    982677      ! ------------------------------------ 
    983       ! handel rim1, do as if rim 0 was land 
     678      ! handle rim1, do as if rim 0 was land 
    984679      ! ------------------------------------ 
    985680 
    … …  
    1000695         DO ii = 2, jpim1 
    1001696            zfmask(ii,ij) = ztmask(ii,ij  ) * ztmask(ii+1,ij  )   & 
    1002            &              * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 
    1003          END DO       
     697               &              * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 
     698         END DO 
    1004699      END DO 
    1005700      CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. ) 
    … …  
    1019714 
    1020715      CALL bdy_rim_treat( zumask, zvmask, zfmask, .false. )   ! compute flagu, flagv, ntreat on rim 1 
    1021  
    1022716 
    1023717      ! 
    … …  
    1053747               IF( iibi == 0     .OR. ii1 == 0     .OR. ii2 == 0     .OR. ii3 == 0     )   lrecv_bdyint(ib_bdy,igrd,1,ir) = .true. 
    1054748               IF( iibi == jpi+1 .OR. ii1 == jpi+1 .OR. ii2 == jpi+1 .OR. ii3 == jpi+1 )   lrecv_bdyint(ib_bdy,igrd,2,ir) = .true.   
    1055                IF( iibe == 0     )   lrecv_bdyext(ib_bdy,igrd,1,ir) = .true. 
    1056                IF( iibe == jpi+1 )   lrecv_bdyext(ib_bdy,igrd,2,ir) = .true.   
     749               IF( iibe == 0                                                           )   lrecv_bdyext(ib_bdy,igrd,1,ir) = .true. 
     750               IF( iibe == jpi+1                                                       )   lrecv_bdyext(ib_bdy,igrd,2,ir) = .true.   
    1057751               ! Check if neighbour has its rim parallel to its mpi subdomain border and located next to its halo 
    1058752               ! :¨¨¨¨¨|¨¨-->    |                                             |    <--¨¨|¨¨¨¨¨:  
    … …  
    1060754               ! :.....|_._:_____|   (1) W neighbour         E neighbour (2)   |_____:_._|.....: 
    1061755               IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) .AND. & 
    1062                  & ( iibi == 3     .OR. ii1 == 3     .OR. ii2 == 3     .OR. ii3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,1,ir) = .true. 
     756                  & ( iibi == 3     .OR. ii1 == 3     .OR. ii2 == 3     .OR. ii3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,1,ir)=.true. 
    1063757               IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. & 
    1064                  & ( iibi == jpi-2 .OR. ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2) )   lsend_bdyint(ib_bdy,igrd,2,ir) = .true. 
    1065                IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) .AND. iibe == 3    )   lsend_bdyext(ib_bdy,igrd,1,ir) = .true. 
    1066                IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. iibe == jpi-2)   lsend_bdyext(ib_bdy,igrd,2,ir) = .true. 
     758                  & ( iibi == jpi-2 .OR. ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2) )   lsend_bdyint(ib_bdy,igrd,2,ir)=.true. 
     759               IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) .AND. iibe == 3     )   lsend_bdyext(ib_bdy,igrd,1,ir)=.true. 
     760               IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. iibe == jpi-2 )   lsend_bdyext(ib_bdy,igrd,2,ir)=.true. 
    1067761               ! 
    1068762               ! search neighbour in the north/south direction    
    … …  
    1073767               IF( ijbi == 0     .OR. ij1 == 0     .OR. ij2 == 0     .OR. ij3 == 0     )   lrecv_bdyint(ib_bdy,igrd,3,ir) = .true. 
    1074768               IF( ijbi == jpj+1 .OR. ij1 == jpj+1 .OR. ij2 == jpj+1 .OR. ij3 == jpj+1 )   lrecv_bdyint(ib_bdy,igrd,4,ir) = .true. 
    1075                IF( ijbe == 0     )   lrecv_bdyext(ib_bdy,igrd,3,ir) = .true. 
    1076                IF( ijbe == jpj+1 )   lrecv_bdyext(ib_bdy,igrd,4,ir) = .true. 
     769               IF( ijbe == 0                                                           )   lrecv_bdyext(ib_bdy,igrd,3,ir) = .true. 
     770               IF( ijbe == jpj+1                                                       )   lrecv_bdyext(ib_bdy,igrd,4,ir) = .true. 
    1077771               ! Check if neighbour has its rim parallel to its mpi subdomain     _________  border and next to its halo 
    1078772               !   ^  |    o    |                                                :         :  
    … …  
    1080774               !      :_________:  (3) S neighbour          N neighbour (4)   v  |    o    |    
    1081775               IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) .AND. & 
    1082                  & ( ijbi == 3     .OR. ij1 == 3     .OR. ij2 == 3     .OR. ij3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,3,ir) = .true. 
     776                  & ( ijbi == 3     .OR. ij1 == 3     .OR. ij2 == 3     .OR. ij3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,3,ir)=.true. 
    1083777               IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. & 
    1084                  & ( ijbi == jpj-2 .OR. ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2) )   lsend_bdyint(ib_bdy,igrd,4,ir) = .true. 
    1085                IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) .AND. ijbe == 3    )   lsend_bdyext(ib_bdy,igrd,3,ir) = .true. 
    1086                IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. ijbe == jpj-2)   lsend_bdyext(ib_bdy,igrd,4,ir) = .true. 
     778                  & ( ijbi == jpj-2 .OR. ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2) )   lsend_bdyint(ib_bdy,igrd,4,ir)=.true. 
     779               IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) .AND. ijbe == 3     )   lsend_bdyext(ib_bdy,igrd,3,ir)=.true. 
     780               IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. ijbe == jpj-2 )   lsend_bdyext(ib_bdy,igrd,4,ir)=.true. 
    1087781            END DO 
    1088782         END DO 
    … …  
    1091785      DO ib_bdy = 1,nb_bdy 
    1092786         IF( cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. & 
    1093            & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. & 
    1094            & cn_tra(ib_bdy)   == 'orlanski' .OR. cn_tra(ib_bdy)   == 'orlanski_npo'      ) THEN 
     787            & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. & 
     788            & cn_tra(ib_bdy)   == 'orlanski' .OR. cn_tra(ib_bdy)   == 'orlanski_npo'      ) THEN 
    1095789            DO igrd = 1, jpbgrd 
    1096790               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    … …  
    1110804      ! Tidy up 
    1111805      !-------- 
    1112       IF( nb_bdy>0 )   DEALLOCATE( nbidta, nbjdta, nbrdta ) 
    1113       ! 
    1114    END SUBROUTINE bdy_segs 
     806      DEALLOCATE( nbidta, nbjdta, nbrdta ) 
     807      ! 
     808   END SUBROUTINE bdy_def 
    1115809 
    1116810 
    … …  
    1120814      !! 
    1121815      !! ** Purpose :   Initialize structures ( flagu, flagv, ntreat ) indicating how rim points 
    1122       !!                  are to be handeled in the boundary condition treatment 
    1123       !! 
    1124       !! ** Method  :   - to handel rim 0 zmasks must indicate ocean points      (set at one on rim 0 and rim 1 and interior) 
     816      !!                  are to be handled in the boundary condition treatment 
     817      !! 
     818      !! ** Method  :   - to handle rim 0 zmasks must indicate ocean points      (set at one on rim 0 and rim 1 and interior) 
    1125819      !!                            and bdymasks must be set at 0 on rim 0       (set at one on rim 1 and interior) 
    1126820      !!                    (as if rim 1 was free ocean) 
    1127       !!                - to handel rim 1 zmasks must be set at 0 on rim 0       (set at one on rim 1 and interior) 
     821      !!                - to handle rim 1 zmasks must be set at 0 on rim 0       (set at one on rim 1 and interior) 
    1128822      !!                            and bdymasks must indicate free ocean points (set at one on interior) 
    1129823      !!                    (as if rim 0 was land) 
    … …  
    13731067      CASE( 16 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij+1   ;   ii3 = ii+1   ;   ij3 = ij 
    13741068      END SELECT 
    1375     END SUBROUTINE find_neib 
     1069   END SUBROUTINE find_neib 
    13761070     
    13771071 
     1072   SUBROUTINE bdy_read_seg( kb_bdy, knblendta )  
     1073      !!---------------------------------------------------------------------- 
     1074      !!                 ***  ROUTINE bdy_coords_seg  *** 
     1075      !! 
     1076      !! ** Purpose :  build bdy coordinates with segments defined in namelist 
     1077      !! 
     1078      !! ** Method  :  read namelist nambdy_index blocks 
     1079      !! 
     1080      !!---------------------------------------------------------------------- 
     1081      INTEGER                   , INTENT (in   ) ::   kb_bdy           ! bdy number 
     1082      INTEGER, DIMENSION(jpbgrd), INTENT (  out) ::   knblendta        ! length of index arrays  
     1083      !! 
     1084      INTEGER          ::   ios                 ! Local integer output status for namelist read 
     1085      INTEGER          ::   nbdyind, nbdybeg, nbdyend 
     1086      CHARACTER(LEN=1) ::   ctypebdy   !     -        -  
     1087      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend 
     1088      !!---------------------------------------------------------------------- 
     1089 
     1090      ! No REWIND here because may need to read more than one nambdy_index namelist. 
     1091      ! Read only namelist_cfg to avoid unseccessfull overwrite  
     1092      ! keep full control of the configuration namelist 
     1093      READ  ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 ) 
     1094904   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_index in configuration namelist', lwp ) 
     1095      IF(lwm) WRITE ( numond, nambdy_index ) 
     1096       
     1097      SELECT CASE ( TRIM(ctypebdy) ) 
     1098      CASE( 'N' ) 
     1099         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
     1100            nbdyind  = jpjglo - 2  ! set boundary to whole side of model domain. 
     1101            nbdybeg  = 2 
     1102            nbdyend  = jpiglo - 1 
     1103         ENDIF 
     1104         nbdysegn = nbdysegn + 1 
     1105         npckgn(nbdysegn) = kb_bdy ! Save bdy package number 
     1106         jpjnob(nbdysegn) = nbdyind 
     1107         jpindt(nbdysegn) = nbdybeg 
     1108         jpinft(nbdysegn) = nbdyend 
     1109         ! 
     1110      CASE( 'S' ) 
     1111         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
     1112            nbdyind  = 2           ! set boundary to whole side of model domain. 
     1113            nbdybeg  = 2 
     1114            nbdyend  = jpiglo - 1 
     1115         ENDIF 
     1116         nbdysegs = nbdysegs + 1 
     1117         npckgs(nbdysegs) = kb_bdy ! Save bdy package number 
     1118         jpjsob(nbdysegs) = nbdyind 
     1119         jpisdt(nbdysegs) = nbdybeg 
     1120         jpisft(nbdysegs) = nbdyend 
     1121         ! 
     1122      CASE( 'E' ) 
     1123         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
     1124            nbdyind  = jpiglo - 2  ! set boundary to whole side of model domain. 
     1125            nbdybeg  = 2 
     1126            nbdyend  = jpjglo - 1 
     1127         ENDIF 
     1128         nbdysege = nbdysege + 1  
     1129         npckge(nbdysege) = kb_bdy ! Save bdy package number 
     1130         jpieob(nbdysege) = nbdyind 
     1131         jpjedt(nbdysege) = nbdybeg 
     1132         jpjeft(nbdysege) = nbdyend 
     1133         ! 
     1134      CASE( 'W' ) 
     1135         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
     1136            nbdyind  = 2           ! set boundary to whole side of model domain. 
     1137            nbdybeg  = 2 
     1138            nbdyend  = jpjglo - 1 
     1139         ENDIF 
     1140         nbdysegw = nbdysegw + 1 
     1141         npckgw(nbdysegw) = kb_bdy ! Save bdy package number 
     1142         jpiwob(nbdysegw) = nbdyind 
     1143         jpjwdt(nbdysegw) = nbdybeg 
     1144         jpjwft(nbdysegw) = nbdyend 
     1145         ! 
     1146      CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' ) 
     1147      END SELECT 
     1148       
     1149      ! For simplicity we assume that in case of straight bdy, arrays have the same length 
     1150      ! (even if it is true that last tangential velocity points 
     1151      ! are useless). This simplifies a little bit boundary data format (and agrees with format 
     1152      ! used so far in obc package) 
     1153       
     1154      knblendta(1:jpbgrd) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(kb_bdy) 
     1155       
     1156   END SUBROUTINE bdy_read_seg 
     1157 
     1158    
    13781159   SUBROUTINE bdy_ctl_seg 
    13791160      !!---------------------------------------------------------------------- 
    … …  
    17191500   END SUBROUTINE bdy_ctl_seg 
    17201501 
    1721  
     1502    
     1503   SUBROUTINE bdy_coords_seg( nbidta, nbjdta, nbrdta )  
     1504      !!---------------------------------------------------------------------- 
     1505      !!                 ***  ROUTINE bdy_coords_seg  *** 
     1506      !! 
     1507      !! ** Purpose :  build nbidta, nbidta, nbrdta for bdy built with segments 
     1508      !! 
     1509      !! ** Method  :   
     1510      !! 
     1511      !!---------------------------------------------------------------------- 
     1512      INTEGER, DIMENSION(:,:,:), intent(  out)  ::   nbidta, nbjdta, nbrdta   ! Index arrays: i and j indices of bdy dta 
     1513      !! 
     1514      INTEGER  ::   ii, ij, ir, iseg 
     1515      INTEGER  ::   igrd         ! grid type (t=1, u=2, v=3) 
     1516      INTEGER  ::   icount       !  
     1517      INTEGER  ::   ib_bdy       ! bdy number 
     1518      !!---------------------------------------------------------------------- 
     1519 
     1520      ! East 
     1521      !----- 
     1522      DO iseg = 1, nbdysege 
     1523         ib_bdy = npckge(iseg) 
     1524         ! 
     1525         ! ------------ T points ------------- 
     1526         igrd=1 
     1527         icount=0 
     1528         DO ir = 1, nn_rimwidth(ib_bdy) 
     1529            DO ij = jpjedt(iseg), jpjeft(iseg) 
     1530               icount = icount + 1 
     1531               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 
     1532               nbjdta(icount, igrd, ib_bdy) = ij 
     1533               nbrdta(icount, igrd, ib_bdy) = ir 
     1534            ENDDO 
     1535         ENDDO 
     1536         ! 
     1537         ! ------------ U points ------------- 
     1538         igrd=2 
     1539         icount=0 
     1540         DO ir = 1, nn_rimwidth(ib_bdy) 
     1541            DO ij = jpjedt(iseg), jpjeft(iseg) 
     1542               icount = icount + 1 
     1543               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir 
     1544               nbjdta(icount, igrd, ib_bdy) = ij 
     1545               nbrdta(icount, igrd, ib_bdy) = ir 
     1546            ENDDO 
     1547         ENDDO 
     1548         ! 
     1549         ! ------------ V points ------------- 
     1550         igrd=3 
     1551         icount=0 
     1552         DO ir = 1, nn_rimwidth(ib_bdy) 
     1553            !            DO ij = jpjedt(iseg), jpjeft(iseg) - 1 
     1554            DO ij = jpjedt(iseg), jpjeft(iseg) 
     1555               icount = icount + 1 
     1556               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 
     1557               nbjdta(icount, igrd, ib_bdy) = ij 
     1558               nbrdta(icount, igrd, ib_bdy) = ir 
     1559            ENDDO 
     1560            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     1561            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     1562         ENDDO 
     1563      ENDDO 
     1564      ! 
     1565      ! West 
     1566      !----- 
     1567      DO iseg = 1, nbdysegw 
     1568         ib_bdy = npckgw(iseg) 
     1569         ! 
     1570         ! ------------ T points ------------- 
     1571         igrd=1 
     1572         icount=0 
     1573         DO ir = 1, nn_rimwidth(ib_bdy) 
     1574            DO ij = jpjwdt(iseg), jpjwft(iseg) 
     1575               icount = icount + 1 
     1576               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     1577               nbjdta(icount, igrd, ib_bdy) = ij 
     1578               nbrdta(icount, igrd, ib_bdy) = ir 
     1579            ENDDO 
     1580         ENDDO 
     1581         ! 
     1582         ! ------------ U points ------------- 
     1583         igrd=2 
     1584         icount=0 
     1585         DO ir = 1, nn_rimwidth(ib_bdy) 
     1586            DO ij = jpjwdt(iseg), jpjwft(iseg) 
     1587               icount = icount + 1 
     1588               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     1589               nbjdta(icount, igrd, ib_bdy) = ij 
     1590               nbrdta(icount, igrd, ib_bdy) = ir 
     1591            ENDDO 
     1592         ENDDO 
     1593         ! 
     1594         ! ------------ V points ------------- 
     1595         igrd=3 
     1596         icount=0 
     1597         DO ir = 1, nn_rimwidth(ib_bdy) 
     1598            !            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 
     1599            DO ij = jpjwdt(iseg), jpjwft(iseg) 
     1600               icount = icount + 1 
     1601               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     1602               nbjdta(icount, igrd, ib_bdy) = ij 
     1603               nbrdta(icount, igrd, ib_bdy) = ir 
     1604            ENDDO 
     1605            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     1606            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     1607         ENDDO 
     1608      ENDDO 
     1609      ! 
     1610      ! North 
     1611      !----- 
     1612      DO iseg = 1, nbdysegn 
     1613         ib_bdy = npckgn(iseg) 
     1614         ! 
     1615         ! ------------ T points ------------- 
     1616         igrd=1 
     1617         icount=0 
     1618         DO ir = 1, nn_rimwidth(ib_bdy) 
     1619            DO ii = jpindt(iseg), jpinft(iseg) 
     1620               icount = icount + 1 
     1621               nbidta(icount, igrd, ib_bdy) = ii 
     1622               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir  
     1623               nbrdta(icount, igrd, ib_bdy) = ir 
     1624            ENDDO 
     1625         ENDDO 
     1626         ! 
     1627         ! ------------ U points ------------- 
     1628         igrd=2 
     1629         icount=0 
     1630         DO ir = 1, nn_rimwidth(ib_bdy) 
     1631            !            DO ii = jpindt(iseg), jpinft(iseg) - 1 
     1632            DO ii = jpindt(iseg), jpinft(iseg) 
     1633               icount = icount + 1 
     1634               nbidta(icount, igrd, ib_bdy) = ii 
     1635               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 
     1636               nbrdta(icount, igrd, ib_bdy) = ir 
     1637            ENDDO 
     1638            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     1639            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     1640         ENDDO 
     1641         ! 
     1642         ! ------------ V points ------------- 
     1643         igrd=3 
     1644         icount=0 
     1645         DO ir = 1, nn_rimwidth(ib_bdy) 
     1646            DO ii = jpindt(iseg), jpinft(iseg) 
     1647               icount = icount + 1 
     1648               nbidta(icount, igrd, ib_bdy) = ii 
     1649               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir 
     1650               nbrdta(icount, igrd, ib_bdy) = ir 
     1651            ENDDO 
     1652         ENDDO 
     1653      ENDDO 
     1654      ! 
     1655      ! South 
     1656      !----- 
     1657      DO iseg = 1, nbdysegs 
     1658         ib_bdy = npckgs(iseg) 
     1659         ! 
     1660         ! ------------ T points ------------- 
     1661         igrd=1 
     1662         icount=0 
     1663         DO ir = 1, nn_rimwidth(ib_bdy) 
     1664            DO ii = jpisdt(iseg), jpisft(iseg) 
     1665               icount = icount + 1 
     1666               nbidta(icount, igrd, ib_bdy) = ii 
     1667               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
     1668               nbrdta(icount, igrd, ib_bdy) = ir 
     1669            ENDDO 
     1670         ENDDO 
     1671         ! 
     1672         ! ------------ U points ------------- 
     1673         igrd=2 
     1674         icount=0 
     1675         DO ir = 1, nn_rimwidth(ib_bdy) 
     1676            !            DO ii = jpisdt(iseg), jpisft(iseg) - 1 
     1677            DO ii = jpisdt(iseg), jpisft(iseg) 
     1678               icount = icount + 1 
     1679               nbidta(icount, igrd, ib_bdy) = ii 
     1680               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
     1681               nbrdta(icount, igrd, ib_bdy) = ir 
     1682            ENDDO 
     1683            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     1684            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     1685         ENDDO 
     1686         ! 
     1687         ! ------------ V points ------------- 
     1688         igrd=3 
     1689         icount=0 
     1690         DO ir = 1, nn_rimwidth(ib_bdy) 
     1691            DO ii = jpisdt(iseg), jpisft(iseg) 
     1692               icount = icount + 1 
     1693               nbidta(icount, igrd, ib_bdy) = ii 
     1694               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
     1695               nbrdta(icount, igrd, ib_bdy) = ir 
     1696            ENDDO 
     1697         ENDDO 
     1698      ENDDO 
     1699 
     1700       
     1701   END SUBROUTINE bdy_coords_seg 
     1702    
     1703    
    17221704   SUBROUTINE bdy_ctl_corn( ib1, ib2 ) 
    17231705      !!---------------------------------------------------------------------- 

  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdytides.F90

    r11048 r11223  
    7070      INTEGER                                   ::   inum, igrd 
    7171      INTEGER, DIMENSION(3)                     ::   ilen0       !: length of boundary data (from OBC arrays) 
    72       INTEGER, POINTER, DIMENSION(:)            ::   nblen, nblenrim     ! short cuts 
    7372      INTEGER                                   ::   ios                 ! Local integer output status for namelist read 
    7473      CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file  
    … …  
    7776      !! 
    7877      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut    
    79       TYPE(MAP_POINTER), DIMENSION(jpbgrd)      ::   ibmap_ptr           !: array of pointers to nbmap 
    8078      !! 
    8179      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj 
    8280      !!---------------------------------------------------------------------- 
    8381      ! 
    84       IF (nb_bdy>0) THEN 
    85          IF(lwp) WRITE(numout,*) 
    86          IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries' 
    87          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    88       ENDIF 
     82      IF(lwp) WRITE(numout,*) 
     83      IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries' 
     84      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    8985 
    9086      REWIND(numnam_cfg) 
    … …  
    9490            ! 
    9591            td => tides(ib_bdy) 
    96             nblen => idx_bdy(ib_bdy)%nblen 
    97             nblenrim => idx_bdy(ib_bdy)%nblenrim 
    9892 
    9993            ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries 
    10094            filtide(:) = '' 
    10195 
    102             ! Don't REWIND here - may need to read more than one of these namelists.  
     96            REWIND( numnam_ref ) 
    10397            READ  ( numnam_ref, nambdy_tide, IOSTAT = ios, ERR = 901) 
    10498901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_tide in reference namelist', lwp ) 
     99            ! Don't REWIND here - may need to read more than one of these namelists.  
    105100            READ  ( numnam_cfg, nambdy_tide, IOSTAT = ios, ERR = 902 ) 
    106101902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_tide in configuration namelist', lwp ) 
    … …  
    125120            ! JC: If FRS scheme is used, we assume that tidal is needed over the whole 
    126121            ! relaxation area       
    127             IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = nblen   (:) 
    128             ELSE                                   ;   ilen0(:) = nblenrim(:) 
     122            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = idx_bdy(ib_bdy)%nblen   (:) 
     123            ELSE                                   ;   ilen0(:) = idx_bdy(ib_bdy)%nblenrim(:) 
    129124            ENDIF 
    130125 
    … …  
    210205               ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) ) 
    211206               ! 
    212                ! Set map structure 
    213                ibmap_ptr(1)%ptr => idx_bdy(ib_bdy)%nbmap(:,1)   ;   ibmap_ptr(1)%ll_unstruc = ln_coords_file(ib_bdy) 
    214                ibmap_ptr(2)%ptr => idx_bdy(ib_bdy)%nbmap(:,2)   ;   ibmap_ptr(2)%ll_unstruc = ln_coords_file(ib_bdy) 
    215                ibmap_ptr(3)%ptr => idx_bdy(ib_bdy)%nbmap(:,3)   ;   ibmap_ptr(3)%ll_unstruc = ln_coords_file(ib_bdy) 
    216  
    217207               ! Open files and read in tidal forcing data 
    218208               ! ----------------------------------------- 
    … …  
    222212                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc' 
    223213                  CALL iom_open( clfile, inum ) 
    224                   CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1,  ibmap_ptr(1) ) 
     214                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
    225215                  td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1) 
    226                   CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1,  ibmap_ptr(1) ) 
     216                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
    227217                  td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1) 
    228218                  CALL iom_close( inum ) 
    … …  
    230220                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc' 
    231221                  CALL iom_open( clfile, inum ) 
    232                   CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, ibmap_ptr(2) ) 
     222                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
    233223                  td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1) 
    234                   CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, ibmap_ptr(2) ) 
     224                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
    235225                  td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1) 
    236226                  CALL iom_close( inum ) 
    … …  
    238228                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc' 
    239229                  CALL iom_open( clfile, inum ) 
    240                   CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, ibmap_ptr(3) ) 
     230                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
    241231                  td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1) 
    242                   CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, ibmap_ptr(3) ) 
     232                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
    243233                  td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1) 
    244234                  CALL iom_close( inum ) 
    … …  
    272262 
    273263 
    274    SUBROUTINE bdytide_update( kt, idx, dta, td, jit, time_offset ) 
     264   SUBROUTINE bdytide_update( kt, idx, dta, td, kit, kt_offset ) 
    275265      !!---------------------------------------------------------------------- 
    276266      !!                 ***  SUBROUTINE bdytide_update  *** 
    … …  
    283273      TYPE(OBC_DATA)   , INTENT(inout) ::   dta         ! OBC external data 
    284274      TYPE(TIDES_DATA) , INTENT(inout) ::   td          ! tidal harmonics data 
    285       INTEGER, OPTIONAL, INTENT(in   ) ::   jit         ! Barotropic timestep counter (for timesplitting option) 
    286       INTEGER, OPTIONAL, INTENT(in   ) ::   time_offset ! time offset in units of timesteps. NB. if jit 
     275      INTEGER, OPTIONAL, INTENT(in   ) ::   kit         ! Barotropic timestep counter (for timesplitting option) 
     276      INTEGER, OPTIONAL, INTENT(in   ) ::   kt_offset   ! time offset in units of timesteps. NB. if kit 
    287277      !                                                 ! is present then units = subcycle timesteps. 
    288       !                                                 ! time_offset = 0  => get data at "now"    time level 
    289       !                                                 ! time_offset = -1 => get data at "before" time level 
    290       !                                                 ! time_offset = +1 => get data at "after"  time level 
     278      !                                                 ! kt_offset = 0  => get data at "now"    time level 
     279      !                                                 ! kt_offset = -1 => get data at "before" time level 
     280      !                                                 ! kt_offset = +1 => get data at "after"  time level 
    291281      !                                                 ! etc. 
    292282      ! 
    … …  
    303293 
    304294      zflag=1 
    305       IF ( PRESENT(jit) ) THEN 
    306         IF ( jit /= 1 ) zflag=0 
     295      IF ( PRESENT(kit) ) THEN 
     296        IF ( kit /= 1 ) zflag=0 
    307297      ENDIF 
    308298 
    … …  
    323313 
    324314      time_add = 0 
    325       IF( PRESENT(time_offset) ) THEN 
    326          time_add = time_offset 
     315      IF( PRESENT(kt_offset) ) THEN 
     316         time_add = kt_offset 
    327317      ENDIF 
    328318          
    329       IF( PRESENT(jit) ) THEN   
    330          z_arg = ((kt-kt_tide) * rdt + (jit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) ) 
     319      IF( PRESENT(kit) ) THEN   
     320         z_arg = ((kt-kt_tide) * rdt + (kit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) ) 
    331321      ELSE                               
    332322         z_arg = ((kt-kt_tide)+time_add) * rdt 
    … …  
    361351 
    362352 
    363    SUBROUTINE bdy_dta_tides( kt, kit, time_offset ) 
     353   SUBROUTINE bdy_dta_tides( kt, kit, kt_offset ) 
    364354      !!---------------------------------------------------------------------- 
    365355      !!                 ***  SUBROUTINE bdy_dta_tides  *** 
    … …  
    370360      INTEGER,           INTENT(in) ::   kt          ! Main timestep counter 
    371361      INTEGER, OPTIONAL, INTENT(in) ::   kit         ! Barotropic timestep counter (for timesplitting option) 
    372       INTEGER, OPTIONAL, INTENT(in) ::   time_offset ! time offset in units of timesteps. NB. if kit 
     362      INTEGER, OPTIONAL, INTENT(in) ::   kt_offset  ! time offset in units of timesteps. NB. if kit 
    373363      !                                              ! is present then units = subcycle timesteps. 
    374       !                                              ! time_offset = 0  => get data at "now"    time level 
    375       !                                              ! time_offset = -1 => get data at "before" time level 
    376       !                                              ! time_offset = +1 => get data at "after"  time level 
     364      !                                              ! kt_offset = 0  => get data at "now"    time level 
     365      !                                              ! kt_offset = -1 => get data at "before" time level 
     366      !                                              ! kt_offset = +1 => get data at "after"  time level 
    377367      !                                              ! etc. 
    378368      ! 
    … …  
    389379 
    390380      time_add = 0 
    391       IF( PRESENT(time_offset) ) THEN 
    392          time_add = time_offset 
     381      IF( PRESENT(kt_offset) ) THEN 
     382         time_add = kt_offset 
    393383      ENDIF 
    394384       
    … …  
    435425            ! If time splitting, initialize arrays from slow varying open boundary data: 
    436426            IF ( PRESENT(kit) ) THEN            
    437                IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) 
    438                IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) 
    439                IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) 
     427               IF ( dta_bdy(ib_bdy)%lneed_ssh  ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) 
     428               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) 
     429               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) 
    440430            ENDIF 
    441431            ! 
    … …  
    447437               z_sist = zramp * SIN( z_sarg ) 
    448438               ! 
    449                IF ( dta_bdy(ib_bdy)%ll_ssh ) THEN 
     439               IF ( dta_bdy(ib_bdy)%lneed_ssh ) THEN 
    450440                  igrd=1                              ! SSH on tracer grid 
    451441                  DO ib = 1, ilen0(igrd) 
    … …  
    456446               ENDIF 
    457447               ! 
    458                IF ( dta_bdy(ib_bdy)%ll_u2d ) THEN 
     448               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) THEN 
    459449                  igrd=2                              ! U grid 
    460450                  DO ib = 1, ilen0(igrd) 
    … …  
    463453                        &                        tides(ib_bdy)%u(ib,itide,2)*z_sist ) 
    464454                  END DO 
    465                ENDIF 
    466                ! 
    467                IF ( dta_bdy(ib_bdy)%ll_v2d ) THEN 
    468455                  igrd=3                              ! V grid 
    469456                  DO ib = 1, ilen0(igrd)  

  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DOM/dommsk.F90

    r10425 r11223  
    101101         &             cn_ice, nn_ice_dta,                                     & 
    102102         &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     & 
    103          &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy 
     103         &             ln_vol, nn_volctl, nn_rimwidth 
    104104      !!--------------------------------------------------------------------- 
    105105      ! 

  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynspg_ts.F90

    r10742 r11223  
    718718         !                                                !  ------------------ 
    719719         ! Update only tidal forcing at open boundaries 
    720          IF( ln_bdy      .AND. ln_tide )   CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 
    721          IF( ln_tide_pot .AND. ln_tide )   CALL upd_tide     ( kt, kit=jn, time_offset= noffset   ) 
     720         IF( ln_bdy      .AND. ln_tide )   CALL bdy_dta_tides( kt, kit=jn, kt_offset= noffset+1 ) 
     721         IF( ln_tide_pot .AND. ln_tide )   CALL upd_tide     ( kt, kit=jn, kt_offset= noffset   ) 
    722722         ! 
    723723         ! Set extrapolation coefficients for predictor step: 

  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/IOM/iom.F90

    r11194 r11223  
    835835 
    836836 
    837    FUNCTION iom_varid ( kiomid, cdvar, kdimsz, kndims, ldstop )   
     837   FUNCTION iom_varid ( kiomid, cdvar, kdimsz, kndims, lduld, ldstop )   
    838838      !!----------------------------------------------------------------------- 
    839839      !!                  ***  FUNCTION  iom_varid  *** 
    … …  
    844844      CHARACTER(len=*)     , INTENT(in   )           ::   cdvar    ! name of the variable 
    845845      INTEGER, DIMENSION(:), INTENT(  out), OPTIONAL ::   kdimsz   ! size of each dimension 
    846       INTEGER,               INTENT(  out), OPTIONAL ::   kndims   ! size of the dimensions 
     846      INTEGER              , INTENT(  out), OPTIONAL ::   kndims   ! number of dimensions 
     847      LOGICAL              , INTENT(  out), OPTIONAL ::   lduld    ! true if the last dimension is unlimited (time) 
    847848      LOGICAL              , INTENT(in   ), OPTIONAL ::   ldstop   ! stop if looking for non-existing variable (default = .TRUE.) 
    848849      ! 
    … …  
    874875               iiv = iiv + 1 
    875876               IF( iiv <= jpmax_vars ) THEN 
    876                   iom_varid = iom_nf90_varid( kiomid, cdvar, iiv, kdimsz, kndims ) 
     877                  iom_varid = iom_nf90_varid( kiomid, cdvar, iiv, kdimsz, kndims, lduld ) 
    877878               ELSE 
    878879                  CALL ctl_stop( trim(clinfo), 'Too many variables in the file '//iom_file(kiomid)%name,   & 
    … …  
    892893               ENDIF 
    893894               IF( PRESENT(kndims) )  kndims = iom_file(kiomid)%ndims(iiv) 
     895               IF( PRESENT( lduld) )  lduld  = iom_file(kiomid)%luld( iiv) 
    894896            ENDIF 
    895897         ENDIF 

  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/IOM/iom_nf90.F90

    r10522 r11223  
    187187 
    188188 
    189    FUNCTION iom_nf90_varid ( kiomid, cdvar, kiv, kdimsz, kndims )   
     189   FUNCTION iom_nf90_varid ( kiomid, cdvar, kiv, kdimsz, kndims, lduld )   
    190190      !!----------------------------------------------------------------------- 
    191191      !!                  ***  FUNCTION  iom_varid  *** 
    … …  
    198198      INTEGER, DIMENSION(:), INTENT(  out), OPTIONAL ::   kdimsz   ! size of the dimensions 
    199199      INTEGER,               INTENT(  out), OPTIONAL ::   kndims   ! size of the dimensions 
     200      LOGICAL              , INTENT(  out), OPTIONAL ::   lduld    ! true if the last dimension is unlimited (time) 
    200201      ! 
    201202      INTEGER                        ::   iom_nf90_varid   ! iom variable Id 
    … …  
    251252         ENDIF 
    252253         IF( PRESENT(kndims) )  kndims = iom_file(kiomid)%ndims(kiv) 
     254         IF( PRESENT( lduld) )  lduld  = iom_file(kiomid)%luld(kiv) 
    253255      ELSE   
    254256         iom_nf90_varid = -1   !   variable not found, return error code: -1 

  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/LBC/mppini.F90

    r11192 r11223  
    167167           &             cn_ice, nn_ice_dta,                                     & 
    168168           &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     & 
    169            &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy 
     169           &             ln_vol, nn_volctl, nn_rimwidth 
    170170      !!---------------------------------------------------------------------- 
    171171 

  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/fldread.F90

    r10425 r11223  
    8080      INTEGER                         ::   nreclast     ! last record to be read in the current file 
    8181      CHARACTER(len = 256)            ::   lsmname      ! current name of the NetCDF mask file acting as a key 
     82      !                                                 ! variables related to BDY 
    8283      INTEGER                         ::   igrd         ! grid type for bdy data 
    8384      INTEGER                         ::   ibdy         ! bdy set id number 
     85      INTEGER, POINTER, DIMENSION(:)  ::   imap         ! Array of integer pointers to 1D arrays 
     86      LOGICAL                         ::   ltotvel      ! total velocity or not (T/F) 
    8487   END TYPE FLD 
    85  
    86    TYPE, PUBLIC ::   MAP_POINTER      !: Map from input data file to local domain 
    87       INTEGER, POINTER, DIMENSION(:)  ::  ptr           ! Array of integer pointers to 1D arrays 
    88       LOGICAL                         ::  ll_unstruc    ! Unstructured (T) or structured (F) boundary data file 
    89    END TYPE MAP_POINTER 
    9088 
    9189!$AGRIF_DO_NOT_TREAT 
    … …  
    129127CONTAINS 
    130128 
    131    SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset, jpk_bdy, fvl ) 
     129   SUBROUTINE fld_read( kt, kn_fsbc, sd, kit, kt_offset ) 
    132130      !!--------------------------------------------------------------------- 
    133131      !!                    ***  ROUTINE fld_read  *** 
    … …  
    144142      INTEGER  , INTENT(in   )               ::   kn_fsbc   ! sbc computation period (in time step)  
    145143      TYPE(FLD), INTENT(inout), DIMENSION(:) ::   sd        ! input field related variables 
    146       TYPE(MAP_POINTER),INTENT(in), OPTIONAL, DIMENSION(:) ::   map   ! global-to-local mapping indices 
    147144      INTEGER  , INTENT(in   ), OPTIONAL     ::   kit       ! subcycle timestep for timesplitting option 
    148145      INTEGER  , INTENT(in   ), OPTIONAL     ::   kt_offset ! provide fields at time other than "now" 
    … …  
    150147      !                                                     !   kt_offset = +1 => fields at "after"  time level 
    151148      !                                                     !   etc. 
    152       INTEGER  , INTENT(in   ), OPTIONAL     ::   jpk_bdy   ! number of vertical levels in the BDY data 
    153       LOGICAL  , INTENT(in   ), OPTIONAL     ::   fvl   ! number of vertical levels in the BDY data 
    154149      !! 
    155150      INTEGER  ::   itmp         ! local variable 
    … …  
    166161      REAL(wp) ::   ztintb       ! ratio applied to before records when doing time interpolation 
    167162      CHARACTER(LEN=1000) ::   clfmt  ! write format 
    168       TYPE(MAP_POINTER)   ::   imap   ! global-to-local mapping indices 
    169163      !!--------------------------------------------------------------------- 
    170164      ll_firstcall = kt == nit000 
    … …  
    175169      ENDIF 
    176170      IF( PRESENT(kt_offset) )   it_offset = kt_offset 
    177  
    178       imap%ptr => NULL() 
    179171 
    180172      ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar  
    … …  
    188180      IF( ll_firstcall ) THEN                      ! initialization 
    189181         DO jf = 1, imf  
    190             IF( PRESENT(map) ) imap = map(jf) 
    191                IF( PRESENT(jpk_bdy) ) THEN 
    192                   CALL fld_init( kn_fsbc, sd(jf), imap, jpk_bdy, fvl )  ! read each before field (put them in after as they will be swapped) 
    193                ELSE 
    194                   CALL fld_init( kn_fsbc, sd(jf), imap )  ! read each before field (put them in after as they will be swapped) 
    195                ENDIF 
     182            IF( TRIM(sd(jf)%clrootname) == 'NOT USED' )   CYCLE 
     183            CALL fld_init( kn_fsbc, sd(jf) )       ! read each before field (put them in after as they will be swapped) 
    196184         END DO 
    197185         IF( lwp ) CALL wgt_print()                ! control print 
    … …  
    202190         ! 
    203191         DO jf = 1, imf                            ! ---   loop over field   --- ! 
    204              
     192 
     193            IF( TRIM(sd(jf)%clrootname) == 'NOT USED' )   CYCLE 
     194                       
    205195            IF( isecsbc > sd(jf)%nrec_a(2) .OR. ll_firstcall ) THEN    ! read/update the after data? 
    206  
    207                IF( PRESENT(map) )   imap = map(jf)   ! temporary definition of map 
    208196 
    209197               sd(jf)%nrec_b(:) = sd(jf)%nrec_a(:)                                  ! swap before record informations 
    … …  
    222210                  itmp = sd(jf)%nrec_a(1)                       ! temporary storage 
    223211                  sd(jf)%nrec_a(1) = sd(jf)%nreclast            ! read the last record of the file currently opened 
    224                   CALL fld_get( sd(jf), imap )                  ! read after data 
     212                  CALL fld_get( sd(jf) )                        ! read after data 
    225213                  sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)   ! re-swap before record field 
    226214                  sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1)           ! update before record informations 
    … …  
    240228                     &                   .AND. sd(jf)%nrec_b(1) /= sd(jf)%nrec_a(1) - 1 ) THEN    
    241229                     sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) - 1       ! move back to before record 
    242                      CALL fld_get( sd(jf), imap )                  ! read after data 
     230                     CALL fld_get( sd(jf) )                        ! read after data 
    243231                     sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)   ! re-swap before record field 
    244232                     sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1)           ! update before record informations 
    … …  
    285273                   
    286274               ! read after data 
    287                IF( PRESENT(jpk_bdy) ) THEN 
    288                   CALL fld_get( sd(jf), imap, jpk_bdy, fvl) 
    289                ELSE 
    290                   CALL fld_get( sd(jf), imap ) 
    291                ENDIF 
     275               CALL fld_get( sd(jf) ) 
     276                
    292277            ENDIF   ! read new data? 
    293278         END DO                                    ! --- end loop over field --- ! 
    … …  
    296281 
    297282         DO jf = 1, imf                            ! ---   loop over field   --- ! 
     283            ! 
     284            IF( TRIM(sd(jf)%clrootname) == 'NOT USED' )   CYCLE 
    298285            ! 
    299286            IF( sd(jf)%ln_tint ) THEN              ! temporal interpolation 
    … …  
    327314 
    328315 
    329    SUBROUTINE fld_init( kn_fsbc, sdjf, map , jpk_bdy, fvl) 
     316   SUBROUTINE fld_init( kn_fsbc, sdjf ) 
    330317      !!--------------------------------------------------------------------- 
    331318      !!                    ***  ROUTINE fld_init  *** 
    … …  
    336323      INTEGER  , INTENT(in   ) ::   kn_fsbc      ! sbc computation period (in time step)  
    337324      TYPE(FLD), INTENT(inout) ::   sdjf         ! input field related variables 
    338       TYPE(MAP_POINTER),INTENT(in) ::   map      ! global-to-local mapping indices 
    339       INTEGER  , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data 
    340       LOGICAL  , INTENT(in), OPTIONAL :: fvl     ! number of vertical levels in the BDY data 
    341325      !! 
    342326      LOGICAL :: llprevyr              ! are we reading previous year  file? 
    … …  
    351335      CHARACTER(LEN=1000) ::   clfmt   ! write format 
    352336      !!--------------------------------------------------------------------- 
     337      ! 
    353338      llprevyr   = .FALSE. 
    354339      llprevmth  = .FALSE. 
    … …  
    433418         ! 
    434419         ! read before data in after arrays(as we will swap it later) 
    435          IF( PRESENT(jpk_bdy) ) THEN 
    436             CALL fld_get( sdjf, map, jpk_bdy, fvl ) 
    437          ELSE 
    438             CALL fld_get( sdjf, map ) 
    439          ENDIF 
     420         CALL fld_get( sdjf ) 
    440421         ! 
    441422         clfmt = "('   fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')" 
    … …  
    613594 
    614595 
    615    SUBROUTINE fld_get( sdjf, map, jpk_bdy, fvl ) 
     596   SUBROUTINE fld_get( sdjf ) 
    616597      !!--------------------------------------------------------------------- 
    617598      !!                    ***  ROUTINE fld_get  *** 
    … …  
    620601      !!---------------------------------------------------------------------- 
    621602      TYPE(FLD)        , INTENT(inout) ::   sdjf   ! input field related variables 
    622       TYPE(MAP_POINTER), INTENT(in   ) ::   map    ! global-to-local mapping indices 
    623       INTEGER  , INTENT(in), OPTIONAL  ::   jpk_bdy ! number of vertical levels in the bdy data 
    624       LOGICAL  , INTENT(in), OPTIONAL  ::   fvl     ! number of vertical levels in the bdy data 
    625603      ! 
    626604      INTEGER ::   ipk      ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
    … …  
    634612      ipk = SIZE( sdjf%fnow, 3 ) 
    635613      ! 
    636       IF( ASSOCIATED(map%ptr) ) THEN 
    637          IF( PRESENT(jpk_bdy) ) THEN 
    638             IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2),                & 
    639                                                         sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl ) 
    640             ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ),                & 
    641                                                         sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl ) 
    642             ENDIF 
    643          ELSE 
    644             IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map ) 
    645             ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1), map ) 
    646             ENDIF 
    647          ENDIF         
     614      IF( ASSOCIATED(sdjf%imap) ) THEN 
     615         IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2),                & 
     616            &                                        sdjf%nrec_a(1), sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel ) 
     617         ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ),                & 
     618            &                                        sdjf%nrec_a(1), sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel ) 
     619         ENDIF 
    648620      ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 
    649621         CALL wgt_list( sdjf, iw ) 
    … …  
    700672   END SUBROUTINE fld_get 
    701673 
    702    SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy, fvl ) 
     674    
     675   SUBROUTINE fld_map( knum, cdvar, pdta, krec, kmap, kgrd, kbdy, ldtotvel ) 
    703676      !!--------------------------------------------------------------------- 
    704677      !!                    ***  ROUTINE fld_map  *** 
    … …  
    707680      !!                using a general mapping (for open boundaries) 
    708681      !!---------------------------------------------------------------------- 
    709  
    710       USE bdy_oce, ONLY: ln_bdy, idx_bdy, dta_global, dta_global_z, dta_global_dz, dta_global2, dta_global2_z, dta_global2_dz                 ! workspace to read in global data arrays 
    711  
    712       INTEGER                   , INTENT(in ) ::   num     ! stream number 
    713       CHARACTER(LEN=*)          , INTENT(in ) ::   clvar   ! variable name 
    714       REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta     ! output field on model grid (2 dimensional) 
    715       INTEGER                   , INTENT(in ) ::   nrec    ! record number to read (ie time slice) 
    716       TYPE(MAP_POINTER)         , INTENT(in ) ::   map     ! global-to-local mapping indices 
    717       INTEGER  , INTENT(in), OPTIONAL         ::   igrd, ibdy, jpk_bdy  ! grid type, set number and number of vertical levels in the bdy data 
    718       LOGICAL  , INTENT(in), OPTIONAL         ::   fvl     ! grid type, set number and number of vertical levels in the bdy data 
    719       INTEGER                                 ::   jpkm1_bdy! number of vertical levels in the bdy data minus 1 
    720       !! 
    721       INTEGER                                 ::   ipi      ! length of boundary data on local process 
    722       INTEGER                                 ::   ipj      ! length of dummy dimension ( = 1 ) 
    723       INTEGER                                 ::   ipk      ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
    724       INTEGER                                 ::   ilendta  ! length of data in file 
    725       INTEGER                                 ::   idvar    ! variable ID 
    726       INTEGER                                 ::   ib, ik, ji, jj   ! loop counters 
    727       INTEGER                                 ::   ierr 
    728       REAL(wp)                                ::   fv          ! fillvalue  
    729       REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read    ! work space for global data 
    730       REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read_z  ! work space for global data 
    731       REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read_dz ! work space for global data 
    732       !!--------------------------------------------------------------------- 
    733       ! 
    734       ipi = SIZE( dta, 1 ) 
    735       ipj = 1 
    736       ipk = SIZE( dta, 3 ) 
    737       ! 
    738       idvar   = iom_varid( num, clvar ) 
    739       ilendta = iom_file(num)%dimsz(1,idvar) 
    740  
    741       IF ( ln_bdy ) THEN 
    742          ipj = iom_file(num)%dimsz(2,idvar) 
    743          IF( map%ll_unstruc) THEN   ! unstructured open boundary data file 
    744             dta_read => dta_global 
    745             IF( PRESENT(jpk_bdy) ) THEN 
    746                IF( jpk_bdy>0 ) THEN 
    747                   dta_read_z => dta_global_z 
    748                   dta_read_dz => dta_global_dz 
    749                   jpkm1_bdy = jpk_bdy-1 
    750                ENDIF 
    751             ENDIF 
    752          ELSE                       ! structured open boundary file 
    753             dta_read => dta_global2 
    754             IF( PRESENT(jpk_bdy) ) THEN 
    755                IF( jpk_bdy>0 ) THEN 
    756                   dta_read_z => dta_global2_z 
    757                   dta_read_dz => dta_global2_dz 
    758                   jpkm1_bdy = jpk_bdy-1 
    759                ENDIF 
    760             ENDIF 
    761          ENDIF 
    762       ENDIF 
    763  
    764       IF(lwp) WRITE(numout,*) 'Dim size for ',        TRIM(clvar),' is ', ilendta 
    765       IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk 
    766       ! 
    767       SELECT CASE( ipk ) 
    768       CASE(1)        ;    
    769       CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1    ), nrec ) 
    770          IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 
    771             DO ib = 1, ipi 
    772               DO ik = 1, ipk 
    773                 dta(ib,1,ik) =  dta_read(map%ptr(ib),1,ik) 
    774               END DO 
    775             END DO 
    776          ELSE ! we assume that this is a structured open boundary file 
    777             DO ib = 1, ipi 
    778                jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    779                ji=map%ptr(ib)-(jj-1)*ilendta 
    780                DO ik = 1, ipk 
    781                   dta(ib,1,ik) =  dta_read(ji,jj,ik) 
    782                END DO 
    783             END DO 
    784          ENDIF 
     682      INTEGER                   , INTENT(in   ) ::   knum         ! stream number 
     683      CHARACTER(LEN=*)          , INTENT(in   ) ::   cdvar        ! variable name 
     684      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta         ! bdy output field on model grid 
     685      INTEGER                   , INTENT(in   ) ::   krec         ! record number to read (ie time slice) 
     686      INTEGER , DIMENSION(:)    , INTENT(in   ) ::   kmap         ! global-to-local bdy mapping indices 
     687      ! optional variables used for vertical interpolation: 
     688      INTEGER, OPTIONAL         , INTENT(in   ) ::   kgrd         ! grid type (t, u, v) 
     689      INTEGER, OPTIONAL         , INTENT(in   ) ::   kbdy         ! bdy number 
     690      LOGICAL, OPTIONAL         , INTENT(in   ) ::   ldtotvel     ! true if toal ( = barotrop + barocline) velocity 
     691      !! 
     692      INTEGER                                   ::   ipi          ! length of boundary data on local process 
     693      INTEGER                                   ::   ipj          ! length of dummy dimension ( = 1 ) 
     694      INTEGER                                   ::   ipk          ! number of vertical levels of pdta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
     695      INTEGER                                   ::   ipkb         ! number of vertical levels in boundary data file 
     696      INTEGER                                   ::   idvar        ! variable ID 
     697      INTEGER                                   ::   indims       ! number of dimensions of the variable 
     698      INTEGER, DIMENSION(4)                     ::   idimsz       ! size of variable dimensions  
     699      REAL(wp)                                  ::   zfv          ! fillvalue  
     700      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zz_read      ! work space for global boundary data 
     701      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zdta_read    ! work space local data requiring vertical interpolation 
     702      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zdta_read_z  ! work space local data requiring vertical interpolation 
     703      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zdta_read_dz ! work space local data requiring vertical interpolation 
     704      CHARACTER(LEN=1),DIMENSION(3)             ::   clgrid 
     705      LOGICAL                                   ::   lluld        ! is the variable using the unlimited dimension 
     706      !!--------------------------------------------------------------------- 
     707      ! 
     708      clgrid = (/'t','u','v'/) 
     709      ! 
     710      ipi = SIZE( pdta, 1 ) 
     711      ipj = SIZE( pdta, 2 )   ! must be equal to 1 
     712      ipk = SIZE( pdta, 3 ) 
     713      ! 
     714      idvar = iom_varid( knum, cdvar, kndims = indims, kdimsz = idimsz, lduld = lluld  ) 
     715      IF( indims == 4 .OR. ( indims == 3 .AND. .NOT. lluld ) ) THEN   ;   ipkb = idimsz(3)   ! xy(zl)t or xy(zl) 
     716      ELSE                                                            ;   ipkb = 1           ! xy or xyt 
     717      ENDIF 
     718      ! 
     719      ALLOCATE( zz_read( idimsz(1), idimsz(2), ipkb ) )  ! ++++++++ !!! this can be very big...          
     720      ! 
     721      IF( ipk == 1 ) THEN 
     722 
     723         IF( ipkb /= 1 ) CALL ctl_stop( 'fld_map : we must have ipkb = 1 to read surface data' ) 
     724         CALL iom_get ( knum, jpdom_unknown, cdvar, zz_read(:,:,1), krec )   ! call iom_get with a 2D file 
     725         CALL fld_map_core( zz_read, kmap, pdta ) 
    785726 
    786727      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    787728      ! Do we include something here to adjust barotropic velocities ! 
    788729      ! in case of a depth difference between bdy files and          ! 
    789       ! bathymetry in the case ln_full_vel = .false. and jpk_bdy>0?  ! 
     730      ! bathymetry in the case ln_totvel = .false. and ipkb>0?       ! 
    790731      ! [as the enveloping and parital cells could change H]         ! 
    791732      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    792733 
    793       CASE DEFAULT   ; 
    794  
    795       IF( PRESENT(jpk_bdy) .AND. jpk_bdy>0 ) THEN       ! boundary data not on model grid: vertical interpolation 
    796          CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar ) 
    797          dta_read(:,:,:) = -ABS(fv) 
    798          dta_read_z(:,:,:) = 0._wp 
    799          dta_read_dz(:,:,:) = 0._wp 
    800          CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:jpk_bdy), nrec ) 
    801          SELECT CASE( igrd )                   
    802             CASE(1) 
    803                CALL iom_get ( num, jpdom_unknown, 'gdept', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 
    804                CALL iom_get ( num, jpdom_unknown, 'e3t',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 
    805             CASE(2)   
    806                CALL iom_get ( num, jpdom_unknown, 'gdepu', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 
    807                CALL iom_get ( num, jpdom_unknown, 'e3u',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 
    808             CASE(3) 
    809                CALL iom_get ( num, jpdom_unknown, 'gdepv', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 
    810                CALL iom_get ( num, jpdom_unknown, 'e3v',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 
    811          END SELECT 
    812  
    813       IF ( ln_bdy ) &  
    814          CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta) 
    815  
    816       ELSE ! boundary data assumed to be on model grid 
    817           
    818          CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec )                     
    819          IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 
    820             DO ib = 1, ipi 
    821               DO ik = 1, ipk 
    822                 dta(ib,1,ik) =  dta_read(map%ptr(ib),1,ik) 
    823               END DO 
     734      ELSE 
     735         ! 
     736         CALL iom_get ( knum, jpdom_unknown, cdvar, zz_read(:,:,:), krec )   ! call iom_get with a 3D file 
     737         ! 
     738         IF( ipkb /= ipk ) THEN   ! boundary data not on model vertical grid : vertical interpolation 
     739            ! 
     740            IF( ipk == jpk .AND. iom_varid(knum,'gdep'//clgrid(kgrd)) /= -1 .AND. iom_varid(knum,'e3'//clgrid(kgrd)) /= -1 ) THEN 
     741                
     742               ALLOCATE( zdta_read(ipi,ipj,ipkb), zdta_read_z(ipi,ipj,ipkb), zdta_read_dz(ipi,ipj,ipkb) ) 
     743                 
     744               CALL fld_map_core( zz_read, kmap, zdta_read ) 
     745               CALL iom_get ( knum, jpdom_unknown, 'gdep'//clgrid(kgrd), zz_read )   ! read only once? Potential temporal evolution? 
     746               CALL fld_map_core( zz_read, kmap, zdta_read_z ) 
     747               CALL iom_get ( knum, jpdom_unknown,   'e3'//clgrid(kgrd), zz_read )   ! read only once? Potential temporal evolution? 
     748               CALL fld_map_core( zz_read, kmap, zdta_read_dz ) 
     749                
     750               CALL iom_getatt(knum, '_FillValue', zfv, cdvar=cdvar ) 
     751               CALL fld_bdy_interp(zdta_read, zdta_read_z, zdta_read_dz, pdta, kgrd, kbdy, zfv, ldtotvel) 
     752               DEALLOCATE( zdta_read, zdta_read_z, zdta_read_dz ) 
     753                
     754            ELSE 
     755               IF( ipk /= jpk ) CALL ctl_stop( 'fld_map : this should be an impossible case...' ) 
     756               WRITE(ctmp1,*) 'fld_map : vertical interpolation for bdy variable '//TRIM(cdvar)//' requires '  
     757               IF( iom_varid(knum, 'gdep'//clgrid(kgrd)) == -1 ) CALL ctl_stop( ctmp1//'gdep'//clgrid(kgrd)//' variable' ) 
     758               IF( iom_varid(knum,   'e3'//clgrid(kgrd)) == -1 ) CALL ctl_stop( ctmp1//  'e3'//clgrid(kgrd)//' variable' ) 
     759 
     760            ENDIF 
     761            ! 
     762         ELSE                            ! bdy data assumed to be the same levels as bdy variables 
     763            ! 
     764            CALL fld_map_core( zz_read, kmap, pdta ) 
     765            ! 
     766         ENDIF   ! ipkb /= ipk 
     767      ENDIF   ! ipk == 1 
     768       
     769      DEALLOCATE( zz_read ) 
     770 
     771   END SUBROUTINE fld_map 
     772 
     773      
     774   SUBROUTINE fld_map_core( pdta_read, kmap, pdta_bdy ) 
     775      !!--------------------------------------------------------------------- 
     776      !!                    ***  ROUTINE fld_map_core  *** 
     777      !! 
     778      !! ** Purpose :  inner core of fld_map 
     779      !!---------------------------------------------------------------------- 
     780      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read    ! global boundary data 
     781      INTEGER,  DIMENSION(:    ), INTENT(in   ) ::   kmap         ! global-to-local bdy mapping indices 
     782      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta_bdy     ! bdy output field on model grid 
     783      !! 
     784      INTEGER,  DIMENSION(3) ::   idim_read,  idim_bdy            ! arrays dimensions 
     785      INTEGER                ::   ji, jj, jk, jb                  ! loop counters 
     786      INTEGER                ::   im1 
     787      !!--------------------------------------------------------------------- 
     788      ! 
     789      idim_read = SHAPE( pdta_read ) 
     790      idim_bdy  = SHAPE( pdta_bdy  ) 
     791      ! 
     792      ! in all cases: idim_bdy(2) == 1 .AND. idim_read(1) * idim_read(2) == idim_bdy(1) 
     793      ! structured BDY with rimwidth > 1                     : idim_read(2) == rimwidth /= 1 
     794      ! structured BDY with rimwidth == 1 or unstructured BDY: idim_read(2) == 1 
     795      ! 
     796      IF( idim_read(2) > 1 ) THEN    ! structured BDY with rimwidth > 1   
     797         DO jk = 1, idim_bdy(3) 
     798            DO jb = 1, idim_bdy(1) 
     799               im1 = kmap(jb) - 1 
     800               jj = im1 / idim_read(1) + 1 
     801               ji = MOD( im1, idim_read(1) ) + 1 
     802               pdta_bdy(jb,1,jk) =  pdta_read(ji,jj,jk) 
    824803            END DO 
    825          ELSE ! we assume that this is a structured open boundary file 
    826             DO ib = 1, ipi 
    827                jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    828                ji=map%ptr(ib)-(jj-1)*ilendta 
    829                DO ik = 1, ipk 
    830                   dta(ib,1,ik) =  dta_read(ji,jj,ik) 
    831                END DO 
     804         END DO 
     805      ELSE 
     806         DO jk = 1, idim_bdy(3) 
     807            DO jb = 1, idim_bdy(1)   ! horizontal remap of bdy data on the local bdy  
     808               pdta_bdy(jb,1,jk) = pdta_read(kmap(jb),1,jk) 
    832809            END DO 
    833          ENDIF 
    834       ENDIF ! PRESENT(jpk_bdy) 
    835       END SELECT 
    836  
    837    END SUBROUTINE fld_map 
     810         END DO 
     811      ENDIF 
     812       
     813   END SUBROUTINE fld_map_core 
    838814    
    839    SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta) 
    840  
     815    
     816   SUBROUTINE fld_bdy_interp(pdta_read, pdta_read_z, pdta_read_dz, pdta, kgrd, kbdy, pfv, ldtotvel) 
    841817      !!--------------------------------------------------------------------- 
    842818      !!                    ***  ROUTINE fld_bdy_interp  *** 
    … …  
    847823      USE bdy_oce, ONLY:  idx_bdy         ! indexing for map <-> ij transformation 
    848824 
    849       REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read      ! work space for global data 
    850       REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_z    ! work space for global data 
    851       REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_dz   ! work space for global data 
    852       REAL(wp) , INTENT(in)                                ::   fv            ! fillvalue and alternative -ABS(fv) 
    853       REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta                        ! output field on model grid (2 dimensional) 
    854       TYPE(MAP_POINTER)         , INTENT(in ) ::   map                        ! global-to-local mapping indices 
    855       LOGICAL  , INTENT(in), OPTIONAL         ::   fvl                        ! grid type, set number and number of vertical levels in the bdy data 
    856       INTEGER  , INTENT(in)                   ::   igrd, ibdy, jpk_bdy        ! number of levels in bdy data 
    857       INTEGER  , INTENT(in)                   ::   ilendta                    ! length of data in file 
    858       !! 
    859       INTEGER                                 ::   ipi                        ! length of boundary data on local process 
    860       INTEGER                                 ::   ipj                        ! length of dummy dimension ( = 1 ) 
    861       INTEGER                                 ::   ipk                        ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
    862       INTEGER                                 ::   jpkm1_bdy                  ! number of levels in bdy data minus 1 
    863       INTEGER                                 ::   ib, ik, ikk                ! loop counters 
    864       INTEGER                                 ::   ji, jj, zij, zjj           ! temporary indices 
    865       REAL(wp)                                ::   zl, zi, zh                 ! tmp variable for current depth and interpolation factor 
    866       REAL(wp)                                ::   fv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(fv) 
    867       CHARACTER (LEN=10)                      ::   ibstr 
     825      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read       ! data read in bdy file 
     826      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdta_read_z     ! depth of the data read in bdy file 
     827      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdta_read_dz    ! thickness of the levels in bdy file 
     828      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta            ! output field on model grid (2 dimensional) 
     829      REAL(wp)                  , INTENT(in   ) ::   pfv             ! fillvalue of the data read in bdy file 
     830      LOGICAL                   , INTENT(in   ) ::   ldtotvel        ! true if toal ( = barotrop + barocline) velocity 
     831      INTEGER                   , INTENT(in   ) ::   kgrd            ! grid type (t, u, v) 
     832      INTEGER                   , INTENT(in   ) ::   kbdy            ! bdy number 
     833      !! 
     834      INTEGER                                   ::   ipi             ! length of boundary data on local process 
     835      INTEGER                                   ::   ipkb            ! number of vertical levels in boundary data file 
     836      INTEGER                                   ::   jb, ji, jj, jk, jkb   ! loop counters 
     837      REAL(wp)                                  ::   zcoef 
     838      REAL(wp)                                  ::   zl, zi, zh      ! tmp variable for current depth and interpolation factor 
     839      REAL(wp)                                  ::   zfv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(pfv) 
     840      REAL(wp), DIMENSION(jpk)                  ::   zdepth, zdhalf  ! level and half-level depth 
    868841      !!--------------------------------------------------------------------- 
    869842      
    870  
    871       ipi       = SIZE( dta, 1 ) 
    872       ipj       = SIZE( dta_read, 2 ) 
    873       ipk       = SIZE( dta, 3 ) 
    874       jpkm1_bdy = jpk_bdy-1 
     843      ipi  = SIZE( pdta, 1 ) 
     844      ipkb = SIZE( pdta_read, 3 ) 
    875845       
    876       fv_alt = -ABS(fv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 
    877       DO ib = 1, ipi 
    878             zij = idx_bdy(ibdy)%nbi(ib,igrd) 
    879             zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    880          IF(narea==2) WRITE(*,*) 'MAPI', ib, igrd, map%ptr(ib), narea-1, zij, zjj 
    881       ENDDO 
    882       ! 
    883       IF ( map%ll_unstruc ) THEN                            ! unstructured open boundary data file 
    884  
    885          DO ib = 1, ipi 
    886             DO ik = 1, jpk_bdy 
    887                IF( ( dta_read(map%ptr(ib),1,ik) == fv ) ) THEN 
    888                   dta_read_z(map%ptr(ib),1,ik)  = fv_alt ! safety: put fillvalue into external depth field so consistent with data 
    889                   dta_read_dz(map%ptr(ib),1,ik) = 0._wp  ! safety: put 0._wp into external thickness factors to ensure transport is correct 
     846      zfv_alt = -ABS(pfv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 
     847      ! 
     848      WHERE( pdta_read == pfv ) 
     849         pdta_read_z  = zfv_alt ! safety: put fillvalue into external depth field so consistent with data 
     850         pdta_read_dz = 0._wp   ! safety: put 0._wp into external thickness factors to ensure transport is correct 
     851      ENDWHERE 
     852       
     853      DO jb = 1, ipi 
     854         ji = idx_bdy(kbdy)%nbi(jb,kgrd) 
     855         jj = idx_bdy(kbdy)%nbj(jb,kgrd) 
     856         zh  = SUM(pdta_read_dz(jb,1,:) ) 
     857         ! 
     858         ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 
     859         SELECT CASE( kgrd )                          
     860         CASE(1) 
     861            IF( ABS( (zh - ht_n(ji,jj)) / ht_n(ji,jj)) * tmask(ji,jj,1) > 0.01_wp ) THEN 
     862               WRITE(ctmp1,"(I10.10)") jb  
     863               CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%') 
     864               !   IF(lwp) WRITE(numout,*) 'DEPTHT', zh, sum(e3t_n(ji,jj,:), mask=tmask(ji,jj,:)==1),  ht_n(ji,jj), jb, jb, ji, jj 
     865            ENDIF 
     866         CASE(2) 
     867            IF( ABS( (zh - hu_n(ji,jj)) * r1_hu_n(ji,jj)) * umask(ji,jj,1) > 0.01_wp ) THEN 
     868               WRITE(ctmp1,"(I10.10)") jb  
     869               CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%') 
     870               !   IF(lwp) WRITE(numout,*) 'DEPTHU', zh, SUM(e3u_n(ji,jj,:), mask=umask(ji,jj,:)==1),  SUM(umask(ji,jj,:)), & 
     871               !      &                hu_n(ji,jj), jb, jb, ji, jj, narea-1, pdta_read(jb,1,:) 
     872            ENDIF 
     873         CASE(3) 
     874            IF( ABS( (zh - hv_n(ji,jj)) * r1_hv_n(ji,jj)) * vmask(ji,jj,1) > 0.01_wp ) THEN 
     875               WRITE(ctmp1,"(I10.10)") jb 
     876               CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%') 
     877            ENDIF 
     878         END SELECT 
     879         ! 
     880         SELECT CASE( kgrd )                          
     881         CASE(1) 
     882            ! depth of T points: 
     883            zdepth(:) = gdept_n(ji,jj,:) 
     884         CASE(2) 
     885            ! depth of U points: we must not use gdept_n as we don't want to do a communication 
     886            !   --> copy what is done for gdept_n in domvvl... 
     887            zdhalf(1) = 0.0_wp 
     888            zdepth(1) = 0.5_wp * e3uw_n(ji,jj,1) 
     889            DO jk = 2, jpk                               ! vertical sum 
     890               !    zcoef = umask - wumask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     891               !                              ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
     892               !                              ! 0.5 where jk = mikt      
     893               !!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ?? 
     894               zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) ) 
     895               zdhalf(jk) = zdhalf(jk-1) + e3u_n(ji,jj,jk-1) 
     896               zdepth(jk) =      zcoef  * ( zdhalf(jk  ) + 0.5 * e3uw_n(ji,jj,jk))  & 
     897                  &         + (1-zcoef) * ( zdepth(jk-1) +       e3uw_n(ji,jj,jk)) 
     898            END DO 
     899         CASE(3) 
     900            ! depth of V points: we must not use gdept_n as we don't want to do a communication 
     901            !   --> copy what is done for gdept_n in domvvl... 
     902            zdhalf(1) = 0.0_wp 
     903            zdepth(1) = 0.5_wp * e3vw_n(ji,jj,1) 
     904            DO jk = 2, jpk                               ! vertical sum 
     905               !    zcoef = vmask - wvmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     906               !                              ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
     907               !                              ! 0.5 where jk = mikt      
     908               !!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ?? 
     909               zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) ) 
     910               zdhalf(jk) = zdhalf(jk-1) + e3v_n(ji,jj,jk-1) 
     911               zdepth(jk) =      zcoef  * ( zdhalf(jk  ) + 0.5 * e3vw_n(ji,jj,jk))  & 
     912                  &         + (1-zcoef) * ( zdepth(jk-1) +       e3vw_n(ji,jj,jk)) 
     913            END DO 
     914         END SELECT 
     915         ! 
     916         DO jk = 1, jpk                       
     917            IF(     zdepth(jk) < pdta_read_z(jb,1,          1) ) THEN                ! above the first level of external data 
     918               pdta(jb,1,jk) =  pdta_read(jb,1,1) 
     919            ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkb) ) THEN                       ! below the last level of external data  
     920               pdta(jb,1,jk) =  pdta_read(jb,1,MAXLOC(pdta_read_z(jb,1,:),1)) 
     921            ELSE                                                             ! inbetween: vertical interpolation between jkb & jkb+1 
     922               DO jkb = 1, ipkb-1                                            ! when  gdept_n(jkb) < zdepth(jk) < gdept_n(jkb+1) 
     923                  IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) <= 0._wp ) & 
     924                     &    .AND. ( pdta_read_z(jb,1,jkb+1) /= zfv_alt) ) THEN   ! linear interpolation between 2 levels 
     925                     zi = ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) / ( pdta_read_z(jb,1,jkb+1) - pdta_read_z(jb,1,jkb) ) 
     926                     pdta(jb,1,jk) = pdta_read(jb,1,jkb) + ( pdta_read  (jb,1,jkb+1) - pdta_read  (jb,1,jkb) ) * zi 
     927                  ENDIF 
     928               END DO 
     929            ENDIF 
     930         END DO   ! jpk 
     931         ! 
     932      END DO   ! ipi 
     933       
     934      IF(kgrd == 2) THEN                                  ! do we need to adjust the transport term? 
     935         DO jb = 1, ipi 
     936            ji = idx_bdy(kbdy)%nbi(jb,kgrd) 
     937            jj = idx_bdy(kbdy)%nbj(jb,kgrd) 
     938            zh  = SUM(pdta_read_dz(jb,1,:) ) 
     939            ztrans = 0._wp 
     940            ztrans_new = 0._wp 
     941            DO jkb = 1, ipkb                              ! calculate transport on input grid 
     942               ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb) 
     943            ENDDO 
     944            DO jk = 1, jpk                                ! calculate transport on model grid 
     945               ztrans_new = ztrans_new +      pdta(jb,1,jk ) *        e3u_n(ji,jj,jk ) * umask(ji,jj,jk) 
     946            ENDDO 
     947            DO jk = 1, jpk                                ! make transport correction 
     948               IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
     949                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu_n(ji,jj) ) * umask(ji,jj,jk) 
     950               ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
     951                  pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hu_n(ji,jj)   * umask(ji,jj,jk) 
    890952               ENDIF 
    891953            ENDDO 
    892          ENDDO  
    893  
    894          DO ib = 1, ipi 
    895             zij = idx_bdy(ibdy)%nbi(ib,igrd) 
    896             zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    897             zh  = SUM(dta_read_dz(map%ptr(ib),1,:) ) 
    898             ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 
    899             SELECT CASE( igrd )                          
    900                CASE(1) 
    901                   IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 
    902                      WRITE(ibstr,"(I10.10)") map%ptr(ib)  
    903                      CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
    904                  !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1),  ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj 
    905                   ENDIF 
    906                CASE(2) 
    907                   IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 
    908                      WRITE(ibstr,"(I10.10)") map%ptr(ib)  
    909                      CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
    910                      IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u_n(zij,zjj,:), mask=umask(zij,zjj,:)==1),  sum(umask(zij,zjj,:)), & 
    911                        &                hu_n(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1  , & 
    912                         &                dta_read(map%ptr(ib),1,:) 
    913                   ENDIF 
    914                CASE(3) 
    915                   IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 
    916                      WRITE(ibstr,"(I10.10)") map%ptr(ib)  
    917                      CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
    918                   ENDIF 
    919             END SELECT 
    920             DO ik = 1, ipk                       
    921                SELECT CASE( igrd )                        
    922                   CASE(1) 
    923                      zl =  gdept_n(zij,zjj,ik)                                          ! if using in step could use fsdept instead of gdept_n? 
    924                   CASE(2) 
    925                      IF(ln_sco) THEN 
    926                         zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n? 
    927                      ELSE 
    928                         zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) )  
    929                      ENDIF 
    930                   CASE(3) 
    931                      IF(ln_sco) THEN 
    932                         zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n? 
    933                      ELSE 
    934                         zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) ) 
    935                      ENDIF 
    936                END SELECT 
    937                IF( zl < dta_read_z(map%ptr(ib),1,1) ) THEN                                         ! above the first level of external data 
    938                   dta(ib,1,ik) =  dta_read(map%ptr(ib),1,1) 
    939                ELSEIF( zl > MAXVAL(dta_read_z(map%ptr(ib),1,:),1) ) THEN                           ! below the last level of external data  
    940                   dta(ib,1,ik) =  dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1)) 
    941                ELSE                                                                                ! inbetween : vertical interpolation between ikk & ikk+1 
    942                   DO ikk = 1, jpkm1_bdy                                                            ! when  gdept_n(ikk) < zl < gdept_n(ikk+1) 
    943                      IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) & 
    944                     &    .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN 
    945                         zi           = ( zl - dta_read_z(map%ptr(ib),1,ikk) ) / & 
    946                        &               ( dta_read_z(map%ptr(ib),1,ikk+1) - dta_read_z(map%ptr(ib),1,ikk) ) 
    947                         dta(ib,1,ik) = dta_read(map%ptr(ib),1,ikk) + & 
    948                        &               ( dta_read(map%ptr(ib),1,ikk+1) - dta_read(map%ptr(ib),1,ikk) ) * zi 
    949                      ENDIF 
    950                   END DO 
    951                ENDIF 
    952             END DO 
    953          END DO 
    954  
    955          IF(igrd == 2) THEN                                 ! do we need to adjust the transport term? 
    956             DO ib = 1, ipi 
    957               zij = idx_bdy(ibdy)%nbi(ib,igrd) 
    958               zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    959               zh  = SUM(dta_read_dz(map%ptr(ib),1,:) ) 
    960               ztrans = 0._wp 
    961               ztrans_new = 0._wp 
    962               DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
    963                   ztrans     = ztrans     + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 
    964               ENDDO 
    965               DO ik = 1, ipk                                ! calculate transport on model grid 
    966                   ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik) 
    967               ENDDO 
    968               DO ik = 1, ipk                                ! make transport correction 
    969                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    970                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 
    971                  ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
    972                     IF( ABS(ztrans * r1_hu_n(zij,zjj)) > 0.01_wp ) & 
    973                    &   CALL ctl_warn('fld_bdy_interp: barotropic component of > 0.01 ms-1 found in baroclinic velocities at') 
    974                     dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) * umask(zij,zjj,ik) 
    975                  ENDIF 
    976               ENDDO 
     954         ENDDO 
     955      ENDIF 
     956       
     957      IF(kgrd == 3) THEN                                  ! do we need to adjust the transport term? 
     958         DO jb = 1, ipi 
     959            ji = idx_bdy(kbdy)%nbi(jb,kgrd) 
     960            jj = idx_bdy(kbdy)%nbj(jb,kgrd) 
     961            zh  = SUM(pdta_read_dz(jb,1,:) ) 
     962            ztrans = 0._wp 
     963            ztrans_new = 0._wp 
     964            DO jkb = 1, ipkb                              ! calculate transport on input grid 
     965               ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb) 
    977966            ENDDO 
    978          ENDIF 
    979  
    980          IF(igrd == 3) THEN                                 ! do we need to adjust the transport term? 
    981             DO ib = 1, ipi 
    982               zij = idx_bdy(ibdy)%nbi(ib,igrd) 
    983               zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    984               zh  = SUM(dta_read_dz(map%ptr(ib),1,:) ) 
    985               ztrans = 0._wp 
    986               ztrans_new = 0._wp 
    987               DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
    988                   ztrans     = ztrans     + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 
    989               ENDDO 
    990               DO ik = 1, ipk                                ! calculate transport on model grid 
    991                   ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik) 
    992               ENDDO 
    993               DO ik = 1, ipk                                ! make transport correction 
    994                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    995                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 
    996                  ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
    997                     dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) * vmask(zij,zjj,ik) 
    998                  ENDIF 
    999               ENDDO 
     967            DO jk = 1, jpk                                ! calculate transport on model grid 
     968               ztrans_new = ztrans_new +      pdta(jb,1,jk ) *        e3v_n(ji,jj,jk ) * vmask(ji,jj,jk) 
    1000969            ENDDO 
    1001          ENDIF 
    1002    
    1003       ELSE ! structured open boundary file 
    1004  
    1005          DO ib = 1, ipi 
    1006             jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    1007             ji=map%ptr(ib)-(jj-1)*ilendta 
    1008             DO ik = 1, jpk_bdy                       
    1009                IF( ( dta_read(ji,jj,ik) == fv ) ) THEN 
    1010                   dta_read_z(ji,jj,ik)  = fv_alt ! safety: put fillvalue into external depth field so consistent with data 
    1011                   dta_read_dz(ji,jj,ik) = 0._wp  ! safety: put 0._wp into external thickness factors to ensure transport is correct 
     970            DO jk = 1, jpk                                ! make transport correction 
     971               IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
     972                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv_n(ji,jj) ) * vmask(ji,jj,jk) 
     973               ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
     974                  pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hv_n(ji,jj)   * vmask(ji,jj,jk) 
    1012975               ENDIF 
    1013976            ENDDO 
    1014          ENDDO  
    1015         
    1016  
    1017          DO ib = 1, ipi 
    1018             jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    1019             ji=map%ptr(ib)-(jj-1)*ilendta 
    1020             zij = idx_bdy(ibdy)%nbi(ib,igrd) 
    1021             zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    1022             zh  = SUM(dta_read_dz(ji,jj,:) ) 
    1023             ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 
    1024             SELECT CASE( igrd )                          
    1025                CASE(1) 
    1026                   IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 
    1027                      WRITE(ibstr,"(I10.10)") map%ptr(ib)  
    1028                      CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
    1029                  !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1),  ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj 
    1030                   ENDIF 
    1031                CASE(2) 
    1032                   IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 
    1033                      WRITE(ibstr,"(I10.10)") map%ptr(ib)  
    1034                      CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
    1035                   ENDIF 
    1036                CASE(3) 
    1037                   IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 
    1038                      WRITE(ibstr,"(I10.10)") map%ptr(ib)  
    1039                      CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 
    1040                   ENDIF 
    1041             END SELECT 
    1042             DO ik = 1, ipk                       
    1043                SELECT CASE( igrd )                          ! coded for sco - need zco and zps option using min 
    1044                   CASE(1) 
    1045                      zl =  gdept_n(zij,zjj,ik)              ! if using in step could use fsdept instead of gdept_n? 
    1046                   CASE(2) 
    1047                      IF(ln_sco) THEN 
    1048                         zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n? 
    1049                      ELSE 
    1050                         zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 
    1051                      ENDIF 
    1052                   CASE(3) 
    1053                      IF(ln_sco) THEN 
    1054                         zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n? 
    1055                      ELSE 
    1056                         zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) ) 
    1057                      ENDIF 
    1058                END SELECT 
    1059                IF( zl < dta_read_z(ji,jj,1) ) THEN                                      ! above the first level of external data 
    1060                   dta(ib,1,ik) =  dta_read(ji,jj,1) 
    1061                ELSEIF( zl > MAXVAL(dta_read_z(ji,jj,:),1) ) THEN                        ! below the last level of external data  
    1062                   dta(ib,1,ik) =  dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1)) 
    1063                ELSE                                                                     ! inbetween : vertical interpolation between ikk & ikk+1 
    1064                   DO ikk = 1, jpkm1_bdy                                                 ! when  gdept_n(ikk) < zl < gdept_n(ikk+1) 
    1065                      IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) & 
    1066                     &    .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN 
    1067                         zi           = ( zl - dta_read_z(ji,jj,ikk) ) / & 
    1068                        &               ( dta_read_z(ji,jj,ikk+1) - dta_read_z(ji,jj,ikk) ) 
    1069                         dta(ib,1,ik) = dta_read(ji,jj,ikk) + & 
    1070                        &               ( dta_read(ji,jj,ikk+1) - dta_read(ji,jj,ikk) ) * zi 
    1071                      ENDIF 
    1072                   END DO 
    1073                ENDIF 
    1074             END DO 
    1075          END DO 
    1076  
    1077          IF(igrd == 2) THEN                                 ! do we need to adjust the transport term? 
    1078             DO ib = 1, ipi 
    1079                jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    1080                ji=map%ptr(ib)-(jj-1)*ilendta 
    1081                zij = idx_bdy(ibdy)%nbi(ib,igrd) 
    1082                zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    1083                zh = SUM(dta_read_dz(ji,jj,:) ) 
    1084                ztrans = 0._wp 
    1085                ztrans_new = 0._wp 
    1086                DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
    1087                   ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 
    1088                ENDDO 
    1089                DO ik = 1, ipk                                ! calculate transport on model grid 
    1090                   ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik) 
    1091                ENDDO 
    1092                DO ik = 1, ipk                                ! make transport correction 
    1093                   IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    1094                      dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 
    1095                   ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
    1096                      dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 
    1097                   ENDIF 
    1098                ENDDO 
    1099             ENDDO 
    1100          ENDIF 
    1101  
    1102          IF(igrd == 3) THEN                                 ! do we need to adjust the transport term? 
    1103             DO ib = 1, ipi 
    1104                jj  = 1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 
    1105                ji  = map%ptr(ib)-(jj-1)*ilendta 
    1106                zij = idx_bdy(ibdy)%nbi(ib,igrd) 
    1107                zjj = idx_bdy(ibdy)%nbj(ib,igrd) 
    1108                zh  = SUM(dta_read_dz(ji,jj,:) ) 
    1109                ztrans = 0._wp 
    1110                ztrans_new = 0._wp 
    1111                DO ik = 1, jpk_bdy                            ! calculate transport on input grid 
    1112                   ztrans     = ztrans     + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 
    1113                ENDDO 
    1114                DO ik = 1, ipk                                ! calculate transport on model grid 
    1115                   ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik) 
    1116                ENDDO 
    1117                DO ik = 1, ipk                                ! make transport correction 
    1118                   IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    1119                      dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 
    1120                   ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
    1121                      dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 
    1122                   ENDIF 
    1123                ENDDO 
    1124             ENDDO 
    1125          ENDIF 
    1126  
    1127       ENDIF ! endif unstructured or structured 
    1128  
     977         ENDDO 
     978      ENDIF 
     979       
    1129980   END SUBROUTINE fld_bdy_interp 
    1130981 
    … …  
    11511002      imf = SIZE( sd ) 
    11521003      DO ju = 1, imf 
     1004         IF( TRIM(sd(ju)%clrootname) == 'NOT USED' )   CYCLE 
    11531005         ill = LEN_TRIM( sd(ju)%vcomp ) 
    11541006         DO jn = 2-COUNT((/sd(ju)%ln_tint/)), 2 
    … …  
    11591011                  iv = -1 
    11601012                  DO jv = 1, imf 
     1013                     IF( TRIM(sd(jv)%clrootname) == 'NOT USED' )   CYCLE 
    11611014                     IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) )   iv = jv 
    11621015                  END DO 
    … …  
    12991152      ! 
    13001153      DO jf = 1, SIZE(sdf) 
    1301          sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname ) 
     1154         sdf(jf)%clrootname = sdf_n(jf)%clname 
     1155         IF( TRIM(sdf_n(jf)%clname) /= 'NOT USED' )   sdf(jf)%clrootname = TRIM( cdir )//sdf(jf)%clrootname 
    13021156         sdf(jf)%clname     = "not yet defined" 
    13031157         sdf(jf)%nfreqh     = sdf_n(jf)%nfreqh 
    … …  
    13081162         sdf(jf)%num        = -1 
    13091163         sdf(jf)%wgtname    = " " 
    1310          IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname ) 
     1164         IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//sdf_n(jf)%wname 
    13111165         sdf(jf)%lsmname = " " 
    1312          IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 )   sdf(jf)%lsmname = TRIM( cdir )//TRIM( sdf_n(jf)%lname ) 
     1166         IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 )   sdf(jf)%lsmname = TRIM( cdir )//sdf_n(jf)%lname 
    13131167         sdf(jf)%vcomp      = sdf_n(jf)%vcomp 
    13141168         sdf(jf)%rotn(:)    = .TRUE.   ! pretend to be rotated -> won't try to rotate data before the first call to fld_get 
    … …  
    13171171         IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim )   & 
    13181172            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.') 
    1319          sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn 
     1173         sdf(jf)%nreclast   = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn 
     1174         sdf(jf)%igrd       = 0 
     1175         sdf(jf)%ibdy       = 0 
     1176         sdf(jf)%imap       => NULL() 
     1177         sdf(jf)%ltotvel    = .FALSE. 
    13201178      END DO 
    13211179      ! 

  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/updtide.F90

    r10068 r11223  
    2727CONTAINS 
    2828 
    29    SUBROUTINE upd_tide( kt, kit, time_offset ) 
     29   SUBROUTINE upd_tide( kt, kit, kt_offset ) 
    3030      !!---------------------------------------------------------------------- 
    3131      !!                 ***  ROUTINE upd_tide  *** 
    … …  
    3939      INTEGER, INTENT(in)           ::   kt      ! ocean time-step index 
    4040      INTEGER, INTENT(in), OPTIONAL ::   kit     ! external mode sub-time-step index (lk_dynspg_ts=T) 
    41       INTEGER, INTENT(in), OPTIONAL ::   time_offset ! time offset in number  
     41      INTEGER, INTENT(in), OPTIONAL ::   kt_offset ! time offset in number  
    4242                                                     ! of internal steps             (lk_dynspg_ts=F) 
    4343                                                     ! of external steps             (lk_dynspg_ts=T) 
    4444      ! 
    45       INTEGER  ::   joffset      ! local integer 
     45      INTEGER  ::   ioffset      ! local integer 
    4646      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    4747      REAL(wp) ::   zt, zramp    ! local scalar 
    … …  
    5252      zt = ( kt - kt_tide ) * rdt 
    5353      ! 
    54       joffset = 0 
    55       IF( PRESENT( time_offset ) )   joffset = time_offset 
     54      ioffset = 0 
     55      IF( PRESENT( kt_offset ) )   ioffset = kt_offset 
    5656      ! 
    5757      IF( PRESENT( kit ) )   THEN 
    58          zt = zt + ( kit +  joffset - 1 ) * rdt / REAL( nn_baro, wp ) 
     58         zt = zt + ( kit +  ioffset - 1 ) * rdt / REAL( nn_baro, wp ) 
    5959      ELSE 
    60          zt = zt + joffset * rdt 
     60         zt = zt + ioffset * rdt 
    6161      ENDIF 
    6262      ! 
    … …  
    7070      IF( ln_tide_ramp ) THEN         ! linear increase if asked 
    7171         zt = ( kt - nit000 ) * rdt 
    72          IF( PRESENT( kit ) )   zt = zt + ( kit + joffset -1) * rdt / REAL( nn_baro, wp ) 
     72         IF( PRESENT( kit ) )   zt = zt + ( kit + ioffset -1) * rdt / REAL( nn_baro, wp ) 
    7373         zramp = MIN(  MAX( zt / (rdttideramp*rday) , 0._wp ) , 1._wp  ) 
    7474         pot_astro(:,:) = zramp * pot_astro(:,:) 

  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/step.F90

    r10364 r11223  
    112112      IF( ln_tide    )   CALL sbc_tide( kstp )                   ! update tide potential 
    113113      IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                   ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)  
    114       IF( ln_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries 
     114      IF( ln_bdy     )   CALL bdy_dta ( kstp, kt_offset = +1 )   ! update dynamic & tracer data at open boundaries 
    115115                         CALL sbc     ( kstp )                   ! Sea Boundary Condition (including sea-ice) 
    116116 
Note: See TracChangeset for help on using the changeset viewer.