source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90 @ 4463

Last change on this file since 4463 was 4463, checked in by trackstand2, 7 years ago

Use of mbkmax in bdy_dta_fra

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