source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydta.F90 @ 11510

Last change on this file since 11510 was 11510, checked in by clem, 14 months ago

add ponds for bdy

  • Property svn:keywords set to Id
File size: 36.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 lib_mpp, ONLY: ctl_stop, ctl_nam
34   USE fldread        ! read input fields
35   USE iom            ! IOM library
36   USE in_out_manager ! I/O logical units
37   USE timing         ! Timing
38
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC   bdy_dta          ! routine called by step.F90 and dynspg_ts.F90
43   PUBLIC   bdy_dta_init     ! routine called by nemogcm.F90
44
45   INTEGER , PARAMETER ::   jpbdyfld  = 14    ! maximum number of files to read
46   INTEGER , PARAMETER ::   jp_bdyssh = 1     !
47   INTEGER , PARAMETER ::   jp_bdyu2d = 2     !
48   INTEGER , PARAMETER ::   jp_bdyv2d = 3     !
49   INTEGER , PARAMETER ::   jp_bdyu3d = 4     !
50   INTEGER , PARAMETER ::   jp_bdyv3d = 5     !
51   INTEGER , PARAMETER ::   jp_bdytem = 6     !
52   INTEGER , PARAMETER ::   jp_bdysal = 7     !
53   INTEGER , PARAMETER ::   jp_bdya_i = 8     !
54   INTEGER , PARAMETER ::   jp_bdyh_i = 9     !
55   INTEGER , PARAMETER ::   jp_bdyh_s = 10    !
56   INTEGER , PARAMETER ::   jp_bdyt_i = 11    !
57   INTEGER , PARAMETER ::   jp_bdyt_s = 12    !
58   INTEGER , PARAMETER ::   jp_bdytsu = 13    !
59   INTEGER , PARAMETER ::   jp_bdys_i = 14    !
60#if ! defined key_si3
61   INTEGER , PARAMETER ::   jpl = 1
62#endif
63
64!$AGRIF_DO_NOT_TREAT
65   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:,:), TARGET ::   bf   ! structure of input fields (file informations, fields read)
66!$AGRIF_END_DO_NOT_TREAT
67
68   !!----------------------------------------------------------------------
69   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
70   !! $Id$
71   !! Software governed by the CeCILL license (see ./LICENSE)
72   !!----------------------------------------------------------------------
73CONTAINS
74
75   SUBROUTINE bdy_dta( kt, kit, kt_offset )
76      !!----------------------------------------------------------------------
77      !!                   ***  SUBROUTINE bdy_dta  ***
78      !!                   
79      !! ** Purpose :   Update external data for open boundary conditions
80      !!
81      !! ** Method  :   Use fldread.F90
82      !!               
83      !!----------------------------------------------------------------------
84      INTEGER, INTENT(in)           ::   kt           ! ocean time-step index
85      INTEGER, INTENT(in), OPTIONAL ::   kit          ! subcycle time-step index (for timesplitting option)
86      INTEGER, INTENT(in), OPTIONAL ::   kt_offset    ! time offset in units of timesteps. NB. if kit
87      !                                               ! is present then units = subcycle timesteps.
88      !                                               ! kt_offset = 0 => get data at "now" time level
89      !                                               ! kt_offset = -1 => get data at "before" time level
90      !                                               ! kt_offset = +1 => get data at "after" time level
91      !                                               ! etc.
92      !
93      INTEGER ::  jbdy, jfld, jstart, jend, ib, jl    ! dummy loop indices
94      INTEGER ::  ii, ij, ik, igrd, ipl               ! local integers
95      INTEGER,   DIMENSION(jpbgrd)     ::   ilen1 
96      INTEGER,   DIMENSION(:), POINTER ::   nblen, nblenrim  ! short cuts
97      TYPE(OBC_DATA)         , POINTER ::   dta_alias        ! short cut
98      TYPE(FLD), DIMENSION(:), POINTER ::   bf_alias
99      !!---------------------------------------------------------------------------
100      !
101      IF( ln_timing )   CALL timing_start('bdy_dta')
102      !
103      ! Initialise data arrays once for all from initial conditions where required
104      !---------------------------------------------------------------------------
105      IF( kt == nit000 .AND. .NOT.PRESENT(kit) ) THEN
106
107         ! Calculate depth-mean currents
108         !-----------------------------
109
110         DO jbdy = 1, nb_bdy
111            !
112            nblen    => idx_bdy(jbdy)%nblen
113            nblenrim => idx_bdy(jbdy)%nblenrim
114            !
115            IF( nn_dyn2d_dta(jbdy) == 0 ) THEN
116               ilen1(:) = nblen(:)
117               IF( dta_bdy(jbdy)%lneed_ssh ) THEN
118                  igrd = 1
119                  DO ib = 1, ilen1(igrd)
120                     ii = idx_bdy(jbdy)%nbi(ib,igrd)
121                     ij = idx_bdy(jbdy)%nbj(ib,igrd)
122                     dta_bdy(jbdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)         
123                  END DO
124               ENDIF
125               IF( dta_bdy(jbdy)%lneed_dyn2d) THEN
126                  igrd = 2
127                  DO ib = 1, ilen1(igrd)
128                     ii = idx_bdy(jbdy)%nbi(ib,igrd)
129                     ij = idx_bdy(jbdy)%nbj(ib,igrd)
130                     dta_bdy(jbdy)%u2d(ib) = un_b(ii,ij) * umask(ii,ij,1)         
131                  END DO
132                  igrd = 3
133                  DO ib = 1, ilen1(igrd)
134                     ii = idx_bdy(jbdy)%nbi(ib,igrd)
135                     ij = idx_bdy(jbdy)%nbj(ib,igrd)
136                     dta_bdy(jbdy)%v2d(ib) = vn_b(ii,ij) * vmask(ii,ij,1)         
137                  END DO
138               ENDIF
139            ENDIF
140            !
141            IF( nn_dyn3d_dta(jbdy) == 0 ) THEN
142               ilen1(:) = nblen(:)
143               IF( dta_bdy(jbdy)%lneed_dyn3d ) THEN
144                  igrd = 2 
145                  DO ib = 1, ilen1(igrd)
146                     DO ik = 1, jpkm1
147                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
148                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
149                        dta_bdy(jbdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - un_b(ii,ij) ) * umask(ii,ij,ik)         
150                     END DO
151                  END DO
152                  igrd = 3 
153                  DO ib = 1, ilen1(igrd)
154                     DO ik = 1, jpkm1
155                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
156                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
157                        dta_bdy(jbdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - vn_b(ii,ij) ) * vmask(ii,ij,ik)         
158                     END DO
159                  END DO
160               ENDIF
161            ENDIF
162
163            IF( nn_tra_dta(jbdy) == 0 ) THEN
164               ilen1(:) = nblen(:)
165               IF( dta_bdy(jbdy)%lneed_tra ) THEN
166                  igrd = 1 
167                  DO ib = 1, ilen1(igrd)
168                     DO ik = 1, jpkm1
169                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
170                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
171                        dta_bdy(jbdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_bdytem) * tmask(ii,ij,ik)         
172                        dta_bdy(jbdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_bdysal) * tmask(ii,ij,ik)         
173                     END DO
174                  END DO
175               ENDIF
176            ENDIF
177
178#if defined key_si3
179            IF( nn_ice_dta(jbdy) == 0 ) THEN    ! set ice to initial values
180               ilen1(:) = nblen(:)
181               IF( dta_bdy(jbdy)%lneed_ice ) THEN
182                  igrd = 1   
183                  DO jl = 1, jpl
184                     DO ib = 1, ilen1(igrd)
185                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
186                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
187                        dta_bdy(jbdy)%a_i(ib,jl) =  a_i (ii,ij,jl) * tmask(ii,ij,1) 
188                        dta_bdy(jbdy)%h_i(ib,jl) =  h_i (ii,ij,jl) * tmask(ii,ij,1) 
189                        dta_bdy(jbdy)%h_s(ib,jl) =  h_s (ii,ij,jl) * tmask(ii,ij,1) 
190                        dta_bdy(jbdy)%t_i(ib,jl) =  SUM(t_i (ii,ij,:,jl)) * r1_nlay_i * tmask(ii,ij,1) 
191                        dta_bdy(jbdy)%t_s(ib,jl) =  SUM(t_s (ii,ij,:,jl)) * r1_nlay_s * tmask(ii,ij,1)
192                        dta_bdy(jbdy)%tsu(ib,jl) =  t_su(ii,ij,jl) * tmask(ii,ij,1) 
193                        dta_bdy(jbdy)%s_i(ib,jl) =  s_i (ii,ij,jl) * tmask(ii,ij,1) 
194                     END DO
195                  END DO
196               ENDIF
197            ENDIF
198#endif
199         END DO ! jbdy
200         !
201      ENDIF ! kt == nit000
202
203      ! update external data from files
204      !--------------------------------
205
206      DO jbdy = 1, nb_bdy
207
208         dta_alias => dta_bdy(jbdy)
209         bf_alias  => bf(:,jbdy)
210
211         ! read/update all bdy data
212         ! ------------------------
213         CALL fld_read( kt, 1, bf_alias, kit = kit, kt_offset = kt_offset )
214
215         ! apply some corrections in some specific cases...
216         ! --------------------------------------------------
217         !
218         ! if runoff condition: change river flow we read (in m3/s) into barotropic velocity (m/s)
219         IF( cn_tra(jbdy) == 'runoff' .AND. TRIM(bf_alias(jp_bdyu2d)%clrootname) /= 'NOT USED' ) THEN   ! runoff and we read u/v2d
220            !
221            igrd = 2                      ! zonal flow (m3/s) to barotropic zonal velocity (m/s)
222            DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
223               ii   = idx_bdy(jbdy)%nbi(ib,igrd)
224               ij   = idx_bdy(jbdy)%nbj(ib,igrd)
225               dta_alias%u2d(ib) = dta_alias%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
226            END DO
227            igrd = 3                      ! meridional flow (m3/s) to barotropic meridional velocity (m/s)
228            DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
229               ii   = idx_bdy(jbdy)%nbi(ib,igrd)
230               ij   = idx_bdy(jbdy)%nbj(ib,igrd)
231               dta_alias%v2d(ib) = dta_alias%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
232            END DO
233         ENDIF
234
235         ! tidal harmonic forcing ONLY: initialise arrays
236         IF( nn_dyn2d_dta(jbdy) == 2 ) THEN   ! we did not read ssh, u/v2d
237            IF( dta_alias%lneed_ssh   ) dta_alias%ssh(:) = 0._wp
238            IF( dta_alias%lneed_dyn2d ) dta_alias%u2d(:) = 0._wp
239            IF( dta_alias%lneed_dyn2d ) dta_alias%v2d(:) = 0._wp
240         ENDIF
241
242         ! If full velocities in boundary data, then split it into barotropic and baroclinic component
243         IF( bf_alias(jp_bdyu3d)%ltotvel ) THEN     ! if we read 3D total velocity (can be true only if u3d was read)
244            !
245            igrd = 2                       ! zonal velocity
246            dta_alias%u2d(:) = 0._wp       ! compute barotrope zonal velocity and put it in u2d
247            DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
248               ii   = idx_bdy(jbdy)%nbi(ib,igrd)
249               ij   = idx_bdy(jbdy)%nbj(ib,igrd)
250               DO ik = 1, jpkm1
251                  dta_alias%u2d(ib) = dta_alias%u2d(ib) + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik)
252               END DO
253               dta_alias%u2d(ib) =  dta_alias%u2d(ib) * r1_hu_n(ii,ij)
254               DO ik = 1, jpkm1            ! compute barocline zonal velocity and put it in u3d
255                  dta_alias%u3d(ib,ik) = dta_alias%u3d(ib,ik) - dta_alias%u2d(ib)
256               END DO
257            END DO
258            igrd = 3                       ! meridional velocity
259            dta_alias%v2d(:) = 0._wp       ! compute barotrope meridional velocity and put it in v2d
260            DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
261               ii   = idx_bdy(jbdy)%nbi(ib,igrd)
262               ij   = idx_bdy(jbdy)%nbj(ib,igrd)
263               DO ik = 1, jpkm1
264                  dta_alias%v2d(ib) = dta_alias%v2d(ib) + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik)
265               END DO
266               dta_alias%v2d(ib) =  dta_alias%v2d(ib) * r1_hv_n(ii,ij)
267               DO ik = 1, jpkm1            ! compute barocline meridional velocity and put it in v3d
268                  dta_alias%v3d(ib,ik) = dta_alias%v3d(ib,ik) - dta_alias%v2d(ib)
269               END DO
270            END DO
271         ENDIF   ! ltotvel
272
273         ! update tidal harmonic forcing
274         IF( PRESENT(kit) .AND. nn_dyn2d_dta(jbdy) .GE. 2 ) THEN
275            CALL bdytide_update( kt = kt, idx = idx_bdy(jbdy), dta = dta_alias, td = tides(jbdy),   & 
276               &                 kit = kit, kt_offset = kt_offset )
277         ENDIF
278
279         !  atm surface pressure : add inverted barometer effect to ssh if it was read
280         IF ( ln_apr_obc .AND. TRIM(bf_alias(jp_bdyssh)%clrootname) /= 'NOT USED' ) THEN
281            igrd = 1
282            DO ib = 1, idx_bdy(jbdy)%nblenrim(igrd)   ! ssh is used only on the rim
283               ii   = idx_bdy(jbdy)%nbi(ib,igrd)
284               ij   = idx_bdy(jbdy)%nbj(ib,igrd)
285               dta_alias%ssh(ib) = dta_alias%ssh(ib) + ssh_ib(ii,ij)
286            END DO
287         ENDIF
288
289#if defined key_si3
290         IF( dta_alias%lneed_ice ) THEN
291            ! fill temperature and salinity arrays
292            IF( TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyt_i)%fnow(:,1,:) = rice_tem(jbdy)
293            IF( TRIM(bf_alias(jp_bdyt_s)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyt_s)%fnow(:,1,:) = rice_tem(jbdy)
294            IF( TRIM(bf_alias(jp_bdytsu)%clrootname) == 'NOT USED' )   bf_alias(jp_bdytsu)%fnow(:,1,:) = rice_tem(jbdy)
295            IF( TRIM(bf_alias(jp_bdys_i)%clrootname) == 'NOT USED' )   bf_alias(jp_bdys_i)%fnow(:,1,:) = rice_sal(jbdy)
296            ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2
297            IF( TRIM(bf_alias(jp_bdytsu)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' ) &
298               &   bf_alias(jp_bdyt_i)%fnow(:,1,:) = 0.5_wp * ( bf_alias(jp_bdytsu)%fnow(:,1,:) + 271.15 )
299            ! if T_su is read and not T_s, set T_s = T_su
300            IF( TRIM(bf_alias(jp_bdytsu)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_s)%clrootname) == 'NOT USED' ) &
301               &   bf_alias(jp_bdyt_s)%fnow(:,1,:) = bf_alias(jp_bdytsu)%fnow(:,1,:)
302            ! if T_s is read and not T_su, set T_su = T_s
303            IF( TRIM(bf_alias(jp_bdyt_s)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdytsu)%clrootname) == 'NOT USED' ) &
304               &   bf_alias(jp_bdytsu)%fnow(:,1,:) = bf_alias(jp_bdyt_s)%fnow(:,1,:)
305            ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2
306            IF( TRIM(bf_alias(jp_bdyt_s)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' ) &
307               &   bf_alias(jp_bdyt_i)%fnow(:,1,:) = 0.5_wp * ( bf_alias(jp_bdyt_s)%fnow(:,1,:) + 271.15 )
308           
309            ! convert N-cat fields (input) into jpl-cat (output)
310            ipl = SIZE(bf_alias(jp_bdya_i)%fnow, 3)           
311            IF( ipl /= jpl ) THEN      ! ice: convert N-cat fields (input) into jpl-cat (output)
312               CALL ice_var_itd( bf_alias(jp_bdyh_i)%fnow(:,1,:), bf_alias(jp_bdyh_s)%fnow(:,1,:), bf_alias(jp_bdya_i)%fnow(:,1,:), &
313                  &              dta_alias%h_i                  , dta_alias%h_s                  , dta_alias%a_i                  , &
314                  &              bf_alias(jp_bdyt_i)%fnow(:,1,:), bf_alias(jp_bdyt_s)%fnow(:,1,:), &
315                  &              bf_alias(jp_bdytsu)%fnow(:,1,:), bf_alias(jp_bdys_i)%fnow(:,1,:), &
316                  &              dta_alias%t_i                  , dta_alias%t_s                  , &
317                  &              dta_alias%tsu                  , dta_alias%s_i )
318            ENDIF
319         ENDIF
320#endif
321      END DO  ! jbdy
322
323      IF ( ln_tide ) THEN
324         IF (ln_dynspg_ts) THEN      ! Fill temporary arrays with slow-varying bdy data                           
325            DO jbdy = 1, nb_bdy      ! Tidal component added in ts loop
326               IF ( nn_dyn2d_dta(jbdy) .GE. 2 ) THEN
327                  nblen => idx_bdy(jbdy)%nblen
328                  nblenrim => idx_bdy(jbdy)%nblenrim
329                  IF( cn_dyn2d(jbdy) == 'frs' ) THEN ; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF
330                     IF ( dta_bdy(jbdy)%lneed_ssh   ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1))
331                     IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2))
332                     IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3))
333                  ENDIF
334               END DO
335            ELSE ! Add tides if not split-explicit free surface else this is done in ts loop
336               !
337               CALL bdy_dta_tides( kt=kt, kt_offset=kt_offset )
338            ENDIF
339         ENDIF
340         !
341         IF( ln_timing )   CALL timing_stop('bdy_dta')
342         !
343      END SUBROUTINE bdy_dta
344
345
346   SUBROUTINE bdy_dta_init
347      !!----------------------------------------------------------------------
348      !!                   ***  SUBROUTINE bdy_dta_init  ***
349      !!                   
350      !! ** Purpose :   Initialise arrays for reading of external data
351      !!                for open boundary conditions
352      !!
353      !! ** Method  :   
354      !!               
355      !!----------------------------------------------------------------------
356      INTEGER ::   jbdy, jfld    ! Local integers
357      INTEGER ::   ierror, ios     !
358      !
359      CHARACTER(len=3)                       ::   cl3           !
360      CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files
361      LOGICAL                                ::   ln_full_vel   ! =T => full velocities in 3D boundary data
362      !                                                         ! =F => baroclinic velocities in 3D boundary data
363      LOGICAL                                ::   ln_zinterp    ! =T => requires a vertical interpolation of the bdydta
364      REAL(wp)                               ::   rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd 
365      INTEGER                                ::   ipk,ipl       !
366      INTEGER                                ::   idvar         ! variable ID
367      INTEGER                                ::   indims        ! number of dimensions of the variable
368      INTEGER                                ::   iszdim        ! number of dimensions of the variable
369      INTEGER, DIMENSION(4)                  ::   i4dimsz       ! size of variable dimensions
370      INTEGER                                ::   igrd          ! index for grid type (1,2,3 = T,U,V)
371      LOGICAL                                ::   lluld         ! is the variable using the unlimited dimension
372      LOGICAL                                ::   llneed        !
373      LOGICAL                                ::   llread        !
374      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! must be an array to be used with fld_fill
375      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read
376      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i       
377      TYPE(FLD_N), DIMENSION(:), POINTER ::   bn_alias                        ! must be an array to be used with fld_fill
378      TYPE(FLD  ), DIMENSION(:), POINTER ::   bf_alias
379      !
380      NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d 
381      NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i
382      NAMELIST/nambdy_dta/ rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd
383      NAMELIST/nambdy_dta/ ln_full_vel, ln_zinterp
384      !!---------------------------------------------------------------------------
385      !
386      IF(lwp) WRITE(numout,*)
387      IF(lwp) WRITE(numout,*) 'bdy_dta_ini : initialization of data at the open boundaries'
388      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
389      IF(lwp) WRITE(numout,*) ''
390
391      ALLOCATE( bf(jpbdyfld,nb_bdy), STAT=ierror )
392      IF( ierror > 0 ) THEN   
393         CALL ctl_stop( 'bdy_dta: unable to allocate bf structure' )   ;   RETURN 
394      ENDIF
395      bf(:,:)%clrootname = 'NOT USED'   ! default definition used as a flag in fld_read to do nothing.
396      bf(:,:)%lzint      = .FALSE.      ! default definition
397      bf(:,:)%ltotvel    = .FALSE.      ! default definition
398 
399      ! Read namelists
400      ! --------------
401      REWIND(numnam_cfg)
402      DO jbdy = 1, nb_bdy
403
404         WRITE(ctmp1, '(a,i2)') 'BDY number ', jbdy
405         WRITE(ctmp2, '(a,i2)') 'block nambdy_dta number ', jbdy
406
407         ! There is only one nambdy_dta block in namelist_ref -> use it for each bdy so we do a rewind
408         REWIND(numnam_ref)
409         READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901)
410901      IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist' )
411
412         !   by-pass nambdy_dta reading if no input data used in this bdy   
413         IF(       ( dta_bdy(jbdy)%lneed_dyn2d .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1 )   &
414            & .OR. ( dta_bdy(jbdy)%lneed_dyn3d .AND.     nn_dyn3d_dta(jbdy)    == 1 )   &
415            & .OR. ( dta_bdy(jbdy)%lneed_tra   .AND.       nn_tra_dta(jbdy)    == 1 )   &
416            & .OR. ( dta_bdy(jbdy)%lneed_ice   .AND.       nn_ice_dta(jbdy)    == 1 )   )   THEN
417            ! WARNING: we don't do a rewind here, each bdy reads its own nambdy_dta block one after another
418            READ  ( numnam_cfg, nambdy_dta, IOSTAT = ios, ERR = 902 )
419902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_dta in configuration namelist' )
420            IF(lwm) WRITE( numond, nambdy_dta )           
421         ENDIF
422
423         ! get the number of ice categories in bdy data file (use a_i information to do this)
424         ipl = jpl   ! default definition
425         IF( dta_bdy(jbdy)%lneed_ice ) THEN    ! if we need ice bdy data
426            IF( nn_ice_dta(jbdy) == 1 ) THEN   ! if we get ice bdy data from netcdf file
427               CALL fld_fill(  bf(jp_bdya_i,jbdy:jbdy), bn_a_i, cn_dir, 'bdy_dta', 'a_i'//' '//ctmp1, ctmp2 )   ! use namelist info
428               CALL fld_clopn( bf(jp_bdya_i,jbdy), nyear, nmonth, nday )   ! not a problem when we call it again after
429               idvar = iom_varid( bf(jp_bdya_i,jbdy)%num, bf(jp_bdya_i,jbdy)%clvar, kndims=indims, kdimsz=i4dimsz, lduld=lluld )
430               IF( indims == 4 .OR. ( indims == 3 .AND. .NOT. lluld ) ) THEN   ;   ipl = i4dimsz(3)   ! xylt or xyl
431               ELSE                                                            ;   ipl = 1            ! xy or xyt
432               ENDIF
433            ENDIF
434         ENDIF
435
436         ! temp, salt, age and ponds of incoming ice
437         rice_tem (jbdy) = rn_ice_tem
438         rice_sal (jbdy) = rn_ice_sal
439         rice_age (jbdy) = rn_ice_age
440         rice_apnd(jbdy) = rn_ice_apnd
441         rice_hpnd(jbdy) = rn_ice_hpnd
442         
443         
444         DO jfld = 1, jpbdyfld
445
446            ! =====================
447            !          ssh
448            ! =====================
449            IF( jfld == jp_bdyssh ) THEN
450               cl3 = 'ssh'
451               igrd = 1                                                    ! T point
452               ipk = 1                                                     ! surface data
453               llneed = dta_bdy(jbdy)%lneed_ssh                            ! dta_bdy(jbdy)%ssh will be needed
454               llread = MOD(nn_dyn2d_dta(jbdy),2) == 1                     ! get data from NetCDF file
455               bf_alias => bf(jp_bdyssh,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
456               bn_alias => bn_ssh                                          ! alias for ssh structure of nambdy_dta
457               iszdim = idx_bdy(jbdy)%nblenrim(igrd)                       ! length of this bdy on this MPI processus : used only on the rim
458            ENDIF
459            ! =====================
460            !         dyn2d
461            ! =====================
462            IF( jfld == jp_bdyu2d ) THEN
463               cl3 = 'u2d'
464               igrd = 2                                                    ! U point
465               ipk = 1                                                     ! surface data
466               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed
467               llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1   ! don't get u2d from u3d and read NetCDF file
468               bf_alias => bf(jp_bdyu2d,jbdy:jbdy)                         ! alias for u2d structure of bdy number jbdy
469               bn_alias => bn_u2d                                          ! alias for u2d structure of nambdy_dta
470               IF( ln_full_vel ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)      ! will be computed from u3d -> need on the full bdy
471               ELSE                    ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)   ! used only on the rim
472               ENDIF
473            ENDIF
474            IF( jfld == jp_bdyv2d ) THEN
475               cl3 = 'v2d'
476               igrd = 3                                                    ! V point
477               ipk = 1                                                     ! surface data
478               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed
479               llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1   ! don't get v2d from v3d and read NetCDF file
480               bf_alias => bf(jp_bdyv2d,jbdy:jbdy)                         ! alias for v2d structure of bdy number jbdy
481               bn_alias => bn_v2d                                          ! alias for v2d structure of nambdy_dta
482               IF( ln_full_vel ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)      ! will be computed from v3d -> need on the full bdy
483               ELSE                    ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)   ! used only on the rim
484               ENDIF
485            ENDIF
486            ! =====================
487            !         dyn3d
488            ! =====================
489            IF( jfld == jp_bdyu3d ) THEN
490               cl3 = 'u3d'
491               igrd = 2                                                    ! U point
492               ipk = jpk                                                   ! 3d data
493               llneed = dta_bdy(jbdy)%lneed_dyn3d .OR.               &     ! dta_bdy(jbdy)%u3d will be needed
494                  &   ( dta_bdy(jbdy)%lneed_dyn2d .AND. ln_full_vel )      !   u3d needed to compute u2d
495               llread = nn_dyn3d_dta(jbdy) == 1                            ! get data from NetCDF file
496               bf_alias => bf(jp_bdyu3d,jbdy:jbdy)                         ! alias for u3d structure of bdy number jbdy
497               bn_alias => bn_u3d                                          ! alias for u3d structure of nambdy_dta
498               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
499           ENDIF
500            IF( jfld == jp_bdyv3d ) THEN
501               cl3 = 'v3d'
502               igrd = 3                                                    ! V point
503               ipk = jpk                                                   ! 3d data
504               llneed = dta_bdy(jbdy)%lneed_dyn3d .OR.               &     ! dta_bdy(jbdy)%v3d will be needed
505                  &   ( dta_bdy(jbdy)%lneed_dyn2d .AND. ln_full_vel )      !   v3d needed to compute v2d
506               llread = nn_dyn3d_dta(jbdy) == 1                            ! get data from NetCDF file
507               bf_alias => bf(jp_bdyv3d,jbdy:jbdy)                         ! alias for v3d structure of bdy number jbdy
508               bn_alias => bn_v3d                                          ! alias for v3d structure of nambdy_dta
509               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
510           ENDIF
511
512            ! =====================
513            !          tra
514            ! =====================
515            IF( jfld == jp_bdytem ) THEN
516               cl3 = 'tem'
517               igrd = 1                                                    ! T point
518               ipk = jpk                                                   ! 3d data
519               llneed = dta_bdy(jbdy)%lneed_tra                            ! dta_bdy(jbdy)%tem will be needed
520               llread = nn_tra_dta(jbdy) == 1                              ! get data from NetCDF file
521               bf_alias => bf(jp_bdytem,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
522               bn_alias => bn_tem                                          ! alias for ssh structure of nambdy_dta
523               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
524            ENDIF
525            IF( jfld == jp_bdysal ) THEN
526               cl3 = 'sal'
527               igrd = 1                                                    ! T point
528               ipk = jpk                                                   ! 3d data
529               llneed = dta_bdy(jbdy)%lneed_tra                            ! dta_bdy(jbdy)%sal will be needed
530               llread = nn_tra_dta(jbdy) == 1                              ! get data from NetCDF file
531               bf_alias => bf(jp_bdysal,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
532               bn_alias => bn_sal                                          ! alias for ssh structure of nambdy_dta
533               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
534            ENDIF
535
536            ! =====================
537            !          ice
538            ! =====================
539            IF(  jfld == jp_bdya_i .OR. jfld == jp_bdyh_i .OR. jfld == jp_bdyh_s .OR. &
540               & jfld == jp_bdyt_i .OR. jfld == jp_bdyt_s .OR. jfld == jp_bdytsu .OR. jfld == jp_bdys_i ) THEN
541               igrd = 1                                                    ! T point
542               ipk = ipl                                                   ! jpl-cat data
543               llneed = dta_bdy(jbdy)%lneed_ice                            ! ice will be needed
544               llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file
545               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
546            ENDIF
547            IF( jfld == jp_bdya_i ) THEN
548               cl3 = 'a_i'
549               bf_alias => bf(jp_bdya_i,jbdy:jbdy)                         ! alias for a_i structure of bdy number jbdy
550               bn_alias => bn_a_i                                          ! alias for a_i structure of nambdy_dta
551            ENDIF
552            IF( jfld == jp_bdyh_i ) THEN
553               cl3 = 'h_i'
554               bf_alias => bf(jp_bdyh_i,jbdy:jbdy)                         ! alias for h_i structure of bdy number jbdy
555               bn_alias => bn_h_i                                          ! alias for h_i structure of nambdy_dta
556            ENDIF
557            IF( jfld == jp_bdyh_s ) THEN
558               cl3 = 'h_s'
559               bf_alias => bf(jp_bdyh_s,jbdy:jbdy)                         ! alias for h_s structure of bdy number jbdy
560               bn_alias => bn_h_s                                          ! alias for h_s structure of nambdy_dta
561            ENDIF
562            IF( jfld == jp_bdyt_i ) THEN
563               cl3 = 't_i'
564               bf_alias => bf(jp_bdyt_i,jbdy:jbdy)                         ! alias for t_i structure of bdy number jbdy
565               bn_alias => bn_t_i                                          ! alias for t_i structure of nambdy_dta
566            ENDIF
567            IF( jfld == jp_bdyt_s ) THEN
568               cl3 = 't_s'
569               bf_alias => bf(jp_bdyt_s,jbdy:jbdy)                         ! alias for t_s structure of bdy number jbdy
570               bn_alias => bn_t_s                                          ! alias for t_s structure of nambdy_dta
571            ENDIF
572            IF( jfld == jp_bdytsu ) THEN
573               cl3 = 'tsu'
574               bf_alias => bf(jp_bdytsu,jbdy:jbdy)                         ! alias for tsu structure of bdy number jbdy
575               bn_alias => bn_tsu                                          ! alias for tsu structure of nambdy_dta
576            ENDIF
577            IF( jfld == jp_bdys_i ) THEN
578               cl3 = 's_i'
579               bf_alias => bf(jp_bdys_i,jbdy:jbdy)                         ! alias for s_i structure of bdy number jbdy
580               bn_alias => bn_s_i                                          ! alias for s_i structure of nambdy_dta
581            ENDIF
582
583            IF( llneed ) THEN                                              ! dta_bdy(jbdy)%xxx will be needed
584               !                                                           !   -> must be associated with an allocated target
585               ALLOCATE( bf_alias(1)%fnow( iszdim, 1, ipk ) )              ! allocate the target
586               !
587               IF( llread ) THEN                                           ! get data from NetCDF file
588                  CALL fld_fill( bf_alias, bn_alias, cn_dir, 'bdy_dta', cl3//' '//ctmp1, ctmp2 )   ! use namelist info
589                  IF( bf_alias(1)%ln_tint ) ALLOCATE( bf_alias(1)%fdta( iszdim, 1, ipk, 2 ) )
590                  bf_alias(1)%imap    => idx_bdy(jbdy)%nbmap(1:iszdim,igrd)   ! associate the mapping used for this bdy
591                  bf_alias(1)%igrd    = igrd                                  ! used only for vertical integration of 3D arrays
592                  bf_alias(1)%ltotvel = ln_full_vel                           ! T if u3d is full velocity
593                  bf_alias(1)%lzint   = ln_zinterp                            ! T if it requires a vertical interpolation
594               ENDIF
595
596               ! associate the pointer and get rid of the dimensions with a size equal to 1
597               IF( jfld == jp_bdyssh )        dta_bdy(jbdy)%ssh => bf_alias(1)%fnow(:,1,1)
598               IF( jfld == jp_bdyu2d )        dta_bdy(jbdy)%u2d => bf_alias(1)%fnow(:,1,1)
599               IF( jfld == jp_bdyv2d )        dta_bdy(jbdy)%v2d => bf_alias(1)%fnow(:,1,1)
600               IF( jfld == jp_bdyu3d )        dta_bdy(jbdy)%u3d => bf_alias(1)%fnow(:,1,:)
601               IF( jfld == jp_bdyv3d )        dta_bdy(jbdy)%v3d => bf_alias(1)%fnow(:,1,:)
602               IF( jfld == jp_bdytem )        dta_bdy(jbdy)%tem => bf_alias(1)%fnow(:,1,:)
603               IF( jfld == jp_bdysal )        dta_bdy(jbdy)%sal => bf_alias(1)%fnow(:,1,:)
604               IF( jfld == jp_bdya_i ) THEN
605                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%a_i => bf_alias(1)%fnow(:,1,:)
606                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%a_i(iszdim,jpl) )
607                  ENDIF
608               ENDIF
609               IF( jfld == jp_bdyh_i ) THEN
610                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%h_i => bf_alias(1)%fnow(:,1,:)
611                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%h_i(iszdim,jpl) )
612                  ENDIF
613               ENDIF
614               IF( jfld == jp_bdyh_s ) THEN
615                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%h_s => bf_alias(1)%fnow(:,1,:)
616                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%h_s(iszdim,jpl) )
617                  ENDIF
618               ENDIF
619               IF( jfld == jp_bdyt_i ) THEN
620                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%t_i => bf_alias(1)%fnow(:,1,:)
621                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%t_i(iszdim,jpl) )
622                  ENDIF
623               ENDIF
624               IF( jfld == jp_bdyt_s ) THEN
625                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%t_s => bf_alias(1)%fnow(:,1,:)
626                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%t_s(iszdim,jpl) )
627                  ENDIF
628               ENDIF
629               IF( jfld == jp_bdytsu ) THEN
630                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%tsu => bf_alias(1)%fnow(:,1,:)
631                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%tsu(iszdim,jpl) )
632                  ENDIF
633               ENDIF
634               IF( jfld == jp_bdys_i ) THEN
635                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%s_i => bf_alias(1)%fnow(:,1,:)
636                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%s_i(iszdim,jpl) )
637                  ENDIF
638               ENDIF
639            ENDIF
640
641         END DO   ! jpbdyfld
642         !
643      END DO ! jbdy
644      !
645   END SUBROUTINE bdy_dta_init
646   
647   !!==============================================================================
648END MODULE bdydta
Note: See TracBrowser for help on using the repository browser.