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 branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90 @ 2287

Last change on this file since 2287 was 2287, checked in by smasson, 13 years ago

update licence of all NEMO files...

  • Property svn:keywords set to Id
File size: 50.1 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_bt
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   !!----------------------------------------------------------------------
13#if defined key_bdy
14   !!----------------------------------------------------------------------
15   !!   'key_bdy'                     Unstructured Open Boundary Conditions
16   !!----------------------------------------------------------------------
17   !!   bdy_dta    : read u, v, t, s data along open boundaries
18   !!   bdy_dta_bt : read depth-mean velocities and elevation along open boundaries       
19   !!----------------------------------------------------------------------
20   USE oce             ! ocean dynamics and tracers
21   USE dom_oce         ! ocean space and time domain
22   USE phycst          ! physical constants
23   USE bdy_oce         ! ocean open boundary conditions
24   USE bdytides        ! tidal forcing at boundaries
25   USE iom
26   USE ioipsl
27   USE in_out_manager  ! I/O logical units
28#if defined key_lim2
29   USE ice_2
30#endif
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   bdy_dta      ! routines called by step.F90
36   PUBLIC   bdy_dta_bt 
37
38   INTEGER ::   numbdyt, numbdyu, numbdyv                      ! logical units for T-, U-, & V-points data file, resp.
39   INTEGER ::   ntimes_bdy                                     ! exact number of time dumps in data files
40   INTEGER ::   nbdy_b, nbdy_a                                 ! record of bdy data file for before and after time step
41   INTEGER ::   numbdyt_bt, numbdyu_bt, numbdyv_bt             ! logical unit for T-, U- & V-points data file, resp.
42   INTEGER ::   ntimes_bdy_bt                                  ! exact number of time dumps in data files
43   INTEGER ::   nbdy_b_bt, nbdy_a_bt                           ! record of bdy data file for before and after time step
44
45   INTEGER, DIMENSION (jpbtime) ::   istep, istep_bt           ! time array in seconds in each data file
46
47   REAL(wp) ::  zoffset                                        ! time offset between time origin in file & start time of model run
48
49   REAL(wp), DIMENSION(jpbdim,jpk,2) ::   tbdydta, sbdydta     ! time interpolated values of T and S bdy data   
50   REAL(wp), DIMENSION(jpbdim,jpk,2) ::   ubdydta, vbdydta     ! time interpolated values of U and V bdy data
51   REAL(wp), DIMENSION(jpbdim,2)     ::   ubtbdydta, vbtbdydta ! Arrays used for time interpolation of bdy data   
52   REAL(wp), DIMENSION(jpbdim,2)     ::   sshbdydta            ! bdy data of ssh
53
54#if defined key_lim2
55   REAL(wp), DIMENSION(jpbdim,2)     ::   frld_bdydta          ! }
56   REAL(wp), DIMENSION(jpbdim,2)     ::   hicif_bdydta         ! } Arrays used for time interp. of ice bdy data
57   REAL(wp), DIMENSION(jpbdim,2)     ::   hsnif_bdydta         ! }
58#endif
59
60   !!----------------------------------------------------------------------
61   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
62   !! $Id$
63   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
64   !!----------------------------------------------------------------------
65CONTAINS
66
67   SUBROUTINE bdy_dta( kt )
68      !!----------------------------------------------------------------------
69      !!                   ***  SUBROUTINE bdy_dta  ***
70      !!                   
71      !! ** Purpose :   Read unstructured boundary data for FRS condition.
72      !!
73      !! ** Method  :   At the first timestep, read in boundary data for two
74      !!                times from the file and time-interpolate. At other
75      !!                timesteps, check to see if we need another time from
76      !!                the file. If so read it in. Time interpolate.
77      !!----------------------------------------------------------------------
78      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index (for timesplitting option, otherwise zero)
79      !!
80      CHARACTER(LEN=80), DIMENSION(3) ::   clfile               ! names of input files
81      CHARACTER(LEN=70 )              ::   clunits              ! units attribute of time coordinate
82      LOGICAL ::   lect                                         ! flag for reading
83      INTEGER ::   it, ib, ik, igrd                             ! dummy loop indices
84      INTEGER ::   igrd_start, igrd_end                         ! start and end of loops on igrd
85      INTEGER ::   idvar                                        ! netcdf var ID
86      INTEGER ::   iman, i15, imois                             ! Time variables for monthly clim forcing
87      INTEGER ::   ntimes_bdyt, ntimes_bdyu, ntimes_bdyv
88      INTEGER ::   itimer, totime
89      INTEGER ::   ii, ij                                       ! array addresses
90      INTEGER ::   ipi, ipj, ipk, inum                          ! local integers (NetCDF read)
91      INTEGER ::   iyear0, imonth0, iday0
92      INTEGER ::   ihours0, iminutes0, isec0
93      INTEGER ::   iyear, imonth, iday, isecs
94      INTEGER, DIMENSION(jpbtime) ::   istept, istepu, istepv   ! time arrays from data files
95      REAL(wp) ::   dayfrac, zxy, zoffsett
96      REAL(wp) ::   zoffsetu, zoffsetv
97      REAL(wp) ::   dayjul0, zdayjulini
98      REAL(wp), DIMENSION(jpbtime)      ::   zstepr             ! REAL time array from data files
99      REAL(wp), DIMENSION(jpbdta,1,jpk) ::   zdta               ! temporary array for data fields
100      !!---------------------------------------------------------------------------
101
102
103      IF( ln_bdy_dyn_frs .OR. ln_bdy_tra_frs    &
104         &               .OR. ln_bdy_ice_frs ) THEN  ! If these are both false then this routine does nothing
105
106      ! -------------------- !
107      !    Initialization    !
108      ! -------------------- !
109
110      lect   = .false.           ! If true, read a time record
111
112      ! Some time variables for monthly climatological forcing:
113      ! *******************************************************
114
115!!gm  here  use directely daymod calendar variables
116 
117      iman = INT( raamo )      ! Number of months in a year
118
119      i15 = INT( 2*REAL( nday, wp ) / ( REAL( nmonth_len(nmonth), wp ) + 0.5 ) )
120      ! i15=0 if the current day is in the first half of the month, else i15=1
121
122      imois = nmonth + i15 - 1            ! imois is the first month record
123      IF( imois == 0 )   imois = iman
124
125      ! Time variable for non-climatological forcing:
126      ! *********************************************
127      itimer = (kt-nit000+1)*rdt      ! current time in seconds for interpolation
128
129
130      !                                                !-------------------!
131      IF( kt == nit000 ) THEN                          !  First call only  !
132         !                                             !-------------------!
133         istep(:) = 0
134         nbdy_b   = 0
135         nbdy_a   = 0
136
137         ! Get time information from bdy data file
138         ! ***************************************
139
140         IF(lwp) WRITE(numout,*)
141         IF(lwp) WRITE(numout,*)    'bdy_dta : Initialize unstructured boundary data'
142         IF(lwp) WRITE(numout,*)    '~~~~~~~' 
143
144         IF     ( nbdy_dta == 0 ) THEN
145            !
146            IF(lwp) WRITE(numout,*) '          Bdy data are taken from initial conditions'
147            !
148         ELSEIF (nbdy_dta == 1) THEN
149            !
150            IF(lwp) WRITE(numout,*) '          Bdy data are read in netcdf files'
151            !
152            dayfrac = adatrj  - REAL( itimer, wp ) / 86400.   ! day fraction at time step kt-1
153            dayfrac = dayfrac - INT ( dayfrac )               !
154            totime  = ( nitend - nit000 + 1 ) * rdt           ! Total time of the run to verify that all the
155            !                                                 ! necessary time dumps in file are included
156            !
157            clfile(1) = filbdy_data_T
158            clfile(2) = filbdy_data_U
159            clfile(3) = filbdy_data_V
160            !                                                 
161            ! how many files are we to read in?
162            igrd_start = 1
163            igrd_end   = 3
164            IF(.NOT. ln_bdy_tra_frs .AND. .NOT. ln_bdy_ice_frs) THEN       ! No T-grid file.
165               igrd_start = 2
166            ELSEIF ( .NOT. ln_bdy_dyn_frs ) THEN                           ! No U-grid or V-grid file.
167               igrd_end   = 1         
168            ENDIF
169
170            DO igrd = igrd_start, igrd_end                     !  loop over T, U & V grid  !
171               !                                               !---------------------------!
172               CALL iom_open( clfile(igrd), inum )
173               CALL iom_gettime( inum, zstepr, kntime=ntimes_bdy, cdunits=clunits ) 
174
175               SELECT CASE( igrd )
176                  CASE (1)   ;   numbdyt = inum
177                  CASE (2)   ;   numbdyu = inum
178                  CASE (3)   ;   numbdyv = inum
179               END SELECT
180
181               ! Calculate time offset
182               READ(clunits,7000) iyear0, imonth0, iday0, ihours0, iminutes0, isec0
183               ! Convert time origin in file to julian days
184               isec0 = isec0 + ihours0*60.*60. + iminutes0*60.
185               CALL ymds2ju(iyear0, imonth0, iday0, REAL(isec0, wp), dayjul0)
186               ! Compute model initialization time
187               iyear  = ndastp / 10000
188               imonth = ( ndastp - iyear * 10000 ) / 100
189               iday   = ndastp - iyear * 10000 - imonth * 100
190               isecs  = dayfrac * 86400
191               CALL ymds2ju(iyear, imonth, iday, REAL(isecs, wp) , zdayjulini)
192               ! offset from initialization date:
193               zoffset = (dayjul0-zdayjulini)*86400
194               !
1957000           FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2)
196
197               !! TO BE DONE... Check consistency between calendar from file
198               !! (available optionally from iom_gettime) and calendar in model
199               !! when calendar in model available outside of IOIPSL.
200
201               IF(lwp) WRITE(numout,*) 'number of times: ',ntimes_bdy
202               IF(lwp) WRITE(numout,*) 'offset: ',zoffset
203               IF(lwp) WRITE(numout,*) 'totime: ',totime
204               IF(lwp) WRITE(numout,*) 'zstepr: ',zstepr(1:ntimes_bdy)
205
206               ! Check that there are not too many times in the file.
207               IF( ntimes_bdy > jpbtime ) THEN
208                  WRITE(ctmp1,*) 'Check file: ', clfile(igrd), 'jpbtime= ', jpbtime, ' ntimes_bdy= ', ntimes_bdy
209                  CALL ctl_stop( 'Number of time dumps in files exceed jpbtime parameter', ctmp1 )
210               ENDIF
211
212               ! Check that time array increases:
213               it = 1
214               DO WHILE( zstepr(it+1) > zstepr(it) .AND. it /= ntimes_bdy - 1 ) 
215                  it = it + 1
216               END DO
217               !
218               IF( it /= ntimes_bdy-1 .AND. ntimes_bdy > 1 ) THEN
219                     WRITE(ctmp1,*) 'Check file: ', clfile(igrd)
220                     CALL ctl_stop( 'Time array in unstructured boundary data files',   &
221                        &           'does not continuously increase.'               , ctmp1 )
222               ENDIF
223               !
224               ! Check that times in file span model run time:
225               IF( zstepr(1) + zoffset > 0 ) THEN
226                     WRITE(ctmp1,*) 'Check file: ', clfile(igrd)
227                     CALL ctl_stop( 'First time dump in bdy file is after model initial time', ctmp1 )
228               END IF
229               IF( zstepr(ntimes_bdy) + zoffset < totime ) THEN
230                     WRITE(ctmp1,*) 'Check file: ', clfile(igrd)
231                     CALL ctl_stop( 'Last time dump in bdy file is before model final time', ctmp1 )
232               END IF
233               !
234               SELECT CASE( igrd )
235                  CASE (1)
236                    ntimes_bdyt = ntimes_bdy
237                    zoffsett = zoffset
238                    istept(:) = INT( zstepr(:) + zoffset )
239                    numbdyt = inum
240                  CASE (2)
241                    ntimes_bdyu = ntimes_bdy
242                    zoffsetu = zoffset
243                    istepu(:) = INT( zstepr(:) + zoffset )
244                    numbdyu = inum
245                  CASE (3)
246                    ntimes_bdyv = ntimes_bdy
247                    zoffsetv = zoffset
248                    istepv(:) = INT( zstepr(:) + zoffset )
249                    numbdyv = inum
250               END SELECT
251               !
252            END DO                                         ! end loop over T, U & V grid
253
254            IF (igrd_start == 1 .and. igrd_end == 3) THEN
255               ! Only test differences if we are reading in 3 files
256               ! Verify time consistency between files 
257               IF( ntimes_bdyu /= ntimes_bdyt .OR. ntimes_bdyv /= ntimes_bdyt ) THEN
258                  CALL ctl_stop( 'Bdy data files must have the same number of time dumps',   &
259                  &           'Multiple time frequencies not implemented yet'  )
260               ENDIF
261               ntimes_bdy = ntimes_bdyt
262               !
263               IF( zoffsetu /= zoffsett .OR. zoffsetv /= zoffsett ) THEN
264                  CALL ctl_stop( 'Bdy data files must have the same time origin',   &
265                  &           'Multiple time frequencies not implemented yet' )
266               ENDIF
267               zoffset = zoffsett
268            ENDIF
269
270            IF( igrd_start == 1 ) THEN   ;   istep(:) = istept(:)
271            ELSE                         ;   istep(:) = istepu(:)
272            ENDIF
273
274            ! Check number of time dumps:             
275            IF( ntimes_bdy == 1 .AND. .NOT. ln_bdy_clim ) THEN
276              CALL ctl_stop( 'There is only one time dump in data files',   &
277                 &           'Choose ln_bdy_clim=.true. in namelist for constant bdy forcing.' )
278            ENDIF
279
280            IF( ln_bdy_clim ) THEN
281              IF( ntimes_bdy /= 1 .AND. ntimes_bdy /= 12 ) THEN
282                 CALL ctl_stop( 'For climatological boundary forcing (ln_bdy_clim=.true.),',   &
283                    &           'bdy data files must contain 1 or 12 time dumps.' )
284              ELSEIF( ntimes_bdy ==  1 ) THEN
285                IF(lwp) WRITE(numout,*)
286                IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files'
287              ELSEIF( ntimes_bdy == 12 ) THEN
288                IF(lwp) WRITE(numout,*)
289                IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files'
290              ENDIF
291            ENDIF
292
293            ! Find index of first record to read (before first model time).
294            it = 1
295            DO WHILE( istep(it+1) <= 0 .AND. it <= ntimes_bdy - 1 )
296               it = it + 1
297            END DO
298            nbdy_b = it
299            !
300            IF(lwp) WRITE(numout,*) 'Time offset is ',zoffset
301            IF(lwp) WRITE(numout,*) 'First record to read is ',nbdy_b
302
303         ENDIF ! endif (nbdy_dta == 1)
304
305
306         ! 1.2  Read first record in file if necessary (ie if nbdy_dta == 1)
307         ! *****************************************************************
308
309         IF( nbdy_dta == 0 ) THEN      ! boundary data arrays are filled with initial conditions
310            !
311            IF (ln_bdy_tra_frs) THEN
312               igrd = 1            ! T-points data
313               DO ib = 1, nblen(igrd)
314                  ii = nbi(ib,igrd)
315                  ij = nbj(ib,igrd)
316                  DO ik = 1, jpkm1
317                     tbdy(ib,ik) = tn(ii,ij,ik)
318                     sbdy(ib,ik) = sn(ii,ij,ik)
319                  END DO
320               END DO
321            ENDIF
322
323            IF(ln_bdy_dyn_frs) THEN
324               igrd = 2            ! U-points data
325               DO ib = 1, nblen(igrd)
326                  ii = nbi(ib,igrd)
327                  ij = nbj(ib,igrd)
328                  DO ik = 1, jpkm1
329                     ubdy(ib,ik) = un(ii, ij, ik)
330                  END DO
331               END DO
332               !
333               igrd = 3            ! V-points data
334               DO ib = 1, nblen(igrd)           
335                  ii = nbi(ib,igrd)
336                  ij = nbj(ib,igrd)
337                  DO ik = 1, jpkm1
338                     vbdy(ib,ik) = vn(ii, ij, ik)
339                  END DO
340               END DO
341            ENDIF
342            !
343#if defined key_lim2
344            IF( ln_bdy_ice_frs ) THEN
345               igrd = 1            ! T-points data
346               DO ib = 1, nblen(igrd)
347                  frld_bdy (ib) =  frld(nbi(ib,igrd), nbj(ib,igrd))
348                  hicif_bdy(ib) = hicif(nbi(ib,igrd), nbj(ib,igrd))
349                  hsnif_bdy(ib) = hsnif(nbi(ib,igrd), nbj(ib,igrd))
350               END DO
351            ENDIF
352#endif
353         ELSEIF( nbdy_dta == 1 ) THEN    ! Set first record in the climatological case:   
354            !
355            IF( ln_bdy_clim .AND. ntimes_bdy == 1 ) THEN
356               nbdy_a = 1
357            ELSEIF( ln_bdy_clim .AND. ntimes_bdy == iman ) THEN
358               nbdy_b = 0
359               nbdy_a = imois
360            ELSE
361               nbdy_a = nbdy_b
362            ENDIF
363   
364            ! Read first record:
365            ipj  = 1
366            ipk  = jpk
367            igrd = 1
368            ipi  = nblendta(igrd)
369
370            IF(ln_bdy_tra_frs) THEN
371               !
372               igrd = 1                                           ! Temperature
373               IF( nblendta(igrd) <=  0 ) THEN
374                  idvar = iom_varid( numbdyt, 'votemper' )
375                  nblendta(igrd) = iom_file(numbdyt)%dimsz(1,idvar)
376               ENDIF
377               IF(lwp) WRITE(numout,*) 'Dim size for votemper is ', nblendta(igrd)
378               ipi = nblendta(igrd)
379               CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a )
380               !
381               DO ib = 1, nblen(igrd)
382                  DO ik = 1, jpkm1
383                     tbdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik)
384                  END DO
385               END DO
386               !
387               igrd = 1                                           ! salinity
388               IF( nblendta(igrd) .le. 0 ) THEN
389                  idvar = iom_varid( numbdyt, 'vosaline' )
390                  nblendta(igrd) = iom_file(numbdyt)%dimsz(1,idvar)
391               ENDIF
392               IF(lwp) WRITE(numout,*) 'Dim size for vosaline is ', nblendta(igrd)
393               ipi = nblendta(igrd)
394               CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a )
395               !
396               DO ib = 1, nblen(igrd)
397                  DO ik = 1, jpkm1
398                     sbdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik)
399                  END DO
400               END DO
401            ENDIF  ! ln_bdy_tra_frs
402 
403            IF( ln_bdy_dyn_frs ) THEN
404               !
405               igrd = 2                                           ! u-velocity
406               IF ( nblendta(igrd) .le. 0 ) THEN
407                 idvar = iom_varid( numbdyu,'vozocrtx' )
408                 nblendta(igrd) = iom_file(numbdyu)%dimsz(1,idvar)
409               ENDIF
410               IF(lwp) WRITE(numout,*) 'Dim size for vozocrtx is ', nblendta(igrd)
411               ipi = nblendta(igrd)
412               CALL iom_get ( numbdyu, jpdom_unknown,'vozocrtx',zdta(1:ipi,1:ipj,1:ipk),nbdy_a )
413               DO ib = 1, nblen(igrd)
414                  DO ik = 1, jpkm1
415                     ubdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik)
416                  END DO
417               END DO
418               !
419               igrd = 3                                           ! v-velocity
420               IF ( nblendta(igrd) .le. 0 ) THEN
421                 idvar = iom_varid( numbdyv,'vomecrty' )
422                 nblendta(igrd) = iom_file(numbdyv)%dimsz(1,idvar)
423               ENDIF
424               IF(lwp) WRITE(numout,*) 'Dim size for vomecrty is ', nblendta(igrd)
425               ipi = nblendta(igrd)
426               CALL iom_get ( numbdyv, jpdom_unknown,'vomecrty',zdta(1:ipi,1:ipj,1:ipk),nbdy_a )
427               DO ib = 1, nblen(igrd)
428                  DO ik = 1, jpkm1
429                     vbdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik)
430                  END DO
431               END DO
432            ENDIF ! ln_bdy_dyn_frs
433
434#if defined key_lim2
435            IF( ln_bdy_ice_frs ) THEN
436              !
437              igrd=1                                              ! leads fraction
438              IF(lwp) WRITE(numout,*) 'Dim size for ildsconc is ',nblendta(igrd)
439              ipi=nblendta(igrd)
440              CALL iom_get ( numbdyt, jpdom_unknown,'ildsconc',zdta(1:ipi,:,1),nbdy_a )
441              DO ib=1, nblen(igrd)
442                frld_bdydta(ib,2) =  zdta(nbmap(ib,igrd),1,1)
443              END DO
444              !
445              igrd=1                                              ! ice thickness
446              IF(lwp) WRITE(numout,*) 'Dim size for iicethic is ',nblendta(igrd)
447              ipi=nblendta(igrd)
448              CALL iom_get ( numbdyt, jpdom_unknown,'iicethic',zdta(1:ipi,:,1),nbdy_a )
449              DO ib=1, nblen(igrd)
450                hicif_bdydta(ib,2) =  zdta(nbmap(ib,igrd),1,1)
451              END DO
452              !
453              igrd=1                                              ! snow thickness
454              IF(lwp) WRITE(numout,*) 'Dim size for isnowthi is ',nblendta(igrd)
455              ipi=nblendta(igrd)
456              CALL iom_get ( numbdyt, jpdom_unknown,'isnowthi',zdta(1:ipi,:,1),nbdy_a )
457              DO ib=1, nblen(igrd)
458                hsnif_bdydta(ib,2) =  zdta(nbmap(ib,igrd),1,1)
459              END DO
460            ENDIF ! just if ln_bdy_ice_frs is set
461#endif
462
463            IF( .NOT.ln_bdy_clim .AND. istep(1) > 0 ) THEN     ! First data time is after start of run
464               nbdy_b = nbdy_a                                 ! Put first value in both time levels
465               IF( ln_bdy_tra_frs ) THEN
466                 tbdydta(:,:,1) = tbdydta(:,:,2)
467                 sbdydta(:,:,1) = sbdydta(:,:,2)
468               ENDIF
469               IF( ln_bdy_dyn_frs ) THEN
470                 ubdydta(:,:,1) = ubdydta(:,:,2)
471                 vbdydta(:,:,1) = vbdydta(:,:,2)
472               ENDIF
473#if defined key_lim2
474               IF( ln_bdy_ice_frs ) THEN
475                  frld_bdydta (:,1) =  frld_bdydta(:,2)
476                  hicif_bdydta(:,1) = hicif_bdydta(:,2)
477                  hsnif_bdydta(:,1) = hsnif_bdydta(:,2)
478               ENDIF
479#endif
480            END IF
481            !
482         END IF   ! nbdy_dta == 0/1
483 
484         ! In the case of constant boundary forcing fill bdy arrays once for all
485         IF( ln_bdy_clim .AND. ntimes_bdy == 1 ) THEN
486            IF( ln_bdy_tra_frs ) THEN
487               tbdy  (:,:) = tbdydta  (:,:,2)
488               sbdy  (:,:) = sbdydta  (:,:,2)
489            ENDIF
490            IF( ln_bdy_dyn_frs) THEN
491               ubdy  (:,:) = ubdydta  (:,:,2)
492               vbdy  (:,:) = vbdydta  (:,:,2)
493            ENDIF
494#if defined key_lim2
495            IF( ln_bdy_ice_frs ) THEN
496               frld_bdy (:) = frld_bdydta (:,2)
497               hicif_bdy(:) = hicif_bdydta(:,2)
498               hsnif_bdy(:) = hsnif_bdydta(:,2)
499            ENDIF
500#endif
501
502            IF( ln_bdy_tra_frs .OR. ln_bdy_ice_frs) CALL iom_close( numbdyt )
503            IF( ln_bdy_dyn_frs                    ) CALL iom_close( numbdyu )
504            IF( ln_bdy_dyn_frs                    ) CALL iom_close( numbdyv )
505         END IF
506         !
507      ENDIF                                            ! End if nit000
508
509
510      !                                                !---------------------!
511      IF( nbdy_dta == 1 .AND. ntimes_bdy > 1 ) THEN    !  at each time step  !
512         !                                             !---------------------!
513         ! Read one more record if necessary
514         !**********************************
515
516         IF( ln_bdy_clim .AND. imois /= nbdy_b ) THEN      ! remember that nbdy_b=0 for kt=nit000
517            nbdy_b = imois
518            nbdy_a = imois + 1
519            nbdy_b = MOD( nbdy_b, iman )   ;   IF( nbdy_b == 0 ) nbdy_b = iman
520            nbdy_a = MOD( nbdy_a, iman )   ;   IF( nbdy_a == 0 ) nbdy_a = iman
521            lect=.true.
522         ELSEIF( .NOT.ln_bdy_clim .AND. itimer >= istep(nbdy_a) ) THEN
523
524            IF( nbdy_a < ntimes_bdy ) THEN
525               nbdy_b = nbdy_a
526               nbdy_a = nbdy_a + 1
527               lect  =.true.
528            ELSE
529               ! We have reached the end of the file
530               ! put the last data time into both time levels
531               nbdy_b = nbdy_a
532               IF(ln_bdy_tra_frs) THEN
533                  tbdydta(:,:,1) =  tbdydta(:,:,2)
534                  sbdydta(:,:,1) =  sbdydta(:,:,2)
535               ENDIF
536               IF(ln_bdy_dyn_frs) THEN
537                  ubdydta(:,:,1) =  ubdydta(:,:,2)
538                  vbdydta(:,:,1) =  vbdydta(:,:,2)
539               ENDIF
540#if defined key_lim2
541               IF(ln_bdy_ice_frs) THEN
542                  frld_bdydta (:,1) =  frld_bdydta (:,2)
543                  hicif_bdydta(:,1) =  hicif_bdydta(:,2)
544                  hsnif_bdydta(:,1) =  hsnif_bdydta(:,2)
545               ENDIF
546#endif
547            END IF ! nbdy_a < ntimes_bdy
548            !
549        END IF
550         
551        IF( lect ) THEN           ! Swap arrays
552           IF( ln_bdy_tra_frs ) THEN
553             tbdydta(:,:,1) =  tbdydta(:,:,2)
554             sbdydta(:,:,1) =  sbdydta(:,:,2)
555           ENDIF
556           IF( ln_bdy_dyn_frs ) THEN
557             ubdydta(:,:,1) =  ubdydta(:,:,2)
558             vbdydta(:,:,1) =  vbdydta(:,:,2)
559           ENDIF
560#if defined key_lim2
561           IF( ln_bdy_ice_frs ) THEN
562             frld_bdydta (:,1) =  frld_bdydta (:,2)
563             hicif_bdydta(:,1) =  hicif_bdydta(:,2)
564             hsnif_bdydta(:,1) =  hsnif_bdydta(:,2)
565           ENDIF
566#endif 
567           ! read another set
568           ipj  = 1
569           ipk  = jpk
570
571           IF( ln_bdy_tra_frs ) THEN
572              !
573              igrd = 1                                   ! temperature
574              ipi  = nblendta(igrd)
575              CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a )
576              DO ib = 1, nblen(igrd)
577                 DO ik = 1, jpkm1
578                    tbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik)
579                 END DO
580              END DO
581              !
582              igrd = 1                                   ! salinity
583              ipi  = nblendta(igrd)
584              CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a )
585              DO ib = 1, nblen(igrd)
586                 DO ik = 1, jpkm1
587                    sbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik)
588                 END DO
589              END DO
590           ENDIF ! ln_bdy_tra_frs
591
592           IF(ln_bdy_dyn_frs) THEN
593              !
594              igrd = 2                                   ! u-velocity
595              ipi  = nblendta(igrd)
596              CALL iom_get ( numbdyu, jpdom_unknown,'vozocrtx',zdta(1:ipi,1:ipj,1:ipk),nbdy_a )
597              DO ib = 1, nblen(igrd)
598                DO ik = 1, jpkm1
599                  ubdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik)
600                END DO
601              END DO
602              !
603              igrd = 3                                   ! v-velocity
604              ipi  = nblendta(igrd)
605              CALL iom_get ( numbdyv, jpdom_unknown,'vomecrty',zdta(1:ipi,1:ipj,1:ipk),nbdy_a )
606              DO ib = 1, nblen(igrd)
607                 DO ik = 1, jpkm1
608                    vbdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik)
609                 END DO
610              END DO
611           ENDIF ! ln_bdy_dyn_frs
612           !
613#if defined key_lim2
614           IF(ln_bdy_ice_frs) THEN
615             !
616             igrd = 1                                    ! ice concentration
617             ipi=nblendta(igrd)
618             CALL iom_get ( numbdyt, jpdom_unknown,'ildsconc',zdta(1:ipi,:,1),nbdy_a )
619             DO ib=1, nblen(igrd)
620               frld_bdydta(ib,2) =  zdta( nbmap(ib,igrd), 1, 1 )
621             END DO
622             !
623             igrd=1                                      ! ice thickness
624             ipi=nblendta(igrd)
625             CALL iom_get ( numbdyt, jpdom_unknown,'iicethic',zdta(1:ipi,:,1),nbdy_a )
626             DO ib=1, nblen(igrd)
627               hicif_bdydta(ib,2) =  zdta( nbmap(ib,igrd), 1, 1 )
628             END DO
629             !
630             igrd=1                                      ! snow thickness
631             ipi=nblendta(igrd)
632             CALL iom_get ( numbdyt, jpdom_unknown,'isnowthi',zdta(1:ipi,:,1),nbdy_a )
633             DO ib=1, nblen(igrd)
634               hsnif_bdydta(ib,2) =  zdta( nbmap(ib,igrd), 1, 1 )
635             END DO
636           ENDIF ! ln_bdy_ice_frs
637#endif
638           !
639           IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy_b ',nbdy_b
640           IF(lwp) WRITE(numout,*) '~~~~~~~~  last  record file used nbdy_a ',nbdy_a
641           IF (.NOT.ln_bdy_clim) THEN
642              IF(lwp) WRITE(numout,*) 'first  record time (s): ', istep(nbdy_b)
643              IF(lwp) WRITE(numout,*) 'model time (s)        : ', itimer
644              IF(lwp) WRITE(numout,*) 'second record time (s): ', istep(nbdy_a)
645           ENDIF
646           !
647       ENDIF ! end lect=.true.
648
649
650       ! Interpolate linearly
651       ! ********************
652       !
653       IF( ln_bdy_clim ) THEN   ;   zxy = REAL( nday                   ) / REAL( nmonth_len(nbdy_b) ) + 0.5 - i15
654       ELSEIF( istep(nbdy_b) == istep(nbdy_a) ) THEN
655                                    zxy = 0.0_wp
656       ELSE                     ;   zxy = REAL( istep(nbdy_b) - itimer ) / REAL( istep(nbdy_b) - istep(nbdy_a) )
657       END IF
658
659          IF(ln_bdy_tra_frs) THEN
660             igrd = 1                                   ! temperature & salinity
661             DO ib = 1, nblen(igrd)
662               DO ik = 1, jpkm1
663                 tbdy(ib,ik) = zxy * tbdydta(ib,ik,2) + (1.-zxy) * tbdydta(ib,ik,1)
664                 sbdy(ib,ik) = zxy * sbdydta(ib,ik,2) + (1.-zxy) * sbdydta(ib,ik,1)
665               END DO
666             END DO
667          ENDIF
668
669          IF(ln_bdy_dyn_frs) THEN
670             igrd = 2                                   ! u-velocity
671             DO ib = 1, nblen(igrd)
672               DO ik = 1, jpkm1
673                 ubdy(ib,ik) = zxy * ubdydta(ib,ik,2) + (1.-zxy) * ubdydta(ib,ik,1)   
674               END DO
675             END DO
676             !
677             igrd = 3                                   ! v-velocity
678             DO ib = 1, nblen(igrd)
679               DO ik = 1, jpkm1
680                 vbdy(ib,ik) = zxy * vbdydta(ib,ik,2) + (1.-zxy) * vbdydta(ib,ik,1)   
681               END DO
682             END DO
683          ENDIF
684
685#if defined key_lim2
686          IF(ln_bdy_ice_frs) THEN
687            igrd=1
688            DO ib=1, nblen(igrd)
689               frld_bdy(ib) = zxy *  frld_bdydta(ib,2) + (1.-zxy) *  frld_bdydta(ib,1)
690              hicif_bdy(ib) = zxy * hicif_bdydta(ib,2) + (1.-zxy) * hicif_bdydta(ib,1)
691              hsnif_bdy(ib) = zxy * hsnif_bdydta(ib,2) + (1.-zxy) * hsnif_bdydta(ib,1)
692            END DO
693          ENDIF ! just if ln_bdy_ice_frs is true
694#endif
695
696      END IF                       !end if ((nbdy_dta==1).AND.(ntimes_bdy>1))
697   
698
699      !                                                !---------------------!
700      !                                                !     last call       !
701      !                                                !---------------------!
702      IF( kt == nitend ) THEN
703          IF(ln_bdy_tra_frs .or. ln_bdy_ice_frs) CALL iom_close( numbdyt )              ! Closing of the 3 files
704          IF(ln_bdy_dyn_frs) CALL iom_close( numbdyu )
705          IF(ln_bdy_dyn_frs) CALL iom_close( numbdyv )
706      ENDIF
707      !
708      ENDIF ! ln_bdy_dyn_frs .OR. ln_bdy_tra_frs
709      !
710   END SUBROUTINE bdy_dta
711
712
713   SUBROUTINE bdy_dta_bt( kt, jit, icycl )
714      !!---------------------------------------------------------------------------
715      !!                      ***  SUBROUTINE bdy_dta_bt  ***
716      !!                   
717      !! ** Purpose :   Read unstructured boundary data for Flather condition
718      !!
719      !! ** Method  :  At the first timestep, read in boundary data for two
720      !!               times from the file and time-interpolate. At other
721      !!               timesteps, check to see if we need another time from
722      !!               the file. If so read it in. Time interpolate.
723      !!---------------------------------------------------------------------------
724!!gm DOCTOR names :   argument integer :  start with "k"
725      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
726      INTEGER, INTENT( in ) ::   jit         ! barotropic time step index
727      INTEGER, INTENT( in ) ::   icycl       ! number of cycles need for final file close
728      !                                      ! (for timesplitting option, otherwise zero)
729      !!
730      LOGICAL ::   lect                      ! flag for reading
731      INTEGER ::   it, ib, igrd              ! dummy loop indices
732      INTEGER ::   idvar                     ! netcdf var ID
733      INTEGER ::   iman, i15, imois          ! Time variables for monthly clim forcing
734      INTEGER ::   ntimes_bdyt, ntimes_bdyu, ntimes_bdyv
735      INTEGER ::   itimer, totime
736      INTEGER ::   ipi, ipj, ipk, inum       ! temporary integers (NetCDF read)
737      INTEGER ::   iyear0, imonth0, iday0
738      INTEGER ::   ihours0, iminutes0, isec0
739      INTEGER ::   iyear, imonth, iday, isecs
740      INTEGER, DIMENSION(jpbtime) ::   istept, istepu, istepv   ! time arrays from data files
741      REAL(wp) ::   dayfrac, zxy, zoffsett
742      REAL(wp) ::   zoffsetu, zoffsetv
743      REAL(wp) ::   dayjul0, zdayjulini
744      REAL(wp) ::   zinterval_s, zinterval_e                    ! First and last interval in time axis
745      REAL(wp), DIMENSION(jpbtime)      ::   zstepr             ! REAL time array from data files
746      REAL(wp), DIMENSION(jpbdta,1)     ::   zdta               ! temporary array for data fields
747      CHARACTER(LEN=80), DIMENSION(6)   ::   clfile
748      CHARACTER(LEN=70 )                ::   clunits            ! units attribute of time coordinate
749      !!---------------------------------------------------------------------------
750
751!!gm   add here the same style as in bdy_dta
752!!gm      clearly bdy_dta_bt and bdy_dta  can be combined...   
753!!gm      too many things duplicated in the read of data...   simplification can be done
754
755      ! -------------------- !
756      !    Initialization    !
757      ! -------------------- !
758
759      lect   = .false.           ! If true, read a time record
760
761      ! Some time variables for monthly climatological forcing:
762      ! *******************************************************
763 !!gm  here  use directely daymod variables
764 
765      iman  = INT( raamo ) ! Number of months in a year
766
767      i15 = INT( 2*REAL( nday, wp ) / ( REAL( nmonth_len(nmonth), wp ) + 0.5 ) )
768      ! i15=0 if the current day is in the first half of the month, else i15=1
769
770      imois = nmonth + i15 - 1            ! imois is the first month record
771      IF( imois == 0 ) imois = iman
772
773      ! Time variable for non-climatological forcing:
774      ! *********************************************
775
776      itimer = ((kt-1)-nit000+1)*rdt                      ! current time in seconds for interpolation
777      itimer = itimer + jit*rdt/REAL(nn_baro,wp)      ! in non-climatological case
778
779      IF ( ln_bdy_tides ) THEN
780
781         ! -------------------------------------!
782         ! Update BDY fields with tidal forcing !
783         ! -------------------------------------! 
784
785         CALL tide_update( kt, jit ) 
786 
787      ENDIF
788
789      IF ( ln_bdy_dyn_fla ) THEN
790
791         ! -------------------------------------!
792         ! Update BDY fields with model data    !
793         ! -------------------------------------! 
794
795      !                                                !-------------------!
796      IF( kt == nit000 .and. jit ==2 ) THEN            !  First call only  !
797         !                                             !-------------------!
798         istep_bt(:) = 0
799         nbdy_b_bt    = 0
800         nbdy_a_bt    = 0
801
802         ! Get time information from bdy data file
803         ! ***************************************
804
805        IF(lwp) WRITE(numout,*)
806        IF(lwp) WRITE(numout,*)    'bdy_dta_bt :Initialize unstructured boundary data for barotropic variables.'
807        IF(lwp) WRITE(numout,*)    '~~~~~~~' 
808
809        IF( nbdy_dta == 0 ) THEN
810          IF(lwp) WRITE(numout,*)  'Bdy data are taken from initial conditions'
811
812        ELSEIF (nbdy_dta == 1) THEN
813          IF(lwp) WRITE(numout,*)  'Bdy data are read in netcdf files'
814
815          dayfrac = adatrj  - REAL(itimer,wp)/86400. ! day fraction at time step kt-1
816          dayfrac = dayfrac - INT (dayfrac)          !
817          totime = (nitend-nit000+1)*rdt             ! Total time of the run to verify that all the
818                                                     ! necessary time dumps in file are included
819
820          clfile(4) = filbdy_data_bt_T
821          clfile(5) = filbdy_data_bt_U
822          clfile(6) = filbdy_data_bt_V
823
824          DO igrd = 4,6
825
826            CALL iom_open( clfile(igrd), inum )
827            CALL iom_gettime( inum, zstepr, kntime=ntimes_bdy_bt, cdunits=clunits ) 
828
829            SELECT CASE( igrd )
830               CASE (4) 
831                  numbdyt_bt = inum
832               CASE (5) 
833                  numbdyu_bt = inum
834               CASE (6) 
835                  numbdyv_bt = inum
836            END SELECT
837
838            ! Calculate time offset
839            READ(clunits,7000) iyear0, imonth0, iday0, ihours0, iminutes0, isec0
840            ! Convert time origin in file to julian days
841            isec0 = isec0 + ihours0*60.*60. + iminutes0*60.
842            CALL ymds2ju(iyear0, imonth0, iday0, REAL(isec0, wp), dayjul0)
843            ! Compute model initialization time
844            iyear  = ndastp / 10000
845            imonth = ( ndastp - iyear * 10000 ) / 100
846            iday   = ndastp - iyear * 10000 - imonth * 100
847            isecs  = dayfrac * 86400
848            CALL ymds2ju(iyear, imonth, iday, REAL(isecs, wp) , zdayjulini)
849            ! zoffset from initialization date:
850            zoffset = (dayjul0-zdayjulini)*86400
851            !
852
8537000 FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2)
854
855            !! TO BE DONE... Check consistency between calendar from file
856            !! (available optionally from iom_gettime) and calendar in model
857            !! when calendar in model available outside of IOIPSL.
858
859            ! Check that there are not too many times in the file.
860            IF (ntimes_bdy_bt > jpbtime) CALL ctl_stop( &
861                 'Number of time dumps in bdy file exceed jpbtime parameter', &
862                 'Check file:' // TRIM(clfile(igrd))  )
863
864            ! Check that time array increases (or interp will fail):
865            DO it = 2, ntimes_bdy_bt
866               IF ( zstepr(it-1) >= zstepr(it) ) THEN
867                  CALL ctl_stop('Time array in unstructured boundary data file', &
868                       'does not continuously increase.',               &
869                       'Check file:' // TRIM(clfile(igrd))  )
870                  EXIT
871               END IF
872            END DO
873
874            IF ( .NOT. ln_bdy_clim ) THEN
875               ! Check that times in file span model run time:
876
877               ! Note: the fields may be time means, so we allow nit000 to be before
878               ! first time in the file, provided that it falls inside the meaning
879               ! period of the first field.  Until we can get the meaning period
880               ! from the file, use the interval between fields as a proxy.
881               ! If nit000 is before the first time, use the value at first time
882               ! instead of extrapolating.  This is done by putting time 1 into
883               ! both time levels.
884               ! The same applies to the last time level: see setting of lect below.
885
886               IF ( ntimes_bdy_bt == 1 ) CALL ctl_stop( &
887                    'There is only one time dump in data files', &
888                    'Set ln_bdy_clim=.true. in namelist for constant bdy forcing.' )
889
890               zinterval_s = zstepr(2) - zstepr(1)
891               zinterval_e = zstepr(ntimes_bdy_bt) - zstepr(ntimes_bdy_bt-1)
892
893               IF( zstepr(1) + zoffset > 0 ) THEN
894                     WRITE(ctmp1,*) 'Check file: ', clfile(igrd)
895                     CALL ctl_stop( 'First time dump in bdy file is after model initial time', ctmp1 )
896               END IF
897               IF( zstepr(ntimes_bdy_bt) + zoffset < totime ) THEN
898                     WRITE(ctmp1,*) 'Check file: ', clfile(igrd)
899                     CALL ctl_stop( 'Last time dump in bdy file is before model final time', ctmp1 )
900               END IF
901            END IF ! .NOT. ln_bdy_clim
902
903            IF ( igrd .EQ. 4) THEN
904              ntimes_bdyt = ntimes_bdy_bt
905              zoffsett = zoffset
906              istept(:) = INT( zstepr(:) + zoffset )
907            ELSE IF (igrd .EQ. 5) THEN
908              ntimes_bdyu = ntimes_bdy_bt
909              zoffsetu = zoffset
910              istepu(:) = INT( zstepr(:) + zoffset )
911            ELSE IF (igrd .EQ. 6) THEN
912              ntimes_bdyv = ntimes_bdy_bt
913              zoffsetv = zoffset
914              istepv(:) = INT( zstepr(:) + zoffset )
915            ENDIF
916
917          ENDDO
918
919      ! Verify time consistency between files 
920
921          IF ( ntimes_bdyu /= ntimes_bdyt .OR. ntimes_bdyv /= ntimes_bdyt ) THEN
922             CALL ctl_stop( &
923             'Time axis lengths differ between bdy data files', &
924             'Multiple time frequencies not implemented yet' )
925          ELSE
926            ntimes_bdy_bt = ntimes_bdyt
927          ENDIF
928
929          IF (zoffsetu.NE.zoffsett .OR. zoffsetv.NE.zoffsett) THEN
930            CALL ctl_stop( & 
931            'Bdy data files must have the same time origin', &
932            'Multiple time frequencies not implemented yet'  )
933          ENDIF
934          zoffset = zoffsett
935
936      !! Check that times are the same in the three files... HERE.
937          istep_bt(:) = istept(:)
938
939      ! Check number of time dumps:             
940          IF (ln_bdy_clim) THEN
941            SELECT CASE ( ntimes_bdy_bt )
942            CASE( 1 )
943              IF(lwp) WRITE(numout,*)
944              IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files'
945              IF(lwp) WRITE(numout,*)             
946            CASE( 12 )
947              IF(lwp) WRITE(numout,*)
948              IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files'
949              IF(lwp) WRITE(numout,*) 
950            CASE DEFAULT
951              CALL ctl_stop( &
952                'For climatological boundary forcing (ln_bdy_clim=.true.),',&
953                'bdy data files must contain 1 or 12 time dumps.' )
954            END SELECT
955          ENDIF
956
957      ! Find index of first record to read (before first model time).
958
959          it=1
960          DO WHILE ( ((istep_bt(it+1)) <= 0 ).AND.(it.LE.(ntimes_bdy_bt-1)))
961            it=it+1
962          END DO
963          nbdy_b_bt = it
964
965          IF(lwp) WRITE(numout,*) 'Time offset is ',zoffset
966          IF(lwp) WRITE(numout,*) 'First record to read is ',nbdy_b_bt
967
968        ENDIF ! endif (nbdy_dta == 1)
969
970      ! 1.2  Read first record in file if necessary (ie if nbdy_dta == 1)
971      ! *****************************************************************
972
973        IF ( nbdy_dta == 0) THEN
974          ! boundary data arrays are filled with initial conditions
975          igrd = 5            ! U-points data
976          DO ib = 1, nblen(igrd)             
977            ubtbdy(ib) = un(nbi(ib,igrd), nbj(ib,igrd), 1)
978          END DO
979
980          igrd = 6            ! V-points data
981          DO ib = 1, nblen(igrd)             
982            vbtbdy(ib) = vn(nbi(ib,igrd), nbj(ib,igrd), 1)
983          END DO
984
985          igrd = 4            ! T-points data
986          DO ib = 1, nblen(igrd)             
987            sshbdy(ib) = sshn(nbi(ib,igrd), nbj(ib,igrd))
988          END DO
989
990        ELSEIF (nbdy_dta == 1) THEN
991 
992        ! Set first record in the climatological case:   
993          IF ((ln_bdy_clim).AND.(ntimes_bdy_bt==1)) THEN
994            nbdy_a_bt = 1
995          ELSEIF ((ln_bdy_clim).AND.(ntimes_bdy_bt==iman)) THEN
996            nbdy_b_bt = 0
997            nbdy_a_bt = imois
998          ELSE
999            nbdy_a_bt = nbdy_b_bt
1000          END IF
1001 
1002         ! Open Netcdf files:
1003
1004          CALL iom_open ( filbdy_data_bt_T, numbdyt_bt )
1005          CALL iom_open ( filbdy_data_bt_U, numbdyu_bt )
1006          CALL iom_open ( filbdy_data_bt_V, numbdyv_bt )
1007
1008         ! Read first record:
1009          ipj=1
1010          igrd=4
1011          ipi=nblendta(igrd)
1012
1013          ! ssh
1014          igrd=4
1015          IF ( nblendta(igrd) .le. 0 ) THEN
1016            idvar = iom_varid( numbdyt_bt,'sossheig' )
1017            nblendta(igrd) = iom_file(numbdyt_bt)%dimsz(1,idvar)
1018          ENDIF
1019          WRITE(numout,*) 'Dim size for sossheig is ',nblendta(igrd)
1020          ipi=nblendta(igrd)
1021
1022          CALL iom_get ( numbdyt_bt, jpdom_unknown,'sossheig',zdta(1:ipi,1:ipj),nbdy_a_bt )
1023
1024          DO ib=1, nblen(igrd)
1025            sshbdydta(ib,2) =  zdta(nbmap(ib,igrd),1)
1026          END DO
1027 
1028          ! u-velocity
1029          igrd=5
1030          IF ( nblendta(igrd) .le. 0 ) THEN
1031            idvar = iom_varid( numbdyu_bt,'vobtcrtx' )
1032            nblendta(igrd) = iom_file(numbdyu_bt)%dimsz(1,idvar)
1033          ENDIF
1034          WRITE(numout,*) 'Dim size for vobtcrtx is ',nblendta(igrd)
1035          ipi=nblendta(igrd)
1036
1037          CALL iom_get ( numbdyu_bt, jpdom_unknown,'vobtcrtx',zdta(1:ipi,1:ipj),nbdy_a_bt )
1038
1039          DO ib=1, nblen(igrd)
1040            ubtbdydta(ib,2) =  zdta(nbmap(ib,igrd),1)
1041          END DO
1042
1043          ! v-velocity
1044          igrd=6
1045          IF ( nblendta(igrd) .le. 0 ) THEN
1046            idvar = iom_varid( numbdyv_bt,'vobtcrty' )
1047            nblendta(igrd) = iom_file(numbdyv_bt)%dimsz(1,idvar)
1048          ENDIF
1049          WRITE(numout,*) 'Dim size for vobtcrty is ',nblendta(igrd)
1050          ipi=nblendta(igrd)
1051
1052          CALL iom_get ( numbdyv_bt, jpdom_unknown,'vobtcrty',zdta(1:ipi,1:ipj),nbdy_a_bt )
1053
1054          DO ib=1, nblen(igrd)
1055            vbtbdydta(ib,2) =  zdta(nbmap(ib,igrd),1)
1056          END DO
1057
1058        END IF
1059 
1060        ! In the case of constant boundary forcing fill bdy arrays once for all
1061        IF ((ln_bdy_clim).AND.(ntimes_bdy_bt==1)) THEN
1062
1063          ubtbdy  (:) = ubtbdydta  (:,2)
1064          vbtbdy  (:) = vbtbdydta  (:,2)
1065          sshbdy  (:) = sshbdydta  (:,2)
1066
1067          CALL iom_close( numbdyt_bt )
1068          CALL iom_close( numbdyu_bt )
1069          CALL iom_close( numbdyv_bt )
1070
1071        END IF
1072
1073      ENDIF ! End if nit000
1074
1075      ! -------------------- !
1076      ! 2. At each time step !
1077      ! -------------------- !
1078
1079      IF ((nbdy_dta==1).AND.(ntimes_bdy_bt>1)) THEN 
1080
1081      ! 2.1 Read one more record if necessary
1082      !**************************************
1083
1084        IF ( (ln_bdy_clim).AND.(imois/=nbdy_b_bt) ) THEN ! remember that nbdy_b_bt=0 for kt=nit000
1085         nbdy_b_bt = imois
1086         nbdy_a_bt = imois+1
1087         nbdy_b_bt = MOD( nbdy_b_bt, iman )
1088         IF( nbdy_b_bt == 0 ) nbdy_b_bt = iman
1089         nbdy_a_bt = MOD( nbdy_a_bt, iman )
1090         IF( nbdy_a_bt == 0 ) nbdy_a_bt = iman
1091         lect=.true.
1092
1093        ELSEIF ((.NOT.ln_bdy_clim).AND.(itimer >= istep_bt(nbdy_a_bt))) THEN
1094          nbdy_b_bt=nbdy_a_bt
1095          nbdy_a_bt=nbdy_a_bt+1
1096          lect=.true.
1097        END IF
1098         
1099        IF (lect) THEN
1100
1101        ! Swap arrays
1102          sshbdydta(:,1) =  sshbdydta(:,2)
1103          ubtbdydta(:,1) =  ubtbdydta(:,2)
1104          vbtbdydta(:,1) =  vbtbdydta(:,2)
1105 
1106        ! read another set
1107
1108          ipj=1
1109          ipk=jpk
1110          igrd=4
1111          ipi=nblendta(igrd)
1112
1113         
1114          ! ssh
1115          igrd=4
1116          ipi=nblendta(igrd)
1117
1118          CALL iom_get ( numbdyt_bt, jpdom_unknown,'sossheig',zdta(1:ipi,1:ipj),nbdy_a_bt )
1119
1120          DO ib=1, nblen(igrd)
1121            sshbdydta(ib,2) =  zdta(nbmap(ib,igrd),1)
1122          END DO
1123
1124          ! u-velocity
1125          igrd=5
1126          ipi=nblendta(igrd)
1127
1128          CALL iom_get ( numbdyu_bt, jpdom_unknown,'vobtcrtx',zdta(1:ipi,1:ipj),nbdy_a_bt )
1129
1130          DO ib=1, nblen(igrd)
1131            ubtbdydta(ib,2) =  zdta(nbmap(ib,igrd),1)
1132          END DO
1133
1134          ! v-velocity
1135          igrd=6
1136          ipi=nblendta(igrd)
1137
1138          CALL iom_get ( numbdyv_bt, jpdom_unknown,'vobtcrty',zdta(1:ipi,1:ipj),nbdy_a_bt )
1139
1140          DO ib=1, nblen(igrd)
1141            vbtbdydta(ib,2) =  zdta(nbmap(ib,igrd),1)
1142          END DO
1143
1144
1145         IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy_b_bt ',nbdy_b_bt
1146         IF(lwp) WRITE(numout,*) '~~~~~~~~  last  record file used nbdy_a_bt ',nbdy_a_bt
1147         IF (.NOT.ln_bdy_clim) THEN
1148           IF(lwp) WRITE(numout,*) 'first  record time (s): ', istep_bt(nbdy_b_bt)
1149           IF(lwp) WRITE(numout,*) 'model time (s)        : ', itimer
1150           IF(lwp) WRITE(numout,*) 'second record time (s): ', istep_bt(nbdy_a_bt)
1151         ENDIF
1152        END IF ! end lect=.true.
1153
1154
1155      ! 2.2   Interpolate linearly:
1156      ! ***************************
1157   
1158        IF (ln_bdy_clim) THEN
1159          zxy = REAL( nday, wp ) / REAL( nmonth_len(nbdy_b_bt), wp ) + 0.5 - i15
1160        ELSE         
1161          zxy = REAL(istep_bt(nbdy_b_bt)-itimer, wp) / REAL(istep_bt(nbdy_b_bt)-istep_bt(nbdy_a_bt), wp)
1162        END IF
1163
1164          igrd=4
1165          DO ib=1, nblen(igrd)
1166            sshbdy(ib) = zxy      * sshbdydta(ib,2) + &
1167                       (1.-zxy) * sshbdydta(ib,1)   
1168          END DO
1169
1170          igrd=5
1171          DO ib=1, nblen(igrd)
1172            ubtbdy(ib) = zxy      * ubtbdydta(ib,2) + &
1173                         (1.-zxy) * ubtbdydta(ib,1)   
1174          END DO
1175
1176          igrd=6
1177          DO ib=1, nblen(igrd)
1178            vbtbdy(ib) = zxy      * vbtbdydta(ib,2) + &
1179                         (1.-zxy) * vbtbdydta(ib,1)   
1180          END DO
1181
1182
1183      END IF !end if ((nbdy_dta==1).AND.(ntimes_bdy_bt>1))
1184   
1185      ! ------------------- !
1186      ! Last call kt=nitend !
1187      ! ------------------- !
1188
1189      ! Closing of the 3 files
1190      IF( kt == nitend   .and. jit == icycl ) THEN
1191          CALL iom_close( numbdyt_bt )
1192          CALL iom_close( numbdyu_bt )
1193          CALL iom_close( numbdyv_bt )
1194      ENDIF
1195
1196      ENDIF ! ln_bdy_dyn_frs
1197
1198      END SUBROUTINE bdy_dta_bt
1199
1200
1201#else
1202   !!----------------------------------------------------------------------
1203   !!   Dummy module                   NO Unstruct Open Boundary Conditions
1204   !!----------------------------------------------------------------------
1205CONTAINS
1206   SUBROUTINE bdy_dta( kt )              ! Empty routine
1207      WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt
1208   END SUBROUTINE bdy_dta
1209   SUBROUTINE bdy_dta_bt( kt, kit, icycle )      ! Empty routine
1210      WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt, kit
1211   END SUBROUTINE bdy_dta_bt
1212#endif
1213
1214   !!==============================================================================
1215END MODULE bdydta
Note: See TracBrowser for help on using the repository browser.