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 NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC – NEMO

source: NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/fldread.F90 @ 9923

Last change on this file since 9923 was 9923, checked in by gm, 6 years ago

#1911 (ENHANCE-04): step I.2: dev_r9838_ENHANCE04_MLF

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