source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/SBC/fldread.F90 @ 11053

Last change on this file since 11053 was 11053, checked in by davestorkey, 19 months ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap : Merge in latest changes from main branch and finish conversion of "h" variables. NB. This version still doesn't work!

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