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/DEV_r1986_BDY_updates/NEMO/OPA_SRC/BDY – NEMO

source: branches/DEV_r1986_BDY_updates/NEMO/OPA_SRC/BDY/bdydta.F90 @ 2094

Last change on this file since 2094 was 2094, checked in by davestorkey, 14 years ago

Correction to bug fix.

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