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

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

add the final touch to the famous gaston's branch. More precisely, add the possibility to have melt ponds as input file when using bdy

  • Property svn:keywords set to Id
File size: 39.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 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_bdytem) * tmask(ii,ij,ik)         
174                        dta_bdy(jbdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_bdysal) * 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            ENDIF
450         ENDIF
451
452         IF( .NOT.ln_pnd ) THEN
453            rn_ice_apnd = 0. ; rn_ice_hpnd = 0.
454            CALL ctl_warn( 'rn_ice_apnd & rn_ice_hpnd = 0 when no ponds' )
455         ENDIF
456
457         ! temp, salt, age and ponds of incoming ice
458         rice_tem (jbdy) = rn_ice_tem
459         rice_sal (jbdy) = rn_ice_sal
460         rice_age (jbdy) = rn_ice_age
461         rice_apnd(jbdy) = rn_ice_apnd
462         rice_hpnd(jbdy) = rn_ice_hpnd
463         
464         
465         DO jfld = 1, jpbdyfld
466
467            ! =====================
468            !          ssh
469            ! =====================
470            IF( jfld == jp_bdyssh ) THEN
471               cl3 = 'ssh'
472               igrd = 1                                                    ! T point
473               ipk = 1                                                     ! surface data
474               llneed = dta_bdy(jbdy)%lneed_ssh                            ! dta_bdy(jbdy)%ssh will be needed
475               llread = MOD(nn_dyn2d_dta(jbdy),2) == 1                     ! get data from NetCDF file
476               bf_alias => bf(jp_bdyssh,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
477               bn_alias => bn_ssh                                          ! alias for ssh structure of nambdy_dta
478               iszdim = idx_bdy(jbdy)%nblenrim(igrd)                       ! length of this bdy on this MPI processus : used only on the rim
479            ENDIF
480            ! =====================
481            !         dyn2d
482            ! =====================
483            IF( jfld == jp_bdyu2d ) THEN
484               cl3 = 'u2d'
485               igrd = 2                                                    ! U point
486               ipk = 1                                                     ! surface data
487               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed
488               llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1   ! don't get u2d from u3d and read NetCDF file
489               bf_alias => bf(jp_bdyu2d,jbdy:jbdy)                         ! alias for u2d structure of bdy number jbdy
490               bn_alias => bn_u2d                                          ! alias for u2d structure of nambdy_dta
491               IF( ln_full_vel ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)      ! will be computed from u3d -> need on the full bdy
492               ELSE                    ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)   ! used only on the rim
493               ENDIF
494            ENDIF
495            IF( jfld == jp_bdyv2d ) THEN
496               cl3 = 'v2d'
497               igrd = 3                                                    ! V point
498               ipk = 1                                                     ! surface data
499               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed
500               llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1   ! don't get v2d from v3d and read NetCDF file
501               bf_alias => bf(jp_bdyv2d,jbdy:jbdy)                         ! alias for v2d structure of bdy number jbdy
502               bn_alias => bn_v2d                                          ! alias for v2d structure of nambdy_dta
503               IF( ln_full_vel ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)      ! will be computed from v3d -> need on the full bdy
504               ELSE                    ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)   ! used only on the rim
505               ENDIF
506            ENDIF
507            ! =====================
508            !         dyn3d
509            ! =====================
510            IF( jfld == jp_bdyu3d ) THEN
511               cl3 = 'u3d'
512               igrd = 2                                                    ! U point
513               ipk = jpk                                                   ! 3d data
514               llneed = dta_bdy(jbdy)%lneed_dyn3d .OR.               &     ! dta_bdy(jbdy)%u3d will be needed
515                  &   ( dta_bdy(jbdy)%lneed_dyn2d .AND. ln_full_vel )      !   u3d needed to compute u2d
516               llread = nn_dyn3d_dta(jbdy) == 1                            ! get data from NetCDF file
517               bf_alias => bf(jp_bdyu3d,jbdy:jbdy)                         ! alias for u3d structure of bdy number jbdy
518               bn_alias => bn_u3d                                          ! alias for u3d structure of nambdy_dta
519               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
520           ENDIF
521            IF( jfld == jp_bdyv3d ) THEN
522               cl3 = 'v3d'
523               igrd = 3                                                    ! V point
524               ipk = jpk                                                   ! 3d data
525               llneed = dta_bdy(jbdy)%lneed_dyn3d .OR.               &     ! dta_bdy(jbdy)%v3d will be needed
526                  &   ( dta_bdy(jbdy)%lneed_dyn2d .AND. ln_full_vel )      !   v3d needed to compute v2d
527               llread = nn_dyn3d_dta(jbdy) == 1                            ! get data from NetCDF file
528               bf_alias => bf(jp_bdyv3d,jbdy:jbdy)                         ! alias for v3d structure of bdy number jbdy
529               bn_alias => bn_v3d                                          ! alias for v3d structure of nambdy_dta
530               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
531           ENDIF
532
533            ! =====================
534            !          tra
535            ! =====================
536            IF( jfld == jp_bdytem ) THEN
537               cl3 = 'tem'
538               igrd = 1                                                    ! T point
539               ipk = jpk                                                   ! 3d data
540               llneed = dta_bdy(jbdy)%lneed_tra                            ! dta_bdy(jbdy)%tem will be needed
541               llread = nn_tra_dta(jbdy) == 1                              ! get data from NetCDF file
542               bf_alias => bf(jp_bdytem,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
543               bn_alias => bn_tem                                          ! alias for ssh structure of nambdy_dta
544               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
545            ENDIF
546            IF( jfld == jp_bdysal ) THEN
547               cl3 = 'sal'
548               igrd = 1                                                    ! T point
549               ipk = jpk                                                   ! 3d data
550               llneed = dta_bdy(jbdy)%lneed_tra                            ! dta_bdy(jbdy)%sal will be needed
551               llread = nn_tra_dta(jbdy) == 1                              ! get data from NetCDF file
552               bf_alias => bf(jp_bdysal,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy
553               bn_alias => bn_sal                                          ! alias for ssh structure of nambdy_dta
554               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
555            ENDIF
556
557            ! =====================
558            !          ice
559            ! =====================
560            IF(  jfld == jp_bdya_i .OR. jfld == jp_bdyh_i .OR. jfld == jp_bdyh_s .OR. &
561               & jfld == jp_bdyt_i .OR. jfld == jp_bdyt_s .OR. jfld == jp_bdytsu .OR. &
562               & jfld == jp_bdys_i .OR. jfld == jp_bdyaip .OR. jfld == jp_bdyhip      ) THEN
563               igrd = 1                                                    ! T point
564               ipk = ipl                                                   ! jpl-cat data
565               llneed = dta_bdy(jbdy)%lneed_ice                            ! ice will be needed
566               llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file
567               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus
568            ENDIF
569            IF( jfld == jp_bdya_i ) THEN
570               cl3 = 'a_i'
571               bf_alias => bf(jp_bdya_i,jbdy:jbdy)                         ! alias for a_i structure of bdy number jbdy
572               bn_alias => bn_a_i                                          ! alias for a_i structure of nambdy_dta
573            ENDIF
574            IF( jfld == jp_bdyh_i ) THEN
575               cl3 = 'h_i'
576               bf_alias => bf(jp_bdyh_i,jbdy:jbdy)                         ! alias for h_i structure of bdy number jbdy
577               bn_alias => bn_h_i                                          ! alias for h_i structure of nambdy_dta
578            ENDIF
579            IF( jfld == jp_bdyh_s ) THEN
580               cl3 = 'h_s'
581               bf_alias => bf(jp_bdyh_s,jbdy:jbdy)                         ! alias for h_s structure of bdy number jbdy
582               bn_alias => bn_h_s                                          ! alias for h_s structure of nambdy_dta
583            ENDIF
584            IF( jfld == jp_bdyt_i ) THEN
585               cl3 = 't_i'
586               bf_alias => bf(jp_bdyt_i,jbdy:jbdy)                         ! alias for t_i structure of bdy number jbdy
587               bn_alias => bn_t_i                                          ! alias for t_i structure of nambdy_dta
588            ENDIF
589            IF( jfld == jp_bdyt_s ) THEN
590               cl3 = 't_s'
591               bf_alias => bf(jp_bdyt_s,jbdy:jbdy)                         ! alias for t_s structure of bdy number jbdy
592               bn_alias => bn_t_s                                          ! alias for t_s structure of nambdy_dta
593            ENDIF
594            IF( jfld == jp_bdytsu ) THEN
595               cl3 = 'tsu'
596               bf_alias => bf(jp_bdytsu,jbdy:jbdy)                         ! alias for tsu structure of bdy number jbdy
597               bn_alias => bn_tsu                                          ! alias for tsu structure of nambdy_dta
598            ENDIF
599            IF( jfld == jp_bdys_i ) THEN
600               cl3 = 's_i'
601               bf_alias => bf(jp_bdys_i,jbdy:jbdy)                         ! alias for s_i structure of bdy number jbdy
602               bn_alias => bn_s_i                                          ! alias for s_i structure of nambdy_dta
603            ENDIF
604            IF( jfld == jp_bdyaip ) THEN
605               cl3 = 'aip'
606               bf_alias => bf(jp_bdyaip,jbdy:jbdy)                         ! alias for aip structure of bdy number jbdy
607               bn_alias => bn_aip                                          ! alias for aip structure of nambdy_dta
608            ENDIF
609            IF( jfld == jp_bdyhip ) THEN
610               cl3 = 'hip'
611               bf_alias => bf(jp_bdyhip,jbdy:jbdy)                         ! alias for hip structure of bdy number jbdy
612               bn_alias => bn_hip                                          ! alias for hip structure of nambdy_dta
613            ENDIF
614
615            IF( llneed ) THEN                                              ! dta_bdy(jbdy)%xxx will be needed
616               !                                                           !   -> must be associated with an allocated target
617               ALLOCATE( bf_alias(1)%fnow( iszdim, 1, ipk ) )              ! allocate the target
618               !
619               IF( llread ) THEN                                           ! get data from NetCDF file
620                  CALL fld_fill( bf_alias, bn_alias, cn_dir, 'bdy_dta', cl3//' '//ctmp1, ctmp2 )   ! use namelist info
621                  IF( bf_alias(1)%ln_tint ) ALLOCATE( bf_alias(1)%fdta( iszdim, 1, ipk, 2 ) )
622                  bf_alias(1)%imap    => idx_bdy(jbdy)%nbmap(1:iszdim,igrd)   ! associate the mapping used for this bdy
623                  bf_alias(1)%igrd    = igrd                                  ! used only for vertical integration of 3D arrays
624                  bf_alias(1)%ltotvel = ln_full_vel                           ! T if u3d is full velocity
625                  bf_alias(1)%lzint   = ln_zinterp                            ! T if it requires a vertical interpolation
626               ENDIF
627
628               ! associate the pointer and get rid of the dimensions with a size equal to 1
629               IF( jfld == jp_bdyssh )        dta_bdy(jbdy)%ssh => bf_alias(1)%fnow(:,1,1)
630               IF( jfld == jp_bdyu2d )        dta_bdy(jbdy)%u2d => bf_alias(1)%fnow(:,1,1)
631               IF( jfld == jp_bdyv2d )        dta_bdy(jbdy)%v2d => bf_alias(1)%fnow(:,1,1)
632               IF( jfld == jp_bdyu3d )        dta_bdy(jbdy)%u3d => bf_alias(1)%fnow(:,1,:)
633               IF( jfld == jp_bdyv3d )        dta_bdy(jbdy)%v3d => bf_alias(1)%fnow(:,1,:)
634               IF( jfld == jp_bdytem )        dta_bdy(jbdy)%tem => bf_alias(1)%fnow(:,1,:)
635               IF( jfld == jp_bdysal )        dta_bdy(jbdy)%sal => bf_alias(1)%fnow(:,1,:)
636               IF( jfld == jp_bdya_i ) THEN
637                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%a_i => bf_alias(1)%fnow(:,1,:)
638                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%a_i(iszdim,jpl) )
639                  ENDIF
640               ENDIF
641               IF( jfld == jp_bdyh_i ) THEN
642                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%h_i => bf_alias(1)%fnow(:,1,:)
643                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%h_i(iszdim,jpl) )
644                  ENDIF
645               ENDIF
646               IF( jfld == jp_bdyh_s ) THEN
647                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%h_s => bf_alias(1)%fnow(:,1,:)
648                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%h_s(iszdim,jpl) )
649                  ENDIF
650               ENDIF
651               IF( jfld == jp_bdyt_i ) THEN
652                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%t_i => bf_alias(1)%fnow(:,1,:)
653                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%t_i(iszdim,jpl) )
654                  ENDIF
655               ENDIF
656               IF( jfld == jp_bdyt_s ) THEN
657                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%t_s => bf_alias(1)%fnow(:,1,:)
658                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%t_s(iszdim,jpl) )
659                  ENDIF
660               ENDIF
661               IF( jfld == jp_bdytsu ) THEN
662                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%tsu => bf_alias(1)%fnow(:,1,:)
663                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%tsu(iszdim,jpl) )
664                  ENDIF
665               ENDIF
666               IF( jfld == jp_bdys_i ) THEN
667                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%s_i => bf_alias(1)%fnow(:,1,:)
668                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%s_i(iszdim,jpl) )
669                  ENDIF
670               ENDIF
671               IF( jfld == jp_bdyaip ) THEN
672                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%aip => bf_alias(1)%fnow(:,1,:)
673                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%aip(iszdim,jpl) )
674                  ENDIF
675               ENDIF
676               IF( jfld == jp_bdyhip ) THEN
677                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%hip => bf_alias(1)%fnow(:,1,:)
678                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%hip(iszdim,jpl) )
679                  ENDIF
680               ENDIF
681            ENDIF
682
683         END DO   ! jpbdyfld
684         !
685      END DO ! jbdy
686      !
687   END SUBROUTINE bdy_dta_init
688   
689   !!==============================================================================
690END MODULE bdydta
Note: See TracBrowser for help on using the repository browser.