source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 10770

Last change on this file since 10770 was 10770, checked in by andmirek, 20 months ago

GMED 450 changes for writing to numout

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