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/branches/UKMO/NEMO_4.0.1_GO8_package/src/OCE/BDY – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.1_GO8_package/src/OCE/BDY/bdydta.F90

Last change on this file was 12088, checked in by deazer, 5 years ago

Updating GO8 Package branch to bring in required BDY bug fixes frouse with CO8
The mirror branch is already updated to have this change, where we merge in the mirror to the package branch

File size: 39.6 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            ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2
305            IF( TRIM(bf_alias(jp_bdytsu)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' ) &
306               &   bf_alias(jp_bdyt_i)%fnow(:,1,:) = 0.5_wp * ( bf_alias(jp_bdytsu)%fnow(:,1,:) + 271.15 )
307            ! if T_su is read and not T_s, set T_s = T_su
308            IF( TRIM(bf_alias(jp_bdytsu)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_s)%clrootname) == 'NOT USED' ) &
309               &   bf_alias(jp_bdyt_s)%fnow(:,1,:) = bf_alias(jp_bdytsu)%fnow(:,1,:)
310            ! if T_s is read and not T_su, set T_su = T_s
311            IF( TRIM(bf_alias(jp_bdyt_s)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdytsu)%clrootname) == 'NOT USED' ) &
312               &   bf_alias(jp_bdytsu)%fnow(:,1,:) = bf_alias(jp_bdyt_s)%fnow(:,1,:)
313            ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2
314            IF( TRIM(bf_alias(jp_bdyt_s)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' ) &
315               &   bf_alias(jp_bdyt_i)%fnow(:,1,:) = 0.5_wp * ( bf_alias(jp_bdyt_s)%fnow(:,1,:) + 271.15 )
316
317            ! make sure ponds = 0 if no ponds scheme
318            IF ( .NOT.ln_pnd ) THEN
319               bf_alias(jp_bdyaip)%fnow(:,1,:) = 0._wp
320               bf_alias(jp_bdyhip)%fnow(:,1,:) = 0._wp
321            ENDIF
322           
323            ! convert N-cat fields (input) into jpl-cat (output)
324            ipl = SIZE(bf_alias(jp_bdya_i)%fnow, 3)           
325            IF( ipl /= jpl ) THEN      ! ice: convert N-cat fields (input) into jpl-cat (output)
326               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,:), &
327                  &              dta_alias%h_i                  , dta_alias%h_s                  , dta_alias%a_i                  , &
328                  &              bf_alias(jp_bdyt_i)%fnow(:,1,:), bf_alias(jp_bdyt_s)%fnow(:,1,:), &
329                  &              bf_alias(jp_bdytsu)%fnow(:,1,:), bf_alias(jp_bdys_i)%fnow(:,1,:), &
330                  &              bf_alias(jp_bdyaip)%fnow(:,1,:), bf_alias(jp_bdyhip)%fnow(:,1,:), &
331                  &              dta_alias%t_i                  , dta_alias%t_s                  , &
332                  &              dta_alias%tsu                  , dta_alias%s_i                  , &
333                  &              dta_alias%aip                  , dta_alias%hip )
334            ENDIF
335         ENDIF
336#endif
337      END DO  ! jbdy
338
339      IF ( ln_tide ) THEN
340         IF (ln_dynspg_ts) THEN      ! Fill temporary arrays with slow-varying bdy data                           
341            DO jbdy = 1, nb_bdy      ! Tidal component added in ts loop
342               IF ( nn_dyn2d_dta(jbdy) .GE. 2 ) THEN
343                  nblen => idx_bdy(jbdy)%nblen
344                  nblenrim => idx_bdy(jbdy)%nblenrim
345                  IF( cn_dyn2d(jbdy) == 'frs' ) THEN ; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF
346                     IF ( dta_bdy(jbdy)%lneed_ssh   ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1))
347                     IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2))
348                     IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3))
349                  ENDIF
350               END DO
351            ELSE ! Add tides if not split-explicit free surface else this is done in ts loop
352               !
353               CALL bdy_dta_tides( kt=kt, kt_offset=kt_offset )
354            ENDIF
355         ENDIF
356         !
357         IF( ln_timing )   CALL timing_stop('bdy_dta')
358         !
359      END SUBROUTINE bdy_dta
360
361
362   SUBROUTINE bdy_dta_init
363      !!----------------------------------------------------------------------
364      !!                   ***  SUBROUTINE bdy_dta_init  ***
365      !!                   
366      !! ** Purpose :   Initialise arrays for reading of external data
367      !!                for open boundary conditions
368      !!
369      !! ** Method  :   
370      !!               
371      !!----------------------------------------------------------------------
372      INTEGER ::   jbdy, jfld    ! Local integers
373      INTEGER ::   ierror, ios     !
374      !
375      CHARACTER(len=3)                       ::   cl3           !
376      CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files
377      LOGICAL                                ::   ln_full_vel   ! =T => full velocities in 3D boundary data
378      !                                                         ! =F => baroclinic velocities in 3D boundary data
379      LOGICAL                                ::   ln_zinterp    ! =T => requires a vertical interpolation of the bdydta
380      REAL(wp)                               ::   rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd 
381      INTEGER                                ::   ipk,ipl       !
382      INTEGER                                ::   idvar         ! variable ID
383      INTEGER                                ::   indims        ! number of dimensions of the variable
384      INTEGER                                ::   iszdim        ! number of dimensions of the variable
385      INTEGER, DIMENSION(4)                  ::   i4dimsz       ! size of variable dimensions
386      INTEGER                                ::   igrd          ! index for grid type (1,2,3 = T,U,V)
387      LOGICAL                                ::   lluld         ! is the variable using the unlimited dimension
388      LOGICAL                                ::   llneed        !
389      LOGICAL                                ::   llread        !
390      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! must be an array to be used with fld_fill
391      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read
392      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       
393      TYPE(FLD_N), DIMENSION(:), POINTER ::   bn_alias                        ! must be an array to be used with fld_fill
394      TYPE(FLD  ), DIMENSION(:), POINTER ::   bf_alias
395      !
396      NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d 
397      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
398      NAMELIST/nambdy_dta/ rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd
399      NAMELIST/nambdy_dta/ ln_full_vel, ln_zinterp
400      !!---------------------------------------------------------------------------
401      !
402      IF(lwp) WRITE(numout,*)
403      IF(lwp) WRITE(numout,*) 'bdy_dta_ini : initialization of data at the open boundaries'
404      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
405      IF(lwp) WRITE(numout,*) ''
406
407      ALLOCATE( bf(jpbdyfld,nb_bdy), STAT=ierror )
408      IF( ierror > 0 ) THEN   
409         CALL ctl_stop( 'bdy_dta: unable to allocate bf structure' )   ;   RETURN 
410      ENDIF
411      bf(:,:)%clrootname = 'NOT USED'   ! default definition used as a flag in fld_read to do nothing.
412      bf(:,:)%lzint      = .FALSE.      ! default definition
413      bf(:,:)%ltotvel    = .FALSE.      ! default definition
414 
415      ! Read namelists
416      ! --------------
417      REWIND(numnam_cfg)
418      DO jbdy = 1, nb_bdy
419
420         WRITE(ctmp1, '(a,i2)') 'BDY number ', jbdy
421         WRITE(ctmp2, '(a,i2)') 'block nambdy_dta number ', jbdy
422
423         ! There is only one nambdy_dta block in namelist_ref -> use it for each bdy so we do a rewind
424         REWIND(numnam_ref)
425         READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901)
426901      IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist' )
427
428         !   by-pass nambdy_dta reading if no input data used in this bdy   
429         IF(       ( dta_bdy(jbdy)%lneed_dyn2d .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1 )   &
430            & .OR. ( dta_bdy(jbdy)%lneed_dyn3d .AND.     nn_dyn3d_dta(jbdy)    == 1 )   &
431            & .OR. ( dta_bdy(jbdy)%lneed_tra   .AND.       nn_tra_dta(jbdy)    == 1 )   &
432            & .OR. ( dta_bdy(jbdy)%lneed_ice   .AND.       nn_ice_dta(jbdy)    == 1 )   )   THEN
433            ! WARNING: we don't do a rewind here, each bdy reads its own nambdy_dta block one after another
434            READ  ( numnam_cfg, nambdy_dta, IOSTAT = ios, ERR = 902 )
435902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_dta in configuration namelist' )
436            IF(lwm) WRITE( numond, nambdy_dta )           
437         ENDIF
438
439         ! get the number of ice categories in bdy data file (use a_i information to do this)
440         ipl = jpl   ! default definition
441         IF( dta_bdy(jbdy)%lneed_ice ) THEN    ! if we need ice bdy data
442            IF( nn_ice_dta(jbdy) == 1 ) THEN   ! if we get ice bdy data from netcdf file
443               CALL fld_fill(  bf(jp_bdya_i,jbdy:jbdy), bn_a_i, cn_dir, 'bdy_dta', 'a_i'//' '//ctmp1, ctmp2 )   ! use namelist info
444               CALL fld_clopn( bf(jp_bdya_i,jbdy), nyear, nmonth, nday )   ! not a problem when we call it again after
445               idvar = iom_varid( bf(jp_bdya_i,jbdy)%num, bf(jp_bdya_i,jbdy)%clvar, kndims=indims, kdimsz=i4dimsz, lduld=lluld )
446               IF( indims == 4 .OR. ( indims == 3 .AND. .NOT. lluld ) ) THEN   ;   ipl = i4dimsz(3)   ! xylt or xyl
447               ELSE                                                            ;   ipl = 1            ! xy or xyt
448               ENDIF
449               bf(jp_bdya_i,jbdy)%clrootname = 'NOT USED'   ! reset to default value as this subdomain may not need to read this bdy
450            ENDIF
451         ENDIF
452
453#if defined key_si3
454         IF( .NOT.ln_pnd ) THEN
455            rn_ice_apnd = 0. ; rn_ice_hpnd = 0.
456            CALL ctl_warn( 'rn_ice_apnd & rn_ice_hpnd = 0 when no ponds' )
457         ENDIF
458#endif
459
460         ! temp, salt, age and ponds of incoming ice
461         rice_tem (jbdy) = rn_ice_tem
462         rice_sal (jbdy) = rn_ice_sal
463         rice_age (jbdy) = rn_ice_age
464         rice_apnd(jbdy) = rn_ice_apnd
465         rice_hpnd(jbdy) = rn_ice_hpnd
466         
467         
468         DO jfld = 1, jpbdyfld
469
470            ! =====================
471            !          ssh
472            ! =====================
473            IF( jfld == jp_bdyssh ) THEN
474               cl3 = 'ssh'
475               igrd = 1                                                    ! T point
476               ipk = 1                                                     ! surface data
477               llneed = dta_bdy(jbdy)%lneed_ssh                            ! dta_bdy(jbdy)%ssh will be needed
478               llread = MOD(nn_dyn2d_dta(jbdy),2) == 1                     ! get data from NetCDF file
479               bf_alias => bf(jp_bdyssh,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
480               bn_alias => bn_ssh                                          ! alias for ssh structure of nambdy_dta
481               iszdim = idx_bdy(jbdy)%nblenrim(igrd)                       ! length of this bdy on this MPI processus : used only on the rim
482            ENDIF
483            ! =====================
484            !         dyn2d
485            ! =====================
486            IF( jfld == jp_bdyu2d ) THEN
487               cl3 = 'u2d'
488               igrd = 2                                                    ! U point
489               ipk = 1                                                     ! surface data
490               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed
491               llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1   ! don't get u2d from u3d and read NetCDF file
492               bf_alias => bf(jp_bdyu2d,jbdy:jbdy)                         ! alias for u2d structure of bdy number jbdy
493               bn_alias => bn_u2d                                          ! alias for u2d structure of nambdy_dta
494               IF( ln_full_vel ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)      ! will be computed from u3d -> need on the full bdy
495               ELSE                    ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)   ! used only on the rim
496               ENDIF
497            ENDIF
498            IF( jfld == jp_bdyv2d ) THEN
499               cl3 = 'v2d'
500               igrd = 3                                                    ! V point
501               ipk = 1                                                     ! surface data
502               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed
503               llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1   ! don't get v2d from v3d and read NetCDF file
504               bf_alias => bf(jp_bdyv2d,jbdy:jbdy)                         ! alias for v2d structure of bdy number jbdy
505               bn_alias => bn_v2d                                          ! alias for v2d structure of nambdy_dta
506               IF( ln_full_vel ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)      ! will be computed from v3d -> need on the full bdy
507               ELSE                    ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)   ! used only on the rim
508               ENDIF
509            ENDIF
510            ! =====================
511            !         dyn3d
512            ! =====================
513            IF( jfld == jp_bdyu3d ) THEN
514               cl3 = 'u3d'
515               igrd = 2                                                    ! U point
516               ipk = jpk                                                   ! 3d data
517               llneed = dta_bdy(jbdy)%lneed_dyn3d .OR.               &     ! dta_bdy(jbdy)%u3d will be needed
518                  &   ( dta_bdy(jbdy)%lneed_dyn2d .AND. ln_full_vel )      !   u3d needed to compute u2d
519               llread = nn_dyn3d_dta(jbdy) == 1                            ! get data from NetCDF file
520               bf_alias => bf(jp_bdyu3d,jbdy:jbdy)                         ! alias for u3d structure of bdy number jbdy
521               bn_alias => bn_u3d                                          ! alias for u3d structure of nambdy_dta
522               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
523           ENDIF
524            IF( jfld == jp_bdyv3d ) THEN
525               cl3 = 'v3d'
526               igrd = 3                                                    ! V point
527               ipk = jpk                                                   ! 3d data
528               llneed = dta_bdy(jbdy)%lneed_dyn3d .OR.               &     ! dta_bdy(jbdy)%v3d will be needed
529                  &   ( dta_bdy(jbdy)%lneed_dyn2d .AND. ln_full_vel )      !   v3d needed to compute v2d
530               llread = nn_dyn3d_dta(jbdy) == 1                            ! get data from NetCDF file
531               bf_alias => bf(jp_bdyv3d,jbdy:jbdy)                         ! alias for v3d structure of bdy number jbdy
532               bn_alias => bn_v3d                                          ! alias for v3d structure of nambdy_dta
533               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
534           ENDIF
535
536            ! =====================
537            !          tra
538            ! =====================
539            IF( jfld == jp_bdytem ) THEN
540               cl3 = 'tem'
541               igrd = 1                                                    ! T point
542               ipk = jpk                                                   ! 3d data
543               llneed = dta_bdy(jbdy)%lneed_tra                            ! dta_bdy(jbdy)%tem will be needed
544               llread = nn_tra_dta(jbdy) == 1                              ! get data from NetCDF file
545               bf_alias => bf(jp_bdytem,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
546               bn_alias => bn_tem                                          ! alias for ssh structure of nambdy_dta
547               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
548            ENDIF
549            IF( jfld == jp_bdysal ) THEN
550               cl3 = 'sal'
551               igrd = 1                                                    ! T point
552               ipk = jpk                                                   ! 3d data
553               llneed = dta_bdy(jbdy)%lneed_tra                            ! dta_bdy(jbdy)%sal will be needed
554               llread = nn_tra_dta(jbdy) == 1                              ! get data from NetCDF file
555               bf_alias => bf(jp_bdysal,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
556               bn_alias => bn_sal                                          ! alias for ssh structure of nambdy_dta
557               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
558            ENDIF
559
560            ! =====================
561            !          ice
562            ! =====================
563            IF(  jfld == jp_bdya_i .OR. jfld == jp_bdyh_i .OR. jfld == jp_bdyh_s .OR. &
564               & jfld == jp_bdyt_i .OR. jfld == jp_bdyt_s .OR. jfld == jp_bdytsu .OR. &
565               & jfld == jp_bdys_i .OR. jfld == jp_bdyaip .OR. jfld == jp_bdyhip      ) THEN
566               igrd = 1                                                    ! T point
567               ipk = ipl                                                   ! jpl-cat data
568               llneed = dta_bdy(jbdy)%lneed_ice                            ! ice will be needed
569               llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file
570               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
571            ENDIF
572            IF( jfld == jp_bdya_i ) THEN
573               cl3 = 'a_i'
574               bf_alias => bf(jp_bdya_i,jbdy:jbdy)                         ! alias for a_i structure of bdy number jbdy
575               bn_alias => bn_a_i                                          ! alias for a_i structure of nambdy_dta
576            ENDIF
577            IF( jfld == jp_bdyh_i ) THEN
578               cl3 = 'h_i'
579               bf_alias => bf(jp_bdyh_i,jbdy:jbdy)                         ! alias for h_i structure of bdy number jbdy
580               bn_alias => bn_h_i                                          ! alias for h_i structure of nambdy_dta
581            ENDIF
582            IF( jfld == jp_bdyh_s ) THEN
583               cl3 = 'h_s'
584               bf_alias => bf(jp_bdyh_s,jbdy:jbdy)                         ! alias for h_s structure of bdy number jbdy
585               bn_alias => bn_h_s                                          ! alias for h_s structure of nambdy_dta
586            ENDIF
587            IF( jfld == jp_bdyt_i ) THEN
588               cl3 = 't_i'
589               bf_alias => bf(jp_bdyt_i,jbdy:jbdy)                         ! alias for t_i structure of bdy number jbdy
590               bn_alias => bn_t_i                                          ! alias for t_i structure of nambdy_dta
591            ENDIF
592            IF( jfld == jp_bdyt_s ) THEN
593               cl3 = 't_s'
594               bf_alias => bf(jp_bdyt_s,jbdy:jbdy)                         ! alias for t_s structure of bdy number jbdy
595               bn_alias => bn_t_s                                          ! alias for t_s structure of nambdy_dta
596            ENDIF
597            IF( jfld == jp_bdytsu ) THEN
598               cl3 = 'tsu'
599               bf_alias => bf(jp_bdytsu,jbdy:jbdy)                         ! alias for tsu structure of bdy number jbdy
600               bn_alias => bn_tsu                                          ! alias for tsu structure of nambdy_dta
601            ENDIF
602            IF( jfld == jp_bdys_i ) THEN
603               cl3 = 's_i'
604               bf_alias => bf(jp_bdys_i,jbdy:jbdy)                         ! alias for s_i structure of bdy number jbdy
605               bn_alias => bn_s_i                                          ! alias for s_i structure of nambdy_dta
606            ENDIF
607            IF( jfld == jp_bdyaip ) THEN
608               cl3 = 'aip'
609               bf_alias => bf(jp_bdyaip,jbdy:jbdy)                         ! alias for aip structure of bdy number jbdy
610               bn_alias => bn_aip                                          ! alias for aip structure of nambdy_dta
611            ENDIF
612            IF( jfld == jp_bdyhip ) THEN
613               cl3 = 'hip'
614               bf_alias => bf(jp_bdyhip,jbdy:jbdy)                         ! alias for hip structure of bdy number jbdy
615               bn_alias => bn_hip                                          ! alias for hip structure of nambdy_dta
616            ENDIF
617
618            IF( llneed .AND. iszdim > 0 ) THEN                             ! dta_bdy(jbdy)%xxx will be needed
619               !                                                           !   -> must be associated with an allocated target
620               ALLOCATE( bf_alias(1)%fnow( iszdim, 1, ipk ) )              ! allocate the target
621               !
622               IF( llread ) THEN                                           ! get data from NetCDF file
623                  CALL fld_fill( bf_alias, bn_alias, cn_dir, 'bdy_dta', cl3//' '//ctmp1, ctmp2 )   ! use namelist info
624                  IF( bf_alias(1)%ln_tint ) ALLOCATE( bf_alias(1)%fdta( iszdim, 1, ipk, 2 ) )
625                  bf_alias(1)%imap    => idx_bdy(jbdy)%nbmap(1:iszdim,igrd)   ! associate the mapping used for this bdy
626                  bf_alias(1)%igrd    = igrd                                  ! used only for vertical integration of 3D arrays
627                  bf_alias(1)%ibdy    = jbdy                                  !  "    "    "     "          "      "  "    "   
628                  bf_alias(1)%ltotvel = ln_full_vel                           ! T if u3d is full velocity
629                  bf_alias(1)%lzint   = ln_zinterp                            ! T if it requires a vertical interpolation
630               ENDIF
631
632               ! associate the pointer and get rid of the dimensions with a size equal to 1
633               IF( jfld == jp_bdyssh )        dta_bdy(jbdy)%ssh => bf_alias(1)%fnow(:,1,1)
634               IF( jfld == jp_bdyu2d )        dta_bdy(jbdy)%u2d => bf_alias(1)%fnow(:,1,1)
635               IF( jfld == jp_bdyv2d )        dta_bdy(jbdy)%v2d => bf_alias(1)%fnow(:,1,1)
636               IF( jfld == jp_bdyu3d )        dta_bdy(jbdy)%u3d => bf_alias(1)%fnow(:,1,:)
637               IF( jfld == jp_bdyv3d )        dta_bdy(jbdy)%v3d => bf_alias(1)%fnow(:,1,:)
638               IF( jfld == jp_bdytem )        dta_bdy(jbdy)%tem => bf_alias(1)%fnow(:,1,:)
639               IF( jfld == jp_bdysal )        dta_bdy(jbdy)%sal => bf_alias(1)%fnow(:,1,:)
640               IF( jfld == jp_bdya_i ) THEN
641                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%a_i => bf_alias(1)%fnow(:,1,:)
642                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%a_i(iszdim,jpl) )
643                  ENDIF
644               ENDIF
645               IF( jfld == jp_bdyh_i ) THEN
646                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%h_i => bf_alias(1)%fnow(:,1,:)
647                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%h_i(iszdim,jpl) )
648                  ENDIF
649               ENDIF
650               IF( jfld == jp_bdyh_s ) THEN
651                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%h_s => bf_alias(1)%fnow(:,1,:)
652                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%h_s(iszdim,jpl) )
653                  ENDIF
654               ENDIF
655               IF( jfld == jp_bdyt_i ) THEN
656                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%t_i => bf_alias(1)%fnow(:,1,:)
657                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%t_i(iszdim,jpl) )
658                  ENDIF
659               ENDIF
660               IF( jfld == jp_bdyt_s ) THEN
661                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%t_s => bf_alias(1)%fnow(:,1,:)
662                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%t_s(iszdim,jpl) )
663                  ENDIF
664               ENDIF
665               IF( jfld == jp_bdytsu ) THEN
666                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%tsu => bf_alias(1)%fnow(:,1,:)
667                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%tsu(iszdim,jpl) )
668                  ENDIF
669               ENDIF
670               IF( jfld == jp_bdys_i ) THEN
671                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%s_i => bf_alias(1)%fnow(:,1,:)
672                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%s_i(iszdim,jpl) )
673                  ENDIF
674               ENDIF
675               IF( jfld == jp_bdyaip ) THEN
676                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%aip => bf_alias(1)%fnow(:,1,:)
677                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%aip(iszdim,jpl) )
678                  ENDIF
679               ENDIF
680               IF( jfld == jp_bdyhip ) THEN
681                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%hip => bf_alias(1)%fnow(:,1,:)
682                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%hip(iszdim,jpl) )
683                  ENDIF
684               ENDIF
685            ENDIF
686
687         END DO   ! jpbdyfld
688         !
689      END DO ! jbdy
690      !
691   END SUBROUTINE bdy_dta_init
692   
693   !!==============================================================================
694END MODULE bdydta
Note: See TracBrowser for help on using the repository browser.