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/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/BDY – NEMO

source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/BDY/bdydta.F90 @ 13257

Last change on this file since 13257 was 13257, checked in by orioltp, 4 years ago

Updated with trunk at r13245 and small change allocating variables in icb_oce.F90.

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