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 branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 5 years ago

GMED 450 add flush after prints

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