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

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

bug fix to fldread to allow code to run with bt vels properly

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