source: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/fldread.F90 @ 11857

Last change on this file since 11857 was 11857, checked in by gsamson, 11 months ago

dev_r11265_ABL: bulks compatibility with 3D atmospheric forcing (by reading only the first level) (#2131)

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