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.
bdydta.F90 in NEMO/trunk/src/OCE/BDY – NEMO

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

Last change on this file since 10068 was 10068, checked in by nicolasmartin, 6 years ago

First part of modifications to have a common default header : fix typos and SVN keywords properties

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