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

Last change on this file since 3851 was 3851, checked in by smasson, 9 years ago

trunk: bugfix in fldread when nn_fsbc * rn_rdt > frcing file period, see #958

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