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

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 6487

Last change on this file since 6487 was 6487, checked in by davestorkey, 9 years ago

Changes from nemo_v3_6_STABLE_copy branch.
Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r6237 cf. r5781 of /branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM@6486

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