Ticket #1222: fldread.F90

File fldread.F90, 84.0 KB (added by rfurner, 8 years ago)

proposed fix to fld_read

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