source: NEMO/trunk/src/OCE/SBC/fldread.F90 @ 12248

Last change on this file since 12248 was 12248, checked in by smasson, 11 months ago

trunk : clean useless warnings in fldread

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