source: NEMO/releases/release-4.0-HEAD/src/OCE/BDY/bdydta.F90 @ 12395

Last change on this file since 12395 was 12395, checked in by clem, 11 months ago

correct heat diffusion in the ice in case of Jules coupling (Met Office), and add a small improvement when reading ice temperature with bdy

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