source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/fldread.F90 @ 12182

Last change on this file since 12182 was 12182, checked in by davestorkey, 22 months ago

2019/dev_r11943_MERGE_2019: Merge in dev_ASINTER-01-05_merge.

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