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

Last change on this file since 11536 was 11536, checked in by smasson, 12 months ago

trunk: merge dev_r10984_HPC-13 into the trunk

  • Property svn:keywords set to Id
File size: 101.3 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         ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary?
868         SELECT CASE( kgrd )                         
869         CASE(1)
870            IF( ABS( (zh - ht_n(ji,jj)) / ht_n(ji,jj)) * tmask(ji,jj,1) > 0.01_wp ) THEN
871               WRITE(ctmp1,"(I10.10)") jb 
872               CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%')
873               !   IF(lwp) WRITE(numout,*) 'DEPTHT', zh, sum(e3t_n(ji,jj,:), mask=tmask(ji,jj,:)==1),  ht_n(ji,jj), jb, jb, ji, jj
874            ENDIF
875         CASE(2)
876            IF( ABS( (zh - hu_n(ji,jj)) * r1_hu_n(ji,jj)) * umask(ji,jj,1) > 0.01_wp ) THEN
877               WRITE(ctmp1,"(I10.10)") jb 
878               CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%')
879               !   IF(lwp) WRITE(numout,*) 'DEPTHU', zh, SUM(e3u_n(ji,jj,:), mask=umask(ji,jj,:)==1),  SUM(umask(ji,jj,:)), &
880               !      &                hu_n(ji,jj), jb, jb, ji, jj, narea-1, pdta_read(jb,1,:)
881            ENDIF
882         CASE(3)
883            IF( ABS( (zh - hv_n(ji,jj)) * r1_hv_n(ji,jj)) * vmask(ji,jj,1) > 0.01_wp ) THEN
884               WRITE(ctmp1,"(I10.10)") jb
885               CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%')
886            ENDIF
887         END SELECT
888         !
889         SELECT CASE( kgrd )                         
890         CASE(1)
891            ! depth of T points:
892            zdepth(:) = gdept_n(ji,jj,:)
893         CASE(2)
894            ! depth of U points: we must not use gdept_n as we don't want to do a communication
895            !   --> copy what is done for gdept_n in domvvl...
896            zdhalf(1) = 0.0_wp
897            zdepth(1) = 0.5_wp * e3uw_n(ji,jj,1)
898            DO jk = 2, jpk                               ! vertical sum
899               !    zcoef = umask - wumask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
900               !                              ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
901               !                              ! 0.5 where jk = mikt     
902               !!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
903               zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) )
904               zdhalf(jk) = zdhalf(jk-1) + e3u_n(ji,jj,jk-1)
905               zdepth(jk) =      zcoef  * ( zdhalf(jk  ) + 0.5 * e3uw_n(ji,jj,jk))  &
906                  &         + (1-zcoef) * ( zdepth(jk-1) +       e3uw_n(ji,jj,jk))
907            END DO
908         CASE(3)
909            ! depth of V points: we must not use gdept_n as we don't want to do a communication
910            !   --> copy what is done for gdept_n in domvvl...
911            zdhalf(1) = 0.0_wp
912            zdepth(1) = 0.5_wp * e3vw_n(ji,jj,1)
913            DO jk = 2, jpk                               ! vertical sum
914               !    zcoef = vmask - wvmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
915               !                              ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
916               !                              ! 0.5 where jk = mikt     
917               !!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
918               zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) )
919               zdhalf(jk) = zdhalf(jk-1) + e3v_n(ji,jj,jk-1)
920               zdepth(jk) =      zcoef  * ( zdhalf(jk  ) + 0.5 * e3vw_n(ji,jj,jk))  &
921                  &         + (1-zcoef) * ( zdepth(jk-1) +       e3vw_n(ji,jj,jk))
922            END DO
923         END SELECT
924         !
925         DO jk = 1, jpk                     
926            IF(     zdepth(jk) < pdta_read_z(jb,1,          1) ) THEN                ! above the first level of external data
927               pdta(jb,1,jk) =  pdta_read(jb,1,1)
928            ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkb) ) THEN                       ! below the last level of external data
929               pdta(jb,1,jk) =  pdta_read(jb,1,MAXLOC(pdta_read_z(jb,1,:),1))
930            ELSE                                                             ! inbetween: vertical interpolation between jkb & jkb+1
931               DO jkb = 1, ipkb-1                                            ! when  gdept_n(jkb) < zdepth(jk) < gdept_n(jkb+1)
932                  IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) <= 0._wp ) &
933                     &    .AND. ( pdta_read_z(jb,1,jkb+1) /= zfv_alt) ) THEN   ! linear interpolation between 2 levels
934                     zi = ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) / ( pdta_read_z(jb,1,jkb+1) - pdta_read_z(jb,1,jkb) )
935                     pdta(jb,1,jk) = pdta_read(jb,1,jkb) + ( pdta_read  (jb,1,jkb+1) - pdta_read  (jb,1,jkb) ) * zi
936                  ENDIF
937               END DO
938            ENDIF
939         END DO   ! jpk
940         !
941      END DO   ! ipi
942     
943      IF(kgrd == 2) THEN                                  ! do we need to adjust the transport term?
944         DO jb = 1, ipi
945            ji = idx_bdy(kbdy)%nbi(jb,kgrd)
946            jj = idx_bdy(kbdy)%nbj(jb,kgrd)
947            zh  = SUM(pdta_read_dz(jb,1,:) )
948            ztrans = 0._wp
949            ztrans_new = 0._wp
950            DO jkb = 1, ipkb                              ! calculate transport on input grid
951               ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)
952            ENDDO
953            DO jk = 1, jpk                                ! calculate transport on model grid
954               ztrans_new = ztrans_new +      pdta(jb,1,jk ) *        e3u_n(ji,jj,jk ) * umask(ji,jj,jk)
955            ENDDO
956            DO jk = 1, jpk                                ! make transport correction
957               IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data
958                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu_n(ji,jj) ) * umask(ji,jj,jk)
959               ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
960                  pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hu_n(ji,jj)   * umask(ji,jj,jk)
961               ENDIF
962            ENDDO
963         ENDDO
964      ENDIF
965     
966      IF(kgrd == 3) THEN                                  ! do we need to adjust the transport term?
967         DO jb = 1, ipi
968            ji = idx_bdy(kbdy)%nbi(jb,kgrd)
969            jj = idx_bdy(kbdy)%nbj(jb,kgrd)
970            zh  = SUM(pdta_read_dz(jb,1,:) )
971            ztrans = 0._wp
972            ztrans_new = 0._wp
973            DO jkb = 1, ipkb                              ! calculate transport on input grid
974               ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)
975            ENDDO
976            DO jk = 1, jpk                                ! calculate transport on model grid
977               ztrans_new = ztrans_new +      pdta(jb,1,jk ) *        e3v_n(ji,jj,jk ) * vmask(ji,jj,jk)
978            ENDDO
979            DO jk = 1, jpk                                ! make transport correction
980               IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data
981                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv_n(ji,jj) ) * vmask(ji,jj,jk)
982               ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
983                  pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hv_n(ji,jj)   * vmask(ji,jj,jk)
984               ENDIF
985            ENDDO
986         ENDDO
987      ENDIF
988     
989   END SUBROUTINE fld_bdy_interp
990
991
992   SUBROUTINE fld_rot( kt, sd )
993      !!---------------------------------------------------------------------
994      !!                    ***  ROUTINE fld_rot  ***
995      !!
996      !! ** Purpose :   Vector fields may need to be rotated onto the local grid direction
997      !!----------------------------------------------------------------------
998      INTEGER                , INTENT(in   ) ::   kt   ! ocean time step
999      TYPE(FLD), DIMENSION(:), INTENT(inout) ::   sd   ! input field related variables
1000      !
1001      INTEGER ::   ju, jv, jk, jn  ! loop indices
1002      INTEGER ::   imf             ! size of the structure sd
1003      INTEGER ::   ill             ! character length
1004      INTEGER ::   iv              ! indice of V component
1005      CHARACTER (LEN=100)          ::   clcomp       ! dummy weight name
1006      REAL(wp), DIMENSION(jpi,jpj) ::   utmp, vtmp   ! temporary arrays for vector rotation
1007      !!---------------------------------------------------------------------
1008      !
1009      !! (sga: following code should be modified so that pairs arent searched for each time
1010      !
1011      imf = SIZE( sd )
1012      DO ju = 1, imf
1013         IF( TRIM(sd(ju)%clrootname) == 'NOT USED' )   CYCLE
1014         ill = LEN_TRIM( sd(ju)%vcomp )
1015         DO jn = 2-COUNT((/sd(ju)%ln_tint/)), 2
1016            IF( ill > 0 .AND. .NOT. sd(ju)%rotn(jn) ) THEN   ! find vector rotations required             
1017               IF( sd(ju)%vcomp(1:1) == 'U' ) THEN      ! east-west component has symbolic name starting with 'U'
1018                  ! look for the north-south component which has same symbolic name but with 'U' replaced with 'V'
1019                  clcomp = 'V' // sd(ju)%vcomp(2:ill)   ! works even if ill == 1
1020                  iv = -1
1021                  DO jv = 1, imf
1022                     IF( TRIM(sd(jv)%clrootname) == 'NOT USED' )   CYCLE
1023                     IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) )   iv = jv
1024                  END DO
1025                  IF( iv > 0 ) THEN   ! fields ju and iv are two components which need to be rotated together
1026                     DO jk = 1, SIZE( sd(ju)%fnow, 3 )
1027                        IF( sd(ju)%ln_tint )THEN
1028                           CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->i', utmp(:,:) )
1029                           CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->j', vtmp(:,:) )
1030                           sd(ju)%fdta(:,:,jk,jn) = utmp(:,:)   ;   sd(iv)%fdta(:,:,jk,jn) = vtmp(:,:)
1031                        ELSE
1032                           CALL rot_rep( sd(ju)%fnow(:,:,jk  ), sd(iv)%fnow(:,:,jk  ), 'T', 'en->i', utmp(:,:) )
1033                           CALL rot_rep( sd(ju)%fnow(:,:,jk  ), sd(iv)%fnow(:,:,jk  ), 'T', 'en->j', vtmp(:,:) )
1034                           sd(ju)%fnow(:,:,jk   ) = utmp(:,:)   ;   sd(iv)%fnow(:,:,jk   ) = vtmp(:,:)
1035                        ENDIF
1036                     END DO
1037                     sd(ju)%rotn(jn) = .TRUE.               ! vector was rotated
1038                     IF( lwp .AND. kt == nit000 )   WRITE(numout,*)   &
1039                        &   'fld_read: vector pair ('//TRIM(sd(ju)%clvar)//', '//TRIM(sd(iv)%clvar)//') rotated on to model grid'
1040                  ENDIF
1041               ENDIF
1042            ENDIF
1043         END DO
1044       END DO
1045      !
1046   END SUBROUTINE fld_rot
1047
1048
1049   SUBROUTINE fld_clopn( sdjf, kyear, kmonth, kday, ldstop )
1050      !!---------------------------------------------------------------------
1051      !!                    ***  ROUTINE fld_clopn  ***
1052      !!
1053      !! ** Purpose :   update the file name and close/open the files
1054      !!----------------------------------------------------------------------
1055      TYPE(FLD)        , INTENT(inout) ::   sdjf     ! input field related variables
1056      INTEGER, OPTIONAL, INTENT(in   ) ::   kyear    ! year value
1057      INTEGER, OPTIONAL, INTENT(in   ) ::   kmonth   ! month value
1058      INTEGER, OPTIONAL, INTENT(in   ) ::   kday     ! day value
1059      LOGICAL, OPTIONAL, INTENT(in   ) ::   ldstop   ! stop if open to read a non-existing file (default = .TRUE.)
1060      !
1061      LOGICAL  :: llprevyr              ! are we reading previous year  file?
1062      LOGICAL  :: llprevmth             ! are we reading previous month file?
1063      INTEGER  :: iyear, imonth, iday   ! first day of the current file in yyyy mm dd
1064      INTEGER  :: isec_week             ! number of seconds since start of the weekly file
1065      INTEGER  :: indexyr               ! year undex (O/1/2: previous/current/next)
1066      REAL(wp) :: zyear_len, zmonth_len ! length (days) of iyear and imonth             !
1067      CHARACTER(len = 256) ::   clname  ! temporary file name
1068      !!----------------------------------------------------------------------
1069      IF( PRESENT(kyear) ) THEN                             ! use given values
1070         iyear = kyear
1071         imonth = kmonth
1072         iday = kday
1073         IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week
1074            isec_week = ksec_week( sdjf%cltype(6:8) )- (86400 * 8 ) 
1075            llprevmth  = isec_week > nsec_month             ! longer time since beginning of the week than the month
1076            llprevyr   = llprevmth .AND. nmonth == 1
1077            iyear  = nyear  - COUNT((/llprevyr /))
1078            imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
1079            iday   = nday   + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
1080         ENDIF
1081      ELSE                                                  ! use current day values
1082         IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week
1083            isec_week  = ksec_week( sdjf%cltype(6:8) )      ! second since the beginning of the week
1084            llprevmth  = isec_week > nsec_month             ! longer time since beginning of the week than the month
1085            llprevyr   = llprevmth .AND. nmonth == 1
1086         ELSE
1087            isec_week  = 0
1088            llprevmth  = .FALSE.
1089            llprevyr   = .FALSE.
1090         ENDIF
1091         iyear  = nyear  - COUNT((/llprevyr /))
1092         imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
1093         iday   = nday   + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
1094      ENDIF
1095
1096      ! build the new filename if not climatological data
1097      clname=TRIM(sdjf%clrootname)
1098      !
1099      ! note that sdjf%ln_clim is is only acting on the presence of the year in the file name
1100      IF( .NOT. sdjf%ln_clim ) THEN   
1101                                         WRITE(clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), iyear    ! add year
1102         IF( sdjf%cltype /= 'yearly' )   WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname          ), imonth   ! add month
1103      ELSE
1104         ! build the new filename if climatological data
1105         IF( sdjf%cltype /= 'yearly' )   WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), imonth   ! add month
1106      ENDIF
1107      IF( sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) &
1108            &                            WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname          ), iday     ! add day
1109      !
1110      IF( TRIM(clname) /= TRIM(sdjf%clname) .OR. sdjf%num == 0 ) THEN   ! new file to be open
1111         !
1112         sdjf%clname = TRIM(clname)
1113         IF( sdjf%num /= 0 )   CALL iom_close( sdjf%num )   ! close file if already open
1114         CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof =  LEN(TRIM(sdjf%wgtname)) > 0 )
1115         !
1116         ! find the last record to be read -> update sdjf%nreclast
1117         indexyr = iyear - nyear + 1
1118         zyear_len = REAL(nyear_len( indexyr ), wp)
1119         SELECT CASE ( indexyr )
1120         CASE ( 0 )   ;   zmonth_len = 31.   ! previous year -> imonth = 12
1121         CASE ( 1 )   ;   zmonth_len = REAL(nmonth_len(imonth), wp)
1122         CASE ( 2 )   ;   zmonth_len = 31.   ! next     year -> imonth = 1
1123         END SELECT
1124         !
1125         ! last record to be read in the current file
1126         IF    ( sdjf%freqh == -12. ) THEN                 ;   sdjf%nreclast = 1    !  yearly mean
1127         ELSEIF( sdjf%freqh ==  -1. ) THEN                                          ! monthly mean
1128            IF(     sdjf%cltype      == 'monthly' ) THEN   ;   sdjf%nreclast = 1
1129            ELSE                                           ;   sdjf%nreclast = 12
1130            ENDIF
1131         ELSE                                                                       ! higher frequency mean (in hours)
1132            IF(     sdjf%cltype      == 'monthly' ) THEN   ;   sdjf%nreclast = NINT( 24. * zmonth_len / sdjf%freqh )
1133            ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   sdjf%nreclast = NINT( 24. * 7.         / sdjf%freqh )
1134            ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   sdjf%nreclast = NINT( 24.              / sdjf%freqh )
1135            ELSE                                           ;   sdjf%nreclast = NINT( 24. * zyear_len  / sdjf%freqh )
1136            ENDIF
1137         ENDIF
1138         !
1139      ENDIF
1140      !
1141   END SUBROUTINE fld_clopn
1142
1143
1144   SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam, knoprint )
1145      !!---------------------------------------------------------------------
1146      !!                    ***  ROUTINE fld_fill  ***
1147      !!
1148      !! ** Purpose :   fill the data structure (sdf) with the associated information
1149      !!              read in namelist (sdf_n) and control print
1150      !!----------------------------------------------------------------------
1151      TYPE(FLD)  , DIMENSION(:)          , INTENT(inout) ::   sdf        ! structure of input fields (file informations, fields read)
1152      TYPE(FLD_N), DIMENSION(:)          , INTENT(in   ) ::   sdf_n      ! array of namelist information structures
1153      CHARACTER(len=*)                   , INTENT(in   ) ::   cdir       ! Root directory for location of flx files
1154      CHARACTER(len=*)                   , INTENT(in   ) ::   cdcaller   ! name of the calling routine
1155      CHARACTER(len=*)                   , INTENT(in   ) ::   cdtitle    ! description of the calling routine
1156      CHARACTER(len=*)                   , INTENT(in   ) ::   cdnam      ! name of the namelist from which sdf_n comes
1157      INTEGER                  , OPTIONAL, INTENT(in   ) ::   knoprint   ! no calling routine information printed
1158      !
1159      INTEGER  ::   jf   ! dummy indices
1160      !!---------------------------------------------------------------------
1161      !
1162      DO jf = 1, SIZE(sdf)
1163         sdf(jf)%clrootname = sdf_n(jf)%clname
1164         IF( TRIM(sdf_n(jf)%clname) /= 'NOT USED' )   sdf(jf)%clrootname = TRIM( cdir )//sdf(jf)%clrootname
1165         sdf(jf)%clname     = "not yet defined"
1166         sdf(jf)%freqh      = sdf_n(jf)%freqh
1167         sdf(jf)%clvar      = sdf_n(jf)%clvar
1168         sdf(jf)%ln_tint    = sdf_n(jf)%ln_tint
1169         sdf(jf)%ln_clim    = sdf_n(jf)%ln_clim
1170         sdf(jf)%cltype     = sdf_n(jf)%cltype
1171         sdf(jf)%num        = -1
1172         sdf(jf)%wgtname    = " "
1173         IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//sdf_n(jf)%wname
1174         sdf(jf)%lsmname = " "
1175         IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 )   sdf(jf)%lsmname = TRIM( cdir )//sdf_n(jf)%lname
1176         sdf(jf)%vcomp      = sdf_n(jf)%vcomp
1177         sdf(jf)%rotn(:)    = .TRUE.   ! pretend to be rotated -> won't try to rotate data before the first call to fld_get
1178         IF( sdf(jf)%cltype(1:4) == 'week' .AND. nn_leapy == 0  )   &
1179            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs nn_leapy = 1')
1180         IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim )   &
1181            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.')
1182         sdf(jf)%nreclast   = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
1183         sdf(jf)%igrd       = 0
1184         sdf(jf)%ibdy       = 0
1185         sdf(jf)%imap       => NULL()
1186         sdf(jf)%ltotvel    = .FALSE.
1187         sdf(jf)%lzint      = .FALSE.
1188      END DO
1189      !
1190      IF(lwp) THEN      ! control print
1191         WRITE(numout,*)
1192         IF( .NOT.PRESENT( knoprint) ) THEN
1193            WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
1194            WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
1195         ENDIF
1196         WRITE(numout,*) '   fld_fill : fill data structure with information from namelist '//TRIM( cdnam )
1197         WRITE(numout,*) '   ~~~~~~~~'
1198         WRITE(numout,*) '      list of files and frequency (>0: in hours ; <0 in months)'
1199         DO jf = 1, SIZE(sdf)
1200            WRITE(numout,*) '      root filename: '  , TRIM( sdf(jf)%clrootname ), '   variable name: ', TRIM( sdf(jf)%clvar )
1201            WRITE(numout,*) '         frequency: '      ,       sdf(jf)%freqh       ,   &
1202               &                  '   time interp: '    ,       sdf(jf)%ln_tint     ,   &
1203               &                  '   climatology: '    ,       sdf(jf)%ln_clim
1204            WRITE(numout,*) '         weights: '        , TRIM( sdf(jf)%wgtname    ),   &
1205               &                  '   pairing: '        , TRIM( sdf(jf)%vcomp      ),   &
1206               &                  '   data type: '      ,       sdf(jf)%cltype      ,   &
1207               &                  '   land/sea mask:'   , TRIM( sdf(jf)%lsmname    )
1208            call flush(numout)
1209         END DO
1210      ENDIF
1211      !
1212   END SUBROUTINE fld_fill
1213
1214
1215   SUBROUTINE wgt_list( sd, kwgt )
1216      !!---------------------------------------------------------------------
1217      !!                    ***  ROUTINE wgt_list  ***
1218      !!
1219      !! ** Purpose :   search array of WGTs and find a weights file entry,
1220      !!                or return a new one adding it to the end if new entry.
1221      !!                the weights data is read in and restructured (fld_weight)
1222      !!----------------------------------------------------------------------
1223      TYPE( FLD ), INTENT(in   ) ::   sd        ! field with name of weights file
1224      INTEGER    , INTENT(inout) ::   kwgt      ! index of weights
1225      !
1226      INTEGER ::   kw, nestid   ! local integer
1227      LOGICAL ::   found        ! local logical
1228      !!----------------------------------------------------------------------
1229      !
1230      !! search down linked list
1231      !! weights filename is either present or we hit the end of the list
1232      found = .FALSE.
1233      !
1234      !! because agrif nest part of filenames are now added in iom_open
1235      !! to distinguish between weights files on the different grids, need to track
1236      !! nest number explicitly
1237      nestid = 0
1238#if defined key_agrif
1239      nestid = Agrif_Fixed()
1240#endif
1241      DO kw = 1, nxt_wgt-1
1242         IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
1243             ref_wgts(kw)%nestid == nestid) THEN
1244            kwgt = kw
1245            found = .TRUE.
1246            EXIT
1247         ENDIF
1248      END DO
1249      IF( .NOT.found ) THEN
1250         kwgt = nxt_wgt
1251         CALL fld_weight( sd )
1252      ENDIF
1253      !
1254   END SUBROUTINE wgt_list
1255
1256
1257   SUBROUTINE wgt_print( )
1258      !!---------------------------------------------------------------------
1259      !!                    ***  ROUTINE wgt_print  ***
1260      !!
1261      !! ** Purpose :   print the list of known weights
1262      !!----------------------------------------------------------------------
1263      INTEGER ::   kw   !
1264      !!----------------------------------------------------------------------
1265      !
1266      DO kw = 1, nxt_wgt-1
1267         WRITE(numout,*) 'weight file:  ',TRIM(ref_wgts(kw)%wgtname)
1268         WRITE(numout,*) '      ddims:  ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
1269         WRITE(numout,*) '     numwgt:  ',ref_wgts(kw)%numwgt
1270         WRITE(numout,*) '     jpiwgt:  ',ref_wgts(kw)%jpiwgt
1271         WRITE(numout,*) '     jpjwgt:  ',ref_wgts(kw)%jpjwgt
1272         WRITE(numout,*) '    botleft:  ',ref_wgts(kw)%botleft
1273         WRITE(numout,*) '   topright:  ',ref_wgts(kw)%topright
1274         IF( ref_wgts(kw)%cyclic ) THEN
1275            WRITE(numout,*) '       cyclical'
1276            IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) '              with overlap of ', ref_wgts(kw)%overlap
1277         ELSE
1278            WRITE(numout,*) '       not cyclical'
1279         ENDIF
1280         IF( ASSOCIATED(ref_wgts(kw)%data_wgt) )  WRITE(numout,*) '       allocated'
1281      END DO
1282      !
1283   END SUBROUTINE wgt_print
1284
1285
1286   SUBROUTINE fld_weight( sd )
1287      !!---------------------------------------------------------------------
1288      !!                    ***  ROUTINE fld_weight  ***
1289      !!
1290      !! ** Purpose :   create a new WGT structure and fill in data from file,
1291      !!              restructuring as required
1292      !!----------------------------------------------------------------------
1293      TYPE( FLD ), INTENT(in) ::   sd   ! field with name of weights file
1294      !!
1295      INTEGER ::   jn         ! dummy loop indices
1296      INTEGER ::   inum       ! local logical unit
1297      INTEGER ::   id         ! local variable id
1298      INTEGER ::   ipk        ! local vertical dimension
1299      INTEGER ::   zwrap      ! local integer
1300      LOGICAL ::   cyclical   !
1301      CHARACTER (len=5) ::   aname   !
1302      INTEGER , DIMENSION(:), ALLOCATABLE ::   ddims
1303      INTEGER,  DIMENSION(jpi,jpj) ::   data_src
1304      REAL(wp), DIMENSION(jpi,jpj) ::   data_tmp
1305      !!----------------------------------------------------------------------
1306      !
1307      IF( nxt_wgt > tot_wgts ) THEN
1308        CALL ctl_stop("fld_weight: weights array size exceeded, increase tot_wgts")
1309      ENDIF
1310      !
1311      !! new weights file entry, add in extra information
1312      !! a weights file represents a 2D grid of a certain shape, so we assume that the current
1313      !! input data file is representative of all other files to be opened and processed with the
1314      !! current weights file
1315
1316      !! open input data file (non-model grid)
1317      CALL iom_open( sd%clname, inum, ldiof =  LEN(TRIM(sd%wgtname)) > 0 )
1318
1319      !! get dimensions
1320      IF ( SIZE(sd%fnow, 3) > 1 ) THEN
1321         ALLOCATE( ddims(4) )
1322      ELSE
1323         ALLOCATE( ddims(3) )
1324      ENDIF
1325      id = iom_varid( inum, sd%clvar, ddims )
1326
1327      !! close it
1328      CALL iom_close( inum )
1329
1330      !! now open the weights file
1331
1332      CALL iom_open ( sd%wgtname, inum )   ! interpolation weights
1333      IF ( inum > 0 ) THEN
1334
1335         !! determine whether we have an east-west cyclic grid
1336         !! from global attribute called "ew_wrap" in the weights file
1337         !! note that if not found, iom_getatt returns -999 and cyclic with no overlap is assumed
1338         !! since this is the most common forcing configuration
1339
1340         CALL iom_getatt(inum, 'ew_wrap', zwrap)
1341         IF( zwrap >= 0 ) THEN
1342            cyclical = .TRUE.
1343         ELSE IF( zwrap == -999 ) THEN
1344            cyclical = .TRUE.
1345            zwrap = 0
1346         ELSE
1347            cyclical = .FALSE.
1348         ENDIF
1349
1350         ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
1351         ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
1352         ref_wgts(nxt_wgt)%wgtname = sd%wgtname
1353         ref_wgts(nxt_wgt)%overlap = zwrap
1354         ref_wgts(nxt_wgt)%cyclic = cyclical
1355         ref_wgts(nxt_wgt)%nestid = 0
1356#if defined key_agrif
1357         ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
1358#endif
1359         !! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
1360         !! for each weight wgtNN there is an integer array srcNN which gives the point in
1361         !! the input data grid which is to be multiplied by the weight
1362         !! they are both arrays on the model grid so the result of the multiplication is
1363         !! added into an output array on the model grid as a running sum
1364
1365         !! two possible cases: bilinear (4 weights) or bicubic (16 weights)
1366         id = iom_varid(inum, 'src05', ldstop=.FALSE.)
1367         IF( id <= 0) THEN
1368            ref_wgts(nxt_wgt)%numwgt = 4
1369         ELSE
1370            ref_wgts(nxt_wgt)%numwgt = 16
1371         ENDIF
1372
1373         ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
1374         ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
1375         ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
1376
1377         DO jn = 1,4
1378            aname = ' '
1379            WRITE(aname,'(a3,i2.2)') 'src',jn
1380            data_tmp(:,:) = 0
1381            CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
1382            data_src(:,:) = INT(data_tmp(:,:))
1383            ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
1384            ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
1385         END DO
1386
1387         DO jn = 1, ref_wgts(nxt_wgt)%numwgt
1388            aname = ' '
1389            WRITE(aname,'(a3,i2.2)') 'wgt',jn
1390            ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
1391            CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
1392         END DO
1393         CALL iom_close (inum)
1394 
1395         ! find min and max indices in grid
1396         ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1397         ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1398         ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1399         ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1400
1401         ! and therefore dimensions of the input box
1402         ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
1403         ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
1404
1405         ! shift indexing of source grid
1406         ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
1407         ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
1408
1409         ! create input grid, give it a halo to allow gradient calculations
1410         ! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
1411         ! a more robust solution will be given in next release
1412         ipk =  SIZE(sd%fnow, 3)
1413         ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
1414         IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
1415         !
1416         nxt_wgt = nxt_wgt + 1
1417         !
1418      ELSE
1419         CALL ctl_stop( '    fld_weight : unable to read the file ' )
1420      ENDIF
1421
1422      DEALLOCATE (ddims )
1423      !
1424   END SUBROUTINE fld_weight
1425
1426
1427   SUBROUTINE apply_seaoverland( clmaskfile, zfieldo, jpi1_lsm, jpi2_lsm, jpj1_lsm,   &
1428      &                          jpj2_lsm, itmpi, itmpj, itmpz, rec1_lsm, recn_lsm )
1429      !!---------------------------------------------------------------------
1430      !!                    ***  ROUTINE apply_seaoverland  ***
1431      !!
1432      !! ** Purpose :   avoid spurious fluxes in coastal or near-coastal areas
1433      !!                due to the wrong usage of "land" values from the coarse
1434      !!                atmospheric model when spatial interpolation is required
1435      !!      D. Delrosso INGV         
1436      !!----------------------------------------------------------------------
1437      INTEGER,                   INTENT(in   ) :: itmpi,itmpj,itmpz                    ! lengths
1438      INTEGER,                   INTENT(in   ) :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm  ! temporary indices
1439      INTEGER, DIMENSION(3),     INTENT(in   ) :: rec1_lsm,recn_lsm                    ! temporary arrays for start and length
1440      REAL(wp),DIMENSION (:,:,:),INTENT(inout) :: zfieldo                              ! input/output array for seaoverland application
1441      CHARACTER (len=100),       INTENT(in   ) :: clmaskfile                           ! land/sea mask file name
1442      !
1443      INTEGER :: inum,jni,jnj,jnz,jc   ! local indices
1444      REAL(wp),DIMENSION (:,:,:),ALLOCATABLE :: zslmec1             ! local array for land point detection
1445      REAL(wp),DIMENSION (:,:),  ALLOCATABLE :: zfieldn   ! array of forcing field with undeff for land points
1446      REAL(wp),DIMENSION (:,:),  ALLOCATABLE :: zfield    ! array of forcing field
1447      !!---------------------------------------------------------------------
1448      !
1449      ALLOCATE ( zslmec1(itmpi,itmpj,itmpz), zfieldn(itmpi,itmpj), zfield(itmpi,itmpj) )
1450      !
1451      ! Retrieve the land sea mask data
1452      CALL iom_open( clmaskfile, inum )
1453      SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1454      CASE(1)
1455         CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm)
1456      CASE DEFAULT
1457         CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm)
1458      END SELECT
1459      CALL iom_close( inum )
1460      !
1461      DO jnz=1,rec1_lsm(3)             !! Loop over k dimension
1462         !
1463         DO jni = 1, itmpi                               !! copy the original field into a tmp array
1464            DO jnj = 1, itmpj                            !! substituting undeff over land points
1465               zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz)
1466               IF( zslmec1(jni,jnj,jnz) == 1. )   zfieldn(jni,jnj) = undeff_lsm
1467            END DO
1468         END DO
1469         !
1470         CALL seaoverland( zfieldn, itmpi, itmpj, zfield )
1471         DO jc = 1, nn_lsm
1472            CALL seaoverland( zfield, itmpi, itmpj, zfield )
1473         END DO
1474         !
1475         !   Check for Undeff and substitute original values
1476         IF( ANY(zfield==undeff_lsm) ) THEN
1477            DO jni = 1, itmpi
1478               DO jnj = 1, itmpj
1479                  IF( zfield(jni,jnj)==undeff_lsm )   zfield(jni,jnj) = zfieldo(jni,jnj,jnz)
1480               END DO
1481            END DO
1482         ENDIF
1483         !
1484         zfieldo(:,:,jnz) = zfield(:,:)
1485         !
1486      END DO                           !! End Loop over k dimension
1487      !
1488      DEALLOCATE ( zslmec1, zfieldn, zfield )
1489      !
1490   END SUBROUTINE apply_seaoverland 
1491
1492
1493   SUBROUTINE seaoverland( zfieldn, ileni, ilenj, zfield )
1494      !!---------------------------------------------------------------------
1495      !!                    ***  ROUTINE seaoverland  ***
1496      !!
1497      !! ** Purpose :   create shifted matrices for seaoverland application 
1498      !!      D. Delrosso INGV
1499      !!----------------------------------------------------------------------
1500      INTEGER                      , INTENT(in   ) :: ileni,ilenj   ! lengths
1501      REAL, DIMENSION (ileni,ilenj), INTENT(in   ) :: zfieldn       ! array of forcing field with undeff for land points
1502      REAL, DIMENSION (ileni,ilenj), INTENT(  out) :: zfield        ! array of forcing field
1503      !
1504      REAL   , DIMENSION (ileni,ilenj)   :: zmat1, zmat2, zmat3, zmat4  ! local arrays
1505      REAL   , DIMENSION (ileni,ilenj)   :: zmat5, zmat6, zmat7, zmat8  !   -     -
1506      REAL   , DIMENSION (ileni,ilenj)   :: zlsm2d                      !   -     -
1507      REAL   , DIMENSION (ileni,ilenj,8) :: zlsm3d                      !   -     -
1508      LOGICAL, DIMENSION (ileni,ilenj,8) :: ll_msknan3d                 ! logical mask for undeff detection
1509      LOGICAL, DIMENSION (ileni,ilenj)   :: ll_msknan2d                 ! logical mask for undeff detection
1510      !!----------------------------------------------------------------------
1511      zmat8 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(:,1)/)     , DIM=2 )
1512      zmat1 = eoshift( zmat8   , SHIFT=-1 , BOUNDARY = (/zmat8(1,:)/)       , DIM=1 )
1513      zmat2 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(1,:)/)     , DIM=1 )
1514      zmat4 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(:,ilenj)/) , DIM=2 )
1515      zmat3 = eoshift( zmat4   , SHIFT=-1 , BOUNDARY = (/zmat4(1,:)/)       , DIM=1 )
1516      zmat5 = eoshift( zmat4   , SHIFT= 1 , BOUNDARY = (/zmat4(ileni,:)/)   , DIM=1 )
1517      zmat6 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(ileni,:)/) , DIM=1 )
1518      zmat7 = eoshift( zmat8   , SHIFT= 1 , BOUNDARY = (/zmat8(ileni,:)/)   , DIM=1 )
1519      !
1520      zlsm3d  = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /))
1521      ll_msknan3d = .NOT.( zlsm3d  == undeff_lsm )
1522      ll_msknan2d = .NOT.( zfieldn == undeff_lsm )  ! FALSE where is Undeff (land)
1523      zlsm2d = SUM( zlsm3d, 3 , ll_msknan3d ) / MAX( 1 , COUNT( ll_msknan3d , 3 ) )
1524      WHERE( COUNT( ll_msknan3d , 3 ) == 0._wp )   zlsm2d = undeff_lsm
1525      zfield = MERGE( zfieldn, zlsm2d, ll_msknan2d )
1526      !
1527   END SUBROUTINE seaoverland
1528
1529
1530   SUBROUTINE fld_interp( num, clvar, kw, kk, dta,  &
1531                          &         nrec, lsmfile)     
1532      !!---------------------------------------------------------------------
1533      !!                    ***  ROUTINE fld_interp  ***
1534      !!
1535      !! ** Purpose :   apply weights to input gridded data to create data
1536      !!                on model grid
1537      !!----------------------------------------------------------------------
1538      INTEGER                   , INTENT(in   ) ::   num     ! stream number
1539      CHARACTER(LEN=*)          , INTENT(in   ) ::   clvar   ! variable name
1540      INTEGER                   , INTENT(in   ) ::   kw      ! weights number
1541      INTEGER                   , INTENT(in   ) ::   kk      ! vertical dimension of kk
1542      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   dta     ! output field on model grid
1543      INTEGER                   , INTENT(in   ) ::   nrec    ! record number to read (ie time slice)
1544      CHARACTER(LEN=*)          , INTENT(in   ) ::   lsmfile ! land sea mask file name
1545      !
1546      INTEGER, DIMENSION(3) ::   rec1, recn           ! temporary arrays for start and length
1547      INTEGER, DIMENSION(3) ::   rec1_lsm, recn_lsm   ! temporary arrays for start and length in case of seaoverland
1548      INTEGER ::   ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2    ! temporary indices
1549      INTEGER ::   jk, jn, jm, jir, jjr               ! loop counters
1550      INTEGER ::   ni, nj                             ! lengths
1551      INTEGER ::   jpimin,jpiwid                      ! temporary indices
1552      INTEGER ::   jpimin_lsm,jpiwid_lsm              ! temporary indices
1553      INTEGER ::   jpjmin,jpjwid                      ! temporary indices
1554      INTEGER ::   jpjmin_lsm,jpjwid_lsm              ! temporary indices
1555      INTEGER ::   jpi1,jpi2,jpj1,jpj2                ! temporary indices
1556      INTEGER ::   jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm   ! temporary indices
1557      INTEGER ::   itmpi,itmpj,itmpz                     ! lengths
1558      REAL(wp),DIMENSION(:,:,:), ALLOCATABLE ::   ztmp_fly_dta                 ! local array of values on input grid     
1559      !!----------------------------------------------------------------------
1560      !
1561      !! for weighted interpolation we have weights at four corners of a box surrounding
1562      !! a model grid point, each weight is multiplied by a grid value (bilinear case)
1563      !! or by a grid value and gradients at the corner point (bicubic case)
1564      !! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
1565
1566      !! sub grid from non-model input grid which encloses all grid points in this nemo process
1567      jpimin = ref_wgts(kw)%botleft(1)
1568      jpjmin = ref_wgts(kw)%botleft(2)
1569      jpiwid = ref_wgts(kw)%jpiwgt
1570      jpjwid = ref_wgts(kw)%jpjwgt
1571
1572      !! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients
1573      rec1(1) = MAX( jpimin-1, 1 )
1574      rec1(2) = MAX( jpjmin-1, 1 )
1575      rec1(3) = 1
1576      recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
1577      recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1578      recn(3) = kk
1579
1580      !! where we need to put it in the non-nemo grid fly_dta
1581      !! note that jpi1 and jpj1 only differ from 1 when jpimin and jpjmin are 1
1582      !! (ie at the extreme west or south of the whole input grid) and similarly for jpi2 and jpj2
1583      jpi1 = 2 + rec1(1) - jpimin
1584      jpj1 = 2 + rec1(2) - jpjmin
1585      jpi2 = jpi1 + recn(1) - 1
1586      jpj2 = jpj1 + recn(2) - 1
1587
1588
1589      IF( LEN( TRIM(lsmfile) ) > 0 ) THEN
1590      !! indeces for ztmp_fly_dta
1591      ! --------------------------
1592         rec1_lsm(1)=MAX(rec1(1)-nn_lsm,1)  ! starting index for enlarged external data, x direction
1593         rec1_lsm(2)=MAX(rec1(2)-nn_lsm,1)  ! starting index for enlarged external data, y direction
1594         rec1_lsm(3) = 1                    ! vertical dimension
1595         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
1596         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
1597         recn_lsm(3) = kk                   ! number of vertical levels in the input file
1598
1599      !  Avoid out of bound
1600         jpimin_lsm = MAX( rec1_lsm(1)+1, 1 )
1601         jpjmin_lsm = MAX( rec1_lsm(2)+1, 1 )
1602         jpiwid_lsm = MIN( recn_lsm(1)-2,ref_wgts(kw)%ddims(1)-rec1(1)+1)
1603         jpjwid_lsm = MIN( recn_lsm(2)-2,ref_wgts(kw)%ddims(2)-rec1(2)+1)
1604
1605         jpi1_lsm = 2+rec1_lsm(1)-jpimin_lsm
1606         jpj1_lsm = 2+rec1_lsm(2)-jpjmin_lsm
1607         jpi2_lsm = jpi1_lsm + recn_lsm(1) - 1
1608         jpj2_lsm = jpj1_lsm + recn_lsm(2) - 1
1609
1610
1611         itmpi=jpi2_lsm-jpi1_lsm+1
1612         itmpj=jpj2_lsm-jpj1_lsm+1
1613         itmpz=kk
1614         ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz))
1615         ztmp_fly_dta(:,:,:) = 0.0
1616         SELECT CASE( SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1617         CASE(1)
1618              CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1),   &
1619                 &                                                                nrec, rec1_lsm, recn_lsm)
1620         CASE DEFAULT
1621              CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),   &
1622                 &                                                                nrec, rec1_lsm, recn_lsm)
1623         END SELECT
1624         CALL apply_seaoverland(lsmfile,ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),                  &
1625                 &                                      jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm,                  &
1626                 &                                      itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
1627
1628
1629         ! Relative indeces for remapping
1630         ii_lsm1 = (rec1(1)-rec1_lsm(1))+1
1631         ii_lsm2 = (ii_lsm1+recn(1))-1
1632         ij_lsm1 = (rec1(2)-rec1_lsm(2))+1
1633         ij_lsm2 = (ij_lsm1+recn(2))-1
1634
1635         ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1636         ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:) = ztmp_fly_dta(ii_lsm1:ii_lsm2,ij_lsm1:ij_lsm2,:)
1637         DEALLOCATE(ztmp_fly_dta)
1638
1639      ELSE
1640         
1641         ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1642         SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) )
1643         CASE(1)
1644              CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn)
1645         CASE DEFAULT
1646              CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn)
1647         END SELECT
1648      ENDIF
1649     
1650
1651      !! first four weights common to both bilinear and bicubic
1652      !! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft
1653      !! note that we have to offset by 1 into fly_dta array because of halo
1654      dta(:,:,:) = 0.0
1655      DO jk = 1,4
1656        DO jn = 1, jpj
1657          DO jm = 1,jpi
1658            ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1659            nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1660            dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,:)
1661          END DO
1662        END DO
1663      END DO
1664
1665      IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
1666
1667        !! fix up halo points that we couldnt read from file
1668        IF( jpi1 == 2 ) THEN
1669           ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:)
1670        ENDIF
1671        IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1672           ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:)
1673        ENDIF
1674        IF( jpj1 == 2 ) THEN
1675           ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:)
1676        ENDIF
1677        IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
1678           ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:)
1679        ENDIF
1680
1681        !! if data grid is cyclic we can do better on east-west edges
1682        !! but have to allow for whether first and last columns are coincident
1683        IF( ref_wgts(kw)%cyclic ) THEN
1684           rec1(2) = MAX( jpjmin-1, 1 )
1685           recn(1) = 1
1686           recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1687           jpj1 = 2 + rec1(2) - jpjmin
1688           jpj2 = jpj1 + recn(2) - 1
1689           IF( jpi1 == 2 ) THEN
1690              rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap
1691              SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1692              CASE(1)
1693                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1694              CASE DEFAULT
1695                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1696              END SELECT     
1697              ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1698           ENDIF
1699           IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1700              rec1(1) = 1 + ref_wgts(kw)%overlap
1701              SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1702              CASE(1)
1703                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1704              CASE DEFAULT
1705                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1706              END SELECT
1707              ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1708           ENDIF
1709        ENDIF
1710
1711        ! gradient in the i direction
1712        DO jk = 1,4
1713          DO jn = 1, jpj
1714            DO jm = 1,jpi
1715              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1716              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1717              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 *         &
1718                               (ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:))
1719            END DO
1720          END DO
1721        END DO
1722
1723        ! gradient in the j direction
1724        DO jk = 1,4
1725          DO jn = 1, jpj
1726            DO jm = 1,jpi
1727              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1728              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1729              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 *         &
1730                               (ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:))
1731            END DO
1732          END DO
1733        END DO
1734
1735         ! gradient in the ij direction
1736         DO jk = 1,4
1737            DO jn = 1, jpj
1738               DO jm = 1,jpi
1739                  ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1740                  nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1741                  dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
1742                               (ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni  ,nj+2,:)) -   &
1743                               (ref_wgts(kw)%fly_dta(ni+2,nj  ,:) - ref_wgts(kw)%fly_dta(ni  ,nj  ,:)))
1744               END DO
1745            END DO
1746         END DO
1747         !
1748      END IF
1749      !
1750   END SUBROUTINE fld_interp
1751
1752
1753   FUNCTION ksec_week( cdday )
1754      !!---------------------------------------------------------------------
1755      !!                    ***  FUNCTION kshift_week ***
1756      !!
1757      !! ** Purpose :   return the first 3 letters of the first day of the weekly file
1758      !!---------------------------------------------------------------------
1759      CHARACTER(len=*), INTENT(in)   ::   cdday   ! first 3 letters of the first day of the weekly file
1760      !!
1761      INTEGER                        ::   ksec_week      ! output variable
1762      INTEGER                        ::   ijul, ishift   ! local integer
1763      CHARACTER(len=3),DIMENSION(7)  ::   cl_week 
1764      !!----------------------------------------------------------------------
1765      cl_week = (/"sun","sat","fri","thu","wed","tue","mon"/)
1766      DO ijul = 1, 7
1767         IF( cl_week(ijul) == TRIM(cdday) ) EXIT
1768      END DO
1769      IF( ijul .GT. 7 )   CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): '//TRIM(cdday) )
1770      !
1771      ishift = ijul * NINT(rday)
1772      !
1773      ksec_week = nsec_week + ishift
1774      ksec_week = MOD( ksec_week, 7*NINT(rday) )
1775      !
1776   END FUNCTION ksec_week
1777
1778   !!======================================================================
1779END MODULE fldread
Note: See TracBrowser for help on using the repository browser.