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

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