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

Last change on this file was 62, checked in by smasson, 12 years ago

bugfix in fldread when rotating field with no time interpolation, see nemo ticket 940

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