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

Last change on this file since 10922 was 10922, checked in by acc, 19 months ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert IOM, LDF, OBS and SBC directories and compatibility changes elsewhere that these changes enforce. Changes pass SETTE and compare with original trunk results. Outstanding issues (currently with work-arounds) in DIU/step_diu.F90 and fld_bdy_interp within SBC/fldread.F90; proper soltions pending

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