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 4292 for branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

Ignore:
Timestamp:
2013-11-20T17:28:04+01:00 (10 years ago)
Author:
cetlod
Message:

dev_MERGE_2013 : 1st step of the merge, see ticket #1185

Location:
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY
Files:
1 deleted
10 edited
2 copied

Legend:

Unmodified
Added
Removed

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r3651 r4292  
    88   !!            3.3  !  2010-09  (D. Storkey) add ice boundary conditions 
    99   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
     10   !!            3.6  !  2012-01  (C. Rousset) add ice boundary conditions for lim3 
    1011   !!---------------------------------------------------------------------- 
    1112#if defined key_bdy  
    … …  
    2728      INTEGER, POINTER, DIMENSION(:,:)   ::  nbr 
    2829      INTEGER, POINTER, DIMENSION(:,:)   ::  nbmap 
    29       REAL   , POINTER, DIMENSION(:,:)   ::  nbw 
    30       REAL   , POINTER, DIMENSION(:,:)   ::  nbd 
    31       REAL   , POINTER, DIMENSION(:)     ::  flagu 
    32       REAL   , POINTER, DIMENSION(:)     ::  flagv 
     30      REAL(wp)   , POINTER, DIMENSION(:,:)   ::  nbw 
     31      REAL(wp)   , POINTER, DIMENSION(:,:)   ::  nbd 
     32      REAL(wp)   , POINTER, DIMENSION(:,:)   ::  nbdout 
     33      REAL(wp)   , POINTER, DIMENSION(:,:)   ::  flagu 
     34      REAL(wp)   , POINTER, DIMENSION(:,:)   ::  flagv 
    3335   END TYPE OBC_INDEX 
    3436 
     37   !! Logicals in OBC_DATA structure are true if the chosen algorithm requires this 
     38   !! field as external data. If true the data can come from external files 
     39   !! or model initial conditions. If false then no "external" data array 
     40   !! is required for this field.  
     41 
    3542   TYPE, PUBLIC ::   OBC_DATA     !: Storage for external data 
    36       REAL, POINTER, DIMENSION(:)     ::  ssh 
    37       REAL, POINTER, DIMENSION(:)     ::  u2d 
    38       REAL, POINTER, DIMENSION(:)     ::  v2d 
    39       REAL, POINTER, DIMENSION(:,:)   ::  u3d 
    40       REAL, POINTER, DIMENSION(:,:)   ::  v3d 
    41       REAL, POINTER, DIMENSION(:,:)   ::  tem 
    42       REAL, POINTER, DIMENSION(:,:)   ::  sal 
     43      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      REAL(wp), POINTER, DIMENSION(:)     ::  ssh 
     52      REAL(wp), POINTER, DIMENSION(:)     ::  u2d 
     53      REAL(wp), POINTER, DIMENSION(:)     ::  v2d 
     54      REAL(wp), POINTER, DIMENSION(:,:)   ::  u3d 
     55      REAL(wp), POINTER, DIMENSION(:,:)   ::  v3d 
     56      REAL(wp), POINTER, DIMENSION(:,:)   ::  tem 
     57      REAL(wp), POINTER, DIMENSION(:,:)   ::  sal 
    4358#if defined key_lim2 
    44       REAL, POINTER, DIMENSION(:)     ::  frld 
    45       REAL, POINTER, DIMENSION(:)     ::  hicif 
    46       REAL, POINTER, DIMENSION(:)     ::  hsnif 
     59      LOGICAL                         ::  ll_frld 
     60      LOGICAL                         ::  ll_hicif 
     61      LOGICAL                         ::  ll_hsnif 
     62      REAL(wp), POINTER, DIMENSION(:)     ::  frld 
     63      REAL(wp), POINTER, DIMENSION(:)     ::  hicif 
     64      REAL(wp), POINTER, DIMENSION(:)     ::  hsnif 
     65#elif defined key_lim3 
     66      LOGICAL                         ::  ll_a_i 
     67      LOGICAL                         ::  ll_ht_i 
     68      LOGICAL                         ::  ll_ht_s 
     69      REAL, POINTER, DIMENSION(:,:)   ::  a_i   !: now ice leads fraction climatology 
     70      REAL, POINTER, DIMENSION(:,:)   ::  ht_i  !: Now ice  thickness climatology 
     71      REAL, POINTER, DIMENSION(:,:)   ::  ht_s  !: now snow thickness 
    4772#endif 
    4873   END TYPE OBC_DATA 
    … …  
    6388   INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P  
    6489   !                                                        !  = 1 the volume will be constant during all the integration. 
    65    INTEGER, DIMENSION(jp_bdy) ::   nn_dyn2d                 ! Choice of boundary condition for barotropic variables (U,V,SSH) 
    66    INTEGER, DIMENSION(jp_bdy) ::   nn_dyn2d_dta           !: = 0 use the initial state as bdy dta ;  
     90   CHARACTER(len=20), DIMENSION(jp_bdy) ::   cn_dyn2d       ! Choice of boundary condition for barotropic variables (U,V,SSH) 
     91   INTEGER, DIMENSION(jp_bdy)           ::   nn_dyn2d_dta   !: = 0 use the initial state as bdy dta ;  
    6792                                                            !: = 1 read it in a NetCDF file 
    6893                                                            !: = 2 read tidal harmonic forcing from a NetCDF file 
    6994                                                            !: = 3 read external data AND tidal harmonic forcing from NetCDF files 
    70    INTEGER, DIMENSION(jp_bdy) ::   nn_dyn3d                 ! Choice of boundary condition for baroclinic velocities  
    71    INTEGER, DIMENSION(jp_bdy) ::   nn_dyn3d_dta           !: = 0 use the initial state as bdy dta ;  
     95   CHARACTER(len=20), DIMENSION(jp_bdy) ::   cn_dyn3d       ! Choice of boundary condition for baroclinic velocities  
     96   INTEGER, DIMENSION(jp_bdy)           ::   nn_dyn3d_dta   !: = 0 use the initial state as bdy dta ;  
    7297                                                            !: = 1 read it in a NetCDF file 
    73    INTEGER, DIMENSION(jp_bdy) ::   nn_tra                   ! Choice of boundary condition for active tracers (T and S) 
    74    INTEGER, DIMENSION(jp_bdy) ::   nn_tra_dta             !: = 0 use the initial state as bdy dta ;  
     98   CHARACTER(len=20), DIMENSION(jp_bdy) ::   cn_tra         ! Choice of boundary condition for active tracers (T and S) 
     99   INTEGER, DIMENSION(jp_bdy)           ::   nn_tra_dta     !: = 0 use the initial state as bdy dta ;  
    75100                                                            !: = 1 read it in a NetCDF file 
    76101   LOGICAL, DIMENSION(jp_bdy) ::   ln_tra_dmp               !: =T Tracer damping 
    77102   LOGICAL, DIMENSION(jp_bdy) ::   ln_dyn3d_dmp             !: =T Baroclinic velocity damping 
    78    REAL,    DIMENSION(jp_bdy) ::   rn_time_dmp              !: Damping time scale in days 
     103   REAL(wp),    DIMENSION(jp_bdy) ::   rn_time_dmp              !: Damping time scale in days 
     104   REAL(wp),    DIMENSION(jp_bdy) ::   rn_time_dmp_out          !: Damping time scale in days at radiation outflow points 
    79105 
    80106#if defined key_lim2 
    81    INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim2              ! Choice of boundary condition for sea ice variables  
    82    INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim2_dta          !: = 0 use the initial state as bdy dta ;  
    83                                                             !: = 1 read it in a NetCDF file 
     107   CHARACTER(len=20), DIMENSION(jp_bdy) ::   nn_ice_lim2      ! Choice of boundary condition for sea ice variables  
     108   INTEGER, DIMENSION(jp_bdy)           ::   nn_ice_lim2_dta  !: = 0 use the initial state as bdy dta ;  
     109                                                              !: = 1 read it in a NetCDF file 
    84110#endif 
    85111   ! 
    … …  
    88114   !! Global variables 
    89115   !!---------------------------------------------------------------------- 
    90    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdytmask   !: Mask defining computational domain at T-points 
    91    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdyumask   !: Mask defining computational domain at U-points 
    92    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdyvmask   !: Mask defining computational domain at V-points 
     116   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   bdytmask   !: Mask defining computational domain at T-points 
     117   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   bdyumask   !: Mask defining computational domain at U-points 
     118   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   bdyvmask   !: Mask defining computational domain at V-points 
    93119 
    94120   REAL(wp)                                    ::   bdysurftot !: Lateral surface of unstructured open boundary 
    95121 
    96    REAL(wp), POINTER, DIMENSION(:,:)           ::   pssh       !:  
    97    REAL(wp), POINTER, DIMENSION(:,:)           ::   phur       !:  
    98    REAL(wp), POINTER, DIMENSION(:,:)           ::   phvr       !: Pointers for barotropic fields  
    99    REAL(wp), POINTER, DIMENSION(:,:)           ::   pu2d       !:  
    100    REAL(wp), POINTER, DIMENSION(:,:)           ::   pv2d       !:  
     122   REAL(wp), POINTER, DIMENSION(:,:)           ::   pssh                  !:  
     123   REAL(wp), POINTER, DIMENSION(:,:)           ::   phur                  !:  
     124   REAL(wp), POINTER, DIMENSION(:,:)           ::   phvr                  !: Pointers for barotropic fields  
     125   REAL(wp), POINTER, DIMENSION(:,:)           ::   pub2d, pun2d, pua2d   !:  
     126   REAL(wp), POINTER, DIMENSION(:,:)           ::   pvb2d, pvn2d, pva2d   !:  
    101127 
    102128   !!---------------------------------------------------------------------- 
    … …  
    109135   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2       !: workspace for reading in global data arrays (struct. bdy) 
    110136   TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process) 
    111    TYPE(OBC_DATA) , DIMENSION(jp_bdy)              ::   dta_bdy           !: bdy external data (local process) 
     137   TYPE(OBC_DATA) , DIMENSION(jp_bdy), TARGET      ::   dta_bdy           !: bdy external data (local process) 
    112138 
    113139   !!---------------------------------------------------------------------- 
    … …  
    125151      !!---------------------------------------------------------------------- 
    126152      ! 
    127       ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj),                    &   
     153      ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj),     &   
    128154         &      STAT=bdy_oce_alloc ) 
    129155         ! 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_par.F90

    r3294 r4292  
    2323# endif 
    2424   INTEGER, PUBLIC, PARAMETER ::   jp_bdy  = 10       !: Maximum number of bdy sets 
    25    INTEGER, PUBLIC, PARAMETER ::   jpbtime = 1000     !: Max number of time dumps per file 
    2625   INTEGER, PUBLIC, PARAMETER ::   jpbgrd  = 3        !: Number of horizontal grid types used  (T, U, V) 
    2726 
    28    !! Flags for choice of schemes 
    29    INTEGER, PUBLIC, PARAMETER ::   jp_none         = 0        !: Flag for no open boundary condition 
    30    INTEGER, PUBLIC, PARAMETER ::   jp_frs          = 1        !: Flag for Flow Relaxation Scheme 
    31    INTEGER, PUBLIC, PARAMETER ::   jp_flather      = 2        !: Flag for Flather 
    3227#else 
    3328   !!---------------------------------------------------------------------- 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r4230 r4292  
    1111   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions 
    1212   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
     13   !!            3.6  !  2012-01  (C. Rousset) add ice boundary conditions for lim3 
    1314   !!---------------------------------------------------------------------- 
    1415#if defined key_bdy 
    … …  
    2930   USE iom             ! IOM library 
    3031   USE in_out_manager  ! I/O logical units 
     32   USE dynspg_oce, ONLY: lk_dynspg_ts ! Split-explicit free surface flag 
    3133#if defined key_lim2 
    3234   USE ice_2 
     35#elif defined key_lim3 
     36   USE par_ice 
     37   USE ice 
     38   USE limcat_1D          ! redistribute ice input into categories 
    3339#endif 
    3440   USE sbcapr 
    … …  
    4955 
    5056   TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:) :: nbmap_ptr   ! array of pointers to nbmap 
     57 
     58#if defined key_lim3 
     59   LOGICAL :: ll_bdylim3                  ! determine whether ice input is lim2 (F) or lim3 (T) type 
     60   INTEGER :: jfld_hti, jfld_hts, jfld_ai ! indices of ice thickness, snow thickness and concentration in bf structure 
     61#endif 
    5162 
    5263#  include "domzgr_substitute.h90" 
    … …  
    7788                                                        ! etc. 
    7889      !! 
    79       INTEGER     ::  ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd  ! local indices 
     90      INTEGER     ::  ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd, jl  ! local indices 
    8091      INTEGER,          DIMENSION(jpbgrd) ::   ilen1  
    8192      INTEGER, POINTER, DIMENSION(:)      ::   nblen, nblenrim  ! short cuts 
     93      TYPE(OBC_DATA), POINTER             ::   dta              ! short cut 
    8294      !! 
    8395      !!--------------------------------------------------------------------------- 
    … …  
    91103         ! Calculate depth-mean currents 
    92104         !----------------------------- 
    93          CALL wrk_alloc(jpi,jpj,pu2d,pv2d)  
    94  
    95          pu2d(:,:) = 0.e0 
    96          pv2d(:,:) = 0.e0 
    97  
     105         CALL wrk_alloc(jpi,jpj,pun2d,pvn2d)  
     106 
     107         pun2d(:,:) = 0.e0 
     108         pvn2d(:,:) = 0.e0 
    98109         DO ik = 1, jpkm1   !! Vertically integrated momentum trends 
    99              pu2d(:,:) = pu2d(:,:) + fse3u(:,:,ik) * umask(:,:,ik) * un(:,:,ik) 
    100              pv2d(:,:) = pv2d(:,:) + fse3v(:,:,ik) * vmask(:,:,ik) * vn(:,:,ik) 
     110             pun2d(:,:) = pun2d(:,:) + fse3u(:,:,ik) * umask(:,:,ik) * un(:,:,ik) 
     111             pvn2d(:,:) = pvn2d(:,:) + fse3v(:,:,ik) * vmask(:,:,ik) * vn(:,:,ik) 
    101112         END DO 
    102          pu2d(:,:) = pu2d(:,:) * hur(:,:) 
    103          pv2d(:,:) = pv2d(:,:) * hvr(:,:) 
     113         pun2d(:,:) = pun2d(:,:) * hur(:,:) 
     114         pvn2d(:,:) = pvn2d(:,:) * hvr(:,:) 
    104115          
    105116         DO ib_bdy = 1, nb_bdy 
    … …  
    107118            nblen => idx_bdy(ib_bdy)%nblen 
    108119            nblenrim => idx_bdy(ib_bdy)%nblenrim 
    109  
    110             IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 0 ) THEN  
     120            dta => dta_bdy(ib_bdy) 
     121 
     122            IF( nn_dyn2d_dta(ib_bdy) .eq. 0 ) THEN  
    111123               ilen1(:) = nblen(:) 
    112                igrd = 1 
    113                DO ib = 1, ilen1(igrd) 
    114                   ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    115                   ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    116                   dta_bdy(ib_bdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)          
    117                END DO  
    118                igrd = 2 
    119                DO ib = 1, ilen1(igrd) 
    120                   ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    121                   ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    122                   dta_bdy(ib_bdy)%u2d(ib) = pu2d(ii,ij) * umask(ii,ij,1)          
    123                END DO  
    124                igrd = 3 
    125                DO ib = 1, ilen1(igrd) 
    126                   ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    127                   ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    128                   dta_bdy(ib_bdy)%v2d(ib) = pv2d(ii,ij) * vmask(ii,ij,1)          
    129                END DO  
    130             ENDIF 
    131  
    132             IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN  
    133                ilen1(:) = nblen(:) 
    134                igrd = 2  
    135                DO ib = 1, ilen1(igrd) 
    136                   DO ik = 1, jpkm1 
     124               IF( dta%ll_ssh ) THEN  
     125                  igrd = 1 
     126                  DO ib = 1, ilen1(igrd) 
    137127                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    138128                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    139                      dta_bdy(ib_bdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - pu2d(ii,ij) ) * umask(ii,ij,ik)          
    140                   END DO 
    141                END DO  
    142                igrd = 3  
    143                DO ib = 1, ilen1(igrd) 
    144                   DO ik = 1, jpkm1 
     129                     dta_bdy(ib_bdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)          
     130                  END DO  
     131               END IF 
     132               IF( dta%ll_u2d ) THEN  
     133                  igrd = 2 
     134                  DO ib = 1, ilen1(igrd) 
    145135                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    146136                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    147                      dta_bdy(ib_bdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - pv2d(ii,ij) ) * vmask(ii,ij,ik)          
    148                      END DO 
    149                END DO  
    150             ENDIF 
    151  
    152             IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 0 ) THEN  
    153                ilen1(:) = nblen(:) 
    154                igrd = 1                       ! Everything is at T-points here 
    155                DO ib = 1, ilen1(igrd) 
    156                   DO ik = 1, jpkm1 
     137                     dta_bdy(ib_bdy)%u2d(ib) = pun2d(ii,ij) * umask(ii,ij,1)          
     138                  END DO  
     139               END IF 
     140               IF( dta%ll_v2d ) THEN  
     141                  igrd = 3 
     142                  DO ib = 1, ilen1(igrd) 
    157143                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    158144                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    159                      dta_bdy(ib_bdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)          
    160                      dta_bdy(ib_bdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)          
     145                     dta_bdy(ib_bdy)%v2d(ib) = pvn2d(ii,ij) * vmask(ii,ij,1)          
     146                  END DO  
     147               END IF 
     148            ENDIF 
     149 
     150            IF( nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN  
     151               ilen1(:) = nblen(:) 
     152               IF( dta%ll_u3d ) THEN  
     153                  igrd = 2  
     154                  DO ib = 1, ilen1(igrd) 
     155                     DO ik = 1, jpkm1 
     156                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     157                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     158                        dta_bdy(ib_bdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - pun2d(ii,ij) ) * umask(ii,ij,ik)          
     159                     END DO 
     160                  END DO  
     161               END IF 
     162               IF( dta%ll_v3d ) THEN  
     163                  igrd = 3  
     164                  DO ib = 1, ilen1(igrd) 
     165                     DO ik = 1, jpkm1 
     166                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     167                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     168                        dta_bdy(ib_bdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - pvn2d(ii,ij) ) * vmask(ii,ij,ik)          
     169                        END DO 
     170                  END DO  
     171               END IF 
     172            ENDIF 
     173 
     174            IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN  
     175               ilen1(:) = nblen(:) 
     176               IF( dta%ll_tem ) THEN 
     177                  igrd = 1  
     178                  DO ib = 1, ilen1(igrd) 
     179                     DO ik = 1, jpkm1 
     180                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     181                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     182                        dta_bdy(ib_bdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)          
     183                     END DO 
     184                  END DO  
     185               END IF 
     186               IF( dta%ll_sal ) THEN 
     187                  igrd = 1  
     188                  DO ib = 1, ilen1(igrd) 
     189                     DO ik = 1, jpkm1 
     190                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     191                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     192                        dta_bdy(ib_bdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)          
     193                     END DO 
     194                  END DO  
     195               END IF 
     196            ENDIF 
     197 
     198#if defined key_lim2 
     199            IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN  
     200               ilen1(:) = nblen(:) 
     201               IF( dta%ll_frld ) THEN 
     202                  igrd = 1  
     203                  DO ib = 1, ilen1(igrd) 
     204                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     205                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     206                     dta_bdy(ib_bdy)%frld(ib) = frld(ii,ij) * tmask(ii,ij,1)          
     207                  END DO  
     208               END IF 
     209               IF( dta%ll_hicif ) THEN 
     210                  igrd = 1  
     211                  DO ib = 1, ilen1(igrd) 
     212                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     213                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     214                     dta_bdy(ib_bdy)%hicif(ib) = hicif(ii,ij) * tmask(ii,ij,1)          
     215                  END DO  
     216               END IF 
     217               IF( dta%ll_hsnif ) THEN 
     218                  igrd = 1  
     219                  DO ib = 1, ilen1(igrd) 
     220                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     221                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     222                     dta_bdy(ib_bdy)%hsnif(ib) = hsnif(ii,ij) * tmask(ii,ij,1)          
     223                  END DO  
     224               END IF 
     225            ENDIF 
     226#elif defined key_lim3 
     227            IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN  
     228               ilen1(:) = nblen(:) 
     229               IF( dta%ll_a_i ) THEN 
     230                  igrd = 1    
     231                  DO jl = 1, jpl 
     232                     DO ib = 1, ilen1(igrd) 
     233                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     234                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     235                        dta_bdy(ib_bdy)%a_i (ib,jl) =  a_i(ii,ij,jl) * tmask(ii,ij,1)  
     236                     END DO 
    161237                  END DO 
    162                END DO  
    163             ENDIF 
    164  
    165 #if defined key_lim2 
    166             IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN  
    167                ilen1(:) = nblen(:) 
    168                igrd = 1                       ! Everything is at T-points here 
    169                DO ib = 1, ilen1(igrd) 
    170                   ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    171                   ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    172                   dta_bdy(ib_bdy)%frld(ib) = frld(ii,ij) * tmask(ii,ij,1)          
    173                   dta_bdy(ib_bdy)%hicif(ib) = hicif(ii,ij) * tmask(ii,ij,1)          
    174                   dta_bdy(ib_bdy)%hsnif(ib) = hsnif(ii,ij) * tmask(ii,ij,1)          
    175                END DO  
     238               ENDIF 
     239               IF( dta%ll_ht_i ) THEN 
     240                  igrd = 1    
     241                  DO jl = 1, jpl 
     242                     DO ib = 1, ilen1(igrd) 
     243                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     244                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     245                        dta_bdy(ib_bdy)%ht_i (ib,jl) =  ht_i(ii,ij,jl) * tmask(ii,ij,1)  
     246                     END DO 
     247                  END DO 
     248               ENDIF 
     249               IF( dta%ll_ht_s ) THEN 
     250                  igrd = 1    
     251                  DO jl = 1, jpl 
     252                     DO ib = 1, ilen1(igrd) 
     253                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     254                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     255                        dta_bdy(ib_bdy)%ht_s (ib,jl) =  ht_s(ii,ij,jl) * tmask(ii,ij,1)  
     256                     END DO 
     257                  END DO 
     258               ENDIF 
    176259            ENDIF 
    177260#endif 
    … …  
    179262         ENDDO ! ib_bdy 
    180263 
    181          CALL wrk_dealloc(jpi,jpj,pu2d,pv2d)  
     264         CALL wrk_dealloc(jpi,jpj,pun2d,pvn2d)  
    182265 
    183266      ENDIF ! kt .eq. nit000 
    … …  
    188271      jstart = 1 
    189272      DO ib_bdy = 1, nb_bdy    
     273         dta => dta_bdy(ib_bdy) 
    190274         IF( nn_dta(ib_bdy) .eq. 1 ) THEN ! skip this bit if no external data required 
    191275       
    … …  
    193277               ! Update barotropic boundary conditions only 
    194278               ! jit is optional argument for fld_read and bdytide_update 
    195                IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN 
     279               IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN 
    196280                  IF( nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
    197                      dta_bdy(ib_bdy)%ssh(:) = 0.0 
    198                      dta_bdy(ib_bdy)%u2d(:) = 0.0 
    199                      dta_bdy(ib_bdy)%v2d(:) = 0.0 
     281                     IF( dta%ll_ssh ) dta%ssh(:) = 0.0 
     282                     IF( dta%ll_u2d ) dta%u2d(:) = 0.0 
     283                     IF( dta%ll_u3d ) dta%v2d(:) = 0.0 
    200284                  ENDIF 
    201                   IF (nn_tra(ib_bdy).ne.4) THEN 
    202                      IF( nn_dyn2d_dta(ib_bdy) .EQ. 1 .OR. nn_dyn2d_dta(ib_bdy) .EQ. 3 .OR.  & 
    203                        & (ln_full_vel_array(ib_bdy) .AND. nn_dyn3d_dta(ib_bdy).eq.1) )THEN 
    204  
    205                         ! For the runoff case, no need to update the forcing (already done in the baroclinic part) 
    206                         jend = nb_bdy_fld(ib_bdy) 
    207                         IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend - 2 
     285                  IF (cn_tra(ib_bdy) /= 'runoff') THEN 
     286                     IF( nn_dyn2d_dta(ib_bdy) .EQ. 1 .OR. nn_dyn2d_dta(ib_bdy) .EQ. 3 ) THEN 
     287 
     288                        jend = jstart + dta%nread(2) - 1 
    208289                        CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  & 
    209290                                     & kit=jit, kt_offset=time_offset ) 
    210                         IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend + 2 
    211  
    212                         ! If full velocities in boundary data then split into barotropic and baroclinic data 
     291 
     292                        ! If full velocities in boundary data then extract barotropic velocities from 3D fields 
    213293                        IF( ln_full_vel_array(ib_bdy) .AND.                                             & 
    214294                          &    ( nn_dyn2d_dta(ib_bdy) .EQ. 1 .OR. nn_dyn2d_dta(ib_bdy) .EQ. 3 .OR.  & 
    … …  
    216296 
    217297                           igrd = 2                      ! zonal velocity 
    218                            dta_bdy(ib_bdy)%u2d(:) = 0.0 
     298                           dta%u2d(:) = 0.0 
    219299                           DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    220300                              ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    221301                              ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    222302                              DO ik = 1, jpkm1 
    223                                  dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) & 
    224                        &                          + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik) 
     303                                 dta%u2d(ib) = dta%u2d(ib) & 
     304                       &                          + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik) 
    225305                              END DO 
    226                               dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 
    227                               DO ik = 1, jpkm1 
    228                                  dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib) 
    229                               END DO 
     306                              dta%u2d(ib) =  dta%u2d(ib) * hur(ii,ij) 
    230307                           END DO 
    231308                           igrd = 3                      ! meridional velocity 
    232                            dta_bdy(ib_bdy)%v2d(:) = 0.0 
     309                           dta%v2d(:) = 0.0 
    233310                           DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    234311                              ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    235312                              ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    236313                              DO ik = 1, jpkm1 
    237                                  dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) & 
    238                        &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik) 
     314                                 dta%v2d(ib) = dta%v2d(ib) & 
     315                       &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik) 
    239316                              END DO 
    240                               dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 
    241                               DO ik = 1, jpkm1 
    242                                  dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib) 
    243                               END DO 
     317                              dta%v2d(ib) =  dta%v2d(ib) * hvr(ii,ij) 
    244318                           END DO 
    245319                        ENDIF                     
    246320                     ENDIF 
    247321                     IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
    248                         CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy),   &  
     322                        CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta, td=tides(ib_bdy),   &  
    249323                          &                 jit=jit, time_offset=time_offset ) 
    250324                     ENDIF 
    … …  
    252326               ENDIF 
    253327            ELSE 
    254                IF (nn_tra(ib_bdy).eq.4) then      ! runoff condition 
     328               IF (cn_tra(ib_bdy) == 'runoff') then      ! runoff condition 
    255329                  jend = nb_bdy_fld(ib_bdy) 
    256330                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend),  & 
    … …  
    261335                     ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    262336                     ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    263                      dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) ) 
     337                     dta%u2d(ib) = dta%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) ) 
    264338                  END DO 
    265339                  ! 
    … …  
    268342                     ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    269343                     ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    270                      dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) ) 
     344                     dta%v2d(ib) = dta%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) ) 
    271345                  END DO 
    272346               ELSE 
    273                   IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
    274                      dta_bdy(ib_bdy)%ssh(:) = 0.0 
    275                      dta_bdy(ib_bdy)%u2d(:) = 0.0 
    276                      dta_bdy(ib_bdy)%v2d(:) = 0.0 
     347                  IF( nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
     348                     IF( dta%ll_ssh ) dta%ssh(:) = 0.0 
     349                     IF( dta%ll_u2d ) dta%u2d(:) = 0.0 
     350                     IF( dta%ll_v2d ) dta%v2d(:) = 0.0 
    277351                  ENDIF 
    278                   IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data 
    279                      jend = nb_bdy_fld(ib_bdy) 
     352                  IF( dta%nread(1) .gt. 0 ) THEN ! update external data 
     353                     jend = jstart + dta%nread(1) - 1 
    280354                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), & 
    281355                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset ) 
    … …  
    286360                    &   nn_dyn3d_dta(ib_bdy) .EQ. 1 ) ) THEN 
    287361                     igrd = 2                      ! zonal velocity 
    288                      dta_bdy(ib_bdy)%u2d(:) = 0.0 
     362                     dta%u2d(:) = 0.0 
    289363                     DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    290364                        ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    291365                        ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    292366                        DO ik = 1, jpkm1 
    293                            dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) & 
    294                  &                       + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik) 
     367                           dta%u2d(ib) = dta%u2d(ib) & 
     368                 &                       + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik) 
    295369                        END DO 
    296                         dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 
     370                        dta%u2d(ib) =  dta%u2d(ib) * hur(ii,ij) 
    297371                        DO ik = 1, jpkm1 
    298                            dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib) 
     372                           dta%u3d(ib,ik) = dta%u3d(ib,ik) - dta%u2d(ib) 
    299373                        END DO 
    300374                     END DO 
    301375                     igrd = 3                      ! meridional velocity 
    302                      dta_bdy(ib_bdy)%v2d(:) = 0.0 
     376                     dta%v2d(:) = 0.0 
    303377                     DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    304378                        ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    305379                        ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    306380                        DO ik = 1, jpkm1 
    307                            dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) & 
    308                  &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik) 
     381                           dta%v2d(ib) = dta%v2d(ib) & 
     382                 &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik) 
    309383                        END DO 
    310                         dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 
     384                        dta%v2d(ib) =  dta%v2d(ib) * hvr(ii,ij) 
    311385                        DO ik = 1, jpkm1 
    312                            dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib) 
     386                           dta%v3d(ib,ik) = dta%v3d(ib,ik) - dta%v2d(ib) 
    313387                        END DO 
    314388                     END DO 
    315389                  ENDIF 
    316                   IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
    317                      CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy),  & 
    318                                         & td=tides(ib_bdy), time_offset=time_offset ) 
    319                   ENDIF 
    320                ENDIF 
    321             ENDIF 
    322             jstart = jend+1 
     390 
     391               ENDIF 
     392#if defined key_lim3 
     393               IF( .NOT. ll_bdylim3 .AND. nn_ice_lim(ib_bdy) > 0 .AND. nn_ice_lim_dta(ib_bdy) == 1 ) THEN ! bdy ice input (case input is lim2 type) 
     394                CALL lim_cat_1D ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), & 
     395                                  & dta_bdy(ib_bdy)%ht_i,     dta_bdy(ib_bdy)%ht_s,     dta_bdy(ib_bdy)%a_i     ) 
     396               ENDIF 
     397#endif 
     398            ENDIF 
     399            jstart = jstart + dta%nread(1) 
    323400         END IF ! nn_dta(ib_bdy) = 1 
    324401      END DO  ! ib_bdy 
    325402 
     403      ! bg jchanut tschanges 
     404#if defined key_tide 
     405      ! Add tides if not split-explicit free surface else this is done in ts loop 
     406      IF (.NOT.lk_dynspg_ts) CALL bdy_dta_tides( kt=kt, time_offset=time_offset ) 
     407#endif 
     408      ! end jchanut tschanges 
     409 
    326410      IF ( ln_apr_obc ) THEN 
    327411         DO ib_bdy = 1, nb_bdy 
    328             IF (nn_tra(ib_bdy).NE.4)THEN 
     412            IF (cn_tra(ib_bdy) /= 'runoff')THEN 
    329413               igrd = 1                      ! meridional velocity 
    330414               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    … …  
    349433      !!                for open boundary conditions 
    350434      !! 
    351       !! ** Method  :   Use fldread.F90 
     435      !! ** Method  :    
    352436      !!                 
    353437      !!---------------------------------------------------------------------- 
    … …  
    362446                                                                ! =F => baroclinic velocities in 3D boundary data 
    363447      INTEGER                                ::   ilen_global   ! Max length required for global bdy dta arrays 
    364       INTEGER,              DIMENSION(jpbgrd) ::  ilen0         ! size of local arrays 
    365448      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ilen1, ilen3  ! size of 1st and 3rd dimensions of local arrays 
    366449      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ibdy           ! bdy set for a particular jfld 
    367450      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   igrid         ! index for grid type (1,2,3 = T,U,V) 
    368451      INTEGER, POINTER, DIMENSION(:)         ::   nblen, nblenrim  ! short cuts 
     452      TYPE(OBC_DATA), POINTER                ::   dta           ! short cut 
     453#if defined key_lim3 
     454      INTEGER, DIMENSION(3) ::   zdimsz   ! number of elements in each of the 4 dimensions (i.e. i,j,t,ice-cat) for an array 
     455      INTEGER               ::   zndims   ! number of dimensions in an array (i.e. 3 = wo ice cat; 4 = w ice cat) 
     456      INTEGER               ::   inum,id1 ! local integer 
     457#endif 
    369458      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   blf_i         !  array of namelist information structures 
    370459      TYPE(FLD_N) ::   bn_tem, bn_sal, bn_u3d, bn_v3d   !  
    … …  
    372461#if defined key_lim2 
    373462      TYPE(FLD_N) ::   bn_frld, bn_hicif, bn_hsnif      ! 
     463#elif defined key_lim3 
     464      TYPE(FLD_N) ::   bn_a_i, bn_ht_i, bn_ht_s       
    374465#endif 
    375466      NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d  
    376467#if defined key_lim2 
    377468      NAMELIST/nambdy_dta/ bn_frld, bn_hicif, bn_hsnif 
     469#elif defined key_lim3 
     470      NAMELIST/nambdy_dta/ bn_a_i, bn_ht_i, bn_ht_s 
    378471#endif 
    379472      NAMELIST/nambdy_dta/ ln_full_vel 
    … …  
    392485                               ,nn_dyn3d_dta(ib_bdy)       & 
    393486                               ,nn_tra_dta(ib_bdy)         & 
    394 #if defined key_lim2 
    395                                ,nn_ice_lim2_dta(ib_bdy)    & 
     487#if ( defined key_lim2 || defined key_lim3 ) 
     488                              ,nn_ice_lim_dta(ib_bdy)    & 
    396489#endif 
    397490                              ) 
    … …  
    404497      nb_bdy_fld(:) = 0 
    405498      DO ib_bdy = 1, nb_bdy          
    406          IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN 
     499         IF( cn_dyn2d(ib_bdy) /= 'none' .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN 
    407500            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3 
    408501         ENDIF 
    409          IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) THEN 
     502         IF( cn_dyn3d(ib_bdy) /= 'none' .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) THEN 
    410503            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2 
    411504         ENDIF 
    412          IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 1  ) THEN 
     505         IF( cn_tra(ib_bdy) /= 'none' .and. nn_tra_dta(ib_bdy) .eq. 1  ) THEN 
    413506            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2 
    414507         ENDIF 
    415 #if defined key_lim2 
    416          IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1  ) THEN 
     508#if ( defined key_lim2 || defined key_lim3 ) 
     509         IF( cn_ice_lim(ib_bdy) /= 'none' .and. nn_ice_lim_dta(ib_bdy) .eq. 1  ) THEN 
    417510            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3 
    418511         ENDIF 
    … …  
    458551            nblen => idx_bdy(ib_bdy)%nblen 
    459552            nblenrim => idx_bdy(ib_bdy)%nblenrim 
     553            dta => dta_bdy(ib_bdy) 
     554            dta%nread(2) = 0 
    460555 
    461556            ! Only read in necessary fields for this set. 
    462557            ! Important that barotropic variables come first. 
    463             IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN  
    464  
    465                IF( nn_dyn2d(ib_bdy) .ne. jp_frs .and. nn_tra(ib_bdy) .ne. 4 ) THEN ! runoff condition : no ssh reading 
     558            IF( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN  
     559 
     560               IF( dta%ll_ssh ) THEN  
     561                  if(lwp) write(numout,*) '++++++ reading in ssh field' 
    466562                  jfld = jfld + 1 
    467563                  blf_i(jfld) = bn_ssh 
    … …  
    470566                  ilen1(jfld) = nblen(igrid(jfld)) 
    471567                  ilen3(jfld) = 1 
    472                ENDIF 
    473  
    474                IF( .not. ln_full_vel_array(ib_bdy) ) THEN 
     568                  dta%nread(2) = dta%nread(2) + 1 
     569               ENDIF 
     570 
     571               IF( dta%ll_u2d .and. .not. ln_full_vel_array(ib_bdy) ) THEN 
     572                  if(lwp) write(numout,*) '++++++ reading in u2d field' 
    475573                  jfld = jfld + 1 
    476574                  blf_i(jfld) = bn_u2d 
    … …  
    479577                  ilen1(jfld) = nblen(igrid(jfld)) 
    480578                  ilen3(jfld) = 1 
    481  
     579                  dta%nread(2) = dta%nread(2) + 1 
     580               ENDIF 
     581 
     582               IF( dta%ll_v2d .and. .not. ln_full_vel_array(ib_bdy) ) THEN 
     583                  if(lwp) write(numout,*) '++++++ reading in v2d field' 
    482584                  jfld = jfld + 1 
    483585                  blf_i(jfld) = bn_v2d 
    … …  
    486588                  ilen1(jfld) = nblen(igrid(jfld)) 
    487589                  ilen3(jfld) = 1 
    488                ENDIF 
    489  
    490             ENDIF 
    491  
    492             ! baroclinic velocities 
    493             IF( ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) .or. & 
    494            &      ( ln_full_vel_array(ib_bdy) .and. nn_dyn2d(ib_bdy) .gt. 0 .and.  & 
    495            &        ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN 
    496  
    497                jfld = jfld + 1 
    498                blf_i(jfld) = bn_u3d 
    499                ibdy(jfld) = ib_bdy 
    500                igrid(jfld) = 2 
    501                ilen1(jfld) = nblen(igrid(jfld)) 
    502                ilen3(jfld) = jpk 
    503  
    504                jfld = jfld + 1 
    505                blf_i(jfld) = bn_v3d 
    506                ibdy(jfld) = ib_bdy 
    507                igrid(jfld) = 3 
    508                ilen1(jfld) = nblen(igrid(jfld)) 
    509                ilen3(jfld) = jpk 
     590                  dta%nread(2) = dta%nread(2) + 1 
     591               ENDIF 
     592 
     593            ENDIF 
     594 
     595            ! read 3D velocities if baroclinic velocities require OR if 
     596            ! barotropic velocities required and ln_full_vel set to .true. 
     597            IF( nn_dyn3d_dta(ib_bdy) .eq. 1 .or. & 
     598           &  ( ln_full_vel_array(ib_bdy) .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN 
     599 
     600               IF( dta%ll_u3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) ) THEN 
     601                  if(lwp) write(numout,*) '++++++ reading in u3d field' 
     602                  jfld = jfld + 1 
     603                  blf_i(jfld) = bn_u3d 
     604                  ibdy(jfld) = ib_bdy 
     605                  igrid(jfld) = 2 
     606                  ilen1(jfld) = nblen(igrid(jfld)) 
     607                  ilen3(jfld) = jpk 
     608                  IF( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) dta%nread(2) = dta%nread(2) + 1 
     609               ENDIF 
     610 
     611               IF( dta%ll_v3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) ) THEN 
     612                  if(lwp) write(numout,*) '++++++ reading in v3d field' 
     613                  jfld = jfld + 1 
     614                  blf_i(jfld) = bn_v3d 
     615                  ibdy(jfld) = ib_bdy 
     616                  igrid(jfld) = 3 
     617                  ilen1(jfld) = nblen(igrid(jfld)) 
     618                  ilen3(jfld) = jpk 
     619                  IF( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) dta%nread(2) = dta%nread(2) + 1 
     620               ENDIF 
    510621 
    511622            ENDIF 
    512623 
    513624            ! temperature and salinity 
    514             IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 1 ) THEN 
    515  
    516                jfld = jfld + 1 
    517                blf_i(jfld) = bn_tem 
    518                ibdy(jfld) = ib_bdy 
    519                igrid(jfld) = 1 
    520                ilen1(jfld) = nblen(igrid(jfld)) 
    521                ilen3(jfld) = jpk 
    522  
    523                jfld = jfld + 1 
    524                blf_i(jfld) = bn_sal 
    525                ibdy(jfld) = ib_bdy 
    526                igrid(jfld) = 1 
    527                ilen1(jfld) = nblen(igrid(jfld)) 
    528                ilen3(jfld) = jpk 
     625            IF( nn_tra_dta(ib_bdy) .eq. 1 ) THEN 
     626 
     627               IF( dta%ll_tem ) THEN 
     628                  if(lwp) write(numout,*) '++++++ reading in tem field' 
     629                  jfld = jfld + 1 
     630                  blf_i(jfld) = bn_tem 
     631                  ibdy(jfld) = ib_bdy 
     632                  igrid(jfld) = 1 
     633                  ilen1(jfld) = nblen(igrid(jfld)) 
     634                  ilen3(jfld) = jpk 
     635               ENDIF 
     636 
     637               IF( dta%ll_sal ) THEN 
     638                  if(lwp) write(numout,*) '++++++ reading in sal field' 
     639                  jfld = jfld + 1 
     640                  blf_i(jfld) = bn_sal 
     641                  ibdy(jfld) = ib_bdy 
     642                  igrid(jfld) = 1 
     643                  ilen1(jfld) = nblen(igrid(jfld)) 
     644                  ilen3(jfld) = jpk 
     645               ENDIF 
    529646 
    530647            ENDIF 
    … …  
    532649#if defined key_lim2 
    533650            ! sea ice 
    534             IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1 ) THEN 
    535  
    536                jfld = jfld + 1 
    537                blf_i(jfld) = bn_frld 
    538                ibdy(jfld) = ib_bdy 
    539                igrid(jfld) = 1 
    540                ilen1(jfld) = nblen(igrid(jfld)) 
    541                ilen3(jfld) = 1 
    542  
    543                jfld = jfld + 1 
    544                blf_i(jfld) = bn_hicif 
    545                ibdy(jfld) = ib_bdy 
    546                igrid(jfld) = 1 
    547                ilen1(jfld) = nblen(igrid(jfld)) 
    548                ilen3(jfld) = 1 
    549  
    550                jfld = jfld + 1 
    551                blf_i(jfld) = bn_hsnif 
    552                ibdy(jfld) = ib_bdy 
    553                igrid(jfld) = 1 
    554                ilen1(jfld) = nblen(igrid(jfld)) 
    555                ilen3(jfld) = 1 
    556  
    557             ENDIF 
     651            IF( nn_ice_lim2_dta(ib_bdy) .eq. 1 ) THEN 
     652 
     653               IF( dta%ll_frld ) THEN 
     654                  jfld = jfld + 1 
     655                  blf_i(jfld) = bn_frld 
     656                  ibdy(jfld) = ib_bdy 
     657                  igrid(jfld) = 1 
     658                  ilen1(jfld) = nblen(igrid(jfld)) 
     659                  ilen3(jfld) = 1 
     660               ENDIF 
     661 
     662               IF( dta%ll_hicif ) THEN 
     663                  jfld = jfld + 1 
     664                  blf_i(jfld) = bn_hicif 
     665                  ibdy(jfld) = ib_bdy 
     666                  igrid(jfld) = 1 
     667                  ilen1(jfld) = nblen(igrid(jfld)) 
     668                  ilen3(jfld) = 1 
     669               ENDIF 
     670 
     671               IF( dta%ll_hsnif ) THEN 
     672                  jfld = jfld + 1 
     673                  blf_i(jfld) = bn_hsnif 
     674                  ibdy(jfld) = ib_bdy 
     675                  igrid(jfld) = 1 
     676                  ilen1(jfld) = nblen(igrid(jfld)) 
     677                  ilen3(jfld) = 1 
     678               ENDIF 
     679 
     680            ENDIF 
     681#elif defined key_lim3 
     682            ! sea ice 
     683            IF( nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN 
     684 
     685               ! Test for types of ice input (lim2 or lim3)  
     686               CALL iom_open ( bn_a_i%clname, inum ) 
     687               id1 = iom_varid ( inum, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. ) 
     688               CALL iom_close ( inum ) 
     689               !CALL fld_clopn ( bn_a_i, nyear, nmonth, nday, ldstop=.TRUE. ) 
     690               !CALL iom_open ( bn_a_i %clname, inum ) 
     691               !id1 = iom_varid ( bn_a_i%num, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. ) 
     692                IF ( zndims == 4 ) THEN 
     693                 ll_bdylim3 = .TRUE.   ! lim3 input 
     694               ELSE 
     695                 ll_bdylim3 = .FALSE.  ! lim2 input       
     696               ENDIF 
     697               ! End test 
     698 
     699               IF( dta%ll_a_i ) THEN 
     700                  jfld = jfld + 1 
     701                  blf_i(jfld) = bn_a_i 
     702                  ibdy(jfld) = ib_bdy 
     703                  igrid(jfld) = 1 
     704                  ilen1(jfld) = nblen(igrid(jfld)) 
     705                  IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF 
     706               ENDIF 
     707 
     708               IF( dta%ll_ht_i ) THEN 
     709                  jfld = jfld + 1 
     710                  blf_i(jfld) = bn_ht_i 
     711                  ibdy(jfld) = ib_bdy 
     712                  igrid(jfld) = 1 
     713                  ilen1(jfld) = nblen(igrid(jfld)) 
     714                  IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF 
     715               ENDIF 
     716 
     717               IF( dta%ll_ht_s ) THEN 
     718                  jfld = jfld + 1 
     719                   blf_i(jfld) = bn_ht_s 
     720                  ibdy(jfld) = ib_bdy 
     721                  igrid(jfld) = 1 
     722                  ilen1(jfld) = nblen(igrid(jfld)) 
     723                  IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF 
     724               ENDIF 
     725 
    558726#endif 
    559727            ! Recalculate field counts 
    … …  
    568736            ENDIF 
    569737 
     738            dta%nread(1) = nb_bdy_fld(ib_bdy) 
     739 
    570740         ENDIF ! nn_dta .eq. 1 
    571741      ENDDO ! ib_bdy 
    … …  
    596766 
    597767         nblen => idx_bdy(ib_bdy)%nblen 
    598          nblenrim => idx_bdy(ib_bdy)%nblenrim 
    599  
    600          IF (nn_dyn2d(ib_bdy) .gt. 0) THEN 
    601             IF( nn_dyn2d_dta(ib_bdy) .eq. 0 .or. nn_dyn2d_dta(ib_bdy) .eq. 2 .or. ln_full_vel_array(ib_bdy) ) THEN 
    602                ilen0(1:3) = nblen(1:3) 
    603                ALLOCATE( dta_bdy(ib_bdy)%u2d(ilen0(2)) ) 
    604                ALLOCATE( dta_bdy(ib_bdy)%v2d(ilen0(3)) ) 
    605                IF ( nn_dyn2d(ib_bdy) .ne. jp_frs .and. (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3) ) THEN 
    606                   jfld = jfld + 1 
    607                   dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1) 
     768         dta => dta_bdy(ib_bdy) 
     769 
     770         if(lwp) then 
     771            write(numout,*) '++++++ dta%ll_ssh = ',dta%ll_ssh 
     772            write(numout,*) '++++++ dta%ll_u2d = ',dta%ll_u2d 
     773            write(numout,*) '++++++ dta%ll_v2d = ',dta%ll_v2d 
     774            write(numout,*) '++++++ dta%ll_u3d = ',dta%ll_u3d 
     775            write(numout,*) '++++++ dta%ll_v3d = ',dta%ll_v3d 
     776            write(numout,*) '++++++ dta%ll_tem = ',dta%ll_tem 
     777            write(numout,*) '++++++ dta%ll_sal = ',dta%ll_sal 
     778         endif 
     779 
     780         IF ( nn_dyn2d_dta(ib_bdy) .eq. 0 .or. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN 
     781            if(lwp) write(numout,*) '++++++ dta%ssh/u2d/u3d allocated space' 
     782            IF( dta%ll_ssh ) ALLOCATE( dta%ssh(nblen(1)) ) 
     783            IF( dta%ll_u2d ) ALLOCATE( dta%u2d(nblen(2)) ) 
     784            IF( dta%ll_v2d ) ALLOCATE( dta%v2d(nblen(3)) ) 
     785         ENDIF 
     786         IF ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN 
     787            IF( dta%ll_ssh ) THEN 
     788               if(lwp) write(numout,*) '++++++ dta%ssh pointing to fnow' 
     789               jfld = jfld + 1 
     790               dta%ssh => bf(jfld)%fnow(:,1,1) 
     791            ENDIF 
     792            IF ( dta%ll_u2d ) THEN 
     793               IF ( ln_full_vel_array(ib_bdy) ) THEN 
     794                  if(lwp) write(numout,*) '++++++ dta%u2d allocated space' 
     795                  ALLOCATE( dta%u2d(nblen(2)) ) 
    608796               ELSE 
    609                   ALLOCATE( dta_bdy(ib_bdy)%ssh(nblen(1)) ) 
    610                ENDIF 
    611             ELSE 
    612                IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN 
    613                   jfld = jfld + 1 
    614                   dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1) 
    615                ENDIF 
     797                  if(lwp) write(numout,*) '++++++ dta%u2d pointing to fnow' 
     798                  jfld = jfld + 1 
     799                  dta%u2d => bf(jfld)%fnow(:,1,1) 
     800               ENDIF 
     801            ENDIF 
     802            IF ( dta%ll_v2d ) THEN 
     803               IF ( ln_full_vel_array(ib_bdy) ) THEN 
     804                  if(lwp) write(numout,*) '++++++ dta%v2d allocated space' 
     805                  ALLOCATE( dta%v2d(nblen(3)) ) 
     806               ELSE 
     807                  if(lwp) write(numout,*) '++++++ dta%v2d pointing to fnow' 
     808                  jfld = jfld + 1 
     809                  dta%v2d => bf(jfld)%fnow(:,1,1) 
     810               ENDIF 
     811            ENDIF 
     812         ENDIF 
     813 
     814         IF ( nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 
     815            if(lwp) write(numout,*) '++++++ dta%u3d/v3d allocated space' 
     816            IF( dta%ll_u3d ) ALLOCATE( dta_bdy(ib_bdy)%u3d(nblen(2),jpk) ) 
     817            IF( dta%ll_v3d ) ALLOCATE( dta_bdy(ib_bdy)%v3d(nblen(3),jpk) ) 
     818         ENDIF 
     819         IF ( nn_dyn3d_dta(ib_bdy) .eq. 1 .or. & 
     820           &  ( ln_full_vel_array(ib_bdy) .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN 
     821            IF ( dta%ll_u3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) ) THEN 
     822               if(lwp) write(numout,*) '++++++ dta%u3d pointing to fnow' 
    616823               jfld = jfld + 1 
    617                dta_bdy(ib_bdy)%u2d => bf(jfld)%fnow(:,1,1) 
     824               dta_bdy(ib_bdy)%u3d => bf(jfld)%fnow(:,1,:) 
     825            ENDIF 
     826            IF ( dta%ll_v3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) ) THEN 
     827               if(lwp) write(numout,*) '++++++ dta%v3d pointing to fnow' 
    618828               jfld = jfld + 1 
    619                dta_bdy(ib_bdy)%v2d => bf(jfld)%fnow(:,1,1) 
    620             ENDIF 
    621          ENDIF 
    622  
    623          IF ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 
    624             ilen0(1:3) = nblen(1:3) 
    625             ALLOCATE( dta_bdy(ib_bdy)%u3d(ilen0(2),jpk) ) 
    626             ALLOCATE( dta_bdy(ib_bdy)%v3d(ilen0(3),jpk) ) 
    627          ENDIF 
    628          IF ( ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ).or. & 
    629            &  ( ln_full_vel_array(ib_bdy) .and. nn_dyn2d(ib_bdy) .gt. 0 .and.   & 
    630            &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN 
    631             jfld = jfld + 1 
    632             dta_bdy(ib_bdy)%u3d => bf(jfld)%fnow(:,1,:) 
    633             jfld = jfld + 1 
    634             dta_bdy(ib_bdy)%v3d => bf(jfld)%fnow(:,1,:) 
    635          ENDIF 
    636  
    637          IF (nn_tra(ib_bdy) .gt. 0) THEN 
    638             IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN 
    639                ilen0(1:3) = nblen(1:3) 
    640                ALLOCATE( dta_bdy(ib_bdy)%tem(ilen0(1),jpk) ) 
    641                ALLOCATE( dta_bdy(ib_bdy)%sal(ilen0(1),jpk) ) 
    642             ELSE 
     829               dta_bdy(ib_bdy)%v3d => bf(jfld)%fnow(:,1,:) 
     830            ENDIF 
     831         ENDIF 
     832 
     833         IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN 
     834            if(lwp) write(numout,*) '++++++ dta%tem/sal allocated space' 
     835            IF( dta%ll_tem ) ALLOCATE( dta_bdy(ib_bdy)%tem(nblen(1),jpk) ) 
     836            IF( dta%ll_sal ) ALLOCATE( dta_bdy(ib_bdy)%sal(nblen(1),jpk) ) 
     837         ELSE 
     838            IF( dta%ll_tem ) THEN 
     839               if(lwp) write(numout,*) '++++++ dta%tem pointing to fnow' 
    643840               jfld = jfld + 1 
    644841               dta_bdy(ib_bdy)%tem => bf(jfld)%fnow(:,1,:) 
     842            ENDIF 
     843            IF( dta%ll_sal ) THEN  
     844               if(lwp) write(numout,*) '++++++ dta%sal pointing to fnow' 
    645845               jfld = jfld + 1 
    646846               dta_bdy(ib_bdy)%sal => bf(jfld)%fnow(:,1,:) 
    … …  
    649849 
    650850#if defined key_lim2 
    651          IF (nn_ice_lim2(ib_bdy) .gt. 0) THEN 
     851         IF (nn_ice_lim(ib_bdy) .gt. 0) THEN 
    652852            IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 
    653                ilen0(1:3) = nblen(1:3) 
    654                ALLOCATE( dta_bdy(ib_bdy)%frld(ilen0(1)) ) 
    655                ALLOCATE( dta_bdy(ib_bdy)%hicif(ilen0(1)) ) 
    656                ALLOCATE( dta_bdy(ib_bdy)%hsnif(ilen0(1)) ) 
     853               ALLOCATE( dta_bdy(ib_bdy)%frld(nblen(1)) ) 
     854               ALLOCATE( dta_bdy(ib_bdy)%hicif(nblen(1)) ) 
     855               ALLOCATE( dta_bdy(ib_bdy)%hsnif(nblen(1)) ) 
    657856            ELSE 
    658857               jfld = jfld + 1 
    … …  
    662861               jfld = jfld + 1 
    663862               dta_bdy(ib_bdy)%hsnif => bf(jfld)%fnow(:,1,1) 
     863            ENDIF 
     864         ENDIF 
     865#elif defined key_lim3 
     866         IF (nn_ice_lim(ib_bdy) .gt. 0) THEN 
     867            IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN 
     868               ALLOCATE( dta_bdy(ib_bdy)%a_i (nblen(1),jpl) ) 
     869               ALLOCATE( dta_bdy(ib_bdy)%ht_i(nblen(1),jpl) ) 
     870               ALLOCATE( dta_bdy(ib_bdy)%ht_s(nblen(1),jpl) ) 
     871            ELSE 
     872               IF ( ll_bdylim3 ) THEN ! case input is lim3 type 
     873                  jfld = jfld + 1 
     874                  dta_bdy(ib_bdy)%a_i  => bf(jfld)%fnow(:,1,:) 
     875                  jfld = jfld + 1 
     876                  dta_bdy(ib_bdy)%ht_i => bf(jfld)%fnow(:,1,:) 
     877                  jfld = jfld + 1 
     878                  dta_bdy(ib_bdy)%ht_s => bf(jfld)%fnow(:,1,:) 
     879               ELSE ! case input is lim2 type 
     880                  jfld_ai  = jfld + 1 
     881                  jfld_hti = jfld + 2 
     882                  jfld_hts = jfld + 3 
     883                  jfld     = jfld + 3 
     884                  ALLOCATE( dta_bdy(ib_bdy)%a_i (nblen(1),jpl) ) 
     885                  ALLOCATE( dta_bdy(ib_bdy)%ht_i(nblen(1),jpl) ) 
     886                  ALLOCATE( dta_bdy(ib_bdy)%ht_s(nblen(1),jpl) ) 
     887                  dta_bdy(ib_bdy)%a_i (:,:) = 0.0 
     888                  dta_bdy(ib_bdy)%ht_i(:,:) = 0.0 
     889                  dta_bdy(ib_bdy)%ht_s(:,:) = 0.0 
     890               ENDIF 
     891 
    664892            ENDIF 
    665893         ENDIF 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90

    r4153 r4292  
    3030   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    3131   USE in_out_manager  ! 
    32    USE domvvl          ! variable volume 
     32   USE domvvl 
    3333 
    3434   IMPLICIT NONE 
    … …  
    5757      LOGICAL, INTENT( in ), OPTIONAL :: dyn3d_only       ! T => only update baroclinic velocities 
    5858      !! 
    59       INTEGER               :: jk,ii,ij,ib,igrd     ! Loop counter 
    60       LOGICAL               :: ll_dyn2d, ll_dyn3d   
    61       !! 
     59      INTEGER               :: jk,ii,ij,ib_bdy,ib,igrd     ! Loop counter 
     60      LOGICAL               :: ll_dyn2d, ll_dyn3d, ll_orlanski 
     61      !! 
     62      REAL(wp), POINTER, DIMENSION(:,:) :: phur1, phvr1     ! inverse depth at u and v points 
    6263 
    6364      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn') 
    … …  
    7071      ENDIF 
    7172 
     73      ll_orlanski = .false. 
     74      DO ib_bdy = 1, nb_bdy 
     75         IF ( cn_dyn2d(ib_bdy) == 'orlanski' .or. cn_dyn2d(ib_bdy) == 'orlanski_npo' & 
     76     &   .or. cn_dyn3d(ib_bdy) == 'orlanski' .or. cn_dyn3d(ib_bdy) == 'orlanski_npo') ll_orlanski = .true. 
     77      ENDDO 
     78 
    7279      !------------------------------------------------------- 
    7380      ! Set pointers 
    … …  
    7784      phur => hur 
    7885      phvr => hvr 
    79       CALL wrk_alloc(jpi,jpj,pu2d,pv2d)  
     86      CALL wrk_alloc(jpi,jpj,pua2d,pva2d)  
     87      IF ( ll_orlanski ) CALL wrk_alloc(jpi,jpj,pub2d,pvb2d,phur1,phvr1)  
    8088 
    8189      !------------------------------------------------------- 
    … …  
    8391      !------------------------------------------------------- 
    8492 
    85       pu2d(:,:) = 0.e0 
    86       pv2d(:,:) = 0.e0 
     93      ! "After" velocities:  
     94 
     95      pua2d(:,:) = 0.e0 
     96      pva2d(:,:) = 0.e0 
     97       
    8798      IF (lk_vvl) THEN 
    88          DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
    89             pu2d(:,:) = pu2d(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
    90             pv2d(:,:) = pv2d(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
    91          END DO 
    92          pu2d(:,:) = pu2d(:,:) / ( hu_0(:,:) + sshu_a(:,:) + 1._wp - umask(:,:,1) ) 
    93          pv2d(:,:) = pv2d(:,:) / ( hv_0(:,:) + sshv_a(:,:) + 1._wp - vmask(:,:,1) ) 
     99         phur1(:,:) = 0. 
     100         phvr1(:,:) = 0. 
     101         DO jk = 1, jpkm1 
     102            phur1(:,:) = phur1(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) 
     103            phvr1(:,:) = phvr1(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) 
     104            pua2d(:,:) = pua2d(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
     105            pva2d(:,:) = pva2d(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
     106         END DO 
     107         phur1(:,:) = umask(:,:,1) / ( phur1(:,:) + 1. - umask(:,:,1) ) 
     108         phvr1(:,:) = vmask(:,:,1) / ( phvr1(:,:) + 1. - vmask(:,:,1) ) 
     109         pua2d(:,:) = pua2d(:,:) * phur1(:,:) 
     110         pva2d(:,:) = pva2d(:,:) * phvr1(:,:) 
    94111      ELSE 
    95          DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
    96             pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
    97             pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
    98          END DO 
    99          pu2d(:,:) = pu2d(:,:) * phur(:,:) 
    100          pv2d(:,:) = pv2d(:,:) * phvr(:,:) 
     112         DO jk = 1, jpkm1 
     113            pua2d(:,:) = pua2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
     114            pva2d(:,:) = pva2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
     115         END DO 
     116         pua2d(:,:) = pua2d(:,:) * phur(:,:) 
     117         pva2d(:,:) = pva2d(:,:) * phvr(:,:) 
    101118      ENDIF 
     119 
    102120      DO jk = 1 , jpkm1 
    103          ua(:,:,jk) = ua(:,:,jk) - pu2d(:,:) * umask(:,:,jk) 
    104          va(:,:,jk) = va(:,:,jk) - pv2d(:,:) * vmask(:,:,jk) 
     121         ua(:,:,jk) = ua(:,:,jk) - pua2d(:,:) 
     122         va(:,:,jk) = va(:,:,jk) - pva2d(:,:) 
    105123      END DO 
     124 
     125      ! "Before" velocities (required for Orlanski condition):  
     126 
     127      IF ( ll_orlanski ) THEN           
     128         pub2d(:,:) = 0.e0 
     129         pvb2d(:,:) = 0.e0 
     130 
     131         IF (lk_vvl) THEN 
     132            phur1(:,:) = 0. 
     133            phvr1(:,:) = 0. 
     134            DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
     135               phur1(:,:) = phur1(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 
     136               phvr1(:,:) = phvr1(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 
     137               pub2d(:,:) = pub2d(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) * ub(:,:,jk) 
     138               pvb2d(:,:) = pvb2d(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) * vb(:,:,jk) 
     139            END DO 
     140            phur1(:,:) = umask(:,:,1) / ( phur1(:,:) + 1. - umask(:,:,1) ) 
     141            phvr1(:,:) = vmask(:,:,1) / ( phvr1(:,:) + 1. - vmask(:,:,1) ) 
     142            pub2d(:,:) = pub2d(:,:) * phur1(:,:) 
     143            pvb2d(:,:) = pvb2d(:,:) * phvr1(:,:) 
     144         ELSE 
     145            DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
     146               pub2d(:,:) = pub2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ub(:,:,jk) 
     147               pvb2d(:,:) = pvb2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * vb(:,:,jk) 
     148            END DO 
     149            pub2d(:,:) = pub2d(:,:) * phur(:,:) 
     150            pvb2d(:,:) = pvb2d(:,:) * phvr(:,:) 
     151         ENDIF 
     152 
     153         DO jk = 1 , jpkm1 
     154            ub(:,:,jk) = ub(:,:,jk) - pub2d(:,:) 
     155            vb(:,:,jk) = vb(:,:,jk) - pvb2d(:,:) 
     156         END DO 
     157      END IF 
    106158 
    107159      !------------------------------------------------------- 
    … …  
    119171 
    120172      DO jk = 1 , jpkm1 
    121          ua(:,:,jk) = ( ua(:,:,jk) + pu2d(:,:) ) * umask(:,:,jk) 
    122          va(:,:,jk) = ( va(:,:,jk) + pv2d(:,:) ) * vmask(:,:,jk) 
     173         ua(:,:,jk) = ( ua(:,:,jk) + pua2d(:,:) ) * umask(:,:,jk) 
     174         va(:,:,jk) = ( va(:,:,jk) + pva2d(:,:) ) * vmask(:,:,jk) 
    123175      END DO 
    124176 
    125       CALL wrk_dealloc(jpi,jpj,pu2d,pv2d)  
     177      IF ( ll_orlanski ) THEN 
     178         DO jk = 1 , jpkm1 
     179            ub(:,:,jk) = ( ub(:,:,jk) + pub2d(:,:) ) * umask(:,:,jk) 
     180            vb(:,:,jk) = ( vb(:,:,jk) + pvb2d(:,:) ) * vmask(:,:,jk) 
     181         END DO 
     182      END IF 
     183 
     184      CALL wrk_dealloc(jpi,jpj,pua2d,pva2d)  
     185      IF ( ll_orlanski ) CALL wrk_dealloc(jpi,jpj,pub2d,pvb2d)  
    126186 
    127187      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn') 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90

    r3680 r4292  
    66   !! History :  3.4  !  2011     (D. Storkey) new module as part of BDY rewrite 
    77   !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Optimization of BDY communications 
     8   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_bdy  
    … …  
    1112   !!   'key_bdy' :                    Unstructured Open Boundary Condition 
    1213   !!---------------------------------------------------------------------- 
    13    !!   bdy_dyn2d      : Apply open boundary conditions to barotropic variables. 
    14    !!   bdy_dyn2d_fla    : Apply Flather condition 
     14   !!   bdy_dyn2d          : Apply open boundary conditions to barotropic variables. 
     15   !!   bdy_dyn2d_frs      : Apply Flow Relaxation Scheme  
     16   !!   bdy_dyn2d_fla      : Apply Flather condition 
     17   !!   bdy_dyn2d_orlanski : Orlanski Radiation 
     18   !!   bdy_ssh            : Duplicate sea level across open boundaries 
    1519   !!---------------------------------------------------------------------- 
    1620   USE timing          ! Timing 
    … …  
    1822   USE dom_oce         ! ocean space and time domain 
    1923   USE bdy_oce         ! ocean open boundary conditions 
     24   USE bdylib          ! BDY library routines 
    2025   USE dynspg_oce      ! for barotropic variables 
    2126   USE phycst          ! physical constants 
    … …  
    2631   PRIVATE 
    2732 
    28    PUBLIC   bdy_dyn2d     ! routine called in dynspg_ts and bdy_dyn 
     33   PUBLIC   bdy_dyn2d   ! routine called in dynspg_ts and bdy_dyn 
     34   PUBLIC   bdy_ssh       ! routine called in dynspg_ts or sshwzv 
    2935 
    3036   !!---------------------------------------------------------------------- 
    … …  
    4854      DO ib_bdy=1, nb_bdy 
    4955 
    50          SELECT CASE( nn_dyn2d(ib_bdy) ) 
    51          CASE(jp_none) 
     56         SELECT CASE( cn_dyn2d(ib_bdy) ) 
     57         CASE('none') 
    5258            CYCLE 
    53          CASE(jp_frs) 
     59         CASE('frs') 
    5460            CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy ) 
    55          CASE(jp_flather) 
     61         CASE('flather') 
    5662            CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy ) 
     63         CASE('orlanski') 
     64            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. ) 
     65         CASE('orlanski_npo') 
     66            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. ) 
    5767         CASE DEFAULT 
    5868            CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' ) 
    … …  
    8999         ij   = idx%nbj(jb,igrd) 
    90100         zwgt = idx%nbw(jb,igrd) 
    91          pu2d(ii,ij) = ( pu2d(ii,ij) + zwgt * ( dta%u2d(jb) - pu2d(ii,ij) ) ) * umask(ii,ij,1) 
     101         pua2d(ii,ij) = ( pua2d(ii,ij) + zwgt * ( dta%u2d(jb) - pua2d(ii,ij) ) ) * umask(ii,ij,1) 
    92102      END DO 
    93103      ! 
    … …  
    97107         ij   = idx%nbj(jb,igrd) 
    98108         zwgt = idx%nbw(jb,igrd) 
    99          pv2d(ii,ij) = ( pv2d(ii,ij) + zwgt * ( dta%v2d(jb) - pv2d(ii,ij) ) ) * vmask(ii,ij,1) 
     109         pva2d(ii,ij) = ( pva2d(ii,ij) + zwgt * ( dta%v2d(jb) - pva2d(ii,ij) ) ) * vmask(ii,ij,1) 
    100110      END DO  
    101       CALL lbc_bdy_lnk( pu2d, 'U', -1., ib_bdy )  
    102       CALL lbc_bdy_lnk( pv2d, 'V', -1., ib_bdy)   ! Boundary points should be updated 
     111      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )  
     112      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy)   ! Boundary points should be updated 
    103113      ! 
    104114      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_frs') 
    … …  
    133143      INTEGER  ::   jb, igrd                         ! dummy loop indices 
    134144      INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses 
     145      REAL(wp), POINTER :: flagu, flagv              ! short cuts 
    135146      REAL(wp) ::   zcorr                            ! Flather correction 
    136147      REAL(wp) ::   zforc                            ! temporary scalar 
     148      REAL(wp) ::   zflag, z1_2                      !    "        " 
    137149      !!---------------------------------------------------------------------- 
    138150 
    139151      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_fla') 
     152 
     153      z1_2 = 0.5_wp 
    140154 
    141155      ! ---------------------------------! 
    … …  
    160174         ii  = idx%nbi(jb,igrd) 
    161175         ij  = idx%nbj(jb,igrd)  
    162          iim1 = ii + MAX( 0, INT( idx%flagu(jb) ) )   ! T pts i-indice inside the boundary 
    163          iip1 = ii - MIN( 0, INT( idx%flagu(jb) ) )   ! T pts i-indice outside the boundary  
     176         flagu => idx%flagu(jb,igrd) 
     177         iim1 = ii + MAX( 0, INT( flagu ) )   ! T pts i-indice inside the boundary 
     178         iip1 = ii - MIN( 0, INT( flagu ) )   ! T pts i-indice outside the boundary  
    164179         ! 
    165          zcorr = - idx%flagu(jb) * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 
    166          zforc = dta%u2d(jb) 
    167          pu2d(ii,ij) = zforc + zcorr * umask(ii,ij,1)  
     180         zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 
     181 
     182         ! jchanut tschanges: Set zflag to 0 below to revert to Flather scheme 
     183         ! Use characteristics method instead 
     184         zflag = ABS(flagu) 
     185         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij) 
     186         pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1)  
    168187      END DO 
    169188      ! 
    … …  
    173192         ii  = idx%nbi(jb,igrd) 
    174193         ij  = idx%nbj(jb,igrd)  
    175          ijm1 = ij + MAX( 0, INT( idx%flagv(jb) ) )   ! T pts j-indice inside the boundary 
    176          ijp1 = ij - MIN( 0, INT( idx%flagv(jb) ) )   ! T pts j-indice outside the boundary  
     194         flagv => idx%flagv(jb,igrd) 
     195         ijm1 = ij + MAX( 0, INT( flagv ) )   ! T pts j-indice inside the boundary 
     196         ijp1 = ij - MIN( 0, INT( flagv ) )   ! T pts j-indice outside the boundary  
    177197         ! 
    178          zcorr = - idx%flagv(jb) * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 
    179          zforc = dta%v2d(jb) 
    180          pv2d(ii,ij) = zforc + zcorr * vmask(ii,ij,1) 
    181       END DO 
    182       CALL lbc_bdy_lnk( pu2d, 'U', -1., ib_bdy )   ! Boundary points should be updated 
    183       CALL lbc_bdy_lnk( pv2d, 'V', -1., ib_bdy )   ! 
     198         zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 
     199 
     200         ! jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme 
     201         ! Use characteristics method instead 
     202         zflag = ABS(flagv) 
     203         zforc  = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1) 
     204         pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1) 
     205      END DO 
     206      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated 
     207      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy )   ! 
    184208      ! 
    185209      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_fla') 
    186210      ! 
    187211   END SUBROUTINE bdy_dyn2d_fla 
     212 
     213 
     214   SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, ll_npo ) 
     215      !!---------------------------------------------------------------------- 
     216      !!                 ***  SUBROUTINE bdy_dyn2d_orlanski  *** 
     217      !!              
     218      !!              - Apply Orlanski radiation condition adaptively: 
     219      !!                  - radiation plus weak nudging at outflow points 
     220      !!                  - no radiation and strong nudging at inflow points 
     221      !!  
     222      !! 
     223      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)     
     224      !!---------------------------------------------------------------------- 
     225      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices 
     226      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data 
     227      INTEGER,                      INTENT(in) ::   ib_bdy  ! number of current open boundary set 
     228      LOGICAL,                      INTENT(in) ::   ll_npo  ! flag for NPO version 
     229 
     230      INTEGER  ::   ib, igrd                               ! dummy loop indices 
     231      INTEGER  ::   ii, ij, iibm1, ijbm1                   ! indices 
     232      !!---------------------------------------------------------------------- 
     233 
     234      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_orlanski') 
     235      ! 
     236      igrd = 2      ! Orlanski bc on u-velocity;  
     237      !             
     238      CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, ll_npo ) 
     239 
     240      igrd = 3      ! Orlanski bc on v-velocity 
     241      !   
     242      CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, ll_npo ) 
     243      ! 
     244      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski') 
     245      ! 
     246      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated 
     247      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy )   ! 
     248      ! 
     249      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski') 
     250      ! 
     251   END SUBROUTINE bdy_dyn2d_orlanski 
     252 
     253   SUBROUTINE bdy_ssh( zssh ) 
     254      !!---------------------------------------------------------------------- 
     255      !!                  ***  SUBROUTINE bdy_ssh  *** 
     256      !! 
     257      !! ** Purpose : Duplicate sea level across open boundaries 
     258      !! 
     259      !!---------------------------------------------------------------------- 
     260      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zssh ! Sea level 
     261      !! 
     262      INTEGER  ::   ib_bdy, ib, igrd                        ! local integers 
     263      INTEGER  ::   ii, ij, zcoef, zcoef1, zcoef2, ip, jp   !   "       " 
     264 
     265      igrd = 1                       ! Everything is at T-points here 
     266 
     267      DO ib_bdy = 1, nb_bdy 
     268         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     269            ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     270            ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     271            ! Set gradient direction: 
     272            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  ) 
     273            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1) 
     274            IF ( zcoef1+zcoef2 == 0 ) THEN 
     275               ! corner 
     276!               zcoef = tmask(ii-1,ij,1) + tmask(ii+1,ij,1) +  tmask(ii,ij-1,1) +  tmask(ii,ij+1,1) 
     277!               zssh(ii,ij) = zssh(ii-1,ij  ) * tmask(ii-1,ij  ,1) + & 
     278!                 &           zssh(ii+1,ij  ) * tmask(ii+1,ij  ,1) + & 
     279!                 &           zssh(ii  ,ij-1) * tmask(ii  ,ij-1,1) + & 
     280!                 &           zssh(ii  ,ij+1) * tmask(ii  ,ij+1,1) 
     281               zcoef = bdytmask(ii-1,ij) + bdytmask(ii+1,ij) +  bdytmask(ii,ij-1) +  bdytmask(ii,ij+1) 
     282               zssh(ii,ij) = zssh(ii-1,ij  ) * bdytmask(ii-1,ij  ) + & 
     283                 &           zssh(ii+1,ij  ) * bdytmask(ii+1,ij  ) + & 
     284                 &           zssh(ii  ,ij-1) * bdytmask(ii  ,ij-1) + & 
     285                 &           zssh(ii  ,ij+1) * bdytmask(ii  ,ij+1) 
     286               zssh(ii,ij) = ( zssh(ii,ij) / MAX( 1, zcoef) ) * tmask(ii,ij,1) 
     287            ELSE 
     288               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  ) 
     289               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1) 
     290               zssh(ii,ij) = zssh(ii+ip,ij+jp) * tmask(ii+ip,ij+jp,1) 
     291            ENDIF 
     292         END DO 
     293 
     294         ! Boundary points should be updated 
     295         CALL lbc_bdy_lnk( zssh(:,:), 'T', 1., ib_bdy ) 
     296      END DO 
     297 
     298   END SUBROUTINE bdy_ssh 
     299 
    188300#else 
    189301   !!---------------------------------------------------------------------- 
    … …  
    192304CONTAINS 
    193305   SUBROUTINE bdy_dyn2d( kt )      ! Empty routine 
    194       WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt 
     306      INTEGER, intent(in) :: kt 
     307      WRITE(*,*) 'bdy_dyn2d: You should not have seen this print! error?', kt 
    195308   END SUBROUTINE bdy_dyn2d 
     309 
    196310#endif 
    197311 
    198312   !!====================================================================== 
    199313END MODULE bdydyn2d 
     314 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90

    r3703 r4292  
    1919   USE dom_oce         ! ocean space and time domain 
    2020   USE bdy_oce         ! ocean open boundary conditions 
     21   USE bdylib          ! for orlanski library routines 
    2122   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2223   USE in_out_manager  ! 
    … …  
    5253      DO ib_bdy=1, nb_bdy 
    5354 
    54 !!$         IF ( using Orlanski radiation conditions ) THEN  
    55 !!$            CALL bdy_rad( kt,  bdyidx(ib_bdy) ) 
    56 !!$         ENDIF 
    57  
    58          SELECT CASE( nn_dyn3d(ib_bdy) ) 
    59          CASE(jp_none) 
     55         SELECT CASE( cn_dyn3d(ib_bdy) ) 
     56         CASE('none') 
    6057            CYCLE 
    61          CASE(jp_frs) 
     58         CASE('frs') 
    6259            CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
    63          CASE(2) 
     60         CASE('specified') 
    6461            CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
    65          CASE(3) 
     62         CASE('zero') 
    6663            CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
     64         CASE('orlanski') 
     65            CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. ) 
     66         CASE('orlanski_npo') 
     67            CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. ) 
    6768         CASE DEFAULT 
    6869            CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' ) 
    … …  
    109110         END DO 
    110111      END DO 
    111       CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ;   CALL lbc_bdy_lnk( va, 'V', -1.,ib_bdy )   ! Boundary points should be updated 
     112      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ! Boundary points should be updated   
     113      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )    
    112114      ! 
    113115      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
    … …  
    204206         END DO 
    205207      END DO  
    206       CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ;   CALL lbc_bdy_lnk( va, 'V', -1.,ib_bdy )   ! Boundary points should be updated 
     208      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )    ! Boundary points should be updated 
     209      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )    
    207210      ! 
    208211      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
    … …  
    211214 
    212215   END SUBROUTINE bdy_dyn3d_frs 
     216 
     217   SUBROUTINE bdy_dyn3d_orlanski( idx, dta, ib_bdy, ll_npo ) 
     218      !!---------------------------------------------------------------------- 
     219      !!                 ***  SUBROUTINE bdy_dyn3d_orlanski  *** 
     220      !!              
     221      !!              - Apply Orlanski radiation to baroclinic velocities.  
     222      !!              - Wrapper routine for bdy_orlanski_3d 
     223      !!  
     224      !! 
     225      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)     
     226      !!---------------------------------------------------------------------- 
     227      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices 
     228      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data 
     229      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index 
     230      LOGICAL,                      INTENT(in) ::   ll_npo  ! switch for NPO version 
     231 
     232      INTEGER  ::   jb, igrd                               ! dummy loop indices 
     233      !!---------------------------------------------------------------------- 
     234 
     235      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_orlanski') 
     236      ! 
     237      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities.  
     238      ! 
     239      igrd = 2      ! Orlanski bc on u-velocity;  
     240      !             
     241      CALL bdy_orlanski_3d( idx, igrd, ub, ua, dta%u3d, ll_npo ) 
     242 
     243      igrd = 3      ! Orlanski bc on v-velocity 
     244      !   
     245      CALL bdy_orlanski_3d( idx, igrd, vb, va, dta%v3d, ll_npo ) 
     246      ! 
     247      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )    ! Boundary points should be updated 
     248      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )    
     249      ! 
     250      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_orlanski') 
     251      ! 
     252   END SUBROUTINE bdy_dyn3d_orlanski 
     253 
    213254 
    214255   SUBROUTINE bdy_dyn3d_dmp( kt ) 
    … …  
    225266      REAL(wp) ::   zwgt           ! boundary weight 
    226267      INTEGER  ::  ib_bdy          ! loop index 
     268      REAL(wp), POINTER, DIMENSION(:,:) :: phur1, phvr1     ! inverse depth at u and v points 
    227269      !!---------------------------------------------------------------------- 
    228270      ! 
    … …  
    232274      ! Remove barotropic part from before velocity 
    233275      !------------------------------------------------------- 
    234       CALL wrk_alloc(jpi,jpj,pu2d,pv2d)  
    235  
    236       pu2d(:,:) = 0.e0 
    237       pv2d(:,:) = 0.e0 
    238  
     276      CALL wrk_alloc(jpi,jpj,pub2d,pvb2d,phur1,phvr1)  
     277 
     278      pub2d(:,:) = 0.e0 
     279      pvb2d(:,:) = 0.e0 
     280 
     281      phur1(:,:) = 0. 
     282      phvr1(:,:) = 0. 
    239283      DO jk = 1, jpkm1 
    240284#if defined key_vvl 
    241          pu2d(:,:) = pu2d(:,:) + fse3u_b(:,:,jk)* ub(:,:,jk)   *umask(:,:,jk)  
    242          pv2d(:,:) = pv2d(:,:) + fse3v_b(:,:,jk)* vb(:,:,jk)   *vmask(:,:,jk) 
     285         phur1(:,:) = phur1(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) 
     286         phvr1(:,:) = phvr1(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) 
     287         pub2d(:,:) = pub2d(:,:) + fse3u_b(:,:,jk)* ub(:,:,jk)   *umask(:,:,jk)  
     288         pvb2d(:,:) = pvb2d(:,:) + fse3v_b(:,:,jk)* vb(:,:,jk)   *vmask(:,:,jk) 
    243289#else 
    244          pu2d(:,:) = pu2d(:,:) + fse3u_0(:,:,jk) * ub(:,:,jk)  * umask(:,:,jk) 
    245          pv2d(:,:) = pv2d(:,:) + fse3v_0(:,:,jk) * vb(:,:,jk)  * vmask(:,:,jk) 
     290         pub2d(:,:) = pub2d(:,:) + fse3u_0(:,:,jk) * ub(:,:,jk)  * umask(:,:,jk) 
     291         pvb2d(:,:) = pvb2d(:,:) + fse3v_0(:,:,jk) * vb(:,:,jk)  * vmask(:,:,jk) 
    246292#endif 
    247293      END DO 
    248294 
    249295      IF( lk_vvl ) THEN 
    250          pu2d(:,:) = pu2d(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 
    251          pv2d(:,:) = pv2d(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 
     296         phur1(:,:) = umask(:,:,1) / ( phur1(:,:) + 1. - umask(:,:,1) ) 
     297         phvr1(:,:) = vmask(:,:,1) / ( phvr1(:,:) + 1. - vmask(:,:,1) ) 
     298         pub2d(:,:) = pub2d(:,:) * umask(:,:,1) * phur1(:,:) 
     299         pvb2d(:,:) = pvb2d(:,:) * vmask(:,:,1) * phvr1(:,:) 
    252300      ELSE 
    253          pu2d(:,:) = pv2d(:,:) * hur(:,:) 
    254          pv2d(:,:) = pu2d(:,:) * hvr(:,:) 
     301         pub2d(:,:) = pvb2d(:,:) * hur(:,:) 
     302         pvb2d(:,:) = pub2d(:,:) * hvr(:,:) 
    255303      ENDIF 
    256304 
    257305      DO ib_bdy=1, nb_bdy 
    258          IF ( ln_dyn3d_dmp(ib_bdy).and.nn_dyn3d(ib_bdy).gt.0 ) THEN 
     306         IF ( ln_dyn3d_dmp(ib_bdy) .and. cn_dyn3d(ib_bdy) /= 'none' ) THEN 
    259307            igrd = 2                      ! Relaxation of zonal velocity 
    260308            DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    … …  
    264312               DO jk = 1, jpkm1 
    265313                  ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - & 
    266                                    ub(ii,ij,jk) + pu2d(ii,ij)) ) * umask(ii,ij,jk) 
     314                                   ub(ii,ij,jk) + pub2d(ii,ij)) ) * umask(ii,ij,jk) 
    267315               END DO 
    268316            END DO 
    … …  
    275323               DO jk = 1, jpkm1 
    276324                  va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) -  & 
    277                                    vb(ii,ij,jk) + pv2d(ii,ij)) ) * vmask(ii,ij,jk) 
     325                                   vb(ii,ij,jk) + pvb2d(ii,ij)) ) * vmask(ii,ij,jk) 
    278326               END DO 
    279327            END DO 
    … …  
    281329      ENDDO 
    282330      ! 
    283       CALL wrk_dealloc(jpi,jpj,pu2d,pv2d)  
     331      CALL wrk_dealloc(jpi,jpj,pub2d,pvb2d)  
    284332      ! 
    285333      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r4148 r4292  
    2121   !!   bdy_init       : Initialization of unstructured open boundaries 
    2222   !!---------------------------------------------------------------------- 
     23   USE wrk_nemo        ! Memory Allocation 
    2324   USE timing          ! Timing 
    2425   USE oce             ! ocean dynamics and tracers variables 
    … …  
    7980      INTEGER  ::   jpbdtau, jpbdtas                       !   -       - 
    8081      INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       - 
     82      INTEGER  ::   i_offset, j_offset                     !   -       - 
    8183      INTEGER, POINTER  ::  nbi, nbj, nbr                  ! short cuts 
    82       REAL   , POINTER  ::  flagu, flagv                   !    -   - 
     84      REAL(wp), POINTER  ::  flagu, flagv                  !    -   - 
     85      REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask    ! pointer to 2D mask fields 
    8386      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars 
    8487      INTEGER, DIMENSION (2)                  ::   kdimsz 
    … …  
    9093      INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b  ! Flags for boundaries receiving 
    9194      INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4)                ! Arrays for neighbours coordinates 
     95      REAL(wp), POINTER, DIMENSION(:,:)       ::   zfmask  ! temporary fmask array excluding coastal boundary condition (shlat) 
    9296 
    9397      !! 
    94       NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,             & 
    95          &             ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta, & 
    96          &             nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta,         &   
    97          &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp,             & 
    98 #if defined key_lim2 
    99          &             nn_ice_lim2, nn_ice_lim2_dta,                       & 
     98      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,                 & 
     99         &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,    & 
     100         &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &   
     101         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, & 
     102#if ( defined key_lim2 || defined key_lim3 ) 
     103         &             cn_ice_lim, nn_ice_lim_dta,                           & 
    100104#endif 
    101105         &             ln_vol, nn_volctl, nn_rimwidth 
    … …  
    156160 
    157161        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  ' 
    158         SELECT CASE( nn_dyn2d(ib_bdy) )                   
    159           CASE(jp_none)         ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    160           CASE(jp_frs)          ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    161           CASE(jp_flather)      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition' 
    162           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' ) 
     162        SELECT CASE( cn_dyn2d(ib_bdy) )                   
     163          CASE('none')          
     164             IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     165             dta_bdy(ib_bdy)%ll_ssh = .false. 
     166             dta_bdy(ib_bdy)%ll_u2d = .false. 
     167             dta_bdy(ib_bdy)%ll_v2d = .false. 
     168          CASE('frs')           
     169             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     170             dta_bdy(ib_bdy)%ll_ssh = .false. 
     171             dta_bdy(ib_bdy)%ll_u2d = .true. 
     172             dta_bdy(ib_bdy)%ll_v2d = .true. 
     173          CASE('flather')       
     174             IF(lwp) WRITE(numout,*) '      Flather radiation condition' 
     175             dta_bdy(ib_bdy)%ll_ssh = .true. 
     176             dta_bdy(ib_bdy)%ll_u2d = .true. 
     177             dta_bdy(ib_bdy)%ll_v2d = .true. 
     178          CASE('orlanski')      
     179             IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging' 
     180             dta_bdy(ib_bdy)%ll_ssh = .false. 
     181             dta_bdy(ib_bdy)%ll_u2d = .true. 
     182             dta_bdy(ib_bdy)%ll_v2d = .true. 
     183          CASE('orlanski_npo')  
     184             IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging' 
     185             dta_bdy(ib_bdy)%ll_ssh = .false. 
     186             dta_bdy(ib_bdy)%ll_u2d = .true. 
     187             dta_bdy(ib_bdy)%ll_v2d = .true. 
     188          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn2d' ) 
    163189        END SELECT 
    164         IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN 
     190        IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN 
    165191           SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !  
    166192              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
    … …  
    177203 
    178204        IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  ' 
    179         SELECT CASE( nn_dyn3d(ib_bdy) )                   
    180           CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    181           CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    182           CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value' 
    183           CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
    184           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn3d' ) 
     205        SELECT CASE( cn_dyn3d(ib_bdy) )                   
     206          CASE('none') 
     207             IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     208             dta_bdy(ib_bdy)%ll_u3d = .false. 
     209             dta_bdy(ib_bdy)%ll_v3d = .false. 
     210          CASE('frs')        
     211             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     212             dta_bdy(ib_bdy)%ll_u3d = .true. 
     213             dta_bdy(ib_bdy)%ll_v3d = .true. 
     214          CASE('specified') 
     215             IF(lwp) WRITE(numout,*) '      Specified value' 
     216             dta_bdy(ib_bdy)%ll_u3d = .true. 
     217             dta_bdy(ib_bdy)%ll_v3d = .true. 
     218          CASE('zero') 
     219             IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
     220             dta_bdy(ib_bdy)%ll_u3d = .false. 
     221             dta_bdy(ib_bdy)%ll_v3d = .false. 
     222          CASE('orlanski') 
     223             IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging' 
     224             dta_bdy(ib_bdy)%ll_u3d = .true. 
     225             dta_bdy(ib_bdy)%ll_v3d = .true. 
     226          CASE('orlanski_npo') 
     227             IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging' 
     228             dta_bdy(ib_bdy)%ll_u3d = .true. 
     229             dta_bdy(ib_bdy)%ll_v3d = .true. 
     230          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn3d' ) 
    185231        END SELECT 
    186         IF( nn_dyn3d(ib_bdy) .gt. 0 ) THEN 
     232        IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 
    187233           SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   !  
    188234              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
    … …  
    193239 
    194240        IF ( ln_dyn3d_dmp(ib_bdy) ) THEN 
    195            IF ( nn_dyn3d(ib_bdy).EQ.0 ) THEN 
     241           IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN 
    196242              IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.' 
    197243              ln_dyn3d_dmp(ib_bdy)=.false. 
    198            ELSEIF ( nn_dyn3d(ib_bdy).EQ.1 ) THEN 
     244           ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN 
    199245              CALL ctl_stop( 'Use FRS OR relaxation' ) 
    200246           ELSE 
    … …  
    202248              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
    203249              IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 
     250              dta_bdy(ib_bdy)%ll_u3d = .true. 
     251              dta_bdy(ib_bdy)%ll_v3d = .true. 
    204252           ENDIF 
    205253        ELSE 
    … …  
    209257 
    210258        IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  ' 
    211         SELECT CASE( nn_tra(ib_bdy) )                   
    212           CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    213           CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    214           CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value' 
    215           CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Neumann conditions' 
    216           CASE( 4 )      ;   IF(lwp) WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity' 
    217           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' ) 
     259        SELECT CASE( cn_tra(ib_bdy) )                   
     260          CASE('none') 
     261             IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     262             dta_bdy(ib_bdy)%ll_tem = .false. 
     263             dta_bdy(ib_bdy)%ll_sal = .false. 
     264          CASE('frs') 
     265             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     266             dta_bdy(ib_bdy)%ll_tem = .true. 
     267             dta_bdy(ib_bdy)%ll_sal = .true. 
     268          CASE('specified') 
     269             IF(lwp) WRITE(numout,*) '      Specified value' 
     270             dta_bdy(ib_bdy)%ll_tem = .true. 
     271             dta_bdy(ib_bdy)%ll_sal = .true. 
     272          CASE('neumann') 
     273             IF(lwp) WRITE(numout,*) '      Neumann conditions' 
     274             dta_bdy(ib_bdy)%ll_tem = .false. 
     275             dta_bdy(ib_bdy)%ll_sal = .false. 
     276          CASE('runoff') 
     277             IF(lwp) WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity' 
     278             dta_bdy(ib_bdy)%ll_tem = .false. 
     279             dta_bdy(ib_bdy)%ll_sal = .false. 
     280          CASE('orlanski') 
     281             IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging' 
     282             dta_bdy(ib_bdy)%ll_tem = .true. 
     283             dta_bdy(ib_bdy)%ll_sal = .true. 
     284          CASE('orlanski_npo') 
     285             IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging' 
     286             dta_bdy(ib_bdy)%ll_tem = .true. 
     287             dta_bdy(ib_bdy)%ll_sal = .true. 
     288          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_tra' ) 
    218289        END SELECT 
    219         IF( nn_tra(ib_bdy) .gt. 0 ) THEN 
     290        IF( cn_tra(ib_bdy) /= 'none' ) THEN 
    220291           SELECT CASE( nn_tra_dta(ib_bdy) )                   !  
    221292              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
    … …  
    226297 
    227298        IF ( ln_tra_dmp(ib_bdy) ) THEN 
    228            IF ( nn_tra(ib_bdy).EQ.0 ) THEN 
     299           IF ( cn_tra(ib_bdy) == 'none' ) THEN 
    229300              IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.' 
    230301              ln_tra_dmp(ib_bdy)=.false. 
    231            ELSEIF ( nn_tra(ib_bdy).EQ.1 ) THEN 
     302           ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN 
    232303              CALL ctl_stop( 'Use FRS OR relaxation' ) 
    233304           ELSE 
    234305              IF(lwp) WRITE(numout,*) '      + T/S relaxation zone' 
    235306              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
     307              IF(lwp) WRITE(numout,*) '      Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days' 
    236308              IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 
     309              dta_bdy(ib_bdy)%ll_tem = .true. 
     310              dta_bdy(ib_bdy)%ll_sal = .true. 
    237311           ENDIF 
    238312        ELSE 
    … …  
    243317#if defined key_lim2 
    244318        IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  ' 
    245         SELECT CASE( nn_ice_lim2(ib_bdy) )                   
    246           CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    247           CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    248           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' ) 
     319        SELECT CASE( cn_ice_lim(ib_bdy) )                   
     320          CASE('none') 
     321             IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     322             dta_bdy(ib_bdy)%ll_frld  = .false. 
     323             dta_bdy(ib_bdy)%ll_hicif = .false. 
     324             dta_bdy(ib_bdy)%ll_hsnif = .false. 
     325          CASE('frs') 
     326             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     327             dta_bdy(ib_bdy)%ll_frld  = .true. 
     328             dta_bdy(ib_bdy)%ll_hicif = .true. 
     329             dta_bdy(ib_bdy)%ll_hsnif = .true. 
     330          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice_lim' ) 
    249331        END SELECT 
    250         IF( nn_ice_lim2(ib_bdy) .gt. 0 ) THEN  
    251            SELECT CASE( nn_ice_lim2_dta(ib_bdy) )                   !  
     332        IF( cn_ice_lim(ib_bdy) /= 'none' ) THEN  
     333           SELECT CASE( nn_ice_lim_dta(ib_bdy) )                   !  
    252334              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
    253335              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
    254               CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim2_dta must be 0 or 1' ) 
     336              CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim_dta must be 0 or 1' ) 
     337           END SELECT 
     338        ENDIF 
     339        IF(lwp) WRITE(numout,*) 
     340#elif defined key_lim3 
     341        IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  ' 
     342        SELECT CASE( cn_ice_lim(ib_bdy) )                   
     343          CASE('none') 
     344             IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     345             dta_bdy(ib_bdy)%ll_a_i  = .false. 
     346             dta_bdy(ib_bdy)%ll_ht_i = .false. 
     347             dta_bdy(ib_bdy)%ll_ht_s = .false. 
     348          CASE('frs') 
     349             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     350             dta_bdy(ib_bdy)%ll_a_i  = .true. 
     351             dta_bdy(ib_bdy)%ll_ht_i = .true. 
     352             dta_bdy(ib_bdy)%ll_ht_s = .true. 
     353          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice_lim' ) 
     354        END SELECT 
     355        IF( cn_ice_lim(ib_bdy) /= 'none' ) THEN  
     356           SELECT CASE( nn_ice_lim_dta(ib_bdy) )                   !  
     357              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
     358              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
     359              CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim_dta must be 0 or 1' ) 
    255360           END SELECT 
    256361        ENDIF 
    … …  
    740845               IF(lwp) THEN         ! Since all procs read global data only need to do this check on one proc... 
    741846                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN 
    742                      CALL ctl_stop('bdy_init : ERROR : boundary data in file  & 
    743                                     must be defined in order of distance from edge nbr.', & 
    744                                    'A utility for re-ordering boundary coordinates and data & 
    745                                     files exists in the TOOLS/OBC directory') 
     847                     CALL ctl_stop('bdy_init : ERROR : boundary data in file must be defined in order of distance from edge nbr.', & 
     848                                   'A utility for re-ordering boundary coordinates and data files exists in the TOOLS/OBC directory') 
    746849                  ENDIF     
    747850               ENDIF 
    … …  
    766869         ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) ) 
    767870         ALLOCATE( idx_bdy(ib_bdy)%nbd(ilen1,jpbgrd) ) 
     871         ALLOCATE( idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) ) 
    768872         ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) ) 
    769873         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) ) 
    770          ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) ) 
    771          ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) ) 
     874         ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1,jpbgrd) ) 
     875         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1,jpbgrd) ) 
    772876 
    773877         ! Dispatch mapping indices and discrete distances on each processor 
    … …  
    9371041            ENDDO 
    9381042         ENDDO  
     1043 
    9391044         ! definition of the i- and j- direction local boundaries arrays 
    9401045         ! used for sending the boudaries 
    … …  
    9901095               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 
    9911096               idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) &  
     1097               & *(FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.   ! quadratic 
     1098               idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) &  
    9921099               & *(FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.   ! quadratic 
    9931100            END DO 
    … …  
    10921199      ENDDO 
    10931200 
     1201      ! For the flagu/flagv calculation below we require a version of fmask without 
     1202      ! the land boundary condition (shlat) included: 
     1203      CALL wrk_alloc(jpi,jpj,zfmask)  
     1204      DO ij = 2, jpjm1 
     1205         DO ii = 2, jpim1 
     1206            zfmask(ii,ij) = tmask(ii,ij  ,1) * tmask(ii+1,ij  ,1)   & 
     1207           &              * tmask(ii,ij+1,1) * tmask(ii+1,ij+1,1) 
     1208         END DO       
     1209      END DO 
     1210 
    10941211      ! Lateral boundary conditions 
     1212      CALL lbc_lnk( zfmask       , 'F', 1. ) 
    10951213      CALL lbc_lnk( fmask        , 'F', 1. )   ;   CALL lbc_lnk( bdytmask(:,:), 'T', 1. ) 
    10961214      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) 
    … …  
    10981216      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components 
    10991217 
    1100          idx_bdy(ib_bdy)%flagu(:) = 0.e0 
    1101          idx_bdy(ib_bdy)%flagv(:) = 0.e0 
     1218         idx_bdy(ib_bdy)%flagu(:,:) = 0.e0 
     1219         idx_bdy(ib_bdy)%flagv(:,:) = 0.e0 
    11021220         icount = 0  
    11031221 
    1104          !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward 
    1105          !flagu =  0 : u is tangential 
    1106          !flagu =  1 : u is normal to the boundary and is direction is inward 
     1222         ! Calculate relationship of U direction to the local orientation of the boundary 
     1223         ! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward 
     1224         ! flagu =  0 : u is tangential 
     1225         ! flagu =  1 : u is normal to the boundary and is direction is inward 
    11071226   
    1108          igrd = 2      ! u-component  
    1109          DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   
    1110             nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
    1111             nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
    1112             zefl = bdytmask(nbi  ,nbj) 
    1113             zwfl = bdytmask(nbi+1,nbj) 
    1114             IF( zefl + zwfl == 2 ) THEN 
    1115                icount = icount + 1 
    1116             ELSE 
    1117                idx_bdy(ib_bdy)%flagu(ib)=-zefl+zwfl 
    1118             ENDIF 
     1227         DO igrd = 1,jpbgrd  
     1228            SELECT CASE( igrd ) 
     1229               CASE( 1 ) 
     1230                  pmask => umask(:,:,1) 
     1231                  i_offset = 0 
     1232               CASE( 2 )  
     1233                  pmask => bdytmask 
     1234                  i_offset = 1 
     1235               CASE( 3 )  
     1236                  pmask => zfmask(:,:) 
     1237                  i_offset = 0 
     1238            END SELECT  
     1239            icount = 0 
     1240            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   
     1241               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     1242               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
     1243               zefl = pmask(nbi+i_offset-1,nbj) 
     1244               zwfl = pmask(nbi+i_offset,nbj) 
     1245               ! This error check only works if you are using the bdyXmask arrays 
     1246               IF( i_offset == 1 .and. zefl + zwfl == 2 ) THEN 
     1247                  icount = icount + 1 
     1248                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj) 
     1249               ELSE 
     1250                  idx_bdy(ib_bdy)%flagu(ib,igrd) = -zefl + zwfl 
     1251               ENDIF 
     1252            END DO 
     1253            IF( icount /= 0 ) THEN 
     1254               IF(lwp) WRITE(numout,*) 
     1255               IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid(igrd),' grid points,',   & 
     1256                  ' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy 
     1257               IF(lwp) WRITE(numout,*) ' ========== ' 
     1258               IF(lwp) WRITE(numout,*) 
     1259               nstop = nstop + 1 
     1260            ENDIF  
    11191261         END DO 
    11201262 
    1121          !flagv = -1 : u component is normal to the dynamical boundary but its direction is outward 
    1122          !flagv =  0 : u is tangential 
    1123          !flagv =  1 : u is normal to the boundary and is direction is inward 
    1124  
    1125          igrd = 3      ! v-component 
    1126          DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   
    1127             nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
    1128             nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
    1129             znfl = bdytmask(nbi,nbj  ) 
    1130             zsfl = bdytmask(nbi,nbj+1) 
    1131             IF( znfl + zsfl == 2 ) THEN 
    1132                icount = icount + 1 
    1133             ELSE 
    1134                idx_bdy(ib_bdy)%flagv(ib) = -znfl + zsfl 
    1135             END IF 
     1263         ! Calculate relationship of V direction to the local orientation of the boundary 
     1264         ! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward 
     1265         ! flagv =  0 : v is tangential 
     1266         ! flagv =  1 : v is normal to the boundary and is direction is inward 
     1267 
     1268         DO igrd = 1,jpbgrd  
     1269            SELECT CASE( igrd ) 
     1270               CASE( 1 ) 
     1271                  pmask => vmask(:,:,1) 
     1272                  j_offset = 0 
     1273               CASE( 2 ) 
     1274                  pmask => zfmask(:,:) 
     1275                  j_offset = 0 
     1276               CASE( 3 ) 
     1277                  pmask => bdytmask 
     1278                  j_offset = 1 
     1279            END SELECT  
     1280            icount = 0 
     1281            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   
     1282               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     1283               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
     1284               znfl = pmask(nbi,nbj+j_offset-1  ) 
     1285               zsfl = pmask(nbi,nbj+j_offset) 
     1286               ! This error check only works if you are using the bdyXmask arrays 
     1287               IF( j_offset == 1 .and. znfl + zsfl == 2 ) THEN 
     1288                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj) 
     1289                  icount = icount + 1 
     1290               ELSE 
     1291                  idx_bdy(ib_bdy)%flagv(ib,igrd) = -znfl + zsfl 
     1292               END IF 
     1293            END DO 
     1294            IF( icount /= 0 ) THEN 
     1295               IF(lwp) WRITE(numout,*) 
     1296               IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid(igrd),' grid points,',   & 
     1297                  ' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy 
     1298               IF(lwp) WRITE(numout,*) ' ========== ' 
     1299               IF(lwp) WRITE(numout,*) 
     1300               nstop = nstop + 1 
     1301            ENDIF  
    11361302         END DO 
    11371303 
    1138          IF( icount /= 0 ) THEN 
    1139             IF(lwp) WRITE(numout,*) 
    1140             IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,',   & 
    1141                ' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_bdy 
    1142             IF(lwp) WRITE(numout,*) ' ========== ' 
    1143             IF(lwp) WRITE(numout,*) 
    1144             nstop = nstop + 1 
    1145          ENDIF  
    1146      
    1147       ENDDO 
     1304      END DO 
    11481305 
    11491306      ! Compute total lateral surface for volume correction: 
    … …  
    11571314               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
    11581315               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
    1159                flagu => idx_bdy(ib_bdy)%flagu(ib) 
     1316               flagu => idx_bdy(ib_bdy)%flagu(ib,igrd) 
    11601317               bdysurftot = bdysurftot + hu     (nbi  , nbj)                           & 
    11611318                  &                    * e2u    (nbi  , nbj) * ABS( flagu ) & 
    … …  
    11701327               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
    11711328               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
    1172                flagv => idx_bdy(ib_bdy)%flagv(ib) 
     1329               flagv => idx_bdy(ib_bdy)%flagv(ib,igrd) 
    11731330               bdysurftot = bdysurftot + hv     (nbi, nbj  )                           & 
    11741331                  &                    * e1v    (nbi, nbj  ) * ABS( flagv ) & 
    … …  
    11861343         DEALLOCATE(nbidta, nbjdta, nbrdta) 
    11871344      ENDIF 
     1345 
     1346      CALL wrk_dealloc(jpi,jpj,zfmask)  
    11881347 
    11891348      IF( nn_timing == 1 ) CALL timing_stop('bdy_init') 
    … …  
    15801739      itest = 0 
    15811740 
    1582       IF (nn_dyn2d(ib1)/=nn_dyn2d(ib2)) itest = itest + 1 
    1583       IF (nn_dyn3d(ib1)/=nn_dyn3d(ib2)) itest = itest + 1 
    1584       IF (nn_tra(ib1)/=nn_tra(ib2)) itest = itest + 1 
     1741      IF (cn_dyn2d(ib1)/=cn_dyn2d(ib2)) itest = itest + 1 
     1742      IF (cn_dyn3d(ib1)/=cn_dyn3d(ib2)) itest = itest + 1 
     1743      IF (cn_tra(ib1)/=cn_tra(ib2)) itest = itest + 1 
    15851744      ! 
    15861745      IF (nn_dyn2d_dta(ib1)/=nn_dyn2d_dta(ib2)) itest = itest + 1 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

    r4147 r4292  
    99   !!            3.3  !  2010-09  (D.Storkey and E.O'Dea)  bug fixes 
    1010   !!            3.4  !  2012-09  (G. Reffray and J. Chanut) New inputs + mods 
     11   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes 
    1112   !!---------------------------------------------------------------------- 
    1213#if defined key_bdy 
    … …  
    3233!   USE tide_mod       ! Useless ?? 
    3334   USE fldread, ONLY: fld_map 
     35   USE dynspg_oce, ONLY: lk_dynspg_ts 
    3436 
    3537   IMPLICIT NONE 
    … …  
    3840   PUBLIC   bdytide_init     ! routine called in bdy_init 
    3941   PUBLIC   bdytide_update   ! routine called in bdy_dta 
     42   PUBLIC   bdy_dta_tides    ! routine called in dyn_spg_ts 
    4043 
    4144   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data 
    … …  
    4952 
    5053   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data 
     54   TYPE(OBC_DATA)  , PRIVATE, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component) 
    5155 
    5256   !!---------------------------------------------------------------------- 
    … …  
    131135            ! JC: If FRS scheme is used, we assume that tidal is needed over the whole 
    132136            ! relaxation area       
    133             IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     137            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN 
    134138               ilen0(:)=nblen(:) 
    135139            ELSE 
    … …  
    146150            ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) ) 
    147151 
    148             td%ssh0(:,:,:) = 0.e0 
    149             td%ssh(:,:,:) = 0.e0 
    150             td%u0(:,:,:) = 0.e0 
    151             td%u(:,:,:) = 0.e0 
    152             td%v0(:,:,:) = 0.e0 
    153             td%v(:,:,:) = 0.e0 
     152            td%ssh0(:,:,:) = 0._wp 
     153            td%ssh (:,:,:) = 0._wp 
     154            td%u0  (:,:,:) = 0._wp 
     155            td%u   (:,:,:) = 0._wp 
     156            td%v0  (:,:,:) = 0._wp 
     157            td%v   (:,:,:) = 0._wp 
    154158 
    155159            IF (ln_bdytide_2ddta) THEN 
    … …  
    255259            ENDIF 
    256260            ! 
     261            IF ( lk_dynspg_ts ) THEN ! Allocate arrays to save slowly varying boundary data during 
     262                                     ! time splitting integration 
     263               ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) ) 
     264               ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) ) 
     265               ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) ) 
     266               dta_bdy_s(ib_bdy)%ssh(:) = 0.e0 
     267               dta_bdy_s(ib_bdy)%u2d(:) = 0.e0 
     268               dta_bdy_s(ib_bdy)%v2d(:) = 0.e0 
     269            ENDIF 
     270            ! 
    257271         ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2 
    258272         ! 
    … …  
    300314      ENDIF 
    301315 
    302       IF ( nsec_day == NINT(0.5 * rdttra(1)) .AND. zflag==1 ) THEN 
     316      IF ( nsec_day == NINT(0.5_wp * rdttra(1)) .AND. zflag==1 ) THEN 
    303317        ! 
    304318        kt_tide = kt 
    … …  
    321335          
    322336      IF( PRESENT(jit) ) THEN   
    323          z_arg = ( ((kt-kt_tide)-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) ) 
     337         z_arg = ((kt-kt_tide) * rdt + (jit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) ) 
    324338      ELSE                               
    325339         z_arg = ((kt-kt_tide)+time_add) * rdt 
    … …  
    327341 
    328342      ! Linear ramp on tidal component at open boundaries  
    329       zramp = 1. 
    330       IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0.),1.) 
     343      zramp = 1._wp 
     344      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0._wp),1._wp) 
    331345 
    332346      DO itide = 1, nb_harmo 
    … …  
    354368      ! 
    355369   END SUBROUTINE bdytide_update 
     370 
     371   SUBROUTINE bdy_dta_tides( kt, kit, time_offset ) 
     372      !!---------------------------------------------------------------------- 
     373      !!                 ***  SUBROUTINE bdy_dta_tides  *** 
     374      !!                 
     375      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.  
     376      !!                 
     377      !!---------------------------------------------------------------------- 
     378      INTEGER, INTENT( in )            ::   kt          ! Main timestep counter 
     379      INTEGER, INTENT( in ),OPTIONAL   ::   kit         ! Barotropic timestep counter (for timesplitting option) 
     380      INTEGER, INTENT( in ),OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if kit 
     381                                                        ! is present then units = subcycle timesteps. 
     382                                                        ! time_offset = 0  => get data at "now"    time level 
     383                                                        ! time_offset = -1 => get data at "before" time level 
     384                                                        ! time_offset = +1 => get data at "after"  time level 
     385                                                        ! etc. 
     386      !! 
     387      LOGICAL  :: lk_first_btstp  ! =.TRUE. if time splitting and first barotropic step 
     388      INTEGER,          DIMENSION(jpbgrd) :: ilen0  
     389      INTEGER, DIMENSION(1:jpbgrd) :: nblen, nblenrim  ! short cuts 
     390      INTEGER  :: itide, ib_bdy, ib, igrd                     ! loop indices 
     391      INTEGER  :: time_add                                    ! time offset in units of timesteps 
     392      REAL(wp) :: z_arg, z_sarg, zramp, zoff, z_cost, z_sist       
     393      !!---------------------------------------------------------------------- 
     394 
     395      IF( nn_timing == 1 ) CALL timing_start('bdy_dta_tides') 
     396 
     397      lk_first_btstp=.TRUE. 
     398      IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF 
     399 
     400      time_add = 0 
     401      IF( PRESENT(time_offset) ) THEN 
     402         time_add = time_offset 
     403      ENDIF 
     404       
     405      ! Absolute time from model initialization:    
     406      IF( PRESENT(kit) ) THEN   
     407         z_arg = ( kt + (kit+0.5_wp*(time_add-1)) / REAL(nn_baro,wp) ) * rdt 
     408      ELSE                               
     409         z_arg = ( kt + time_add ) * rdt 
     410      ENDIF 
     411 
     412      ! Linear ramp on tidal component at open boundaries  
     413      zramp = 1. 
     414      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.) 
     415 
     416      DO ib_bdy = 1,nb_bdy 
     417 
     418         ! line below should be simplified (runoff case) 
     419!! CHANUT: TO BE SORTED OUT 
     420!!         IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(nn_tra(ib_bdy).NE.4)) THEN 
     421         IF ( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN 
     422 
     423            nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd) 
     424            nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd) 
     425 
     426            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN 
     427               ilen0(:)=nblen(:) 
     428            ELSE 
     429               ilen0(:)=nblenrim(:) 
     430            ENDIF      
     431 
     432            ! We refresh nodal factors every day below 
     433            ! This should be done somewhere else 
     434            IF ( nsec_day == NINT(0.5_wp * rdttra(1)) .AND. lk_first_btstp ) THEN 
     435               ! 
     436               kt_tide = kt                
     437               ! 
     438               IF(lwp) THEN 
     439               WRITE(numout,*) 
     440               WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt 
     441               WRITE(numout,*) '~~~~~~~~~~~~~~ ' 
     442               ENDIF 
     443               ! 
     444               CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) ) 
     445               CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) ) 
     446               ! 
     447            ENDIF 
     448            zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time 
     449            ! 
     450            ! If time splitting, save data at first barotropic iteration 
     451            IF ( PRESENT(kit) ) THEN 
     452               IF ( lk_first_btstp ) THEN ! Save slow varying open boundary data: 
     453                  dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy(ib_bdy)%ssh(1:ilen0(1)) 
     454                  dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy(ib_bdy)%u2d(1:ilen0(2)) 
     455                  dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy(ib_bdy)%v2d(1:ilen0(3)) 
     456 
     457               ELSE ! Initialize arrays from slow varying open boundary data:             
     458                  dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) 
     459                  dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) 
     460                  dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) 
     461               ENDIF 
     462            ENDIF 
     463            ! 
     464            ! Update open boundary data arrays: 
     465            DO itide = 1, nb_harmo 
     466               ! 
     467               z_sarg = (z_arg + zoff) * omega_tide(itide) 
     468               z_cost = zramp * COS( z_sarg ) 
     469               z_sist = zramp * SIN( z_sarg ) 
     470               ! 
     471               igrd=1                              ! SSH on tracer grid 
     472               DO ib = 1, ilen0(igrd) 
     473                  dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + & 
     474                     &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + & 
     475                     &                        tides(ib_bdy)%ssh(ib,itide,2)*z_sist ) 
     476               END DO 
     477               ! 
     478               igrd=2                              ! U grid 
     479               DO ib = 1, ilen0(igrd) 
     480                  dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + & 
     481                     &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + & 
     482                     &                        tides(ib_bdy)%u(ib,itide,2)*z_sist ) 
     483               END DO 
     484               ! 
     485               igrd=3                              ! V grid 
     486               DO ib = 1, ilen0(igrd)  
     487                  dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + & 
     488                     &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + & 
     489                     &                        tides(ib_bdy)%v(ib,itide,2)*z_sist ) 
     490               END DO 
     491            END DO 
     492         END IF 
     493      END DO 
     494      ! 
     495      IF( nn_timing == 1 ) CALL timing_stop('bdy_dta_tides') 
     496      ! 
     497   END SUBROUTINE bdy_dta_tides 
    356498 
    357499   SUBROUTINE tide_init_elevation( idx, td ) 
    … …  
    460602      WRITE(*,*) 'bdytide_update: You should not have seen this print! error?', kt, jit 
    461603   END SUBROUTINE bdytide_update 
     604   SUBROUTINE bdy_dta_tides( kt, kit, time_offset )     ! Empty routine 
     605      INTEGER, INTENT( in )            ::   kt          ! Dummy argument empty routine       
     606      INTEGER, INTENT( in ),OPTIONAL   ::   kit         ! Dummy argument empty routine 
     607      INTEGER, INTENT( in ),OPTIONAL   ::   time_offset ! Dummy argument empty routine 
     608      WRITE(*,*) 'bdy_dta_tides: You should not have seen this print! error?', kt, jit 
     609   END SUBROUTINE bdy_dta_tides 
    462610#endif 
    463611 
    464612   !!====================================================================== 
    465613END MODULE bdytides 
     614 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90

    r3777 r4292  
    2020   USE dom_oce         ! ocean space and time domain variables  
    2121   USE bdy_oce         ! ocean open boundary conditions 
     22   USE bdylib          ! for orlanski library routines 
    2223   USE bdydta, ONLY:   bf 
    2324   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    … …  
    5152      DO ib_bdy=1, nb_bdy 
    5253 
    53          SELECT CASE( nn_tra(ib_bdy) ) 
    54          CASE(jp_none) 
     54         SELECT CASE( cn_tra(ib_bdy) ) 
     55         CASE('none') 
    5556            CYCLE 
    56          CASE(jp_frs) 
     57         CASE('frs') 
    5758            CALL bdy_tra_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
    58          CASE(2) 
     59         CASE('specified') 
    5960            CALL bdy_tra_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
    60          CASE(3) 
     61         CASE('neumann') 
    6162            CALL bdy_tra_nmn( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
    62          CASE(4) 
     63         CASE('orlanski') 
     64            CALL bdy_tra_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ll_npo=.false. ) 
     65         CASE('orlanski_npo') 
     66            CALL bdy_tra_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ll_npo=.true. ) 
     67         CASE('runoff') 
    6368            CALL bdy_tra_rnf( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
    6469         CASE DEFAULT 
    … …  
    196201      ! 
    197202   END SUBROUTINE bdy_tra_nmn 
     203  
     204 
     205   SUBROUTINE bdy_tra_orlanski( idx, dta, ll_npo ) 
     206      !!---------------------------------------------------------------------- 
     207      !!                 ***  SUBROUTINE bdy_tra_orlanski  *** 
     208      !!              
     209      !!              - Apply Orlanski radiation to temperature and salinity.  
     210      !!              - Wrapper routine for bdy_orlanski_3d 
     211      !!  
     212      !! 
     213      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)     
     214      !!---------------------------------------------------------------------- 
     215      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices 
     216      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data 
     217      LOGICAL,                      INTENT(in) ::   ll_npo  ! switch for NPO version 
     218 
     219      INTEGER  ::   igrd                                    ! grid index 
     220      !!---------------------------------------------------------------------- 
     221 
     222      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_orlanski') 
     223      ! 
     224      igrd = 1      ! Orlanski bc on temperature;  
     225      !             
     226      CALL bdy_orlanski_3d( idx, igrd, tsb(:,:,:,jp_tem), tsa(:,:,:,jp_tem), dta%tem, ll_npo ) 
     227 
     228      igrd = 1      ! Orlanski bc on salinity; 
     229      !   
     230      CALL bdy_orlanski_3d( idx, igrd, tsb(:,:,:,jp_sal), tsa(:,:,:,jp_sal), dta%sal, ll_npo ) 
     231      ! 
     232      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_orlanski') 
     233      ! 
     234 
     235   END SUBROUTINE bdy_tra_orlanski 
     236 
    198237 
    199238   SUBROUTINE bdy_tra_rnf( idx, dta, kt ) 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90

    r3294 r4292  
    104104               ii = idx%nbi(jb,jgrd) 
    105105               ij = idx%nbj(jb,jgrd) 
    106                zubtpecor = zubtpecor + idx%flagu(jb) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 
     106               zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 
    107107            END DO 
    108108         END DO 
    … …  
    112112               ii = idx%nbi(jb,jgrd) 
    113113               ij = idx%nbj(jb,jgrd) 
    114                zubtpecor = zubtpecor + idx%flagv(jb) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk)  
     114               zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk)  
    115115            END DO 
    116116         END DO 
    … …  
    136136               ii = idx%nbi(jb,jgrd) 
    137137               ij = idx%nbj(jb,jgrd) 
    138                ua(ii,ij,jk) = ua(ii,ij,jk) - idx%flagu(jb) * zubtpecor * umask(ii,ij,jk) 
    139                ztranst = ztranst + idx%flagu(jb) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 
     138               ua(ii,ij,jk) = ua(ii,ij,jk) - idx%flagu(jb,jgrd) * zubtpecor * umask(ii,ij,jk) 
     139               ztranst = ztranst + idx%flagu(jb,jgrd) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 
    140140            END DO 
    141141         END DO 
    … …  
    145145               ii = idx%nbi(jb,jgrd) 
    146146               ij = idx%nbj(jb,jgrd) 
    147                va(ii,ij,jk) = va(ii,ij,jk) -idx%flagv(jb) * zubtpecor * vmask(ii,ij,jk) 
    148                ztranst = ztranst + idx%flagv(jb) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 
     147               va(ii,ij,jk) = va(ii,ij,jk) -idx%flagv(jb,jgrd) * zubtpecor * vmask(ii,ij,jk) 
     148               ztranst = ztranst + idx%flagv(jb,jgrd) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 
    149149            END DO 
    150150         END DO 
Note: See TracChangeset for help on using the changeset viewer.