source: NEMO/trunk/src/OCE/BDY/bdydta.F90 @ 9629

Last change on this file since 9629 was 9629, checked in by jamesharle, 3 years ago

Inserted missing code for bdy interp, #2088.

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