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

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

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 5883

Last change on this file since 5883 was 5883, checked in by gm, 8 years ago

#1613: vvl by default: TRA/TRC remove optimization associated with linear free surface

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