New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
fldread.F90 in branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 11080

Last change on this file since 11080 was 11080, checked in by jcastill, 6 years ago

Latest changes from Jason

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