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

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

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

Last change on this file since 12367 was 12367, checked in by clem, 2 years ago

solve ticket #2380

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