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.
fldread.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 5749

Last change on this file since 5749 was 5749, checked in by jpaul, 9 years ago

bug fix: see ticket #1598

  • Property svn:keywords set to Id
File size: 85.4 KB
Line 
1MODULE fldread
2   !!======================================================================
3   !!                       ***  MODULE  fldread  ***
4   !! Ocean forcing:  read input field for surface boundary condition
5   !!=====================================================================
6   !! History :  2.0  !  06-2006  (S. Masson, G. Madec) Original code
7   !!                 !  05-2008  (S. Alderson) Modified for Interpolation in memory
8   !!                 !                         from input grid to model grid
9   !!                 !  10-2013  (D. Delrosso, P. Oddo) implement suppression of
10   !!                 !                         land point prior to interpolation
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   fld_read      : read input fields used for the computation of the
15   !!                   surface boundary condition
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and tracers
18   USE dom_oce         ! ocean space and time domain
19   USE phycst          ! ???
20   USE in_out_manager  ! I/O manager
21   USE iom             ! I/O manager library
22   USE geo2ocean       ! for vector rotation on to model grid
23   USE lib_mpp         ! MPP library
24   USE wrk_nemo        ! work arrays
25   USE lbclnk          ! ocean lateral boundary conditions (C1D case)
26   USE ioipsl, ONLY :   ymds2ju, ju2ymds   ! for calendar
27   USE sbc_oce
28   
29   IMPLICIT NONE
30   PRIVATE   
31 
32   PUBLIC   fld_map    ! routine called by tides_init
33   PUBLIC   fld_read, fld_fill   ! called by sbc... modules
34
35   TYPE, PUBLIC ::   FLD_N      !: Namelist field informations
36      CHARACTER(len = 256) ::   clname      ! generic name of the NetCDF flux file
37      REAL(wp)             ::   nfreqh      ! frequency of each flux file
38      CHARACTER(len = 34)  ::   clvar       ! generic name of the variable in the NetCDF flux file
39      LOGICAL              ::   ln_tint     ! time interpolation or not (T/F)
40      LOGICAL              ::   ln_clim     ! climatology or not (T/F)
41      CHARACTER(len = 8)   ::   cltype      ! type of data file 'daily', 'monthly' or yearly'
42      CHARACTER(len = 256) ::   wname       ! generic name of a NetCDF weights file to be used, blank if not
43      CHARACTER(len = 34)  ::   vcomp       ! symbolic component name if a vector that needs rotation
44      !                                     ! a string starting with "U" or "V" for each component   
45      !                                     ! chars 2 onwards identify which components go together 
46      CHARACTER(len = 34)  ::   lname       ! generic name of a NetCDF land/sea mask file to be used, blank if not
47      !                                     ! 0=sea 1=land
48   END TYPE FLD_N
49
50   TYPE, PUBLIC ::   FLD        !: Input field related variables
51      CHARACTER(len = 256)            ::   clrootname   ! generic name of the NetCDF file
52      CHARACTER(len = 256)            ::   clname       ! current name of the NetCDF file
53      REAL(wp)                        ::   nfreqh       ! frequency of each flux file
54      CHARACTER(len = 34)             ::   clvar        ! generic name of the variable in the NetCDF flux file
55      LOGICAL                         ::   ln_tint      ! time interpolation or not (T/F)
56      LOGICAL                         ::   ln_clim      ! climatology or not (T/F)
57      CHARACTER(len = 8)              ::   cltype       ! type of data file 'daily', 'monthly' or yearly'
58      INTEGER                         ::   num          ! iom id of the jpfld files to be read
59      INTEGER , DIMENSION(2)          ::   nrec_b       ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year)
60      INTEGER , DIMENSION(2)          ::   nrec_a       ! after  record (1: index, 2: second since Jan. 1st 00h of nit000 year)
61      REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:  ) ::   fnow       ! input fields interpolated to now time step
62      REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) ::   fdta       ! 2 consecutive record of input fields
63      CHARACTER(len = 256)            ::   wgtname      ! current name of the NetCDF weight file acting as a key
64      !                                                 ! into the WGTLIST structure
65      CHARACTER(len = 34)             ::   vcomp        ! symbolic name for a vector component that needs rotation
66      LOGICAL, DIMENSION(2)           ::   rotn         ! flag to indicate whether before/after field has been rotated
67      INTEGER                         ::   nreclast     ! last record to be read in the current file
68      CHARACTER(len = 256)            ::   lsmname      ! current name of the NetCDF mask file acting as a key
69   END TYPE FLD
70
71   TYPE, PUBLIC ::   MAP_POINTER      !: Map from input data file to local domain
72      INTEGER, POINTER, DIMENSION(:)  ::  ptr           ! Array of integer pointers to 1D arrays
73      LOGICAL                         ::  ll_unstruc    ! Unstructured (T) or structured (F) boundary data file
74   END TYPE MAP_POINTER
75
76!$AGRIF_DO_NOT_TREAT
77
78   !! keep list of all weights variables so they're only read in once
79   !! need to add AGRIF directives not to process this structure
80   !! also need to force wgtname to include AGRIF nest number
81   TYPE         ::   WGT        !: Input weights related variables
82      CHARACTER(len = 256)                    ::   wgtname      ! current name of the NetCDF weight file
83      INTEGER , DIMENSION(2)                  ::   ddims        ! shape of input grid
84      INTEGER , DIMENSION(2)                  ::   botleft      ! top left corner of box in input grid containing
85      !                                                         ! current processor grid
86      INTEGER , DIMENSION(2)                  ::   topright     ! top right corner of box
87      INTEGER                                 ::   jpiwgt       ! width of box on input grid
88      INTEGER                                 ::   jpjwgt       ! height of box on input grid
89      INTEGER                                 ::   numwgt       ! number of weights (4=bilinear, 16=bicubic)
90      INTEGER                                 ::   nestid       ! for agrif, keep track of nest we're in
91      INTEGER                                 ::   overlap      ! =0 when cyclic grid has no overlapping EW columns
92      !                                                         ! =>1 when they have one or more overlapping columns     
93      !                                                         ! =-1 not cyclic
94      LOGICAL                                 ::   cyclic       ! east-west cyclic or not
95      INTEGER,  DIMENSION(:,:,:), POINTER     ::   data_jpi     ! array of source integers
96      INTEGER,  DIMENSION(:,:,:), POINTER     ::   data_jpj     ! array of source integers
97      REAL(wp), DIMENSION(:,:,:), POINTER     ::   data_wgt     ! array of weights on model grid
98      REAL(wp), DIMENSION(:,:,:), POINTER     ::   fly_dta      ! array of values on input grid
99      REAL(wp), DIMENSION(:,:,:), POINTER     ::   col          ! temporary array for reading in columns
100   END TYPE WGT
101
102   INTEGER,     PARAMETER             ::   tot_wgts = 10
103   TYPE( WGT ), DIMENSION(tot_wgts)   ::   ref_wgts     ! array of wgts
104   INTEGER                            ::   nxt_wgt = 1  ! point to next available space in ref_wgts array
105   REAL(wp), PARAMETER                ::   undeff_lsm = -999.00_wp
106
107!$AGRIF_END_DO_NOT_TREAT
108
109   !!----------------------------------------------------------------------
110   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
111   !! $Id$
112   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
113   !!----------------------------------------------------------------------
114CONTAINS
115
116   SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset )
117      !!---------------------------------------------------------------------
118      !!                    ***  ROUTINE fld_read  ***
119      !!                   
120      !! ** Purpose :   provide at each time step the surface ocean fluxes
121      !!                (momentum, heat, freshwater and runoff)
122      !!
123      !! ** Method  :   READ each input fields in NetCDF files using IOM
124      !!      and intepolate it to the model time-step.
125      !!         Several assumptions are made on the input file:
126      !!      blahblahblah....
127      !!----------------------------------------------------------------------
128      INTEGER  , INTENT(in   )               ::   kt        ! ocean time step
129      INTEGER  , INTENT(in   )               ::   kn_fsbc   ! sbc computation period (in time step)
130      TYPE(FLD), INTENT(inout), DIMENSION(:) ::   sd        ! input field related variables
131      TYPE(MAP_POINTER),INTENT(in), OPTIONAL, DIMENSION(:) ::   map   ! global-to-local mapping indices
132      INTEGER  , INTENT(in   ), OPTIONAL     ::   kit       ! subcycle timestep for timesplitting option
133      INTEGER  , INTENT(in   ), OPTIONAL     ::   kt_offset ! provide fields at time other than "now"
134                                                            !   kt_offset = -1 => fields at "before" time level
135                                                            !   kt_offset = +1 => fields at "after"  time level
136                                                            !   etc.
137      !!
138      INTEGER  ::   itmp       ! temporary variable
139      INTEGER  ::   imf        ! size of the structure sd
140      INTEGER  ::   jf         ! dummy indices
141      INTEGER  ::   isecend    ! number of second since Jan. 1st 00h of nit000 year at nitend
142      INTEGER  ::   isecsbc    ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step
143      INTEGER  ::   it_offset  ! local time offset variable
144      LOGICAL  ::   llnxtyr    ! open next year  file?
145      LOGICAL  ::   llnxtmth   ! open next month file?
146      LOGICAL  ::   llstop     ! stop is the file does not exist
147      LOGICAL  ::   ll_firstcall ! true if this is the first call to fld_read for this set of fields
148      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation
149      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation
150      CHARACTER(LEN=1000) ::   clfmt   ! write format
151      TYPE(MAP_POINTER) ::   imap   ! global-to-local mapping indices
152      !!---------------------------------------------------------------------
153      ll_firstcall = kt == nit000
154      IF( PRESENT(kit) )   ll_firstcall = ll_firstcall .and. kit == 1
155
156      IF ( nn_components == jp_iam_sas ) THEN   ;   it_offset = nn_fsbc
157      ELSE                                      ;   it_offset = 0
158      ENDIF
159      IF( PRESENT(kt_offset) )   it_offset = kt_offset
160
161      imap%ptr => NULL()
162
163      ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
164      IF( present(kit) ) THEN   ! ignore kn_fsbc in this case
165         isecsbc = nsec_year + nsec1jan000 + (kit+it_offset)*NINT( rdt/REAL(nn_baro,wp) )
166      ELSE                      ! middle of sbc time step
167         isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdttra(1)) + it_offset * NINT(rdttra(1))
168      ENDIF
169      imf = SIZE( sd )
170      !
171      IF( ll_firstcall ) THEN                      ! initialization
172         DO jf = 1, imf 
173            IF( PRESENT(map) ) imap = map(jf)
174            CALL fld_init( kn_fsbc, sd(jf), imap )  ! read each before field (put them in after as they will be swapped)
175         END DO
176         IF( lwp ) CALL wgt_print()                ! control print
177      ENDIF
178      !                                            ! ====================================== !
179      IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! update field at each kn_fsbc time-step !
180         !                                         ! ====================================== !
181         !
182         DO jf = 1, imf                            ! ---   loop over field   --- !
183           
184            IF( isecsbc > sd(jf)%nrec_a(2) .OR. ll_firstcall ) THEN    ! read/update the after data?
185
186               IF( PRESENT(map) )   imap = map(jf)   ! temporary definition of map
187
188               sd(jf)%nrec_b(:) = sd(jf)%nrec_a(:)                                  ! swap before record informations
189               sd(jf)%rotn(1) = sd(jf)%rotn(2)                                      ! swap before rotate informations
190               IF( sd(jf)%ln_tint )   sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)   ! swap before record field
191
192               CALL fld_rec( kn_fsbc, sd(jf), kt_offset = it_offset, kit = kit )    ! update after record informations
193
194               ! if kn_fsbc*rdttra is larger than nfreqh (which is kind of odd),
195               ! it is possible that the before value is no more the good one... we have to re-read it
196               ! if before is not the last record of the file currently opened and after is the first record to be read
197               ! in a new file which means after = 1 (the file to be opened corresponds to the current time)
198               ! or after = nreclast + 1 (the file to be opened corresponds to a future time step)
199               IF( .NOT. ll_firstcall .AND. sd(jf)%ln_tint .AND. sd(jf)%nrec_b(1) /= sd(jf)%nreclast &
200                  &                   .AND. MOD( sd(jf)%nrec_a(1), sd(jf)%nreclast ) == 1 ) THEN
201                  itmp = sd(jf)%nrec_a(1)                       ! temporary storage
202                  sd(jf)%nrec_a(1) = sd(jf)%nreclast            ! read the last record of the file currently opened
203                  CALL fld_get( sd(jf), imap )                  ! read after data
204                  sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)   ! re-swap before record field
205                  sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1)           ! update before record informations
206                  sd(jf)%nrec_b(2) = sd(jf)%nrec_a(2) - NINT( sd(jf)%nfreqh * 3600 )  ! assume freq to be in hours in this case
207                  sd(jf)%rotn(1)   = sd(jf)%rotn(2)             ! update before rotate informations
208                  sd(jf)%nrec_a(1) = itmp                       ! move back to after record
209               ENDIF
210
211               CALL fld_clopn( sd(jf) )   ! Do we need to open a new year/month/week/day file?
212               
213               IF( sd(jf)%ln_tint ) THEN
214                 
215                  ! if kn_fsbc*rdttra is larger than nfreqh (which is kind of odd),
216                  ! it is possible that the before value is no more the good one... we have to re-read it
217                  ! if before record is not just just before the after record...
218                  IF( .NOT. ll_firstcall .AND. MOD( sd(jf)%nrec_a(1), sd(jf)%nreclast ) /= 1 &
219                     &                   .AND. sd(jf)%nrec_b(1) /= sd(jf)%nrec_a(1) - 1 ) THEN   
220                     sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) - 1       ! move back to before record
221                     CALL fld_get( sd(jf), imap )                  ! read after data
222                     sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)   ! re-swap before record field
223                     sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1)           ! update before record informations
224                     sd(jf)%nrec_b(2) = sd(jf)%nrec_a(2) - NINT( sd(jf)%nfreqh * 3600 )  ! assume freq to be in hours in this case
225                     sd(jf)%rotn(1)   = sd(jf)%rotn(2)             ! update before rotate informations
226                     sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) + 1       ! move back to after record
227                  ENDIF
228
229                  ! do we have to change the year/month/week/day of the forcing field??
230                  ! if we do time interpolation we will need to open next year/month/week/day file before the end of the current
231                  ! one. If so, we are still before the end of the year/month/week/day when calling fld_rec so sd(jf)%nrec_a(1)
232                  ! will be larger than the record number that should be read for current year/month/week/day
233                  ! do we need next file data?
234                  IF( sd(jf)%nrec_a(1) > sd(jf)%nreclast ) THEN
235                     
236                     sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) - sd(jf)%nreclast   !
237                     
238                     IF( .NOT. ( sd(jf)%ln_clim .AND. sd(jf)%cltype == 'yearly' ) ) THEN   ! close/open the current/new file
239                       
240                        llnxtmth = sd(jf)%cltype == 'monthly' .OR. nday == nmonth_len(nmonth)      ! open next month file?
241                        llnxtyr  = sd(jf)%cltype == 'yearly'  .OR. (nmonth == 12 .AND. llnxtmth)   ! open next year  file?
242
243                        ! if the run finishes at the end of the current year/month/week/day, we will allow next
244                        ! year/month/week/day file to be not present. If the run continue further than the current
245                        ! year/month/week/day, next year/month/week/day file must exist
246                        isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdttra(1))   ! second at the end of the run
247                        llstop = isecend > sd(jf)%nrec_a(2)                                   ! read more than 1 record of next year
248                        ! we suppose that the date of next file is next day (should be ok even for weekly files...)
249                        CALL fld_clopn( sd(jf), nyear  + COUNT((/llnxtyr /))                                           ,         &
250                           &                    nmonth + COUNT((/llnxtmth/)) - 12                 * COUNT((/llnxtyr /)),         &
251                           &                    nday   + 1                   - nmonth_len(nmonth) * COUNT((/llnxtmth/)), llstop )
252
253                        IF( sd(jf)%num <= 0 .AND. .NOT. llstop ) THEN    ! next year file does not exist
254                           CALL ctl_warn('next year/month/week/day file: '//TRIM(sd(jf)%clname)//     &
255                              &     ' not present -> back to current year/month/day')
256                           CALL fld_clopn( sd(jf) )       ! back to the current year/month/day
257                           sd(jf)%nrec_a(1) = sd(jf)%nreclast     ! force to read the last record in the current year file
258                        ENDIF
259                       
260                     ENDIF
261                  ENDIF   ! open need next file?
262                 
263               ENDIF   ! temporal interpolation?
264
265               ! read after data
266               CALL fld_get( sd(jf), imap )
267
268            ENDIF   ! read new data?
269         END DO                                    ! --- end loop over field --- !
270
271         CALL fld_rot( kt, sd )                    ! rotate vector before/now/after fields if needed
272
273         DO jf = 1, imf                            ! ---   loop over field   --- !
274            !
275            IF( sd(jf)%ln_tint ) THEN              ! temporal interpolation
276               IF(lwp .AND. kt - nit000 <= 100 ) THEN
277                  clfmt = "('fld_read: var ', a, ' kt = ', i8, ' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
278                     &    "', records b/a: ', i6.4, '/', i6.4, ' (days ', f9.4,'/', f9.4, ')')"
279                  WRITE(numout, clfmt)  TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,   &           
280                     & sd(jf)%nrec_b(1), sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday
281                  WRITE(numout, *) 'it_offset is : ',it_offset
282               ENDIF
283               ! temporal interpolation weights
284               ztinta =  REAL( isecsbc - sd(jf)%nrec_b(2), wp ) / REAL( sd(jf)%nrec_a(2) - sd(jf)%nrec_b(2), wp )
285               ztintb =  1. - ztinta
286!CDIR COLLAPSE
287               sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:,1) + ztinta * sd(jf)%fdta(:,:,:,2)
288            ELSE   ! nothing to do...
289               IF(lwp .AND. kt - nit000 <= 100 ) THEN
290                  clfmt = "('fld_read: var ', a, ' kt = ', i8,' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
291                     &    "', record: ', i6.4, ' (days ', f9.4, ' <-> ', f9.4, ')')"
292                  WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,    &
293                     &                 sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday
294               ENDIF
295            ENDIF
296            !
297            IF( kt == nitend - kn_fsbc + 1 )   CALL iom_close( sd(jf)%num )   ! Close the input files
298
299         END DO                                    ! --- end loop over field --- !
300         !
301         !                                         ! ====================================== !
302      ENDIF                                        ! update field at each kn_fsbc time-step !
303      !                                            ! ====================================== !
304      !
305   END SUBROUTINE fld_read
306
307
308   SUBROUTINE fld_init( kn_fsbc, sdjf, map )
309      !!---------------------------------------------------------------------
310      !!                    ***  ROUTINE fld_init  ***
311      !!
312      !! ** Purpose :  - first call to fld_rec to define before values
313      !!               - if time interpolation, read before data
314      !!----------------------------------------------------------------------
315      INTEGER  , INTENT(in   ) ::   kn_fsbc   ! sbc computation period (in time step)
316      TYPE(FLD), INTENT(inout) ::   sdjf      ! input field related variables
317      TYPE(MAP_POINTER),INTENT(in) ::   map   ! global-to-local mapping indices
318      !!
319      LOGICAL :: llprevyr              ! are we reading previous year  file?
320      LOGICAL :: llprevmth             ! are we reading previous month file?
321      LOGICAL :: llprevweek            ! are we reading previous week  file?
322      LOGICAL :: llprevday             ! are we reading previous day   file?
323      LOGICAL :: llprev                ! llprevyr .OR. llprevmth .OR. llprevweek .OR. llprevday
324      INTEGER :: idvar                 ! variable id
325      INTEGER :: inrec                 ! number of record existing for this variable
326      INTEGER :: iyear, imonth, iday   ! first day of the current file in yyyy mm dd
327      INTEGER :: isec_week             ! number of seconds since start of the weekly file
328      CHARACTER(LEN=1000) ::   clfmt   ! write format
329      !!---------------------------------------------------------------------
330      llprevyr   = .FALSE.
331      llprevmth  = .FALSE.
332      llprevweek = .FALSE.
333      llprevday  = .FALSE.
334      isec_week  = 0
335           
336      ! define record informations
337      CALL fld_rec( kn_fsbc, sdjf, ldbefore = .TRUE. )  ! return before values in sdjf%nrec_a (as we will swap it later)
338
339      ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
340
341      IF( sdjf%ln_tint ) THEN ! we need to read the previous record and we will put it in the current record structure
342
343         IF( sdjf%nrec_a(1) == 0  ) THEN   ! we redefine record sdjf%nrec_a(1) with the last record of previous year file
344            IF    ( sdjf%nfreqh == -12 ) THEN   ! yearly mean
345               IF( sdjf%cltype == 'yearly' ) THEN             ! yearly file
346                  sdjf%nrec_a(1) = 1                                                       ! force to read the unique record
347                  llprevyr  = .NOT. sdjf%ln_clim                                           ! use previous year  file?
348               ELSE
349                  CALL ctl_stop( "fld_init: yearly mean file must be in a yearly type of file: "//TRIM(sdjf%clrootname) )
350               ENDIF
351            ELSEIF( sdjf%nfreqh ==  -1 ) THEN   ! monthly mean
352               IF( sdjf%cltype == 'monthly' ) THEN            ! monthly file
353                  sdjf%nrec_a(1) = 1                                                       ! force to read the unique record
354                  llprevmth = .TRUE.                                                       ! use previous month file?
355                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file?
356               ELSE                                           ! yearly file
357                  sdjf%nrec_a(1) = 12                                                      ! force to read december mean
358                  llprevyr = .NOT. sdjf%ln_clim                                            ! use previous year  file?
359               ENDIF
360            ELSE                                ! higher frequency mean (in hours)
361               IF    ( sdjf%cltype      == 'monthly' ) THEN   ! monthly file
362                  sdjf%nrec_a(1) = NINT( 24 * nmonth_len(nmonth-1) / sdjf%nfreqh )         ! last record of previous month
363                  llprevmth = .TRUE.                                                       ! use previous month file?
364                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file?
365               ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ! weekly file
366                  llprevweek = .TRUE.                                                      ! use previous week  file?
367                  sdjf%nrec_a(1) = NINT( 24 * 7 / sdjf%nfreqh )                            ! last record of previous week
368                  isec_week = NINT(rday) * 7                                               ! add a shift toward previous week
369               ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ! daily file
370                  sdjf%nrec_a(1) = NINT( 24 / sdjf%nfreqh )                                ! last record of previous day
371                  llprevday = .TRUE.                                                       ! use previous day   file?
372                  llprevmth = llprevday .AND. nday   == 1                                  ! use previous month file?
373                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file?
374               ELSE                                           ! yearly file
375                  sdjf%nrec_a(1) = NINT( 24 * nyear_len(0) / sdjf%nfreqh )                 ! last record of previous year
376                  llprevyr = .NOT. sdjf%ln_clim                                            ! use previous year  file?
377               ENDIF
378            ENDIF
379         ENDIF
380         !
381         IF ( sdjf%cltype(1:4) == 'week' ) THEN
382            isec_week = isec_week + ksec_week( sdjf%cltype(6:8) )   ! second since the beginning of the week
383            llprevmth = isec_week > nsec_month                      ! longer time since the beginning of the week than the month
384            llprevyr  = llprevmth .AND. nmonth == 1
385         ENDIF
386         llprev = llprevyr .OR. llprevmth .OR. llprevweek .OR. llprevday
387         !
388         iyear  = nyear  - COUNT((/llprevyr /))
389         imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
390         iday   = nday   - COUNT((/llprevday/)) + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
391         !
392         CALL fld_clopn( sdjf, iyear, imonth, iday, .NOT. llprev )
393
394         ! if previous year/month/day file does not exist, we switch to the current year/month/day
395         IF( llprev .AND. sdjf%num <= 0 ) THEN
396            CALL ctl_warn( 'previous year/month/week/day file: '//TRIM(sdjf%clrootname)//   &
397               &           ' not present -> back to current year/month/week/day' )
398            ! we force to read the first record of the current year/month/day instead of last record of previous year/month/day
399            llprev = .FALSE.
400            sdjf%nrec_a(1) = 1
401            CALL fld_clopn( sdjf )
402         ENDIF
403         
404         IF( llprev ) THEN   ! check if the record sdjf%nrec_a(1) exists in the file
405            idvar = iom_varid( sdjf%num, sdjf%clvar )                                        ! id of the variable sdjf%clvar
406            IF( idvar <= 0 )   RETURN
407            inrec = iom_file( sdjf%num )%dimsz( iom_file( sdjf%num )%ndims(idvar), idvar )   ! size of the last dim of idvar
408            sdjf%nrec_a(1) = MIN( sdjf%nrec_a(1), inrec )   ! make sure we select an existing record
409         ENDIF
410
411         ! read before data in after arrays(as we will swap it later)
412         CALL fld_get( sdjf, map )
413
414         clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')"
415         IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_a(1), REAL(sdjf%nrec_a(2),wp)/rday
416
417      ENDIF
418      !
419   END SUBROUTINE fld_init
420
421
422   SUBROUTINE fld_rec( kn_fsbc, sdjf, ldbefore, kit, kt_offset )
423      !!---------------------------------------------------------------------
424      !!                    ***  ROUTINE fld_rec  ***
425      !!
426      !! ** Purpose : Compute
427      !!              if sdjf%ln_tint = .TRUE.
428      !!                  nrec_a: record number and its time (nrec_b is obtained from nrec_a when swapping)
429      !!              if sdjf%ln_tint = .FALSE.
430      !!                  nrec_a(1): record number
431      !!                  nrec_b(2) and nrec_a(2): time of the beginning and end of the record (for print only)
432      !!----------------------------------------------------------------------
433      INTEGER  , INTENT(in   )           ::   kn_fsbc   ! sbc computation period (in time step)
434      TYPE(FLD), INTENT(inout)           ::   sdjf      ! input field related variables
435      LOGICAL  , INTENT(in   ), OPTIONAL ::   ldbefore  ! sent back before record values (default = .FALSE.)
436      INTEGER  , INTENT(in   ), OPTIONAL ::   kit       ! index of barotropic subcycle
437                                                        ! used only if sdjf%ln_tint = .TRUE.
438      INTEGER  , INTENT(in   ), OPTIONAL ::   kt_offset ! Offset of required time level compared to "now"
439                                                        !   time level in units of time steps.
440      !!
441      LOGICAL  ::   llbefore    ! local definition of ldbefore
442      INTEGER  ::   iendrec     ! end of this record (in seconds)
443      INTEGER  ::   imth        ! month number
444      INTEGER  ::   ifreq_sec   ! frequency mean (in seconds)
445      INTEGER  ::   isec_week   ! number of seconds since the start of the weekly file
446      INTEGER  ::   it_offset   ! local time offset variable
447      REAL(wp) ::   ztmp        ! temporary variable
448      !!----------------------------------------------------------------------
449      !
450      ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
451      !
452      IF( PRESENT(ldbefore) ) THEN   ;   llbefore = ldbefore .AND. sdjf%ln_tint   ! needed only if sdjf%ln_tint = .TRUE.
453      ELSE                           ;   llbefore = .FALSE.
454      ENDIF
455      !
456      IF ( nn_components == jp_iam_sas ) THEN   ;   it_offset = nn_fsbc
457      ELSE                                      ;   it_offset = 0
458      ENDIF
459      IF( PRESENT(kt_offset) )   it_offset = kt_offset
460      IF( PRESENT(kit) ) THEN   ;   it_offset = ( kit + it_offset ) * NINT( rdt/REAL(nn_baro,wp) )
461      ELSE                      ;   it_offset =         it_offset   * NINT(       rdttra(1)      )
462      ENDIF
463      !
464      !                                      ! =========== !
465      IF    ( sdjf%nfreqh == -12 ) THEN      ! yearly mean
466         !                                   ! =========== !
467         !
468         IF( sdjf%ln_tint ) THEN                 ! time interpolation, shift by 1/2 record
469            !
470            !                  INT( ztmp )
471            !                     /|\
472            !                    1 |    *----
473            !                    0 |----(             
474            !                      |----+----|--> time
475            !                      0   /|\   1   (nday/nyear_len(1))
476            !                           |   
477            !                           |   
478            !       forcing record :    1
479            !                           
480            ztmp = REAL( nsec_year, wp ) / ( REAL( nyear_len(1), wp ) * rday ) + 0.5 &
481           &       + REAL( it_offset, wp ) / ( REAL( nyear_len(1), wp ) * rday )
482            sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
483            ! swap at the middle of the year
484            IF( llbefore ) THEN   ;   sdjf%nrec_a(2) = nsec1jan000 - (1 - INT(ztmp)) * NINT(0.5 * rday) * nyear_len(0) + &
485                                    & INT(ztmp) * NINT( 0.5 * rday) * nyear_len(1) 
486            ELSE                  ;   sdjf%nrec_a(2) = nsec1jan000 + (1 - INT(ztmp)) * NINT(0.5 * rday) * nyear_len(1) + &
487                                    & INT(ztmp) * INT(rday) * nyear_len(1) + INT(ztmp) * NINT( 0.5 * rday) * nyear_len(2) 
488            ENDIF
489         ELSE                                    ! no time interpolation
490            sdjf%nrec_a(1) = 1
491            sdjf%nrec_a(2) = NINT(rday) * nyear_len(1) + nsec1jan000   ! swap at the end    of the year
492            sdjf%nrec_b(2) = nsec1jan000                               ! beginning of the year (only for print)
493         ENDIF
494         !
495         !                                   ! ============ !
496      ELSEIF( sdjf%nfreqh ==  -1 ) THEN      ! monthly mean !
497         !                                   ! ============ !
498         !
499         IF( sdjf%ln_tint ) THEN                 ! time interpolation, shift by 1/2 record
500            !
501            !                  INT( ztmp )
502            !                     /|\
503            !                    1 |    *----
504            !                    0 |----(             
505            !                      |----+----|--> time
506            !                      0   /|\   1   (nday/nmonth_len(nmonth))
507            !                           |   
508            !                           |   
509            !       forcing record :  nmonth
510            !                           
511            ztmp = REAL( nsec_month, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) + 0.5 &
512           &       + REAL( it_offset, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday )
513            imth = nmonth + INT( ztmp ) - COUNT((/llbefore/))
514            IF( sdjf%cltype == 'monthly' ) THEN   ;   sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
515            ELSE                                  ;   sdjf%nrec_a(1) = imth
516            ENDIF
517            sdjf%nrec_a(2) = nmonth_half(   imth ) + nsec1jan000   ! swap at the middle of the month
518         ELSE                                    ! no time interpolation
519            IF( sdjf%cltype == 'monthly' ) THEN   ;   sdjf%nrec_a(1) = 1
520            ELSE                                  ;   sdjf%nrec_a(1) = nmonth
521            ENDIF
522            sdjf%nrec_a(2) =  nmonth_end(nmonth  ) + nsec1jan000   ! swap at the end    of the month
523            sdjf%nrec_b(2) =  nmonth_end(nmonth-1) + nsec1jan000   ! beginning of the month (only for print)
524         ENDIF
525         !
526         !                                   ! ================================ !
527      ELSE                                   ! higher frequency mean (in hours)
528         !                                   ! ================================ !
529         !
530         ifreq_sec = NINT( sdjf%nfreqh * 3600 )                                         ! frequency mean (in seconds)
531         IF( sdjf%cltype(1:4) == 'week' )   isec_week = ksec_week( sdjf%cltype(6:8) )   ! since the first day of the current week
532         ! number of second since the beginning of the file
533         IF(     sdjf%cltype      == 'monthly' ) THEN   ;   ztmp = REAL(nsec_month,wp)  ! since the first day of the current month
534         ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   ztmp = REAL(isec_week ,wp)  ! since the first day of the current week
535         ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   ztmp = REAL(nsec_day  ,wp)  ! since 00h of the current day
536         ELSE                                           ;   ztmp = REAL(nsec_year ,wp)  ! since 00h on Jan 1 of the current year
537         ENDIF
538         ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdttra(1) + REAL( it_offset, wp )  ! centrered in the middle of sbc time step
539         ztmp = ztmp + 0.01 * rdttra(1)                                                 ! avoid truncation error
540         IF( sdjf%ln_tint ) THEN                ! time interpolation, shift by 1/2 record
541            !
542            !          INT( ztmp/ifreq_sec + 0.5 )
543            !                     /|\
544            !                    2 |        *-----(
545            !                    1 |  *-----(
546            !                    0 |--(             
547            !                      |--+--|--+--|--+--|--> time
548            !                      0 /|\ 1 /|\ 2 /|\ 3    (ztmp/ifreq_sec)
549            !                         |     |     |
550            !                         |     |     |
551            !       forcing record :  1     2     3
552            !                   
553            ztmp= ztmp / REAL(ifreq_sec, wp) + 0.5
554         ELSE                                   ! no time interpolation
555            !
556            !           INT( ztmp/ifreq_sec )
557            !                     /|\
558            !                    2 |           *-----(
559            !                    1 |     *-----(
560            !                    0 |-----(             
561            !                      |--+--|--+--|--+--|--> time
562            !                      0 /|\ 1 /|\ 2 /|\ 3    (ztmp/ifreq_sec)
563            !                         |     |     |
564            !                         |     |     |
565            !       forcing record :  1     2     3
566            !                           
567            ztmp= ztmp / REAL(ifreq_sec, wp)
568         ENDIF
569         sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))   ! record number to be read
570
571         iendrec = ifreq_sec * sdjf%nrec_a(1) + nsec1jan000       ! end of this record (in second)
572         ! add the number of seconds between 00h Jan 1 and the end of previous month/week/day (ok if nmonth=1)
573         IF( sdjf%cltype      == 'monthly' )   iendrec = iendrec + NINT(rday) * SUM(nmonth_len(1:nmonth -1))
574         IF( sdjf%cltype(1:4) == 'week'    )   iendrec = iendrec + ( nsec_year - isec_week )
575         IF( sdjf%cltype      == 'daily'   )   iendrec = iendrec + NINT(rday) * ( nday_year - 1 )
576         IF( sdjf%ln_tint ) THEN
577             sdjf%nrec_a(2) = iendrec - ifreq_sec / 2        ! swap at the middle of the record
578         ELSE
579             sdjf%nrec_a(2) = iendrec                        ! swap at the end    of the record
580             sdjf%nrec_b(2) = iendrec - ifreq_sec            ! beginning of the record (only for print)
581         ENDIF
582         !
583      ENDIF
584      !
585   END SUBROUTINE fld_rec
586
587
588   SUBROUTINE fld_get( sdjf, map )
589      !!---------------------------------------------------------------------
590      !!                    ***  ROUTINE fld_get  ***
591      !!
592      !! ** Purpose :   read the data
593      !!----------------------------------------------------------------------
594      TYPE(FLD), INTENT(inout) ::   sdjf   ! input field related variables
595      TYPE(MAP_POINTER),INTENT(in) ::   map   ! global-to-local mapping indices
596      !!
597      INTEGER                  ::   ipk    ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
598      INTEGER                  ::   iw     ! index into wgts array
599      INTEGER                  ::   ipdom  ! index of the domain
600      INTEGER                  ::   idvar  ! variable ID
601      INTEGER                  ::   idmspc ! number of spatial dimensions
602      LOGICAL                  ::   lmoor  ! C1D case: point data
603      !!---------------------------------------------------------------------
604      !
605      ipk = SIZE( sdjf%fnow, 3 )
606      !
607      IF( ASSOCIATED(map%ptr) ) THEN
608         IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map )
609         ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1), map )
610         ENDIF
611      ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
612         CALL wgt_list( sdjf, iw )
613         IF( sdjf%ln_tint ) THEN   ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk  , sdjf%fdta(:,:,:,2), & 
614              & sdjf%nrec_a(1), sdjf%lsmname )
615         ELSE                      ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk  , sdjf%fnow(:,:,:  ), &
616              & sdjf%nrec_a(1), sdjf%lsmname )
617         ENDIF
618      ELSE
619         IF( SIZE(sdjf%fnow, 1) == jpi ) THEN  ;  ipdom = jpdom_data
620         ELSE                                  ;  ipdom = jpdom_unknown
621         ENDIF
622         ! C1D case: If product of spatial dimensions == ipk, then x,y are of
623         ! size 1 (point/mooring data): this must be read onto the central grid point
624         idvar  = iom_varid( sdjf%num, sdjf%clvar )
625         idmspc = iom_file( sdjf%num )%ndims( idvar )
626         IF( iom_file( sdjf%num )%luld( idvar ) )   idmspc = idmspc - 1
627         lmoor  = (idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk)
628         !
629         SELECT CASE( ipk )
630         CASE(1)
631            IF( lk_c1d .AND. lmoor ) THEN
632               IF( sdjf%ln_tint ) THEN
633                  CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fdta(2,2,1,2), sdjf%nrec_a(1) )
634                  CALL lbc_lnk( sdjf%fdta(:,:,1,2),'Z',1. )
635               ELSE
636                  CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fnow(2,2,1  ), sdjf%nrec_a(1) )
637                  CALL lbc_lnk( sdjf%fnow(:,:,1  ),'Z',1. )
638               ENDIF
639            ELSE
640               IF( sdjf%ln_tint ) THEN   ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,1,2), sdjf%nrec_a(1) )
641               ELSE                      ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,1  ), sdjf%nrec_a(1) )
642               ENDIF
643            ENDIF
644         CASE DEFAULT
645            IF (lk_c1d .AND. lmoor ) THEN
646               IF( sdjf%ln_tint ) THEN
647                  CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fdta(2,2,:,2), sdjf%nrec_a(1) )
648                  CALL lbc_lnk( sdjf%fdta(:,:,:,2),'Z',1. )
649               ELSE
650                  CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fnow(2,2,:  ), sdjf%nrec_a(1) )
651                  CALL lbc_lnk( sdjf%fnow(:,:,:  ),'Z',1. )
652               ENDIF
653            ELSE
654               IF( sdjf%ln_tint ) THEN   ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) )
655               ELSE                      ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1) )
656               ENDIF
657            ENDIF
658         END SELECT
659      ENDIF
660      !
661      sdjf%rotn(2) = .false.   ! vector not yet rotated
662
663   END SUBROUTINE fld_get
664
665   SUBROUTINE fld_map( num, clvar, dta, nrec, map )
666      !!---------------------------------------------------------------------
667      !!                    ***  ROUTINE fld_map  ***
668      !!
669      !! ** Purpose :   read global data from file and map onto local data
670      !!                using a general mapping (for open boundaries)
671      !!----------------------------------------------------------------------
672#if defined key_bdy
673      USE bdy_oce, ONLY:  dta_global, dta_global2         ! workspace to read in global data arrays
674#endif
675      INTEGER                   , INTENT(in ) ::   num     ! stream number
676      CHARACTER(LEN=*)          , INTENT(in ) ::   clvar   ! variable name
677      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta   ! output field on model grid (2 dimensional)
678      INTEGER                   , INTENT(in ) ::   nrec    ! record number to read (ie time slice)
679      TYPE(MAP_POINTER)         , INTENT(in ) ::   map     ! global-to-local mapping indices
680      !!
681      INTEGER                                 ::   ipi      ! length of boundary data on local process
682      INTEGER                                 ::   ipj      ! length of dummy dimension ( = 1 )
683      INTEGER                                 ::   ipk      ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk )
684      INTEGER                                 ::   ilendta  ! length of data in file
685      INTEGER                                 ::   idvar    ! variable ID
686      INTEGER                                 ::   ib, ik, ji, jj   ! loop counters
687      INTEGER                                 ::   ierr
688      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read  ! work space for global data
689      !!---------------------------------------------------------------------
690           
691      ipi = SIZE( dta, 1 )
692      ipj = 1
693      ipk = SIZE( dta, 3 )
694
695      idvar   = iom_varid( num, clvar )
696      ilendta = iom_file(num)%dimsz(1,idvar)
697
698#if defined key_bdy
699      ipj = iom_file(num)%dimsz(2,idvar)
700      IF ( map%ll_unstruc) THEN ! unstructured open boundary data file
701         dta_read => dta_global
702      ELSE                      ! structured open boundary data file
703         dta_read => dta_global2
704      ENDIF
705#endif
706
707      IF(lwp) WRITE(numout,*) 'Dim size for ',TRIM(clvar),' is ', ilendta
708      IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk
709
710      SELECT CASE( ipk )
711      CASE(1)        ;   CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1    ), nrec )
712      CASE DEFAULT   ;   CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec )
713      END SELECT
714      !
715      IF ( map%ll_unstruc ) THEN ! unstructured open boundary data file
716         DO ib = 1, ipi
717            DO ik = 1, ipk
718               dta(ib,1,ik) =  dta_read(map%ptr(ib),1,ik)
719            END DO
720         END DO
721      ELSE                       ! structured open boundary data file
722         DO ib = 1, ipi
723            jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
724            ji=map%ptr(ib)-(jj-1)*ilendta
725            DO ik = 1, ipk
726               dta(ib,1,ik) =  dta_read(ji,jj,ik)
727            END DO
728         END DO
729      ENDIF
730
731   END SUBROUTINE fld_map
732
733
734   SUBROUTINE fld_rot( kt, sd )
735      !!---------------------------------------------------------------------
736      !!                    ***  ROUTINE fld_rot  ***
737      !!
738      !! ** Purpose :   Vector fields may need to be rotated onto the local grid direction
739      !!----------------------------------------------------------------------
740      INTEGER  , INTENT(in   )               ::   kt        ! ocean time step
741      TYPE(FLD), INTENT(inout), DIMENSION(:) ::   sd        ! input field related variables
742      !!
743      INTEGER                           ::   ju,jv,jk,jn  ! loop indices
744      INTEGER                           ::   imf          ! size of the structure sd
745      INTEGER                           ::   ill          ! character length
746      INTEGER                           ::   iv           ! indice of V component
747      REAL(wp), POINTER, DIMENSION(:,:) ::   utmp, vtmp   ! temporary arrays for vector rotation
748      CHARACTER (LEN=100)               ::   clcomp       ! dummy weight name
749      !!---------------------------------------------------------------------
750
751      CALL wrk_alloc( jpi,jpj, utmp, vtmp )
752
753      !! (sga: following code should be modified so that pairs arent searched for each time
754      !
755      imf = SIZE( sd )
756      DO ju = 1, imf
757         ill = LEN_TRIM( sd(ju)%vcomp )
758         DO jn = 2-COUNT((/sd(ju)%ln_tint/)), 2
759            IF( ill > 0 .AND. .NOT. sd(ju)%rotn(jn) ) THEN   ! find vector rotations required             
760               IF( sd(ju)%vcomp(1:1) == 'U' ) THEN      ! east-west component has symbolic name starting with 'U'
761                  ! look for the north-south component which has same symbolic name but with 'U' replaced with 'V'
762                  clcomp = 'V' // sd(ju)%vcomp(2:ill)   ! works even if ill == 1
763                  iv = -1
764                  DO jv = 1, imf
765                     IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) )   iv = jv
766                  END DO
767                  IF( iv > 0 ) THEN   ! fields ju and iv are two components which need to be rotated together
768                     DO jk = 1, SIZE( sd(ju)%fnow, 3 )
769                        IF( sd(ju)%ln_tint )THEN
770                           CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->i', utmp(:,:) )
771                           CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->j', vtmp(:,:) )
772                           sd(ju)%fdta(:,:,jk,jn) = utmp(:,:)   ;   sd(iv)%fdta(:,:,jk,jn) = vtmp(:,:)
773                        ELSE
774                           CALL rot_rep( sd(ju)%fnow(:,:,jk  ), sd(iv)%fnow(:,:,jk  ), 'T', 'en->i', utmp(:,:) )
775                           CALL rot_rep( sd(ju)%fnow(:,:,jk  ), sd(iv)%fnow(:,:,jk  ), 'T', 'en->j', vtmp(:,:) )
776                           sd(ju)%fnow(:,:,jk   ) = utmp(:,:)   ;   sd(iv)%fnow(:,:,jk   ) = vtmp(:,:)
777                        ENDIF
778                     END DO
779                     sd(ju)%rotn(jn) = .TRUE.               ! vector was rotated
780                     IF( lwp .AND. kt == nit000 )   WRITE(numout,*)   &
781                        &   'fld_read: vector pair ('//TRIM(sd(ju)%clvar)//', '//TRIM(sd(iv)%clvar)//') rotated on to model grid'
782                  ENDIF
783               ENDIF
784            ENDIF
785         END DO
786       END DO
787      !
788      CALL wrk_dealloc( jpi,jpj, utmp, vtmp )
789      !
790   END SUBROUTINE fld_rot
791
792
793   SUBROUTINE fld_clopn( sdjf, kyear, kmonth, kday, ldstop )
794      !!---------------------------------------------------------------------
795      !!                    ***  ROUTINE fld_clopn  ***
796      !!
797      !! ** Purpose :   update the file name and open the file
798      !!----------------------------------------------------------------------
799      TYPE(FLD)        , INTENT(inout) ::   sdjf     ! input field related variables
800      INTEGER, OPTIONAL, INTENT(in   ) ::   kyear    ! year value
801      INTEGER, OPTIONAL, INTENT(in   ) ::   kmonth   ! month value
802      INTEGER, OPTIONAL, INTENT(in   ) ::   kday     ! day value
803      LOGICAL, OPTIONAL, INTENT(in   ) ::   ldstop   ! stop if open to read a non-existing file (default = .TRUE.)
804      !!
805      LOGICAL :: llprevyr              ! are we reading previous year  file?
806      LOGICAL :: llprevmth             ! are we reading previous month file?
807      INTEGER :: iyear, imonth, iday   ! first day of the current file in yyyy mm dd
808      INTEGER :: isec_week             ! number of seconds since start of the weekly file
809      INTEGER :: indexyr               ! year undex (O/1/2: previous/current/next)
810      INTEGER :: iyear_len, imonth_len ! length (days) of iyear and imonth             !
811      CHARACTER(len = 256)::   clname  ! temporary file name
812      !!----------------------------------------------------------------------
813      IF( PRESENT(kyear) ) THEN                             ! use given values
814         iyear = kyear
815         imonth = kmonth
816         iday = kday
817         IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of
818         the beginning of the week
819            isec_week = ksec_week( sdjf%cltype(6:8) )- (86400 * 8 ) 
820            llprevmth  = isec_week > nsec_month             ! longer time since
821            beginning of the week than the month
822            llprevyr   = llprevmth .AND. nmonth == 1
823            iyear  = nyear  - COUNT((/llprevyr /))
824            imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
825            iday   = nday   + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) -
826            isec_week / NINT(rday)
827         ENDIF
828      ELSE                                                  ! use current day values
829         IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week
830            isec_week  = ksec_week( sdjf%cltype(6:8) )      ! second since the beginning of the week
831            llprevmth  = isec_week > nsec_month             ! longer time since beginning of the week than the month
832            llprevyr   = llprevmth .AND. nmonth == 1
833         ELSE
834            isec_week  = 0
835            llprevmth  = .FALSE.
836            llprevyr   = .FALSE.
837         ENDIF
838         iyear  = nyear  - COUNT((/llprevyr /))
839         imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
840         iday   = nday   + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
841      ENDIF
842
843      ! build the new filename if not climatological data
844      clname=TRIM(sdjf%clrootname)
845      !
846      ! note that sdjf%ln_clim is is only acting on the presence of the year in the file name
847      IF( .NOT. sdjf%ln_clim ) THEN   
848                                         WRITE(clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), iyear    ! add year
849         IF( sdjf%cltype /= 'yearly' )   WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname          ), imonth   ! add month
850      ELSE
851         ! build the new filename if climatological data
852         IF( sdjf%cltype /= 'yearly' )   WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), imonth   ! add month
853      ENDIF
854      IF( sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) &
855            &                            WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname          ), iday     ! add day
856      !
857      IF( TRIM(clname) /= TRIM(sdjf%clname) .OR. sdjf%num == 0 ) THEN   ! new file to be open
858
859         sdjf%clname = TRIM(clname)
860         IF( sdjf%num /= 0 )   CALL iom_close( sdjf%num )   ! close file if already open
861         CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof =  LEN(TRIM(sdjf%wgtname)) > 0 )
862
863         ! find the last record to be read -> update sdjf%nreclast
864         indexyr = iyear - nyear + 1
865         iyear_len = nyear_len( indexyr )
866         SELECT CASE ( indexyr )
867         CASE ( 0 )   ;   imonth_len = 31   ! previous year -> imonth = 12
868         CASE ( 1 )   ;   imonth_len = nmonth_len(imonth) 
869         CASE ( 2 )   ;   imonth_len = 31   ! next     year -> imonth = 1
870         END SELECT
871         
872         ! last record to be read in the current file
873         IF    ( sdjf%nfreqh == -12 ) THEN                 ;   sdjf%nreclast = 1    !  yearly mean
874         ELSEIF( sdjf%nfreqh ==  -1 ) THEN                                          ! monthly mean
875            IF(     sdjf%cltype      == 'monthly' ) THEN   ;   sdjf%nreclast = 1
876            ELSE                                           ;   sdjf%nreclast = 12
877            ENDIF
878         ELSE                                                                       ! higher frequency mean (in hours)
879            IF(     sdjf%cltype      == 'monthly' ) THEN   ;   sdjf%nreclast = NINT( 24 * imonth_len / sdjf%nfreqh )
880            ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   sdjf%nreclast = NINT( 24 * 7          / sdjf%nfreqh )
881            ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   sdjf%nreclast = NINT( 24              / sdjf%nfreqh )
882            ELSE                                           ;   sdjf%nreclast = NINT( 24 * iyear_len  / sdjf%nfreqh )
883            ENDIF
884         ENDIF
885         
886      ENDIF
887      !
888   END SUBROUTINE fld_clopn
889
890
891   SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )
892      !!---------------------------------------------------------------------
893      !!                    ***  ROUTINE fld_fill  ***
894      !!
895      !! ** Purpose :   fill sdf with sdf_n and control print
896      !!----------------------------------------------------------------------
897      TYPE(FLD)  , DIMENSION(:), INTENT(inout) ::   sdf        ! structure of input fields (file informations, fields read)
898      TYPE(FLD_N), DIMENSION(:), INTENT(in   ) ::   sdf_n      ! array of namelist information structures
899      CHARACTER(len=*)         , INTENT(in   ) ::   cdir       ! Root directory for location of flx files
900      CHARACTER(len=*)         , INTENT(in   ) ::   cdcaller   !
901      CHARACTER(len=*)         , INTENT(in   ) ::   cdtitle    !
902      CHARACTER(len=*)         , INTENT(in   ) ::   cdnam      !
903      !
904      INTEGER  ::   jf       ! dummy indices
905      !!---------------------------------------------------------------------
906
907      DO jf = 1, SIZE(sdf)
908         sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
909         sdf(jf)%clname     = "not yet defined"
910         sdf(jf)%nfreqh     = sdf_n(jf)%nfreqh
911         sdf(jf)%clvar      = sdf_n(jf)%clvar
912         sdf(jf)%ln_tint    = sdf_n(jf)%ln_tint
913         sdf(jf)%ln_clim    = sdf_n(jf)%ln_clim
914         sdf(jf)%cltype     = sdf_n(jf)%cltype
915         sdf(jf)%num        = -1
916         sdf(jf)%wgtname    = " "
917         IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
918         sdf(jf)%lsmname = " "
919         IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 )   sdf(jf)%lsmname = TRIM( cdir )//TRIM( sdf_n(jf)%lname )
920         sdf(jf)%vcomp      = sdf_n(jf)%vcomp
921         sdf(jf)%rotn(:)    = .TRUE.   ! pretend to be rotated -> won't try to rotate data before the first call to fld_get
922         IF( sdf(jf)%cltype(1:4) == 'week' .AND. nn_leapy == 0  )   &
923            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs nn_leapy = 1')
924         IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim )   &
925            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.')
926         sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
927      END DO
928
929      IF(lwp) THEN      ! control print
930         WRITE(numout,*)
931         WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
932         WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
933         WRITE(numout,*) '          '//TRIM( cdnam )//' Namelist'
934         WRITE(numout,*) '          list of files and frequency (>0: in hours ; <0 in months)'
935         DO jf = 1, SIZE(sdf)
936            WRITE(numout,*) '               root filename: '  , TRIM( sdf(jf)%clrootname ),   &
937               &                          ' variable name: '  , TRIM( sdf(jf)%clvar      )
938            WRITE(numout,*) '               frequency: '      ,       sdf(jf)%nfreqh      ,   &
939               &                          ' time interp: '    ,       sdf(jf)%ln_tint     ,   &
940               &                          ' climatology: '    ,       sdf(jf)%ln_clim     ,   &
941               &                          ' weights    : '    , TRIM( sdf(jf)%wgtname    ),   &
942               &                          ' pairing    : '    , TRIM( sdf(jf)%vcomp      ),   &
943               &                          ' data type: '      ,       sdf(jf)%cltype      ,   &
944               &                          ' land/sea mask:'   , TRIM( sdf(jf)%lsmname    )
945            call flush(numout)
946         END DO
947      ENDIF
948     
949   END SUBROUTINE fld_fill
950
951
952   SUBROUTINE wgt_list( sd, kwgt )
953      !!---------------------------------------------------------------------
954      !!                    ***  ROUTINE wgt_list  ***
955      !!
956      !! ** Purpose :   search array of WGTs and find a weights file
957      !!                entry, or return a new one adding it to the end
958      !!                if it is a new entry, the weights data is read in and
959      !!                restructured (fld_weight)
960      !!----------------------------------------------------------------------
961      TYPE( FLD ), INTENT(in   ) ::   sd        ! field with name of weights file
962      INTEGER    , INTENT(inout) ::   kwgt      ! index of weights
963      !!
964      INTEGER ::   kw, nestid   ! local integer
965      LOGICAL ::   found        ! local logical
966      !!----------------------------------------------------------------------
967      !
968      !! search down linked list
969      !! weights filename is either present or we hit the end of the list
970      found = .FALSE.
971
972      !! because agrif nest part of filenames are now added in iom_open
973      !! to distinguish between weights files on the different grids, need to track
974      !! nest number explicitly
975      nestid = 0
976#if defined key_agrif
977      nestid = Agrif_Fixed()
978#endif
979      DO kw = 1, nxt_wgt-1
980         IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
981             ref_wgts(kw)%nestid == nestid) THEN
982            kwgt = kw
983            found = .TRUE.
984            EXIT
985         ENDIF
986      END DO
987      IF( .NOT.found ) THEN
988         kwgt = nxt_wgt
989         CALL fld_weight( sd )
990      ENDIF
991      !
992   END SUBROUTINE wgt_list
993
994
995   SUBROUTINE wgt_print( )
996      !!---------------------------------------------------------------------
997      !!                    ***  ROUTINE wgt_print  ***
998      !!
999      !! ** Purpose :   print the list of known weights
1000      !!----------------------------------------------------------------------
1001      INTEGER ::   kw   !
1002      !!----------------------------------------------------------------------
1003      !
1004      DO kw = 1, nxt_wgt-1
1005         WRITE(numout,*) 'weight file:  ',TRIM(ref_wgts(kw)%wgtname)
1006         WRITE(numout,*) '      ddims:  ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
1007         WRITE(numout,*) '     numwgt:  ',ref_wgts(kw)%numwgt
1008         WRITE(numout,*) '     jpiwgt:  ',ref_wgts(kw)%jpiwgt
1009         WRITE(numout,*) '     jpjwgt:  ',ref_wgts(kw)%jpjwgt
1010         WRITE(numout,*) '    botleft:  ',ref_wgts(kw)%botleft
1011         WRITE(numout,*) '   topright:  ',ref_wgts(kw)%topright
1012         IF( ref_wgts(kw)%cyclic ) THEN
1013            WRITE(numout,*) '       cyclical'
1014            IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) '              with overlap of ', ref_wgts(kw)%overlap
1015         ELSE
1016            WRITE(numout,*) '       not cyclical'
1017         ENDIF
1018         IF( ASSOCIATED(ref_wgts(kw)%data_wgt) )  WRITE(numout,*) '       allocated'
1019      END DO
1020      !
1021   END SUBROUTINE wgt_print
1022
1023
1024   SUBROUTINE fld_weight( sd )
1025      !!---------------------------------------------------------------------
1026      !!                    ***  ROUTINE fld_weight  ***
1027      !!
1028      !! ** Purpose :   create a new WGT structure and fill in data from 
1029      !!                file, restructuring as required
1030      !!----------------------------------------------------------------------
1031      TYPE( FLD ), INTENT(in) ::   sd   ! field with name of weights file
1032      !!
1033      INTEGER                           ::   jn            ! dummy loop indices
1034      INTEGER                           ::   inum          ! temporary logical unit
1035      INTEGER                           ::   id            ! temporary variable id
1036      INTEGER                           ::   ipk           ! temporary vertical dimension
1037      CHARACTER (len=5)                 ::   aname
1038      INTEGER , DIMENSION(:), ALLOCATABLE ::   ddims
1039      INTEGER , POINTER, DIMENSION(:,:) ::   data_src
1040      REAL(wp), POINTER, DIMENSION(:,:) ::   data_tmp
1041      LOGICAL                           ::   cyclical
1042      INTEGER                           ::   zwrap      ! local integer
1043      !!----------------------------------------------------------------------
1044      !
1045      CALL wrk_alloc( jpi,jpj, data_src )   ! integer
1046      CALL wrk_alloc( jpi,jpj, data_tmp )
1047      !
1048      IF( nxt_wgt > tot_wgts ) THEN
1049        CALL ctl_stop("fld_weight: weights array size exceeded, increase tot_wgts")
1050      ENDIF
1051      !
1052      !! new weights file entry, add in extra information
1053      !! a weights file represents a 2D grid of a certain shape, so we assume that the current
1054      !! input data file is representative of all other files to be opened and processed with the
1055      !! current weights file
1056
1057      !! open input data file (non-model grid)
1058      CALL iom_open( sd%clname, inum, ldiof =  LEN(TRIM(sd%wgtname)) > 0 )
1059
1060      !! get dimensions
1061      IF ( SIZE(sd%fnow, 3) > 1 ) THEN
1062         ALLOCATE( ddims(4) )
1063      ELSE
1064         ALLOCATE( ddims(3) )
1065      ENDIF
1066      id = iom_varid( inum, sd%clvar, ddims )
1067
1068      !! close it
1069      CALL iom_close( inum )
1070
1071      !! now open the weights file
1072
1073      CALL iom_open ( sd%wgtname, inum )   ! interpolation weights
1074      IF ( inum > 0 ) THEN
1075
1076         !! determine whether we have an east-west cyclic grid
1077         !! from global attribute called "ew_wrap" in the weights file
1078         !! note that if not found, iom_getatt returns -999 and cyclic with no overlap is assumed
1079         !! since this is the most common forcing configuration
1080
1081         CALL iom_getatt(inum, 'ew_wrap', zwrap)
1082         IF( zwrap >= 0 ) THEN
1083            cyclical = .TRUE.
1084         ELSE IF( zwrap == -999 ) THEN
1085            cyclical = .TRUE.
1086            zwrap = 0
1087         ELSE
1088            cyclical = .FALSE.
1089         ENDIF
1090
1091         ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
1092         ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
1093         ref_wgts(nxt_wgt)%wgtname = sd%wgtname
1094         ref_wgts(nxt_wgt)%overlap = zwrap
1095         ref_wgts(nxt_wgt)%cyclic = cyclical
1096         ref_wgts(nxt_wgt)%nestid = 0
1097#if defined key_agrif
1098         ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
1099#endif
1100         !! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
1101         !! for each weight wgtNN there is an integer array srcNN which gives the point in
1102         !! the input data grid which is to be multiplied by the weight
1103         !! they are both arrays on the model grid so the result of the multiplication is
1104         !! added into an output array on the model grid as a running sum
1105
1106         !! two possible cases: bilinear (4 weights) or bicubic (16 weights)
1107         id = iom_varid(inum, 'src05', ldstop=.FALSE.)
1108         IF( id <= 0) THEN
1109            ref_wgts(nxt_wgt)%numwgt = 4
1110         ELSE
1111            ref_wgts(nxt_wgt)%numwgt = 16
1112         ENDIF
1113
1114         ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
1115         ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
1116         ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
1117
1118         DO jn = 1,4
1119            aname = ' '
1120            WRITE(aname,'(a3,i2.2)') 'src',jn
1121            data_tmp(:,:) = 0
1122            CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
1123            data_src(:,:) = INT(data_tmp(:,:))
1124            ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
1125            ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
1126         END DO
1127
1128         DO jn = 1, ref_wgts(nxt_wgt)%numwgt
1129            aname = ' '
1130            WRITE(aname,'(a3,i2.2)') 'wgt',jn
1131            ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
1132            CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
1133         END DO
1134         CALL iom_close (inum)
1135 
1136         ! find min and max indices in grid
1137         ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1138         ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1139         ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1140         ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1141
1142         ! and therefore dimensions of the input box
1143         ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
1144         ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
1145
1146         ! shift indexing of source grid
1147         ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
1148         ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
1149
1150         ! create input grid, give it a halo to allow gradient calculations
1151         ! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
1152         ! a more robust solution will be given in next release
1153         ipk =  SIZE(sd%fnow, 3)
1154         ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
1155         IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
1156
1157         nxt_wgt = nxt_wgt + 1
1158
1159      ELSE
1160         CALL ctl_stop( '    fld_weight : unable to read the file ' )
1161      ENDIF
1162
1163      DEALLOCATE (ddims )
1164
1165      CALL wrk_dealloc( jpi,jpj, data_src )   ! integer
1166      CALL wrk_dealloc( jpi,jpj, data_tmp )
1167      !
1168   END SUBROUTINE fld_weight
1169
1170
1171   SUBROUTINE apply_seaoverland(clmaskfile,zfieldo,jpi1_lsm,jpi2_lsm,jpj1_lsm, &
1172                          &      jpj2_lsm,itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
1173      !!---------------------------------------------------------------------
1174      !!                    ***  ROUTINE apply_seaoverland  ***
1175      !!
1176      !! ** Purpose :   avoid spurious fluxes in coastal or near-coastal areas
1177      !!                due to the wrong usage of "land" values from the coarse
1178      !!                atmospheric model when spatial interpolation is required
1179      !!      D. Delrosso INGV         
1180      !!----------------------------------------------------------------------
1181      INTEGER                                   :: inum,jni,jnj,jnz,jc                  ! temporary indices
1182      INTEGER,                   INTENT(in)     :: itmpi,itmpj,itmpz                    ! lengths
1183      INTEGER,                   INTENT(in)     :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm  ! temporary indices
1184      INTEGER, DIMENSION(3),     INTENT(in)     :: rec1_lsm,recn_lsm                    ! temporary arrays for start and length
1185      REAL(wp),DIMENSION (:,:,:),INTENT(inout)  :: zfieldo                              ! input/output array for seaoverland application
1186      REAL(wp),DIMENSION (:,:,:),ALLOCATABLE    :: zslmec1                              ! temporary array for land point detection
1187      REAL(wp),DIMENSION (:,:),  ALLOCATABLE    :: zfieldn                              ! array of forcing field with undeff for land points
1188      REAL(wp),DIMENSION (:,:),  ALLOCATABLE    :: zfield                               ! array of forcing field
1189      CHARACTER (len=100),       INTENT(in)     :: clmaskfile                           ! land/sea mask file name
1190      !!---------------------------------------------------------------------
1191      ALLOCATE ( zslmec1(itmpi,itmpj,itmpz) )
1192      ALLOCATE ( zfieldn(itmpi,itmpj) )
1193      ALLOCATE ( zfield(itmpi,itmpj) )
1194
1195      ! Retrieve the land sea mask data
1196      CALL iom_open( clmaskfile, inum )
1197      SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1198      CASE(1)
1199           CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm)
1200      CASE DEFAULT
1201           CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm)
1202      END SELECT
1203      CALL iom_close( inum )
1204
1205      DO jnz=1,rec1_lsm(3)                            !! Loop over k dimension
1206
1207         DO jni=1,itmpi                               !! copy the original field into a tmp array
1208            DO jnj=1,itmpj                            !! substituting undeff over land points
1209            zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz)
1210               IF ( zslmec1(jni,jnj,jnz) == 1. ) THEN
1211                  zfieldn(jni,jnj) = undeff_lsm
1212               ENDIF
1213            END DO
1214         END DO
1215 
1216      CALL seaoverland(zfieldn,itmpi,itmpj,zfield)
1217      DO jc=1,nn_lsm
1218         CALL seaoverland(zfield,itmpi,itmpj,zfield)
1219      END DO
1220
1221      !   Check for Undeff and substitute original values
1222      IF(ANY(zfield==undeff_lsm)) THEN
1223         DO jni=1,itmpi
1224            DO jnj=1,itmpj
1225               IF (zfield(jni,jnj)==undeff_lsm) THEN
1226                  zfield(jni,jnj) = zfieldo(jni,jnj,jnz)
1227               ENDIF
1228            ENDDO
1229         ENDDO
1230      ENDIF
1231
1232      zfieldo(:,:,jnz)=zfield(:,:)
1233
1234      END DO                          !! End Loop over k dimension
1235
1236      DEALLOCATE ( zslmec1 )
1237      DEALLOCATE ( zfieldn )
1238      DEALLOCATE ( zfield )
1239
1240   END SUBROUTINE apply_seaoverland 
1241
1242
1243   SUBROUTINE seaoverland(zfieldn,ileni,ilenj,zfield)
1244      !!---------------------------------------------------------------------
1245      !!                    ***  ROUTINE seaoverland  ***
1246      !!
1247      !! ** Purpose :   create shifted matrices for seaoverland application 
1248      !!      D. Delrosso INGV
1249      !!----------------------------------------------------------------------
1250      INTEGER,INTENT(in)                       :: ileni,ilenj              ! lengths
1251      REAL,DIMENSION (ileni,ilenj),INTENT(in)  :: zfieldn                  ! array of forcing field with undeff for land points
1252      REAL,DIMENSION (ileni,ilenj),INTENT(out) :: zfield                   ! array of forcing field
1253      REAL,DIMENSION (ileni,ilenj)             :: zmat1,zmat2,zmat3,zmat4  ! temporary arrays for seaoverland application
1254      REAL,DIMENSION (ileni,ilenj)             :: zmat5,zmat6,zmat7,zmat8  ! temporary arrays for seaoverland application
1255      REAL,DIMENSION (ileni,ilenj)             :: zlsm2d                   ! temporary arrays for seaoverland application
1256      REAL,DIMENSION (ileni,ilenj,8)           :: zlsm3d                   ! temporary arrays for seaoverland application
1257      LOGICAL,DIMENSION (ileni,ilenj,8)        :: ll_msknan3d              ! logical mask for undeff detection
1258      LOGICAL,DIMENSION (ileni,ilenj)          :: ll_msknan2d              ! logical mask for undeff detection
1259      !!----------------------------------------------------------------------
1260      zmat8 = eoshift(zfieldn   ,  SHIFT=-1, BOUNDARY = (/zfieldn(:,1)/)    ,DIM=2)
1261      zmat1 = eoshift(zmat8     ,  SHIFT=-1, BOUNDARY = (/zmat8(1,:)/)      ,DIM=1)
1262      zmat2 = eoshift(zfieldn   ,  SHIFT=-1, BOUNDARY = (/zfieldn(1,:)/)    ,DIM=1)
1263      zmat4 = eoshift(zfieldn   ,  SHIFT= 1, BOUNDARY = (/zfieldn(:,ilenj)/),DIM=2)
1264      zmat3 = eoshift(zmat4     ,  SHIFT=-1, BOUNDARY = (/zmat4(1,:)/)      ,DIM=1)
1265      zmat5 = eoshift(zmat4     ,  SHIFT= 1, BOUNDARY = (/zmat4(ileni,:)/)  ,DIM=1)
1266      zmat6 = eoshift(zfieldn   ,  SHIFT= 1, BOUNDARY = (/zfieldn(ileni,:)/),DIM=1)
1267      zmat7 = eoshift(zmat8     ,  SHIFT= 1, BOUNDARY = (/zmat8(ileni,:)/)  ,DIM=1)
1268
1269      zlsm3d  = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /))
1270      ll_msknan3d = .not.(zlsm3d==undeff_lsm)
1271      ll_msknan2d = .not.(zfieldn==undeff_lsm)  ! FALSE where is Undeff (land)
1272      zlsm2d = (SUM ( zlsm3d, 3 , ll_msknan3d ) )/(MAX(1,(COUNT( ll_msknan3d , 3 ))   ))
1273      WHERE ((COUNT( ll_msknan3d , 3 )) == 0.0_wp)  zlsm2d = undeff_lsm
1274      zfield = MERGE (zfieldn,zlsm2d,ll_msknan2d)
1275   END SUBROUTINE seaoverland
1276
1277
1278   SUBROUTINE fld_interp( num, clvar, kw, kk, dta,  &
1279                          &         nrec, lsmfile)     
1280      !!---------------------------------------------------------------------
1281      !!                    ***  ROUTINE fld_interp  ***
1282      !!
1283      !! ** Purpose :   apply weights to input gridded data to create data
1284      !!                on model grid
1285      !!----------------------------------------------------------------------
1286      INTEGER                   , INTENT(in   ) ::   num     ! stream number
1287      CHARACTER(LEN=*)          , INTENT(in   ) ::   clvar   ! variable name
1288      INTEGER                   , INTENT(in   ) ::   kw      ! weights number
1289      INTEGER                   , INTENT(in   ) ::   kk      ! vertical dimension of kk
1290      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   dta     ! output field on model grid
1291      INTEGER                   , INTENT(in   ) ::   nrec    ! record number to read (ie time slice)
1292      CHARACTER(LEN=*)          , INTENT(in   ) ::   lsmfile ! land sea mask file name
1293      !!
1294      REAL(wp),DIMENSION(:,:,:),ALLOCATABLE     ::   ztmp_fly_dta,zfieldo                  ! temporary array of values on input grid
1295      INTEGER, DIMENSION(3)                     ::   rec1,recn                             ! temporary arrays for start and length
1296      INTEGER, DIMENSION(3)                     ::   rec1_lsm,recn_lsm                     ! temporary arrays for start and length in case of seaoverland
1297      INTEGER                                   ::   ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2       ! temporary indices
1298      INTEGER                                   ::   jk, jn, jm, jir, jjr                  ! loop counters
1299      INTEGER                                   ::   ni, nj                                ! lengths
1300      INTEGER                                   ::   jpimin,jpiwid                         ! temporary indices
1301      INTEGER                                   ::   jpimin_lsm,jpiwid_lsm                 ! temporary indices
1302      INTEGER                                   ::   jpjmin,jpjwid                         ! temporary indices
1303      INTEGER                                   ::   jpjmin_lsm,jpjwid_lsm                 ! temporary indices
1304      INTEGER                                   ::   jpi1,jpi2,jpj1,jpj2                   ! temporary indices
1305      INTEGER                                   ::   jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm   ! temporary indices
1306      INTEGER                                   ::   itmpi,itmpj,itmpz                     ! lengths
1307     
1308      !!----------------------------------------------------------------------
1309      !
1310      !! for weighted interpolation we have weights at four corners of a box surrounding
1311      !! a model grid point, each weight is multiplied by a grid value (bilinear case)
1312      !! or by a grid value and gradients at the corner point (bicubic case)
1313      !! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
1314
1315      !! sub grid from non-model input grid which encloses all grid points in this nemo process
1316      jpimin = ref_wgts(kw)%botleft(1)
1317      jpjmin = ref_wgts(kw)%botleft(2)
1318      jpiwid = ref_wgts(kw)%jpiwgt
1319      jpjwid = ref_wgts(kw)%jpjwgt
1320
1321      !! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients
1322      rec1(1) = MAX( jpimin-1, 1 )
1323      rec1(2) = MAX( jpjmin-1, 1 )
1324      rec1(3) = 1
1325      recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
1326      recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1327      recn(3) = kk
1328
1329      !! where we need to put it in the non-nemo grid fly_dta
1330      !! note that jpi1 and jpj1 only differ from 1 when jpimin and jpjmin are 1
1331      !! (ie at the extreme west or south of the whole input grid) and similarly for jpi2 and jpj2
1332      jpi1 = 2 + rec1(1) - jpimin
1333      jpj1 = 2 + rec1(2) - jpjmin
1334      jpi2 = jpi1 + recn(1) - 1
1335      jpj2 = jpj1 + recn(2) - 1
1336
1337
1338      IF( LEN( TRIM(lsmfile) ) > 0 ) THEN
1339      !! indeces for ztmp_fly_dta
1340      ! --------------------------
1341         rec1_lsm(1)=MAX(rec1(1)-nn_lsm,1)  ! starting index for enlarged external data, x direction
1342         rec1_lsm(2)=MAX(rec1(2)-nn_lsm,1)  ! starting index for enlarged external data, y direction
1343         rec1_lsm(3) = 1                    ! vertical dimension
1344         recn_lsm(1)=MIN(rec1(1)-rec1_lsm(1)+recn(1)+nn_lsm,ref_wgts(kw)%ddims(1)-rec1_lsm(1)) ! n points in x direction
1345         recn_lsm(2)=MIN(rec1(2)-rec1_lsm(2)+recn(2)+nn_lsm,ref_wgts(kw)%ddims(2)-rec1_lsm(2)) ! n points in y direction
1346         recn_lsm(3) = kk                   ! number of vertical levels in the input file
1347
1348      !  Avoid out of bound
1349         jpimin_lsm = MAX( rec1_lsm(1)+1, 1 )
1350         jpjmin_lsm = MAX( rec1_lsm(2)+1, 1 )
1351         jpiwid_lsm = MIN( recn_lsm(1)-2,ref_wgts(kw)%ddims(1)-rec1(1)+1)
1352         jpjwid_lsm = MIN( recn_lsm(2)-2,ref_wgts(kw)%ddims(2)-rec1(2)+1)
1353
1354         jpi1_lsm = 2+rec1_lsm(1)-jpimin_lsm
1355         jpj1_lsm = 2+rec1_lsm(2)-jpjmin_lsm
1356         jpi2_lsm = jpi1_lsm + recn_lsm(1) - 1
1357         jpj2_lsm = jpj1_lsm + recn_lsm(2) - 1
1358
1359
1360         itmpi=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),1)
1361         itmpj=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),2)
1362         itmpz=kk
1363         ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz))
1364         ztmp_fly_dta(:,:,:) = 0.0
1365         SELECT CASE( SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1366         CASE(1)
1367              CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1),   &
1368                 &                                                                nrec, rec1_lsm, recn_lsm)
1369         CASE DEFAULT
1370              CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),   &
1371                 &                                                                nrec, rec1_lsm, recn_lsm)
1372         END SELECT
1373         CALL apply_seaoverland(lsmfile,ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),                  &
1374                 &                                      jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm,                  &
1375                 &                                      itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
1376
1377
1378         ! Relative indeces for remapping
1379         ii_lsm1 = (rec1(1)-rec1_lsm(1))+1
1380         ii_lsm2 = (ii_lsm1+recn(1))-1
1381         ij_lsm1 = (rec1(2)-rec1_lsm(2))+1
1382         ij_lsm2 = (ij_lsm1+recn(2))-1
1383
1384         ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1385         ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:) = ztmp_fly_dta(ii_lsm1:ii_lsm2,ij_lsm1:ij_lsm2,:)
1386         DEALLOCATE(ztmp_fly_dta)
1387
1388      ELSE
1389         
1390         ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1391         SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) )
1392         CASE(1)
1393              CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn)
1394         CASE DEFAULT
1395              CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn)
1396         END SELECT
1397      ENDIF
1398     
1399
1400      !! first four weights common to both bilinear and bicubic
1401      !! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft
1402      !! note that we have to offset by 1 into fly_dta array because of halo
1403      dta(:,:,:) = 0.0
1404      DO jk = 1,4
1405        DO jn = 1, jpj
1406          DO jm = 1,jpi
1407            ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1408            nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1409            dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,:)
1410          END DO
1411        END DO
1412      END DO
1413
1414      IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
1415
1416        !! fix up halo points that we couldnt read from file
1417        IF( jpi1 == 2 ) THEN
1418           ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:)
1419        ENDIF
1420        IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1421           ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:)
1422        ENDIF
1423        IF( jpj1 == 2 ) THEN
1424           ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:)
1425        ENDIF
1426        IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
1427           ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:)
1428        ENDIF
1429
1430        !! if data grid is cyclic we can do better on east-west edges
1431        !! but have to allow for whether first and last columns are coincident
1432        IF( ref_wgts(kw)%cyclic ) THEN
1433           rec1(2) = MAX( jpjmin-1, 1 )
1434           recn(1) = 1
1435           recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1436           jpj1 = 2 + rec1(2) - jpjmin
1437           jpj2 = jpj1 + recn(2) - 1
1438           IF( jpi1 == 2 ) THEN
1439              rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap
1440              SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1441              CASE(1)
1442                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1443              CASE DEFAULT
1444                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1445              END SELECT     
1446              ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1447           ENDIF
1448           IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1449              rec1(1) = 1 + ref_wgts(kw)%overlap
1450              SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1451              CASE(1)
1452                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1453              CASE DEFAULT
1454                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1455              END SELECT
1456              ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1457           ENDIF
1458        ENDIF
1459
1460        ! gradient in the i direction
1461        DO jk = 1,4
1462          DO jn = 1, jpj
1463            DO jm = 1,jpi
1464              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1465              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1466              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 *         &
1467                               (ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:))
1468            END DO
1469          END DO
1470        END DO
1471
1472        ! gradient in the j direction
1473        DO jk = 1,4
1474          DO jn = 1, jpj
1475            DO jm = 1,jpi
1476              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1477              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1478              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 *         &
1479                               (ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:))
1480            END DO
1481          END DO
1482        END DO
1483
1484         ! gradient in the ij direction
1485         DO jk = 1,4
1486            DO jn = 1, jpj
1487               DO jm = 1,jpi
1488                  ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1489                  nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1490                  dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
1491                               (ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni  ,nj+2,:)) -   &
1492                               (ref_wgts(kw)%fly_dta(ni+2,nj  ,:) - ref_wgts(kw)%fly_dta(ni  ,nj  ,:)))
1493               END DO
1494            END DO
1495         END DO
1496         !
1497      END IF
1498      !
1499   END SUBROUTINE fld_interp
1500
1501
1502   FUNCTION ksec_week( cdday )
1503      !!---------------------------------------------------------------------
1504      !!                    ***  FUNCTION kshift_week ***
1505      !!
1506      !! ** Purpose : 
1507      !!---------------------------------------------------------------------
1508      CHARACTER(len=*), INTENT(in)   ::   cdday   !3 first letters of the first day of the weekly file
1509      !!
1510      INTEGER                        ::   ksec_week  ! output variable
1511      INTEGER                        ::   ijul       !temp variable
1512      INTEGER                        ::   ishift     !temp variable
1513      CHARACTER(len=3),DIMENSION(7)  ::   cl_week 
1514      !!----------------------------------------------------------------------
1515      cl_week = (/"sun","sat","fri","thu","wed","tue","mon"/)
1516      DO ijul = 1, 7
1517         IF( cl_week(ijul) == TRIM(cdday) ) EXIT
1518      END DO
1519      IF( ijul .GT. 7 )   CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): '//TRIM(cdday) )
1520      !
1521      ishift = ijul * NINT(rday)
1522      !
1523      ksec_week = nsec_week + ishift
1524      ksec_week = MOD( ksec_week, 7*NINT(rday) )
1525      !
1526   END FUNCTION ksec_week
1527
1528   !!======================================================================
1529END MODULE fldread
Note: See TracBrowser for help on using the repository browser.