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

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

source: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 4694

Last change on this file since 4694 was 4694, checked in by jamesharle, 10 years ago

Update of fldread to handle depth information in BDY files and addition of an interpolation routine. Updated BDY code to handle T/S BDY interpolation on the fly. Conservative remapping of U/V still to be coded. Not compiled or test yet.

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