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 @ 5035

Last change on this file since 5035 was 5035, checked in by jamesharle, 9 years ago

Update to reading in thickness variables in BDY read for bt velocity
adjustment when using interpolation on the fly.

  • Property svn:keywords set to Id
File size: 94.6 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 = 256) ::   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( nsec_year, wp ) / ( REAL( nyear_len(1), wp ) * rday ) + 0.5 &
492           &       + REAL( it_offset, wp ) / ( REAL( nyear_len(1), wp ) * rday )
493            sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
494            ! swap at the middle of the year
495            IF( llbefore ) THEN   ;   sdjf%nrec_a(2) = nsec1jan000 - (1 - INT(ztmp)) * NINT(0.5 * rday) * nyear_len(0) + &
496                                    & INT(ztmp) * NINT( 0.5 * rday) * nyear_len(1) 
497            ELSE                  ;   sdjf%nrec_a(2) = nsec1jan000 + (1 - INT(ztmp)) * NINT(0.5 * rday) * nyear_len(1) + &
498                                    & INT(ztmp) * INT(rday) * nyear_len(1) + INT(ztmp) * NINT( 0.5 * rday) * nyear_len(2) 
499            ENDIF
500         ELSE                                    ! no time interpolation
501            sdjf%nrec_a(1) = 1
502            sdjf%nrec_a(2) = NINT(rday) * nyear_len(1) + nsec1jan000   ! swap at the end    of the year
503            sdjf%nrec_b(2) = nsec1jan000                               ! beginning of the year (only for print)
504         ENDIF
505         !
506         !                                   ! ============ !
507      ELSEIF( sdjf%nfreqh ==  -1 ) THEN      ! monthly mean !
508         !                                   ! ============ !
509         !
510         IF( sdjf%ln_tint ) THEN                 ! time interpolation, shift by 1/2 record
511            !
512            !                  INT( ztmp )
513            !                     /|\
514            !                    1 |    *----
515            !                    0 |----(             
516            !                      |----+----|--> time
517            !                      0   /|\   1   (nday/nmonth_len(nmonth))
518            !                           |   
519            !                           |   
520            !       forcing record :  nmonth
521            !                           
522            ztmp = REAL( nsec_month, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) + 0.5 &
523           &       + REAL( it_offset, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday )
524            imth = nmonth + INT( ztmp ) - COUNT((/llbefore/))
525            IF( sdjf%cltype == 'monthly' ) THEN   ;   sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
526            ELSE                                  ;   sdjf%nrec_a(1) = imth
527            ENDIF
528            sdjf%nrec_a(2) = nmonth_half(   imth ) + nsec1jan000   ! swap at the middle of the month
529         ELSE                                    ! no time interpolation
530            IF( sdjf%cltype == 'monthly' ) THEN   ;   sdjf%nrec_a(1) = 1
531            ELSE                                  ;   sdjf%nrec_a(1) = nmonth
532            ENDIF
533            sdjf%nrec_a(2) =  nmonth_end(nmonth  ) + nsec1jan000   ! swap at the end    of the month
534            sdjf%nrec_b(2) =  nmonth_end(nmonth-1) + nsec1jan000   ! beginning of the month (only for print)
535         ENDIF
536         !
537         !                                   ! ================================ !
538      ELSE                                   ! higher frequency mean (in hours)
539         !                                   ! ================================ !
540         !
541         ifreq_sec = NINT( sdjf%nfreqh * 3600 )                                         ! frequency mean (in seconds)
542         IF( sdjf%cltype(1:4) == 'week' )   isec_week = ksec_week( sdjf%cltype(6:8) )   ! since the first day of the current week
543         ! number of second since the beginning of the file
544         IF(     sdjf%cltype      == 'monthly' ) THEN   ;   ztmp = REAL(nsec_month,wp)  ! since the first day of the current month
545         ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   ztmp = REAL(isec_week ,wp)  ! since the first day of the current week
546         ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   ztmp = REAL(nsec_day  ,wp)  ! since 00h of the current day
547         ELSE                                           ;   ztmp = REAL(nsec_year ,wp)  ! since 00h on Jan 1 of the current year
548         ENDIF
549         ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdttra(1) + REAL( it_offset, wp )  ! centrered in the middle of sbc time step
550         ztmp = ztmp + 0.01 * rdttra(1)                                                 ! avoid truncation error
551         IF( sdjf%ln_tint ) THEN                ! time interpolation, shift by 1/2 record
552            !
553            !          INT( ztmp/ifreq_sec + 0.5 )
554            !                     /|\
555            !                    2 |        *-----(
556            !                    1 |  *-----(
557            !                    0 |--(             
558            !                      |--+--|--+--|--+--|--> time
559            !                      0 /|\ 1 /|\ 2 /|\ 3    (ztmp/ifreq_sec)
560            !                         |     |     |
561            !                         |     |     |
562            !       forcing record :  1     2     3
563            !                   
564            ztmp= ztmp / REAL(ifreq_sec, wp) + 0.5
565         ELSE                                   ! no time interpolation
566            !
567            !           INT( ztmp/ifreq_sec )
568            !                     /|\
569            !                    2 |           *-----(
570            !                    1 |     *-----(
571            !                    0 |-----(             
572            !                      |--+--|--+--|--+--|--> time
573            !                      0 /|\ 1 /|\ 2 /|\ 3    (ztmp/ifreq_sec)
574            !                         |     |     |
575            !                         |     |     |
576            !       forcing record :  1     2     3
577            !                           
578            ztmp= ztmp / REAL(ifreq_sec, wp)
579         ENDIF
580         sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))   ! record number to be read
581
582         iendrec = ifreq_sec * sdjf%nrec_a(1) + nsec1jan000       ! end of this record (in second)
583         ! add the number of seconds between 00h Jan 1 and the end of previous month/week/day (ok if nmonth=1)
584         IF( sdjf%cltype      == 'monthly' )   iendrec = iendrec + NINT(rday) * SUM(nmonth_len(1:nmonth -1))
585         IF( sdjf%cltype(1:4) == 'week'    )   iendrec = iendrec + ( nsec_year - isec_week )
586         IF( sdjf%cltype      == 'daily'   )   iendrec = iendrec + NINT(rday) * ( nday_year - 1 )
587         IF( sdjf%ln_tint ) THEN
588             sdjf%nrec_a(2) = iendrec - ifreq_sec / 2        ! swap at the middle of the record
589         ELSE
590             sdjf%nrec_a(2) = iendrec                        ! swap at the end    of the record
591             sdjf%nrec_b(2) = iendrec - ifreq_sec            ! beginning of the record (only for print)
592         ENDIF
593         !
594      ENDIF
595      !
596   END SUBROUTINE fld_rec
597
598
599   SUBROUTINE fld_get( sdjf, map, jpk_bdy )
600      !!---------------------------------------------------------------------
601      !!                    ***  ROUTINE fld_get  ***
602      !!
603      !! ** Purpose :   read the data
604      !!----------------------------------------------------------------------
605      TYPE(FLD), INTENT(inout)        ::   sdjf    ! input field related variables
606      TYPE(MAP_POINTER),INTENT(in)    ::   map     ! global-to-local mapping indices
607      INTEGER  , INTENT(in), OPTIONAL ::   jpk_bdy ! number of vertical levels in the bdy data
608      !!
609      INTEGER                  ::   ipk    ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
610      INTEGER                  ::   iw     ! index into wgts array
611      INTEGER                  ::   ipdom  ! index of the domain
612      INTEGER                  ::   idvar  ! variable ID
613      INTEGER                  ::   idmspc ! number of spatial dimensions
614      LOGICAL                  ::   lmoor  ! C1D case: point data
615      !!---------------------------------------------------------------------
616      !
617      ipk = SIZE( sdjf%fnow, 3 )
618      !
619      IF( ASSOCIATED(map%ptr) ) THEN
620         IF( PRESENT(jpk_bdy) ) THEN
621            IF( sdjf%ln_tint ) THEN   ;   
622               CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map%ptr, sdjf%igrd, sdjf%ibdy, jpk_bdy )
623            ELSE                      ;   
624               CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1), map%ptr, sdjf%igrd, sdjf%ibdy, jpk_bdy )
625            ENDIF
626         ELSE
627            IF( sdjf%ln_tint ) THEN   ;   
628               CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map%ptr )
629            ELSE                      ;   
630               CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1), map%ptr )
631            ENDIF
632         ENDIF
633      ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
634         CALL wgt_list( sdjf, iw )
635         IF( sdjf%ln_tint ) THEN   ;   
636            CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk  , sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), sdjf%lsmname )
637         ELSE                      ;   
638            CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk  , sdjf%fnow(:,:,:  ), sdjf%nrec_a(1), sdjf%lsmname )
639         ENDIF
640      ELSE
641         IF( SIZE(sdjf%fnow, 1) == jpi ) THEN  ;  ipdom = jpdom_data
642         ELSE                                  ;  ipdom = jpdom_unknown
643         ENDIF
644         ! C1D case: If product of spatial dimensions == ipk, then x,y are of
645         ! size 1 (point/mooring data): this must be read onto the central grid point
646         idvar  = iom_varid( sdjf%num, sdjf%clvar )
647         idmspc = iom_file( sdjf%num )%ndims( idvar )
648         IF( iom_file( sdjf%num )%luld( idvar ) )   idmspc = idmspc - 1
649         lmoor  = (idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk)
650         !
651         SELECT CASE( ipk )
652         CASE(1)
653            IF( lk_c1d .AND. lmoor ) THEN
654               IF( sdjf%ln_tint ) THEN
655                  CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fdta(2,2,1,2), sdjf%nrec_a(1) )
656                  CALL lbc_lnk( sdjf%fdta(:,:,1,2),'Z',1. )
657               ELSE
658                  CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fnow(2,2,1  ), sdjf%nrec_a(1) )
659                  CALL lbc_lnk( sdjf%fnow(:,:,1  ),'Z',1. )
660               ENDIF
661            ELSE
662               IF( sdjf%ln_tint ) THEN   ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,1,2), sdjf%nrec_a(1) )
663               ELSE                      ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,1  ), sdjf%nrec_a(1) )
664               ENDIF
665            ENDIF
666         CASE DEFAULT
667            IF (lk_c1d .AND. lmoor ) THEN
668               IF( sdjf%ln_tint ) THEN
669                  CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fdta(2,2,:,2), sdjf%nrec_a(1) )
670                  CALL lbc_lnk( sdjf%fdta(:,:,:,2),'Z',1. )
671               ELSE
672                  CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fnow(2,2,:  ), sdjf%nrec_a(1) )
673                  CALL lbc_lnk( sdjf%fnow(:,:,:  ),'Z',1. )
674               ENDIF
675            ELSE
676               IF( sdjf%ln_tint ) THEN   ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) )
677               ELSE                      ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1) )
678               ENDIF
679            ENDIF
680         END SELECT
681      ENDIF
682      !
683      sdjf%rotn(2) = .false.   ! vector not yet rotated
684
685   END SUBROUTINE fld_get
686
687   SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy )
688      !!---------------------------------------------------------------------
689      !!                    ***  ROUTINE fld_map  ***
690      !!
691      !! ** Purpose :   read global data from file and map onto local data
692      !!                using a general mapping (for open boundaries)
693      !!----------------------------------------------------------------------
694#if defined key_bdy
695      USE bdy_oce, ONLY:  dta_global, dta_global_z, dta_global2, dta_global2_z         ! workspace to read in global data arrays
696#endif
697      INTEGER                   , INTENT(in ) ::   num     ! stream number
698      CHARACTER(LEN=*)          , INTENT(in ) ::   clvar   ! variable name
699      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta   ! output field on model grid (2 dimensional)
700      INTEGER                   , INTENT(in ) ::   nrec    ! record number to read (ie time slice)
701      INTEGER,  DIMENSION(:)    , INTENT(in ) ::   map     ! global-to-local mapping indices
702      INTEGER  , INTENT(in), OPTIONAL         ::   igrd, ibdy, jpk_bdy  ! grid type, set number and number of vertical levels in the bdy data
703      INTEGER                                 ::   jpkm1_bdy! number of vertical levels in the bdy data minus 1
704      !!
705      INTEGER                                 ::   ipi      ! length of boundary data on local process
706      INTEGER                                 ::   ipj      ! length of dummy dimension ( = 1 )
707      INTEGER                                 ::   ipk      ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk )
708      INTEGER                                 ::   ilendta  ! length of data in file
709      INTEGER                                 ::   idvar    ! variable ID
710      INTEGER                                 ::   ib, ik, ji, jj   ! loop counters
711      INTEGER                                 ::   ierr
712      REAL(wp)                                ::   fv ! fillvalue and alternative -ABS(fv)
713      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read    ! work space for global data
714      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read_z  ! work space for global data
715      !!---------------------------------------------------------------------
716           
717      ipi = SIZE( dta, 1 )
718      ipj = 1
719      ipk = SIZE( dta, 3 )
720
721      idvar   = iom_varid( num, clvar )
722      ilendta = iom_file(num)%dimsz(1,idvar)
723
724#if defined key_bdy
725      ipj = iom_file(num)%dimsz(2,idvar)
726      IF (ipj == 1) THEN
727         dta_read => dta_global
728         IF( PRESENT(jpk_bdy) ) THEN
729            IF( jpk_bdy>0 ) THEN
730               dta_read_z => dta_global_z
731               jpkm1_bdy = jpk_bdy-1
732            ENDIF
733         ENDIF
734      ELSE ! we assume that this is a structured open boundary file
735         dta_read => dta_global2
736         IF( PRESENT(jpk_bdy) ) THEN
737            IF( jpk_bdy>0 ) THEN
738               dta_read_z => dta_global2_z
739               jpkm1_bdy = jpk_bdy-1
740            ENDIF
741         ENDIF
742      ENDIF
743#endif
744
745      IF(lwp) WRITE(numout,*) 'Dim size for ',TRIM(clvar),' is ', ilendta
746      IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk
747
748      SELECT CASE( ipk )
749      CASE(1)        ;   CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1    ), nrec )
750      CASE DEFAULT   ;   
751 
752      IF( PRESENT(jpk_bdy) .AND. jpk_bdy>0 ) THEN       ! boundary data not on model grid: vertical interpolation
753         CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:jpk_bdy), nrec )
754         SELECT CASE( igrd )                 
755            CASE(1)
756               CALL iom_get ( num, jpdom_unknown, 'gdept', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
757               CALL iom_get ( num, jpdom_unknown, 'e3t', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
758            CASE(2) 
759               CALL iom_get ( num, jpdom_unknown, 'gdepu', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
760               CALL iom_get ( num, jpdom_unknown, 'e3u', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
761            CASE(3)
762               CALL iom_get ( num, jpdom_unknown, 'gdepv', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
763               CALL iom_get ( num, jpdom_unknown, 'e3v', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
764         END SELECT
765         CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar )
766#if defined key_bdy
767         CALL fld_bdy_interp(dta_read, dta_read_z, map, jpk_bdy, igrd, ibdy, fv, dta)
768#endif
769      ELSE ! boundary data assumed to be on model grid
770         CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec )                   
771         IF (ipj==1) THEN ! we assume that this is an un-structured open boundary file
772            DO ib = 1, ipi
773              DO ik = 1, ipk
774                dta(ib,1,ik) =  dta_read(map(ib),1,ik)
775              END DO
776            END DO
777         ELSE ! we assume that this is a structured open boundary file
778            DO ib = 1, ipi
779               jj=1+floor(REAL(map(ib)-1)/REAL(ilendta))
780               ji=map(ib)-(jj-1)*ilendta
781               DO ik = 1, ipk
782                  dta(ib,1,ik) =  dta_read(ji,jj,ik)
783               END DO
784            END DO
785         ENDIF
786      ENDIF ! PRESENT(jpk_bdy)
787      END SELECT
788 
789   END SUBROUTINE fld_map
790   
791#if defined key_bdy
792   SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, map, jpk_bdy, igrd, ibdy, fv, dta)
793
794      !!---------------------------------------------------------------------
795      !!                    ***  ROUTINE fld_bdy_interp  ***
796      !!
797      !! ** Purpose :   on the fly vertical interpolation to allow the use of
798      !!                boundary data from non-native vertical grid
799      !!----------------------------------------------------------------------
800      USE bdy_oce, ONLY:  idx_bdy         ! indexing for map <-> ij transformation
801
802      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read    ! work space for global data
803      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_z  ! work space for global data
804      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta        ! output field on model grid (2 dimensional)
805      INTEGER,  DIMENSION(:)    , INTENT(in ) ::   map        ! global-to-local mapping indices
806      INTEGER  , INTENT(in)                   ::   igrd, ibdy, jpk_bdy      ! number of levels in bdy data
807      INTEGER                                 ::   jpkm1_bdy    ! number of levels in bdy data minus 1
808      REAL(wp) , INTENT(in)                                ::   fv ! fillvalue and alternative -ABS(fv)
809      !!
810      INTEGER                                 ::   ipi        ! length of boundary data on local process
811      INTEGER                                 ::   ipj        ! length of dummy dimension ( = 1 )
812      INTEGER                                 ::   ipk, ipkm1 ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk )
813      INTEGER                                 ::   ilendta    ! length of data in file
814      INTEGER                                 ::   ib, ik, ikk! loop counters
815      INTEGER                                 ::   ji, jj ! loop counters
816      REAL(wp)                                ::   zl, zi     ! tmp variable for current depth and interpolation factor
817      REAL(wp)                                ::   fv_alt ! fillvalue and alternative -ABS(fv)
818      !!---------------------------------------------------------------------
819
820      ipi       = SIZE( dta, 1 )
821      ipj       = SIZE( dta_read, 2 )
822      ipk       = SIZE( dta, 3 )
823      ilendta   = SIZE( dta_read, 1 )
824      jpkm1_bdy = jpk_bdy-1
825 
826      fv_alt = -ABS(fv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later
827      !
828      IF (ipj==1) THEN ! we assume that this is an un-structured open boundary file
829         DO ib = 1, ipi
830            DO ik = 1, jpk_bdy
831               IF( ( dta_read(map(ib),1,ik) == fv ) ) THEN
832                  dta_read_z(map(ib),1,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data
833                  dta_read_dz(map(ib),1,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct
834               ENDIF
835      !           dta(ib,1,ik) = fv_alt    ! put fillvalue into new field as if all goes well all wet points will be replaced
836            ENDDO
837         ENDDO 
838      !
839         DO ib = 1, ipi
840            DO ik = 1, ipk                     
841               zl =  gdept_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)   ! if using in step could use fsdept instead of gdept_0?
842               IF( zl < dta_read_z(map(ib),1,1) ) THEN                                         ! above the first level of external data
843                  dta(ib,1,ik) =  dta_read(map(ib),1,1)
844               ELSEIF( zl > MAXVAL(dta_read_z(map(ib),1,:),1) ) THEN                           ! below the last level of external data
845                  dta(ib,1,ik) =  dta_read(map(ib),1,MAXLOC(dta_read_z(map(ib),1,:),1))
846               ELSE                                                                          ! inbetween : vertical interpolation between ikk & ikk+1
847                  DO ikk = 1, jpkm1_bdy                                                          ! when  gdept_0(ikk) < zl < gdept_0(ikk+1)
848                     IF( ( (zl-dta_read_z(map(ib),1,ikk)) * (zl-dta_read_z(map(ib),1,ikk+1)) <= 0._wp)   &
849                    &    .AND. (dta_read_z(map(ib),1,ikk+1) /= fv_alt)) THEN
850                        zi = ( zl - dta_read_z(map(ib),1,ikk) ) / (dta_read_z(map(ib),1,ikk+1)-dta_read_z(map(ib),1,ikk))
851                        dta(ib,1,ik) = dta_read(map(ib),1,ikk) + &
852                      &                ( dta_read(map(ib),1,ikk+1) -  dta_read(map(ib),1,ikk) ) * zi
853                     ENDIF
854                  END DO
855               ENDIF
856            END DO
857         END DO
858      ELSE ! we assume that this is a structured open boundary file
859         DO ib = 1, ipi
860            jj=1+floor(REAL(map(ib)-1)/REAL(ilendta))
861            ji=map(ib)-(jj-1)*ilendta
862            DO ik = 1, jpk_bdy                     
863               IF( ( dta_read(ji,jj,ik) == fv ) ) THEN
864                  dta_read_z(map(ib),1,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data
865               ENDIF
866     !            dta(ib,1,ik) = fv_alt    ! put fillvalue into new field as if all goes well all wet points will be replaced
867            ENDDO
868         ENDDO 
869      !
870         DO ib = 1, ipi
871            jj=1+floor(REAL(map(ib)-1)/REAL(ilendta))
872            ji=map(ib)-(jj-1)*ilendta
873            DO ik = 1, ipk                     
874               zl =  gdept_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)   ! if using in step could use fsdept instead of gdept_0?
875               IF( zl < dta_read_z(ji,jj,1) ) THEN                                         ! above the first level of external data
876                  dta(ib,1,ik) =  dta_read(ji,jj,1)
877               ELSEIF( zl > MAXVAL(dta_read_z(ji,ji,:),1) ) THEN                           ! below the last level of external data
878                  dta(ib,1,ik) =  dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1))
879               ELSE                                                                          ! inbetween : vertical interpolation between ikk & ikk+1
880                  DO ikk = 1, jpkm1_bdy                                                          ! when  gdept_0(ikk) < zl < gdept_0(ikk+1)
881                     IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp)   &
882                    &    .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN
883                        zi = ( zl - dta_read_z(ji,jj,ikk) ) / (dta_read_z(ji,jj,ikk+1)-dta_read_z(ji,jj,ikk))
884                        dta(ib,1,ik) = dta_read(ji,jj,ikk) + &
885                      &                ( dta_read(ji,jj,ikk+1) -  dta_read(ji,jj,ikk) ) * zi
886                     ENDIF
887                  END DO
888               ENDIF
889            END DO
890         END DO
891      ENDIF
892
893   END SUBROUTINE fld_bdy_interp
894   SUBROUTINE fld_bdy_conserve(dta_read, dta_read_z, map, jpk_bdy, igrd, ibdy, fv, dta)
895
896   END SUBROUTINE fld_bdy_conserve
897
898#endif
899
900   SUBROUTINE fld_rot( kt, sd )
901      !!---------------------------------------------------------------------
902      !!                    ***  ROUTINE fld_rot  ***
903      !!
904      !! ** Purpose :   Vector fields may need to be rotated onto the local grid direction
905      !!----------------------------------------------------------------------
906      INTEGER  , INTENT(in   )               ::   kt        ! ocean time step
907      TYPE(FLD), INTENT(inout), DIMENSION(:) ::   sd        ! input field related variables
908      !!
909      INTEGER                           ::   ju,jv,jk,jn  ! loop indices
910      INTEGER                           ::   imf          ! size of the structure sd
911      INTEGER                           ::   ill          ! character length
912      INTEGER                           ::   iv           ! indice of V component
913      REAL(wp), POINTER, DIMENSION(:,:) ::   utmp, vtmp   ! temporary arrays for vector rotation
914      CHARACTER (LEN=100)               ::   clcomp       ! dummy weight name
915      !!---------------------------------------------------------------------
916
917      CALL wrk_alloc( jpi,jpj, utmp, vtmp )
918
919      !! (sga: following code should be modified so that pairs arent searched for each time
920      !
921      imf = SIZE( sd )
922      DO ju = 1, imf
923         ill = LEN_TRIM( sd(ju)%vcomp )
924         DO jn = 2-COUNT((/sd(ju)%ln_tint/)), 2
925            IF( ill > 0 .AND. .NOT. sd(ju)%rotn(jn) ) THEN   ! find vector rotations required             
926               IF( sd(ju)%vcomp(1:1) == 'U' ) THEN      ! east-west component has symbolic name starting with 'U'
927                  ! look for the north-south component which has same symbolic name but with 'U' replaced with 'V'
928                  clcomp = 'V' // sd(ju)%vcomp(2:ill)   ! works even if ill == 1
929                  iv = -1
930                  DO jv = 1, imf
931                     IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) )   iv = jv
932                  END DO
933                  IF( iv > 0 ) THEN   ! fields ju and iv are two components which need to be rotated together
934                     DO jk = 1, SIZE( sd(ju)%fnow, 3 )
935                        IF( sd(ju)%ln_tint )THEN
936                           CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->i', utmp(:,:) )
937                           CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->j', vtmp(:,:) )
938                           sd(ju)%fdta(:,:,jk,jn) = utmp(:,:)   ;   sd(iv)%fdta(:,:,jk,jn) = vtmp(:,:)
939                        ELSE
940                           CALL rot_rep( sd(ju)%fnow(:,:,jk  ), sd(iv)%fnow(:,:,jk  ), 'T', 'en->i', utmp(:,:) )
941                           CALL rot_rep( sd(ju)%fnow(:,:,jk  ), sd(iv)%fnow(:,:,jk  ), 'T', 'en->j', vtmp(:,:) )
942                           sd(ju)%fnow(:,:,jk   ) = utmp(:,:)   ;   sd(iv)%fnow(:,:,jk   ) = vtmp(:,:)
943                        ENDIF
944                     END DO
945                     sd(ju)%rotn(jn) = .TRUE.               ! vector was rotated
946                     IF( lwp .AND. kt == nit000 )   WRITE(numout,*)   &
947                        &   'fld_read: vector pair ('//TRIM(sd(ju)%clvar)//', '//TRIM(sd(iv)%clvar)//') rotated on to model grid'
948                  ENDIF
949               ENDIF
950            ENDIF
951         END DO
952       END DO
953      !
954      CALL wrk_dealloc( jpi,jpj, utmp, vtmp )
955      !
956   END SUBROUTINE fld_rot
957
958
959   SUBROUTINE fld_clopn( sdjf, kyear, kmonth, kday, ldstop )
960      !!---------------------------------------------------------------------
961      !!                    ***  ROUTINE fld_clopn  ***
962      !!
963      !! ** Purpose :   update the file name and open the file
964      !!----------------------------------------------------------------------
965      TYPE(FLD)        , INTENT(inout) ::   sdjf     ! input field related variables
966      INTEGER, OPTIONAL, INTENT(in   ) ::   kyear    ! year value
967      INTEGER, OPTIONAL, INTENT(in   ) ::   kmonth   ! month value
968      INTEGER, OPTIONAL, INTENT(in   ) ::   kday     ! day value
969      LOGICAL, OPTIONAL, INTENT(in   ) ::   ldstop   ! stop if open to read a non-existing file (default = .TRUE.)
970      !!
971      LOGICAL :: llprevyr              ! are we reading previous year  file?
972      LOGICAL :: llprevmth             ! are we reading previous month file?
973      INTEGER :: iyear, imonth, iday   ! first day of the current file in yyyy mm dd
974      INTEGER :: isec_week             ! number of seconds since start of the weekly file
975      INTEGER :: indexyr               ! year undex (O/1/2: previous/current/next)
976      INTEGER :: iyear_len, imonth_len ! length (days) of iyear and imonth             !
977      CHARACTER(len = 256)::   clname  ! temporary file name
978      !!----------------------------------------------------------------------
979      IF( PRESENT(kyear) ) THEN                             ! use given values
980         iyear = kyear
981         imonth = kmonth
982         iday = kday
983      ELSE                                                  ! use current day values
984         IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week
985            isec_week  = ksec_week( sdjf%cltype(6:8) )      ! second since the beginning of the week
986            llprevmth  = isec_week > nsec_month             ! longer time since beginning of the week than the month
987            llprevyr   = llprevmth .AND. nmonth == 1
988         ELSE
989            isec_week  = 0
990            llprevmth  = .FALSE.
991            llprevyr   = .FALSE.
992         ENDIF
993         iyear  = nyear  - COUNT((/llprevyr /))
994         imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
995         iday   = nday   + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
996      ENDIF
997
998      ! build the new filename if not climatological data
999      clname=TRIM(sdjf%clrootname)
1000      !
1001      ! note that sdjf%ln_clim is is only acting on the presence of the year in the file name
1002      IF( .NOT. sdjf%ln_clim ) THEN   
1003                                         WRITE(clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), iyear    ! add year
1004         IF( sdjf%cltype /= 'yearly' )   WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname          ), imonth   ! add month
1005      ELSE
1006         ! build the new filename if climatological data
1007         IF( sdjf%cltype /= 'yearly' )   WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), imonth   ! add month
1008      ENDIF
1009      IF( sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) &
1010            &                            WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname          ), iday     ! add day
1011      !
1012      IF( TRIM(clname) /= TRIM(sdjf%clname) .OR. sdjf%num == 0 ) THEN   ! new file to be open
1013
1014         sdjf%clname = TRIM(clname)
1015         IF( sdjf%num /= 0 )   CALL iom_close( sdjf%num )   ! close file if already open
1016         CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof =  LEN(TRIM(sdjf%wgtname)) > 0 )
1017
1018         ! find the last record to be read -> update sdjf%nreclast
1019         indexyr = iyear - nyear + 1
1020         iyear_len = nyear_len( indexyr )
1021         SELECT CASE ( indexyr )
1022         CASE ( 0 )   ;   imonth_len = 31   ! previous year -> imonth = 12
1023         CASE ( 1 )   ;   imonth_len = nmonth_len(imonth) 
1024         CASE ( 2 )   ;   imonth_len = 31   ! next     year -> imonth = 1
1025         END SELECT
1026         
1027         ! last record to be read in the current file
1028         IF    ( sdjf%nfreqh == -12 ) THEN                 ;   sdjf%nreclast = 1    !  yearly mean
1029         ELSEIF( sdjf%nfreqh ==  -1 ) THEN                                          ! monthly mean
1030            IF(     sdjf%cltype      == 'monthly' ) THEN   ;   sdjf%nreclast = 1
1031            ELSE                                           ;   sdjf%nreclast = 12
1032            ENDIF
1033         ELSE                                                                       ! higher frequency mean (in hours)
1034            IF(     sdjf%cltype      == 'monthly' ) THEN   ;   sdjf%nreclast = NINT( 24 * imonth_len / sdjf%nfreqh )
1035            ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   sdjf%nreclast = NINT( 24 * 7          / sdjf%nfreqh )
1036            ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   sdjf%nreclast = NINT( 24              / sdjf%nfreqh )
1037            ELSE                                           ;   sdjf%nreclast = NINT( 24 * iyear_len  / sdjf%nfreqh )
1038            ENDIF
1039         ENDIF
1040         
1041      ENDIF
1042      !
1043   END SUBROUTINE fld_clopn
1044
1045
1046   SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )
1047      !!---------------------------------------------------------------------
1048      !!                    ***  ROUTINE fld_fill  ***
1049      !!
1050      !! ** Purpose :   fill sdf with sdf_n and control print
1051      !!----------------------------------------------------------------------
1052      TYPE(FLD)  , DIMENSION(:), INTENT(inout) ::   sdf        ! structure of input fields (file informations, fields read)
1053      TYPE(FLD_N), DIMENSION(:), INTENT(in   ) ::   sdf_n      ! array of namelist information structures
1054      CHARACTER(len=*)         , INTENT(in   ) ::   cdir       ! Root directory for location of flx files
1055      CHARACTER(len=*)         , INTENT(in   ) ::   cdcaller   !
1056      CHARACTER(len=*)         , INTENT(in   ) ::   cdtitle    !
1057      CHARACTER(len=*)         , INTENT(in   ) ::   cdnam      !
1058      !
1059      INTEGER  ::   jf       ! dummy indices
1060      !!---------------------------------------------------------------------
1061
1062      DO jf = 1, SIZE(sdf)
1063         sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
1064         sdf(jf)%clname     = "not yet defined"
1065         sdf(jf)%nfreqh     = sdf_n(jf)%nfreqh
1066         sdf(jf)%clvar      = sdf_n(jf)%clvar
1067         sdf(jf)%ln_tint    = sdf_n(jf)%ln_tint
1068         sdf(jf)%ln_clim    = sdf_n(jf)%ln_clim
1069         sdf(jf)%cltype     = sdf_n(jf)%cltype
1070         sdf(jf)%num        = -1
1071         sdf(jf)%wgtname    = " "
1072         IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
1073         sdf(jf)%lsmname = " "
1074         IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 )   sdf(jf)%lsmname = TRIM( cdir )//TRIM( sdf_n(jf)%lname )
1075         sdf(jf)%vcomp      = sdf_n(jf)%vcomp
1076         sdf(jf)%rotn(:)    = .TRUE.   ! pretend to be rotated -> won't try to rotate data before the first call to fld_get
1077         IF( sdf(jf)%cltype(1:4) == 'week' .AND. nn_leapy == 0  )   &
1078            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs nn_leapy = 1')
1079         IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim )   &
1080            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.')
1081         sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
1082         sdf(jf)%igrd     = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during bdy_dta_init
1083         sdf(jf)%ibdy     = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during bdy_dta_init
1084      END DO
1085
1086      IF(lwp) THEN      ! control print
1087         WRITE(numout,*)
1088         WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
1089         WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
1090         WRITE(numout,*) '          '//TRIM( cdnam )//' Namelist'
1091         WRITE(numout,*) '          list of files and frequency (>0: in hours ; <0 in months)'
1092         DO jf = 1, SIZE(sdf)
1093            WRITE(numout,*) '               root filename: '  , TRIM( sdf(jf)%clrootname ),   &
1094               &                          ' variable name: '  , TRIM( sdf(jf)%clvar      )
1095            WRITE(numout,*) '               frequency: '      ,       sdf(jf)%nfreqh      ,   &
1096               &                          ' time interp: '    ,       sdf(jf)%ln_tint     ,   &
1097               &                          ' climatology: '    ,       sdf(jf)%ln_clim     ,   &
1098               &                          ' weights    : '    , TRIM( sdf(jf)%wgtname    ),   &
1099               &                          ' pairing    : '    , TRIM( sdf(jf)%vcomp      ),   &
1100               &                          ' data type: '      ,       sdf(jf)%cltype      ,   &
1101               &                          ' land/sea mask:'   , TRIM( sdf(jf)%lsmname    )
1102            call flush(numout)
1103         END DO
1104      ENDIF
1105     
1106   END SUBROUTINE fld_fill
1107
1108
1109   SUBROUTINE wgt_list( sd, kwgt )
1110      !!---------------------------------------------------------------------
1111      !!                    ***  ROUTINE wgt_list  ***
1112      !!
1113      !! ** Purpose :   search array of WGTs and find a weights file
1114      !!                entry, or return a new one adding it to the end
1115      !!                if it is a new entry, the weights data is read in and
1116      !!                restructured (fld_weight)
1117      !!----------------------------------------------------------------------
1118      TYPE( FLD ), INTENT(in   ) ::   sd        ! field with name of weights file
1119      INTEGER    , INTENT(inout) ::   kwgt      ! index of weights
1120      !!
1121      INTEGER ::   kw, nestid   ! local integer
1122      LOGICAL ::   found        ! local logical
1123      !!----------------------------------------------------------------------
1124      !
1125      !! search down linked list
1126      !! weights filename is either present or we hit the end of the list
1127      found = .FALSE.
1128
1129      !! because agrif nest part of filenames are now added in iom_open
1130      !! to distinguish between weights files on the different grids, need to track
1131      !! nest number explicitly
1132      nestid = 0
1133#if defined key_agrif
1134      nestid = Agrif_Fixed()
1135#endif
1136      DO kw = 1, nxt_wgt-1
1137         IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
1138             ref_wgts(kw)%nestid == nestid) THEN
1139            kwgt = kw
1140            found = .TRUE.
1141            EXIT
1142         ENDIF
1143      END DO
1144      IF( .NOT.found ) THEN
1145         kwgt = nxt_wgt
1146         CALL fld_weight( sd )
1147      ENDIF
1148      !
1149   END SUBROUTINE wgt_list
1150
1151
1152   SUBROUTINE wgt_print( )
1153      !!---------------------------------------------------------------------
1154      !!                    ***  ROUTINE wgt_print  ***
1155      !!
1156      !! ** Purpose :   print the list of known weights
1157      !!----------------------------------------------------------------------
1158      INTEGER ::   kw   !
1159      !!----------------------------------------------------------------------
1160      !
1161      DO kw = 1, nxt_wgt-1
1162         WRITE(numout,*) 'weight file:  ',TRIM(ref_wgts(kw)%wgtname)
1163         WRITE(numout,*) '      ddims:  ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
1164         WRITE(numout,*) '     numwgt:  ',ref_wgts(kw)%numwgt
1165         WRITE(numout,*) '     jpiwgt:  ',ref_wgts(kw)%jpiwgt
1166         WRITE(numout,*) '     jpjwgt:  ',ref_wgts(kw)%jpjwgt
1167         WRITE(numout,*) '    botleft:  ',ref_wgts(kw)%botleft
1168         WRITE(numout,*) '   topright:  ',ref_wgts(kw)%topright
1169         IF( ref_wgts(kw)%cyclic ) THEN
1170            WRITE(numout,*) '       cyclical'
1171            IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) '              with overlap of ', ref_wgts(kw)%overlap
1172         ELSE
1173            WRITE(numout,*) '       not cyclical'
1174         ENDIF
1175         IF( ASSOCIATED(ref_wgts(kw)%data_wgt) )  WRITE(numout,*) '       allocated'
1176      END DO
1177      !
1178   END SUBROUTINE wgt_print
1179
1180
1181   SUBROUTINE fld_weight( sd )
1182      !!---------------------------------------------------------------------
1183      !!                    ***  ROUTINE fld_weight  ***
1184      !!
1185      !! ** Purpose :   create a new WGT structure and fill in data from 
1186      !!                file, restructuring as required
1187      !!----------------------------------------------------------------------
1188      TYPE( FLD ), INTENT(in) ::   sd   ! field with name of weights file
1189      !!
1190      INTEGER                           ::   jn            ! dummy loop indices
1191      INTEGER                           ::   inum          ! temporary logical unit
1192      INTEGER                           ::   id            ! temporary variable id
1193      INTEGER                           ::   ipk           ! temporary vertical dimension
1194      CHARACTER (len=5)                 ::   aname
1195      INTEGER , DIMENSION(3)            ::   ddims
1196      INTEGER , POINTER, DIMENSION(:,:) ::   data_src
1197      REAL(wp), POINTER, DIMENSION(:,:) ::   data_tmp
1198      LOGICAL                           ::   cyclical
1199      INTEGER                           ::   zwrap      ! local integer
1200      !!----------------------------------------------------------------------
1201      !
1202      CALL wrk_alloc( jpi,jpj, data_src )   ! integer
1203      CALL wrk_alloc( jpi,jpj, data_tmp )
1204      !
1205      IF( nxt_wgt > tot_wgts ) THEN
1206        CALL ctl_stop("fld_weight: weights array size exceeded, increase tot_wgts")
1207      ENDIF
1208      !
1209      !! new weights file entry, add in extra information
1210      !! a weights file represents a 2D grid of a certain shape, so we assume that the current
1211      !! input data file is representative of all other files to be opened and processed with the
1212      !! current weights file
1213
1214      !! open input data file (non-model grid)
1215      CALL iom_open( sd%clname, inum, ldiof =  LEN(TRIM(sd%wgtname)) > 0 )
1216
1217      !! get dimensions
1218      id = iom_varid( inum, sd%clvar, ddims )
1219
1220      !! close it
1221      CALL iom_close( inum )
1222
1223      !! now open the weights file
1224
1225      CALL iom_open ( sd%wgtname, inum )   ! interpolation weights
1226      IF ( inum > 0 ) THEN
1227
1228         !! determine whether we have an east-west cyclic grid
1229         !! from global attribute called "ew_wrap" in the weights file
1230         !! note that if not found, iom_getatt returns -999 and cyclic with no overlap is assumed
1231         !! since this is the most common forcing configuration
1232
1233         CALL iom_getatt(inum, 'ew_wrap', zwrap)
1234         IF( zwrap >= 0 ) THEN
1235            cyclical = .TRUE.
1236         ELSE IF( zwrap == -999 ) THEN
1237            cyclical = .TRUE.
1238            zwrap = 0
1239         ELSE
1240            cyclical = .FALSE.
1241         ENDIF
1242
1243         ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
1244         ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
1245         ref_wgts(nxt_wgt)%wgtname = sd%wgtname
1246         ref_wgts(nxt_wgt)%overlap = zwrap
1247         ref_wgts(nxt_wgt)%cyclic = cyclical
1248         ref_wgts(nxt_wgt)%nestid = 0
1249#if defined key_agrif
1250         ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
1251#endif
1252         !! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
1253         !! for each weight wgtNN there is an integer array srcNN which gives the point in
1254         !! the input data grid which is to be multiplied by the weight
1255         !! they are both arrays on the model grid so the result of the multiplication is
1256         !! added into an output array on the model grid as a running sum
1257
1258         !! two possible cases: bilinear (4 weights) or bicubic (16 weights)
1259         id = iom_varid(inum, 'src05', ldstop=.FALSE.)
1260         IF( id <= 0) THEN
1261            ref_wgts(nxt_wgt)%numwgt = 4
1262         ELSE
1263            ref_wgts(nxt_wgt)%numwgt = 16
1264         ENDIF
1265
1266         ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
1267         ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
1268         ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
1269
1270         DO jn = 1,4
1271            aname = ' '
1272            WRITE(aname,'(a3,i2.2)') 'src',jn
1273            data_tmp(:,:) = 0
1274            CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
1275            data_src(:,:) = INT(data_tmp(:,:))
1276            ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
1277            ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
1278         END DO
1279
1280         DO jn = 1, ref_wgts(nxt_wgt)%numwgt
1281            aname = ' '
1282            WRITE(aname,'(a3,i2.2)') 'wgt',jn
1283            ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
1284            CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
1285         END DO
1286         CALL iom_close (inum)
1287 
1288         ! find min and max indices in grid
1289         ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1290         ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1291         ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1292         ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1293
1294         ! and therefore dimensions of the input box
1295         ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
1296         ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
1297
1298         ! shift indexing of source grid
1299         ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
1300         ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
1301
1302         ! create input grid, give it a halo to allow gradient calculations
1303         ! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
1304         ! a more robust solution will be given in next release
1305         ipk =  SIZE(sd%fnow, 3)
1306         ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
1307         IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
1308
1309         nxt_wgt = nxt_wgt + 1
1310
1311      ELSE
1312         CALL ctl_stop( '    fld_weight : unable to read the file ' )
1313      ENDIF
1314
1315      CALL wrk_dealloc( jpi,jpj, data_src )   ! integer
1316      CALL wrk_dealloc( jpi,jpj, data_tmp )
1317      !
1318   END SUBROUTINE fld_weight
1319
1320
1321   SUBROUTINE apply_seaoverland(clmaskfile,zfieldo,jpi1_lsm,jpi2_lsm,jpj1_lsm, &
1322                          &      jpj2_lsm,itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
1323      !!---------------------------------------------------------------------
1324      !!                    ***  ROUTINE apply_seaoverland  ***
1325      !!
1326      !! ** Purpose :   avoid spurious fluxes in coastal or near-coastal areas
1327      !!                due to the wrong usage of "land" values from the coarse
1328      !!                atmospheric model when spatial interpolation is required
1329      !!      D. Delrosso INGV         
1330      !!----------------------------------------------------------------------
1331      INTEGER                                   :: inum,jni,jnj,jnz,jc                  ! temporary indices
1332      INTEGER,                   INTENT(in)     :: itmpi,itmpj,itmpz                    ! lengths
1333      INTEGER,                   INTENT(in)     :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm  ! temporary indices
1334      INTEGER, DIMENSION(3),     INTENT(in)     :: rec1_lsm,recn_lsm                    ! temporary arrays for start and length
1335      REAL(wp),DIMENSION (:,:,:),INTENT(inout)  :: zfieldo                              ! input/output array for seaoverland application
1336      REAL(wp),DIMENSION (:,:,:),ALLOCATABLE    :: zslmec1                              ! temporary array for land point detection
1337      REAL(wp),DIMENSION (:,:),  ALLOCATABLE    :: zfieldn                              ! array of forcing field with undeff for land points
1338      REAL(wp),DIMENSION (:,:),  ALLOCATABLE    :: zfield                               ! array of forcing field
1339      CHARACTER (len=100),       INTENT(in)     :: clmaskfile                           ! land/sea mask file name
1340      !!---------------------------------------------------------------------
1341      ALLOCATE ( zslmec1(itmpi,itmpj,itmpz) )
1342      ALLOCATE ( zfieldn(itmpi,itmpj) )
1343      ALLOCATE ( zfield(itmpi,itmpj) )
1344
1345      ! Retrieve the land sea mask data
1346      CALL iom_open( clmaskfile, inum )
1347      SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1348      CASE(1)
1349           CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm)
1350      CASE DEFAULT
1351           CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm)
1352      END SELECT
1353      CALL iom_close( inum )
1354
1355      DO jnz=1,rec1_lsm(3)                            !! Loop over k dimension
1356
1357         DO jni=1,itmpi                               !! copy the original field into a tmp array
1358            DO jnj=1,itmpj                            !! substituting undeff over land points
1359            zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz)
1360               IF ( zslmec1(jni,jnj,jnz) == 1. ) THEN
1361                  zfieldn(jni,jnj) = undeff_lsm
1362               ENDIF
1363            END DO
1364         END DO
1365 
1366      CALL seaoverland(zfieldn,itmpi,itmpj,zfield)
1367      DO jc=1,nn_lsm
1368         CALL seaoverland(zfield,itmpi,itmpj,zfield)
1369      END DO
1370
1371      !   Check for Undeff and substitute original values
1372      IF(ANY(zfield==undeff_lsm)) THEN
1373         DO jni=1,itmpi
1374            DO jnj=1,itmpj
1375               IF (zfield(jni,jnj)==undeff_lsm) THEN
1376                  zfield(jni,jnj) = zfieldo(jni,jnj,jnz)
1377               ENDIF
1378            ENDDO
1379         ENDDO
1380      ENDIF
1381
1382      zfieldo(:,:,jnz)=zfield(:,:)
1383
1384      END DO                          !! End Loop over k dimension
1385
1386      DEALLOCATE ( zslmec1 )
1387      DEALLOCATE ( zfieldn )
1388      DEALLOCATE ( zfield )
1389
1390   END SUBROUTINE apply_seaoverland 
1391
1392
1393   SUBROUTINE seaoverland(zfieldn,ileni,ilenj,zfield)
1394      !!---------------------------------------------------------------------
1395      !!                    ***  ROUTINE seaoverland  ***
1396      !!
1397      !! ** Purpose :   create shifted matrices for seaoverland application 
1398      !!      D. Delrosso INGV
1399      !!----------------------------------------------------------------------
1400      INTEGER,INTENT(in)                       :: ileni,ilenj              ! lengths
1401      REAL,DIMENSION (ileni,ilenj),INTENT(in)  :: zfieldn                  ! array of forcing field with undeff for land points
1402      REAL,DIMENSION (ileni,ilenj),INTENT(out) :: zfield                   ! array of forcing field
1403      REAL,DIMENSION (ileni,ilenj)             :: zmat1,zmat2,zmat3,zmat4  ! temporary arrays for seaoverland application
1404      REAL,DIMENSION (ileni,ilenj)             :: zmat5,zmat6,zmat7,zmat8  ! temporary arrays for seaoverland application
1405      REAL,DIMENSION (ileni,ilenj)             :: zlsm2d                   ! temporary arrays for seaoverland application
1406      REAL,DIMENSION (ileni,ilenj,8)           :: zlsm3d                   ! temporary arrays for seaoverland application
1407      LOGICAL,DIMENSION (ileni,ilenj,8)        :: ll_msknan3d              ! logical mask for undeff detection
1408      LOGICAL,DIMENSION (ileni,ilenj)          :: ll_msknan2d              ! logical mask for undeff detection
1409      !!----------------------------------------------------------------------
1410      zmat8 = eoshift(zfieldn   ,  SHIFT=-1, BOUNDARY = (/zfieldn(:,1)/)    ,DIM=2)
1411      zmat1 = eoshift(zmat8     ,  SHIFT=-1, BOUNDARY = (/zmat8(1,:)/)      ,DIM=1)
1412      zmat2 = eoshift(zfieldn   ,  SHIFT=-1, BOUNDARY = (/zfieldn(1,:)/)    ,DIM=1)
1413      zmat4 = eoshift(zfieldn   ,  SHIFT= 1, BOUNDARY = (/zfieldn(:,ilenj)/),DIM=2)
1414      zmat3 = eoshift(zmat4     ,  SHIFT=-1, BOUNDARY = (/zmat4(1,:)/)      ,DIM=1)
1415      zmat5 = eoshift(zmat4     ,  SHIFT= 1, BOUNDARY = (/zmat4(ileni,:)/)  ,DIM=1)
1416      zmat6 = eoshift(zfieldn   ,  SHIFT= 1, BOUNDARY = (/zfieldn(ileni,:)/),DIM=1)
1417      zmat7 = eoshift(zmat8     ,  SHIFT= 1, BOUNDARY = (/zmat8(ileni,:)/)  ,DIM=1)
1418
1419      zlsm3d  = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /))
1420      ll_msknan3d = .not.(zlsm3d==undeff_lsm)
1421      ll_msknan2d = .not.(zfieldn==undeff_lsm)  ! FALSE where is Undeff (land)
1422      zlsm2d = (SUM ( zlsm3d, 3 , ll_msknan3d ) )/(MAX(1,(COUNT( ll_msknan3d , 3 ))   ))
1423      WHERE ((COUNT( ll_msknan3d , 3 )) == 0.0_wp)  zlsm2d = undeff_lsm
1424      zfield = MERGE (zfieldn,zlsm2d,ll_msknan2d)
1425   END SUBROUTINE seaoverland
1426
1427
1428   SUBROUTINE fld_interp( num, clvar, kw, kk, dta,  &
1429                          &         nrec, lsmfile)     
1430      !!---------------------------------------------------------------------
1431      !!                    ***  ROUTINE fld_interp  ***
1432      !!
1433      !! ** Purpose :   apply weights to input gridded data to create data
1434      !!                on model grid
1435      !!----------------------------------------------------------------------
1436      INTEGER                   , INTENT(in   ) ::   num     ! stream number
1437      CHARACTER(LEN=*)          , INTENT(in   ) ::   clvar   ! variable name
1438      INTEGER                   , INTENT(in   ) ::   kw      ! weights number
1439      INTEGER                   , INTENT(in   ) ::   kk      ! vertical dimension of kk
1440      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   dta     ! output field on model grid
1441      INTEGER                   , INTENT(in   ) ::   nrec    ! record number to read (ie time slice)
1442      CHARACTER(LEN=*)          , INTENT(in   ) ::   lsmfile ! land sea mask file name
1443      !!
1444      REAL(wp),DIMENSION(:,:,:),ALLOCATABLE     ::   ztmp_fly_dta,zfieldo                  ! temporary array of values on input grid
1445      INTEGER, DIMENSION(3)                     ::   rec1,recn                             ! temporary arrays for start and length
1446      INTEGER, DIMENSION(3)                     ::   rec1_lsm,recn_lsm                     ! temporary arrays for start and length in case of seaoverland
1447      INTEGER                                   ::   ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2       ! temporary indices
1448      INTEGER                                   ::   jk, jn, jm, jir, jjr                  ! loop counters
1449      INTEGER                                   ::   ni, nj                                ! lengths
1450      INTEGER                                   ::   jpimin,jpiwid                         ! temporary indices
1451      INTEGER                                   ::   jpimin_lsm,jpiwid_lsm                 ! temporary indices
1452      INTEGER                                   ::   jpjmin,jpjwid                         ! temporary indices
1453      INTEGER                                   ::   jpjmin_lsm,jpjwid_lsm                 ! temporary indices
1454      INTEGER                                   ::   jpi1,jpi2,jpj1,jpj2                   ! temporary indices
1455      INTEGER                                   ::   jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm   ! temporary indices
1456      INTEGER                                   ::   itmpi,itmpj,itmpz                     ! lengths
1457     
1458      !!----------------------------------------------------------------------
1459      !
1460      !! for weighted interpolation we have weights at four corners of a box surrounding
1461      !! a model grid point, each weight is multiplied by a grid value (bilinear case)
1462      !! or by a grid value and gradients at the corner point (bicubic case)
1463      !! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
1464
1465      !! sub grid from non-model input grid which encloses all grid points in this nemo process
1466      jpimin = ref_wgts(kw)%botleft(1)
1467      jpjmin = ref_wgts(kw)%botleft(2)
1468      jpiwid = ref_wgts(kw)%jpiwgt
1469      jpjwid = ref_wgts(kw)%jpjwgt
1470
1471      !! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients
1472      rec1(1) = MAX( jpimin-1, 1 )
1473      rec1(2) = MAX( jpjmin-1, 1 )
1474      rec1(3) = 1
1475      recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
1476      recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1477      recn(3) = kk
1478
1479      !! where we need to put it in the non-nemo grid fly_dta
1480      !! note that jpi1 and jpj1 only differ from 1 when jpimin and jpjmin are 1
1481      !! (ie at the extreme west or south of the whole input grid) and similarly for jpi2 and jpj2
1482      jpi1 = 2 + rec1(1) - jpimin
1483      jpj1 = 2 + rec1(2) - jpjmin
1484      jpi2 = jpi1 + recn(1) - 1
1485      jpj2 = jpj1 + recn(2) - 1
1486
1487
1488      IF( LEN( TRIM(lsmfile) ) > 0 ) THEN
1489      !! indeces for ztmp_fly_dta
1490      ! --------------------------
1491         rec1_lsm(1)=MAX(rec1(1)-nn_lsm,1)  ! starting index for enlarged external data, x direction
1492         rec1_lsm(2)=MAX(rec1(2)-nn_lsm,1)  ! starting index for enlarged external data, y direction
1493         rec1_lsm(3) = 1                    ! vertical dimension
1494         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
1495         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
1496         recn_lsm(3) = kk                   ! number of vertical levels in the input file
1497
1498      !  Avoid out of bound
1499         jpimin_lsm = MAX( rec1_lsm(1)+1, 1 )
1500         jpjmin_lsm = MAX( rec1_lsm(2)+1, 1 )
1501         jpiwid_lsm = MIN( recn_lsm(1)-2,ref_wgts(kw)%ddims(1)-rec1(1)+1)
1502         jpjwid_lsm = MIN( recn_lsm(2)-2,ref_wgts(kw)%ddims(2)-rec1(2)+1)
1503
1504         jpi1_lsm = 2+rec1_lsm(1)-jpimin_lsm
1505         jpj1_lsm = 2+rec1_lsm(2)-jpjmin_lsm
1506         jpi2_lsm = jpi1_lsm + recn_lsm(1) - 1
1507         jpj2_lsm = jpj1_lsm + recn_lsm(2) - 1
1508
1509
1510         itmpi=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),1)
1511         itmpj=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),2)
1512         itmpz=kk
1513         ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz))
1514         ztmp_fly_dta(:,:,:) = 0.0
1515         SELECT CASE( SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1516         CASE(1)
1517              CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1),   &
1518                 &                                                                nrec, rec1_lsm, recn_lsm)
1519         CASE DEFAULT
1520              CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),   &
1521                 &                                                                nrec, rec1_lsm, recn_lsm)
1522         END SELECT
1523         CALL apply_seaoverland(lsmfile,ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),                  &
1524                 &                                      jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm,                  &
1525                 &                                      itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
1526
1527
1528         ! Relative indeces for remapping
1529         ii_lsm1 = (rec1(1)-rec1_lsm(1))+1
1530         ii_lsm2 = (ii_lsm1+recn(1))-1
1531         ij_lsm1 = (rec1(2)-rec1_lsm(2))+1
1532         ij_lsm2 = (ij_lsm1+recn(2))-1
1533
1534         ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1535         ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:) = ztmp_fly_dta(ii_lsm1:ii_lsm2,ij_lsm1:ij_lsm2,:)
1536         DEALLOCATE(ztmp_fly_dta)
1537
1538      ELSE
1539         
1540         ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1541         SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) )
1542         CASE(1)
1543              CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn)
1544         CASE DEFAULT
1545              CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn)
1546         END SELECT
1547      ENDIF
1548     
1549
1550      !! first four weights common to both bilinear and bicubic
1551      !! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft
1552      !! note that we have to offset by 1 into fly_dta array because of halo
1553      dta(:,:,:) = 0.0
1554      DO jk = 1,4
1555        DO jn = 1, jpj
1556          DO jm = 1,jpi
1557            ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1558            nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1559            dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,:)
1560          END DO
1561        END DO
1562      END DO
1563
1564      IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
1565
1566        !! fix up halo points that we couldnt read from file
1567        IF( jpi1 == 2 ) THEN
1568           ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:)
1569        ENDIF
1570        IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1571           ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:)
1572        ENDIF
1573        IF( jpj1 == 2 ) THEN
1574           ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:)
1575        ENDIF
1576        IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
1577           ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:)
1578        ENDIF
1579
1580        !! if data grid is cyclic we can do better on east-west edges
1581        !! but have to allow for whether first and last columns are coincident
1582        IF( ref_wgts(kw)%cyclic ) THEN
1583           rec1(2) = MAX( jpjmin-1, 1 )
1584           recn(1) = 1
1585           recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1586           jpj1 = 2 + rec1(2) - jpjmin
1587           jpj2 = jpj1 + recn(2) - 1
1588           IF( jpi1 == 2 ) THEN
1589              rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap
1590              SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1591              CASE(1)
1592                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1593              CASE DEFAULT
1594                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1595              END SELECT     
1596              ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1597           ENDIF
1598           IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1599              rec1(1) = 1 + ref_wgts(kw)%overlap
1600              SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1601              CASE(1)
1602                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1603              CASE DEFAULT
1604                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1605              END SELECT
1606              ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1607           ENDIF
1608        ENDIF
1609
1610        ! gradient in the i direction
1611        DO jk = 1,4
1612          DO jn = 1, jpj
1613            DO jm = 1,jpi
1614              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1615              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1616              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 *         &
1617                               (ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:))
1618            END DO
1619          END DO
1620        END DO
1621
1622        ! gradient in the j direction
1623        DO jk = 1,4
1624          DO jn = 1, jpj
1625            DO jm = 1,jpi
1626              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1627              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1628              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 *         &
1629                               (ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:))
1630            END DO
1631          END DO
1632        END DO
1633
1634         ! gradient in the ij direction
1635         DO jk = 1,4
1636            DO jn = 1, jpj
1637               DO jm = 1,jpi
1638                  ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1639                  nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1640                  dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
1641                               (ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni  ,nj+2,:)) -   &
1642                               (ref_wgts(kw)%fly_dta(ni+2,nj  ,:) - ref_wgts(kw)%fly_dta(ni  ,nj  ,:)))
1643               END DO
1644            END DO
1645         END DO
1646         !
1647      END IF
1648      !
1649   END SUBROUTINE fld_interp
1650
1651
1652   FUNCTION ksec_week( cdday )
1653      !!---------------------------------------------------------------------
1654      !!                    ***  FUNCTION kshift_week ***
1655      !!
1656      !! ** Purpose : 
1657      !!---------------------------------------------------------------------
1658      CHARACTER(len=*), INTENT(in)   ::   cdday   !3 first letters of the first day of the weekly file
1659      !!
1660      INTEGER                        ::   ksec_week  ! output variable
1661      INTEGER                        ::   ijul       !temp variable
1662      INTEGER                        ::   ishift     !temp variable
1663      CHARACTER(len=3),DIMENSION(7)  ::   cl_week 
1664      !!----------------------------------------------------------------------
1665      cl_week = (/"sun","sat","fri","thu","wed","tue","mon"/)
1666      DO ijul = 1, 7
1667         IF( cl_week(ijul) == TRIM(cdday) ) EXIT
1668      END DO
1669      IF( ijul .GT. 7 )   CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): '//TRIM(cdday) )
1670      !
1671      ishift = ijul * NINT(rday)
1672      !
1673      ksec_week = nsec_week + ishift
1674      ksec_week = MOD( ksec_week, 7*NINT(rday) )
1675      !
1676   END FUNCTION ksec_week
1677
1678   !!======================================================================
1679END MODULE fldread
Note: See TracBrowser for help on using the repository browser.