Ticket #2063: bdydta.F90

File bdydta.F90, 37.6 KB (added by molines, 3 years ago)

bdydta.F90 fix for dev_merge_2017

Line 
1MODULE bdydta
2   !!======================================================================
3   !!                       ***  MODULE bdydta  ***
4   !! Open boundary data : read the data for the unstructured open boundaries.
5   !!======================================================================
6   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code
7   !!             -   !  2007-01  (D. Storkey) Update to use IOM module
8   !!             -   !  2007-07  (D. Storkey) add bdy_dta_fla
9   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
10   !!            3.3  !  2010-09  (E.O'Dea) modifications for Shelf configurations
11   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions
12   !!            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
14   !!----------------------------------------------------------------------
15
16   !!----------------------------------------------------------------------
17   !!    bdy_dta      : read external data along open boundaries from file
18   !!    bdy_dta_init : initialise arrays etc for reading of external data
19   !!----------------------------------------------------------------------
20   USE oce            ! ocean dynamics and tracers
21   USE dom_oce        ! ocean space and time domain
22   USE phycst         ! physical constants
23   USE sbcapr         ! atmospheric pressure forcing
24   USE sbctide        ! Tidal forcing or not
25   USE bdy_oce        ! ocean open boundary conditions 
26   USE bdytides       ! tidal forcing at boundaries
27#if defined key_lim3
28   USE ice            ! sea-ice variables
29   USE icevar         ! redistribute ice input into categories
30#endif
31   !
32   USE fldread        ! read input fields
33   USE iom            ! IOM library
34   USE in_out_manager ! I/O logical units
35   USE timing         ! Timing
36
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC   bdy_dta          ! routine called by step.F90 and dynspg_ts.F90
41   PUBLIC   bdy_dta_init     ! routine called by nemogcm.F90
42
43   INTEGER, ALLOCATABLE, DIMENSION(:)   ::   nb_bdy_fld        ! Number of fields to update for each boundary set.
44   INTEGER                              ::   nb_bdy_fld_sum    ! Total number of fields to update for all boundary sets.
45   LOGICAL,           DIMENSION(jp_bdy) ::   ln_full_vel_array ! =T => full velocities in 3D boundary conditions
46                                                               ! =F => baroclinic velocities in 3D boundary conditions
47!$AGRIF_DO_NOT_TREAT
48   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:), TARGET ::   bf        ! structure of input fields (file informations, fields read)
49!$AGRIF_END_DO_NOT_TREAT
50   TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:) :: nbmap_ptr   ! array of pointers to nbmap
51
52#if defined key_lim3
53   INTEGER ::   nice_cat                      ! number of categories in the input file
54   INTEGER ::   jfld_hti, jfld_hts, jfld_ai   ! indices of ice thickness, snow thickness and concentration in bf structure
55   INTEGER, DIMENSION(jp_bdy) :: jfld_htit, jfld_htst, jfld_ait
56#endif
57
58   !!----------------------------------------------------------------------
59   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
60   !! $Id: bdydta.F90 9168 2017-12-23 12:27:17Z gm $
61   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
62   !!----------------------------------------------------------------------
63CONTAINS
64
65      SUBROUTINE bdy_dta( kt, jit, time_offset )
66      !!----------------------------------------------------------------------
67      !!                   ***  SUBROUTINE bdy_dta  ***
68      !!                   
69      !! ** Purpose :   Update external data for open boundary conditions
70      !!
71      !! ** Method  :   Use fldread.F90
72      !!               
73      !!----------------------------------------------------------------------
74      INTEGER, INTENT(in)           ::   kt           ! ocean time-step index
75      INTEGER, INTENT(in), OPTIONAL ::   jit          ! subcycle time-step index (for timesplitting option)
76      INTEGER, INTENT(in), OPTIONAL ::   time_offset  ! time offset in units of timesteps. NB. if jit
77      !                                               ! is present then units = subcycle timesteps.
78      !                                               ! time_offset = 0 => get data at "now" time level
79      !                                               ! time_offset = -1 => get data at "before" time level
80      !                                               ! time_offset = +1 => get data at "after" time level
81      !                                               ! etc.
82      !
83      INTEGER ::  jbdy, jfld, jstart, jend, ib, jl  ! dummy loop indices
84      INTEGER ::  ii, ij, ik, igrd                  ! local integers
85      INTEGER,          DIMENSION(jpbgrd) ::   ilen1 
86      INTEGER, POINTER, DIMENSION(:)      ::   nblen, nblenrim  ! short cuts
87      TYPE(OBC_DATA), POINTER             ::   dta              ! short cut
88      !!---------------------------------------------------------------------------
89      !
90      IF( ln_timing )   CALL timing_start('bdy_dta')
91      !
92      ! Initialise data arrays once for all from initial conditions where required
93      !---------------------------------------------------------------------------
94      IF( kt == nit000 .AND. .NOT.PRESENT(jit) ) THEN
95
96         ! Calculate depth-mean currents
97         !-----------------------------
98         
99         DO jbdy = 1, nb_bdy
100            !
101            nblen    => idx_bdy(jbdy)%nblen
102            nblenrim => idx_bdy(jbdy)%nblenrim
103            dta      => dta_bdy(jbdy)
104            !
105            IF( nn_dyn2d_dta(jbdy) == 0 ) THEN
106               ilen1(:) = nblen(:)
107               IF( dta%ll_ssh ) THEN
108                  igrd = 1
109                  DO ib = 1, ilen1(igrd)
110                     ii = idx_bdy(jbdy)%nbi(ib,igrd)
111                     ij = idx_bdy(jbdy)%nbj(ib,igrd)
112                     dta_bdy(jbdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)         
113                  END DO
114               ENDIF
115               IF( dta%ll_u2d ) THEN
116                  igrd = 2
117                  DO ib = 1, ilen1(igrd)
118                     ii = idx_bdy(jbdy)%nbi(ib,igrd)
119                     ij = idx_bdy(jbdy)%nbj(ib,igrd)
120                     dta_bdy(jbdy)%u2d(ib) = un_b(ii,ij) * umask(ii,ij,1)         
121                  END DO
122               ENDIF
123               IF( dta%ll_v2d ) THEN
124                  igrd = 3
125                  DO ib = 1, ilen1(igrd)
126                     ii = idx_bdy(jbdy)%nbi(ib,igrd)
127                     ij = idx_bdy(jbdy)%nbj(ib,igrd)
128                     dta_bdy(jbdy)%v2d(ib) = vn_b(ii,ij) * vmask(ii,ij,1)         
129                  END DO
130               ENDIF
131            ENDIF
132            !
133            IF( nn_dyn3d_dta(jbdy) == 0 ) THEN
134               ilen1(:) = nblen(:)
135               IF( dta%ll_u3d ) THEN
136                  igrd = 2 
137                  DO ib = 1, ilen1(igrd)
138                     DO ik = 1, jpkm1
139                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
140                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
141                        dta_bdy(jbdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - un_b(ii,ij) ) * umask(ii,ij,ik)         
142                     END DO
143                  END DO
144               ENDIF
145               IF( dta%ll_v3d ) THEN
146                  igrd = 3 
147                  DO ib = 1, ilen1(igrd)
148                     DO ik = 1, jpkm1
149                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
150                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
151                        dta_bdy(jbdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - vn_b(ii,ij) ) * vmask(ii,ij,ik)         
152                        END DO
153                  END DO
154               ENDIF
155            ENDIF
156
157            IF( nn_tra_dta(jbdy) == 0 ) THEN
158               ilen1(:) = nblen(:)
159               IF( dta%ll_tem ) THEN
160                  igrd = 1 
161                  DO ib = 1, ilen1(igrd)
162                     DO ik = 1, jpkm1
163                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
164                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
165                        dta_bdy(jbdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)         
166                     END DO
167                  END DO
168               ENDIF
169               IF( dta%ll_sal ) THEN
170                  igrd = 1 
171                  DO ib = 1, ilen1(igrd)
172                     DO ik = 1, jpkm1
173                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
174                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
175                        dta_bdy(jbdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)         
176                     END DO
177                  END DO
178               ENDIF
179            ENDIF
180
181#if defined key_lim3
182            IF( nn_ice_lim_dta(jbdy) == 0 ) THEN    ! set ice to initial values
183               ilen1(:) = nblen(:)
184               IF( dta%ll_a_i ) THEN
185                  igrd = 1   
186                  DO jl = 1, jpl
187                     DO ib = 1, ilen1(igrd)
188                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
189                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
190                        dta_bdy(jbdy)%a_i (ib,jl) =  a_i(ii,ij,jl) * tmask(ii,ij,1) 
191                     END DO
192                  END DO
193               ENDIF
194               IF( dta%ll_h_i ) THEN
195                  igrd = 1   
196                  DO jl = 1, jpl
197                     DO ib = 1, ilen1(igrd)
198                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
199                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
200                        dta_bdy(jbdy)%h_i (ib,jl) =  h_i(ii,ij,jl) * tmask(ii,ij,1) 
201                     END DO
202                  END DO
203               ENDIF
204               IF( dta%ll_h_s ) THEN
205                  igrd = 1   
206                  DO jl = 1, jpl
207                     DO ib = 1, ilen1(igrd)
208                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
209                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
210                        dta_bdy(jbdy)%h_s (ib,jl) =  h_s(ii,ij,jl) * tmask(ii,ij,1) 
211                     END DO
212                  END DO
213               ENDIF
214            ENDIF
215#endif
216         END DO ! jbdy
217         !
218      ENDIF ! kt == nit000
219
220      ! update external data from files
221      !--------------------------------
222     
223      jstart = 1
224      DO jbdy = 1, nb_bdy   
225         dta => dta_bdy(jbdy)
226         IF( nn_dta(jbdy) == 1 ) THEN ! skip this bit if no external data required
227     
228            IF( PRESENT(jit) ) THEN
229               ! Update barotropic boundary conditions only
230               ! jit is optional argument for fld_read and bdytide_update
231               IF( cn_dyn2d(jbdy) /= 'none' ) THEN
232                  IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
233                     IF( dta%ll_ssh ) dta%ssh(:) = 0._wp
234                     IF( dta%ll_u2d ) dta%u2d(:) = 0._wp
235                     IF( dta%ll_u3d ) dta%v2d(:) = 0._wp
236                  ENDIF
237                  IF (cn_tra(jbdy) /= 'runoff') THEN
238                     IF( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 ) THEN
239
240                        jend = jstart + dta%nread(2) - 1
241                        IF( ln_full_vel_array(jbdy) ) THEN
242                           CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  &
243                                     & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy,   &
244                                     & fvl=ln_full_vel_array(jbdy)  )
245                        ELSE
246                           CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  &
247                                     & kit=jit, kt_offset=time_offset  )
248                        ENDIF
249
250                        ! If full velocities in boundary data then extract barotropic velocities from 3D fields
251                        IF( ln_full_vel_array(jbdy) .AND.                                             &
252                          &    ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR.  &
253                          &      nn_dyn3d_dta(jbdy) == 1 ) )THEN
254
255                           igrd = 2                      ! zonal velocity
256                           dta%u2d(:) = 0._wp
257                           DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
258                              ii   = idx_bdy(jbdy)%nbi(ib,igrd)
259                              ij   = idx_bdy(jbdy)%nbj(ib,igrd)
260                              DO ik = 1, jpkm1
261                                 dta%u2d(ib) = dta%u2d(ib) &
262                       &                          + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik)
263                              END DO
264                              dta%u2d(ib) =  dta%u2d(ib) * r1_hu_n(ii,ij)
265                           END DO
266                           igrd = 3                      ! meridional velocity
267                           dta%v2d(:) = 0._wp
268                           DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
269                              ii   = idx_bdy(jbdy)%nbi(ib,igrd)
270                              ij   = idx_bdy(jbdy)%nbj(ib,igrd)
271                              DO ik = 1, jpkm1
272                                 dta%v2d(ib) = dta%v2d(ib) &
273                       &                       + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik)
274                              END DO
275                              dta%v2d(ib) =  dta%v2d(ib) * r1_hv_n(ii,ij)
276                           END DO
277                        ENDIF                   
278                     ENDIF
279                     IF( nn_dyn2d_dta(jbdy) .ge. 2 ) THEN ! update tidal harmonic forcing
280                        CALL bdytide_update( kt=kt, idx=idx_bdy(jbdy), dta=dta, td=tides(jbdy),   & 
281                          &                 jit=jit, time_offset=time_offset )
282                     ENDIF
283                  ENDIF
284               ENDIF
285            ELSE
286               IF (cn_tra(jbdy) == 'runoff') then      ! runoff condition
287                  jend = nb_bdy_fld(jbdy)
288                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend),  &
289                               & map=nbmap_ptr(jstart:jend), kt_offset=time_offset )
290                  !
291                  igrd = 2                      ! zonal velocity
292                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
293                     ii   = idx_bdy(jbdy)%nbi(ib,igrd)
294                     ij   = idx_bdy(jbdy)%nbj(ib,igrd)
295                     dta%u2d(ib) = dta%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
296                  END DO
297                  !
298                  igrd = 3                      ! meridional velocity
299                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
300                     ii   = idx_bdy(jbdy)%nbi(ib,igrd)
301                     ij   = idx_bdy(jbdy)%nbj(ib,igrd)
302                     dta%v2d(ib) = dta%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
303                  END DO
304               ELSE
305                  IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
306                     IF( dta%ll_ssh ) dta%ssh(:) = 0._wp
307                     IF( dta%ll_u2d ) dta%u2d(:) = 0._wp
308                     IF( dta%ll_v2d ) dta%v2d(:) = 0._wp
309                  ENDIF
310                  IF( dta%nread(1) .gt. 0 ) THEN ! update external data
311                     jend = jstart + dta%nread(1) - 1
312                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), &
313                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy,   &
314                                  & fvl=ln_full_vel_array(jbdy) )
315                  ENDIF
316                  ! If full velocities in boundary data then split into barotropic and baroclinic data
317                  IF( ln_full_vel_array(jbdy) .and.                                             &
318                    & ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR. &
319                    &   nn_dyn3d_dta(jbdy) == 1 ) ) THEN
320                     igrd = 2                      ! zonal velocity
321                     dta%u2d(:) = 0._wp
322                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
323                        ii   = idx_bdy(jbdy)%nbi(ib,igrd)
324                        ij   = idx_bdy(jbdy)%nbj(ib,igrd)
325                        DO ik = 1, jpkm1
326                           dta%u2d(ib) = dta%u2d(ib) &
327                 &                       + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik)
328                        END DO
329                        dta%u2d(ib) =  dta%u2d(ib) * r1_hu_n(ii,ij)
330                        DO ik = 1, jpkm1
331                           dta%u3d(ib,ik) = dta%u3d(ib,ik) - dta%u2d(ib)
332                        END DO
333                     END DO
334                     igrd = 3                      ! meridional velocity
335                     dta%v2d(:) = 0._wp
336                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
337                        ii   = idx_bdy(jbdy)%nbi(ib,igrd)
338                        ij   = idx_bdy(jbdy)%nbj(ib,igrd)
339                        DO ik = 1, jpkm1
340                           dta%v2d(ib) = dta%v2d(ib) &
341                 &                       + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik)
342                        END DO
343                        dta%v2d(ib) =  dta%v2d(ib) * r1_hv_n(ii,ij)
344                        DO ik = 1, jpkm1
345                           dta%v3d(ib,ik) = dta%v3d(ib,ik) - dta%v2d(ib)
346                        END DO
347                     END DO
348                  ENDIF
349
350               ENDIF
351#if defined key_lim3
352               IF( cn_ice_lim(jbdy) /= 'none' .AND. nn_ice_lim_dta(jbdy) == 1 ) THEN
353                  jfld_hti = jfld_htit(jbdy)
354                  jfld_hts = jfld_htst(jbdy)
355                  jfld_ai  = jfld_ait(jbdy)
356                  IF( nice_cat == 1 ) THEN ! case input cat = 1
357                     CALL ice_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), &
358                        &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    )
359                  ELSEIF( nice_cat /= 1 .AND. nice_cat /= jpl ) THEN ! case input cat /=1 and /=jpl
360                     CALL ice_var_itd2( bf(jfld_hti)%fnow(:,1,:), bf(jfld_hts)%fnow(:,1,:), bf(jfld_ai)%fnow(:,1,:), &
361                        &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    )
362                  ENDIF
363               ENDIF
364#endif
365            ENDIF
366            jstart = jstart + dta%nread(1)
367         ENDIF    ! nn_dta(jbdy) = 1
368      END DO  ! jbdy
369
370      IF ( ln_tide ) THEN
371         IF (ln_dynspg_ts) THEN      ! Fill temporary arrays with slow-varying bdy data                           
372            DO jbdy = 1, nb_bdy    ! Tidal component added in ts loop
373               IF ( nn_dyn2d_dta(jbdy) .ge. 2 ) THEN
374                  nblen => idx_bdy(jbdy)%nblen
375                  nblenrim => idx_bdy(jbdy)%nblenrim
376                  IF( cn_dyn2d(jbdy) == 'frs' ) THEN; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF
377                  IF ( dta_bdy(jbdy)%ll_ssh ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1))
378                  IF ( dta_bdy(jbdy)%ll_u2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2))
379                  IF ( dta_bdy(jbdy)%ll_v2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3))
380               ENDIF
381            END DO
382         ELSE ! Add tides if not split-explicit free surface else this is done in ts loop
383            !
384            CALL bdy_dta_tides( kt=kt, time_offset=time_offset )
385         ENDIF
386      ENDIF
387
388      IF ( ln_apr_obc ) THEN
389         DO jbdy = 1, nb_bdy
390            IF (cn_tra(jbdy) /= 'runoff')THEN
391               igrd = 1                      ! meridional velocity
392               DO ib = 1, idx_bdy(jbdy)%nblenrim(igrd)
393                  ii   = idx_bdy(jbdy)%nbi(ib,igrd)
394                  ij   = idx_bdy(jbdy)%nbj(ib,igrd)
395                  dta_bdy(jbdy)%ssh(ib) = dta_bdy(jbdy)%ssh(ib) + ssh_ib(ii,ij)
396               END DO
397            ENDIF
398         END DO
399      ENDIF
400      !
401      IF( ln_timing )   CALL timing_stop('bdy_dta')
402      !
403   END SUBROUTINE bdy_dta
404
405
406   SUBROUTINE bdy_dta_init
407      !!----------------------------------------------------------------------
408      !!                   ***  SUBROUTINE bdy_dta_init  ***
409      !!                   
410      !! ** Purpose :   Initialise arrays for reading of external data
411      !!                for open boundary conditions
412      !!
413      !! ** Method  :   
414      !!               
415      !!----------------------------------------------------------------------
416      INTEGER ::   jbdy, jfld, jstart, jend, ierror, ios     ! Local integers
417      !
418      CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files
419      CHARACTER(len=100), DIMENSION(nb_bdy)  ::   cn_dir_array  ! Root directory for location of data files
420      CHARACTER(len = 256)::   clname                           ! temporary file name
421      LOGICAL                                ::   ln_full_vel   ! =T => full velocities in 3D boundary data
422      !                                                         ! =F => baroclinic velocities in 3D boundary data
423      INTEGER                                ::   ilen_global   ! Max length required for global bdy dta arrays
424      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ilen1, ilen3  ! size of 1st and 3rd dimensions of local arrays
425      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ibdy           ! bdy set for a particular jfld
426      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   igrid         ! index for grid type (1,2,3 = T,U,V)
427      INTEGER, POINTER, DIMENSION(:)         ::   nblen, nblenrim  ! short cuts
428      TYPE(OBC_DATA), POINTER                ::   dta           ! short cut
429#if defined key_lim3
430      INTEGER               ::   kndims   ! number of dimensions in an array (i.e. 3 = wo ice cat; 4 = w ice cat)
431      INTEGER, DIMENSION(4) ::   kdimsz   ! size   of dimensions
432      INTEGER               ::   inum,id1 ! local integer
433#endif
434      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   blf_i         !  array of namelist information structures
435      TYPE(FLD_N) ::   bn_tem, bn_sal, bn_u3d, bn_v3d   !
436      TYPE(FLD_N) ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read
437#if defined key_lim3
438      TYPE(FLD_N) ::   bn_a_i, bn_h_i, bn_h_s     
439#endif
440      NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d 
441#if defined key_lim3
442      NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s
443#endif
444      NAMELIST/nambdy_dta/ ln_full_vel, nb_jpk_bdy
445      !!---------------------------------------------------------------------------
446      !
447      IF(lwp) WRITE(numout,*)
448      IF(lwp) WRITE(numout,*) 'bdy_dta_ini : initialization of data at the open boundaries'
449      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
450      IF(lwp) WRITE(numout,*) ''
451
452      ! Set nn_dta
453      DO jbdy = 1, nb_bdy
454         nn_dta(jbdy) = MAX(   nn_dyn2d_dta  (jbdy)    &
455            &                , nn_dyn3d_dta  (jbdy)    &
456            &                , nn_tra_dta    (jbdy)    &
457#if defined key_lim3
458            &                , nn_ice_lim_dta(jbdy)    &
459#endif
460                              )
461         IF(nn_dta(jbdy) > 1)   nn_dta(jbdy) = 1
462      END DO
463
464      ! Work out upper bound of how many fields there are to read in and allocate arrays
465      ! ---------------------------------------------------------------------------
466      ALLOCATE( nb_bdy_fld(nb_bdy) )
467      nb_bdy_fld(:) = 0
468      DO jbdy = 1, nb_bdy         
469         IF( cn_dyn2d(jbdy) /= 'none' .AND. ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) ) THEN
470            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 3
471         ENDIF
472         IF( cn_dyn3d(jbdy) /= 'none' .AND. nn_dyn3d_dta(jbdy) == 1 ) THEN
473            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 2
474         ENDIF
475         IF( cn_tra(jbdy) /= 'none' .AND. nn_tra_dta(jbdy) == 1  ) THEN
476            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 2
477         ENDIF
478#if defined key_lim3
479         IF( cn_ice_lim(jbdy) /= 'none' .AND. nn_ice_lim_dta(jbdy) == 1  ) THEN
480            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 3
481         ENDIF
482#endif               
483         IF(lwp) WRITE(numout,*) 'Maximum number of files to open =', nb_bdy_fld(jbdy)
484      END DO           
485
486      nb_bdy_fld_sum = SUM( nb_bdy_fld )
487
488      ALLOCATE( bf(nb_bdy_fld_sum), STAT=ierror )
489      IF( ierror > 0 ) THEN   
490         CALL ctl_stop( 'bdy_dta: unable to allocate bf structure' )   ;   RETURN 
491      ENDIF
492      ALLOCATE( blf_i(nb_bdy_fld_sum), STAT=ierror )
493      IF( ierror > 0 ) THEN   
494         CALL ctl_stop( 'bdy_dta: unable to allocate blf_i structure' )   ;   RETURN 
495      ENDIF
496      ALLOCATE( nbmap_ptr(nb_bdy_fld_sum), STAT=ierror )
497      IF( ierror > 0 ) THEN   
498         CALL ctl_stop( 'bdy_dta: unable to allocate nbmap_ptr structure' )   ;   RETURN 
499      ENDIF
500      ALLOCATE( ilen1(nb_bdy_fld_sum), ilen3(nb_bdy_fld_sum) ) 
501      ALLOCATE( ibdy(nb_bdy_fld_sum) ) 
502      ALLOCATE( igrid(nb_bdy_fld_sum) ) 
503
504      ! Read namelists
505      ! --------------
506      REWIND(numnam_ref)
507      REWIND(numnam_cfg)
508      jfld = 0 
509      DO jbdy = 1, nb_bdy         
510         IF( nn_dta(jbdy) == 1 ) THEN
511            READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901)
512901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist', lwp )
513            READ  ( numnam_cfg, nambdy_dta, IOSTAT = ios, ERR = 902 )
514902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_dta in configuration namelist', lwp )
515            IF(lwm) WRITE( numond, nambdy_dta )
516
517            cn_dir_array(jbdy) = cn_dir
518            ln_full_vel_array(jbdy) = ln_full_vel
519
520            nblen => idx_bdy(jbdy)%nblen
521            nblenrim => idx_bdy(jbdy)%nblenrim
522            dta => dta_bdy(jbdy)
523            dta%nread(2) = 0
524
525            ! Only read in necessary fields for this set.
526            ! Important that barotropic variables come first.
527            IF( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) THEN
528
529               IF( dta%ll_ssh ) THEN
530                  if(lwp) write(numout,*) '++++++ reading in ssh field'
531                  jfld = jfld + 1
532                  blf_i(jfld) = bn_ssh
533                  ibdy(jfld) = jbdy
534                  igrid(jfld) = 1
535                  ilen1(jfld) = nblen(igrid(jfld))
536                  ilen3(jfld) = 1
537                  dta%nread(2) = dta%nread(2) + 1
538               ENDIF
539
540               IF( dta%ll_u2d .and. .not. ln_full_vel_array(jbdy) ) THEN
541                  if(lwp) write(numout,*) '++++++ reading in u2d field'
542                  jfld = jfld + 1
543                  blf_i(jfld) = bn_u2d
544                  ibdy(jfld) = jbdy
545                  igrid(jfld) = 2
546                  ilen1(jfld) = nblen(igrid(jfld))
547                  ilen3(jfld) = 1
548                  dta%nread(2) = dta%nread(2) + 1
549               ENDIF
550
551               IF( dta%ll_v2d .and. .not. ln_full_vel_array(jbdy) ) THEN
552                  if(lwp) write(numout,*) '++++++ reading in v2d field'
553                  jfld = jfld + 1
554                  blf_i(jfld) = bn_v2d
555                  ibdy(jfld) = jbdy
556                  igrid(jfld) = 3
557                  ilen1(jfld) = nblen(igrid(jfld))
558                  ilen3(jfld) = 1
559                  dta%nread(2) = dta%nread(2) + 1
560               ENDIF
561
562            ENDIF
563
564            ! read 3D velocities if baroclinic velocities require OR if
565            ! barotropic velocities required and ln_full_vel set to .true.
566            IF( nn_dyn3d_dta(jbdy) == 1 .OR. &
567           &  ( ln_full_vel_array(jbdy) .AND. ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 ) ) ) THEN
568
569               IF( dta%ll_u3d .OR. ( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) ) THEN
570                  if(lwp) write(numout,*) '++++++ reading in u3d field'
571                  jfld = jfld + 1
572                  blf_i(jfld) = bn_u3d
573                  ibdy(jfld) = jbdy
574                  igrid(jfld) = 2
575                  ilen1(jfld) = nblen(igrid(jfld))
576                  ilen3(jfld) = jpk
577                  IF( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) dta%nread(2) = dta%nread(2) + 1
578               ENDIF
579
580               IF( dta%ll_v3d .OR. ( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) ) THEN
581                  if(lwp) write(numout,*) '++++++ reading in v3d field'
582                  jfld = jfld + 1
583                  blf_i(jfld) = bn_v3d
584                  ibdy(jfld) = jbdy
585                  igrid(jfld) = 3
586                  ilen1(jfld) = nblen(igrid(jfld))
587                  ilen3(jfld) = jpk
588                  IF( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) dta%nread(2) = dta%nread(2) + 1
589               ENDIF
590
591            ENDIF
592
593            ! temperature and salinity
594            IF( nn_tra_dta(jbdy) == 1 ) THEN
595
596               IF( dta%ll_tem ) THEN
597                  if(lwp) write(numout,*) '++++++ reading in tem field'
598                  jfld = jfld + 1
599                  blf_i(jfld) = bn_tem
600                  ibdy(jfld) = jbdy
601                  igrid(jfld) = 1
602                  ilen1(jfld) = nblen(igrid(jfld))
603                  ilen3(jfld) = jpk
604               ENDIF
605
606               IF( dta%ll_sal ) THEN
607                  if(lwp) write(numout,*) '++++++ reading in sal field'
608                  jfld = jfld + 1
609                  blf_i(jfld) = bn_sal
610                  ibdy(jfld) = jbdy
611                  igrid(jfld) = 1
612                  ilen1(jfld) = nblen(igrid(jfld))
613                  ilen3(jfld) = jpk
614               ENDIF
615
616            ENDIF
617
618#if defined key_lim3
619            ! sea ice
620            IF( nn_ice_lim_dta(jbdy) == 1 ) THEN
621               ! Test for types of ice input (1cat or Xcat)
622               ! Build file name to find dimensions
623               clname=TRIM( cn_dir )//TRIM(bn_a_i%clname)
624               IF( .NOT. bn_a_i%ln_clim ) THEN   
625                                                  WRITE(clname, '(a,"_y",i4.4)' ) TRIM( clname ), nyear    ! add year
626                  IF( bn_a_i%cltype /= 'yearly' ) WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname ), nmonth   ! add month
627               ELSE
628                  IF( bn_a_i%cltype /= 'yearly' ) WRITE(clname, '(a,"_m",i2.2)' ) TRIM( clname ), nmonth   ! add month
629               ENDIF
630               IF( bn_a_i%cltype == 'daily' .OR. bn_a_i%cltype(1:4) == 'week' ) &
631               &                                  WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname ), nday     ! add day
632               !
633               CALL iom_open  ( clname, inum )
634               id1 = iom_varid( inum, bn_a_i%clvar, kdimsz=kdimsz, kndims=kndims, ldstop = .FALSE. )
635               CALL iom_close ( inum )
636
637                IF ( kndims == 4 ) THEN
638                 nice_cat = kdimsz(4)   ! Xcat input
639               ELSE
640                 nice_cat = 1           ! 1cat input     
641               ENDIF
642               ! End test
643
644               IF( dta%ll_a_i ) THEN
645                  jfld = jfld + 1
646                  blf_i(jfld) = bn_a_i
647                  ibdy(jfld)  = jbdy
648                  igrid(jfld) = 1
649                  ilen1(jfld) = nblen(igrid(jfld))
650                  ilen3(jfld) = nice_cat
651               ENDIF
652
653               IF( dta%ll_h_i ) THEN
654                  jfld = jfld + 1
655                  blf_i(jfld) = bn_h_i
656                  ibdy(jfld)  = jbdy
657                  igrid(jfld) = 1
658                  ilen1(jfld) = nblen(igrid(jfld))
659                  ilen3(jfld) = nice_cat
660               ENDIF
661
662               IF( dta%ll_h_s ) THEN
663                  jfld = jfld + 1
664                  blf_i(jfld) = bn_h_s
665                  ibdy(jfld)  = jbdy
666                  igrid(jfld) = 1
667                  ilen1(jfld) = nblen(igrid(jfld))
668                  ilen3(jfld) = nice_cat
669               ENDIF
670
671            ENDIF
672#endif
673            ! Recalculate field counts
674            !-------------------------
675            IF( jbdy == 1 ) THEN
676               nb_bdy_fld_sum = 0
677               nb_bdy_fld(jbdy) = jfld
678               nb_bdy_fld_sum     = jfld             
679            ELSE
680               nb_bdy_fld(jbdy) = jfld - nb_bdy_fld_sum
681               nb_bdy_fld_sum = nb_bdy_fld_sum + nb_bdy_fld(jbdy)
682            ENDIF
683
684            dta%nread(1) = nb_bdy_fld(jbdy)
685
686         ENDIF ! nn_dta == 1
687      ENDDO ! jbdy
688
689      DO jfld = 1, nb_bdy_fld_sum
690         ALLOCATE( bf(jfld)%fnow(ilen1(jfld),1,ilen3(jfld)) )
691         IF( blf_i(jfld)%ln_tint ) ALLOCATE( bf(jfld)%fdta(ilen1(jfld),1,ilen3(jfld),2) )
692         nbmap_ptr(jfld)%ptr => idx_bdy(ibdy(jfld))%nbmap(:,igrid(jfld))
693         nbmap_ptr(jfld)%ll_unstruc = ln_coords_file(ibdy(jfld))
694      ENDDO
695
696      ! fill bf with blf_i and control print
697      !-------------------------------------
698      jstart = 1
699      DO jbdy = 1, nb_bdy
700         jend = jstart - 1 + nb_bdy_fld(jbdy) 
701         CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(jbdy), 'bdy_dta',   &
702         &              'open boundary conditions', 'nambdy_dta' )
703         jstart = jend + 1
704      ENDDO
705
706      ! Initialise local boundary data arrays
707      ! nn_xxx_dta=0 : allocate space - will be filled from initial conditions later
708      ! nn_xxx_dta=1 : point to "fnow" arrays
709      !-------------------------------------
710
711      jfld = 0
712      DO jbdy=1, nb_bdy
713
714         nblen => idx_bdy(jbdy)%nblen
715         dta => dta_bdy(jbdy)
716
717         if(lwp) then
718            write(numout,*) '++++++ dta%ll_ssh = ',dta%ll_ssh
719            write(numout,*) '++++++ dta%ll_u2d = ',dta%ll_u2d
720            write(numout,*) '++++++ dta%ll_v2d = ',dta%ll_v2d
721            write(numout,*) '++++++ dta%ll_u3d = ',dta%ll_u3d
722            write(numout,*) '++++++ dta%ll_v3d = ',dta%ll_v3d
723            write(numout,*) '++++++ dta%ll_tem = ',dta%ll_tem
724            write(numout,*) '++++++ dta%ll_sal = ',dta%ll_sal
725         endif
726
727         IF ( nn_dyn2d_dta(jbdy) == 0 .or. nn_dyn2d_dta(jbdy) == 2 ) THEN
728            if(lwp) write(numout,*) '++++++ dta%ssh/u2d/u3d allocated space'
729            IF( dta%ll_ssh ) ALLOCATE( dta%ssh(nblen(1)) )
730            IF( dta%ll_u2d ) ALLOCATE( dta%u2d(nblen(2)) )
731            IF( dta%ll_v2d ) ALLOCATE( dta%v2d(nblen(3)) )
732         ENDIF
733         IF ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) THEN
734            IF( dta%ll_ssh ) THEN
735               if(lwp) write(numout,*) '++++++ dta%ssh pointing to fnow'
736               jfld = jfld + 1
737               dta%ssh => bf(jfld)%fnow(:,1,1)
738            ENDIF
739            IF ( dta%ll_u2d ) THEN
740               IF ( ln_full_vel_array(jbdy) ) THEN
741                  if(lwp) write(numout,*) '++++++ dta%u2d allocated space'
742                  ALLOCATE( dta%u2d(nblen(2)) )
743               ELSE
744                  if(lwp) write(numout,*) '++++++ dta%u2d pointing to fnow'
745                  jfld = jfld + 1
746                  dta%u2d => bf(jfld)%fnow(:,1,1)
747               ENDIF
748            ENDIF
749            IF ( dta%ll_v2d ) THEN
750               IF ( ln_full_vel_array(jbdy) ) THEN
751                  if(lwp) write(numout,*) '++++++ dta%v2d allocated space'
752                  ALLOCATE( dta%v2d(nblen(3)) )
753               ELSE
754                  if(lwp) write(numout,*) '++++++ dta%v2d pointing to fnow'
755                  jfld = jfld + 1
756                  dta%v2d => bf(jfld)%fnow(:,1,1)
757               ENDIF
758            ENDIF
759         ENDIF
760
761         IF ( nn_dyn3d_dta(jbdy) == 0 ) THEN
762            if(lwp) write(numout,*) '++++++ dta%u3d/v3d allocated space'
763            IF( dta%ll_u3d ) ALLOCATE( dta_bdy(jbdy)%u3d(nblen(2),jpk) )
764            IF( dta%ll_v3d ) ALLOCATE( dta_bdy(jbdy)%v3d(nblen(3),jpk) )
765         ENDIF
766         IF ( nn_dyn3d_dta(jbdy) == 1 .or. &
767           &  ( ln_full_vel_array(jbdy) .and. ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) ) ) THEN
768            IF ( dta%ll_u3d .or. ( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) ) THEN
769               if(lwp) write(numout,*) '++++++ dta%u3d pointing to fnow'
770               jfld = jfld + 1
771               dta_bdy(jbdy)%u3d => bf(jfld)%fnow(:,1,:)
772            ENDIF
773            IF ( dta%ll_v3d .or. ( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) ) THEN
774               if(lwp) write(numout,*) '++++++ dta%v3d pointing to fnow'
775               jfld = jfld + 1
776               dta_bdy(jbdy)%v3d => bf(jfld)%fnow(:,1,:)
777            ENDIF
778         ENDIF
779
780         IF( nn_tra_dta(jbdy) == 0 ) THEN
781            if(lwp) write(numout,*) '++++++ dta%tem/sal allocated space'
782            IF( dta%ll_tem ) ALLOCATE( dta_bdy(jbdy)%tem(nblen(1),jpk) )
783            IF( dta%ll_sal ) ALLOCATE( dta_bdy(jbdy)%sal(nblen(1),jpk) )
784         ELSE
785            IF( dta%ll_tem ) THEN
786               if(lwp) write(numout,*) '++++++ dta%tem pointing to fnow'
787               jfld = jfld + 1
788               dta_bdy(jbdy)%tem => bf(jfld)%fnow(:,1,:)
789            ENDIF
790            IF( dta%ll_sal ) THEN
791               if(lwp) write(numout,*) '++++++ dta%sal pointing to fnow'
792               jfld = jfld + 1
793               dta_bdy(jbdy)%sal => bf(jfld)%fnow(:,1,:)
794            ENDIF
795         ENDIF
796
797#if defined key_lim3
798         IF (cn_ice_lim(jbdy) /= 'none') THEN
799            IF( nn_ice_lim_dta(jbdy) == 0 ) THEN
800               ALLOCATE( dta_bdy(jbdy)%a_i(nblen(1),jpl) )
801               ALLOCATE( dta_bdy(jbdy)%h_i(nblen(1),jpl) )
802               ALLOCATE( dta_bdy(jbdy)%h_s(nblen(1),jpl) )
803            ELSE
804               IF ( nice_cat == jpl ) THEN ! case input cat = jpl
805                  jfld = jfld + 1
806                  dta_bdy(jbdy)%a_i => bf(jfld)%fnow(:,1,:)
807                  jfld = jfld + 1
808                  dta_bdy(jbdy)%h_i => bf(jfld)%fnow(:,1,:)
809                  jfld = jfld + 1
810                  dta_bdy(jbdy)%h_s => bf(jfld)%fnow(:,1,:)
811               ELSE                        ! case input cat = 1 OR (/=1 and /=jpl)
812                  jfld_ait(jbdy)  = jfld + 1
813                  jfld_htit(jbdy) = jfld + 2
814                  jfld_htst(jbdy) = jfld + 3
815                  jfld     = jfld + 3
816                  ALLOCATE( dta_bdy(jbdy)%a_i(nblen(1),jpl) )
817                  ALLOCATE( dta_bdy(jbdy)%h_i(nblen(1),jpl) )
818                  ALLOCATE( dta_bdy(jbdy)%h_s(nblen(1),jpl) )
819                  dta_bdy(jbdy)%a_i(:,:) = 0._wp
820                  dta_bdy(jbdy)%h_i(:,:) = 0._wp
821                  dta_bdy(jbdy)%h_s(:,:) = 0._wp
822               ENDIF
823
824            ENDIF
825         ENDIF
826#endif
827         !
828      END DO ! jbdy
829      !
830   END SUBROUTINE bdy_dta_init
831
832   !!==============================================================================
833END MODULE bdydta