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

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

source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 5132

Last change on this file since 5132 was 5132, checked in by jchanut, 8 years ago

Fixes obc/bdy data files compatibility. Ticket #1377

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