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/UKMO/dev_r5518_obs_oper_update_ersem/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r5518_obs_oper_update_ersem/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 10260

Last change on this file since 10260 was 10260, checked in by dford, 5 years ago

Fixes for running with ERSEM (diaobs.F90) and merging with AMM15_v3_6_STABLE_package (fldread.F90 and trcnxt.F90).

File size: 84.8 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_clopn
33   PUBLIC   fld_map    ! routine called by tides_init
34   PUBLIC   fld_read, fld_fill   ! called by sbc... modules
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. kt - nit000 <= 100 ) 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. kt - nit000 <= 100 ) 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) 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      ELSE                                                  ! use current day values
819         IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week
820            isec_week  = ksec_week( sdjf%cltype(6:8) )      ! second since the beginning of the week
821            llprevmth  = isec_week > nsec_month             ! longer time since beginning of the week than the month
822            llprevyr   = llprevmth .AND. nmonth == 1
823         ELSE
824            isec_week  = 0
825            llprevmth  = .FALSE.
826            llprevyr   = .FALSE.
827         ENDIF
828         iyear  = nyear  - COUNT((/llprevyr /))
829         imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
830         iday   = nday   + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
831      ENDIF
832
833      ! build the new filename if not climatological data
834      clname=TRIM(sdjf%clrootname)
835      !
836      ! note that sdjf%ln_clim is is only acting on the presence of the year in the file name
837      IF( .NOT. sdjf%ln_clim ) THEN   
838                                         WRITE(clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), iyear    ! add year
839         IF( sdjf%cltype /= 'yearly' )   WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname          ), imonth   ! add month
840      ELSE
841         ! build the new filename if climatological data
842         IF( sdjf%cltype /= 'yearly' )   WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), imonth   ! add month
843      ENDIF
844      IF( sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) &
845            &                            WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname          ), iday     ! add day
846      !
847      IF( TRIM(clname) /= TRIM(sdjf%clname) .OR. sdjf%num == 0 ) THEN   ! new file to be open
848
849         sdjf%clname = TRIM(clname)
850         IF( sdjf%num /= 0 )   CALL iom_close( sdjf%num )   ! close file if already open
851         CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof =  LEN(TRIM(sdjf%wgtname)) > 0 )
852
853         ! find the last record to be read -> update sdjf%nreclast
854         indexyr = iyear - nyear + 1
855         iyear_len = nyear_len( indexyr )
856         SELECT CASE ( indexyr )
857         CASE ( 0 )   ;   imonth_len = 31   ! previous year -> imonth = 12
858         CASE ( 1 )   ;   imonth_len = nmonth_len(imonth) 
859         CASE ( 2 )   ;   imonth_len = 31   ! next     year -> imonth = 1
860         END SELECT
861         
862         ! last record to be read in the current file
863         IF    ( sdjf%nfreqh == -12 ) THEN                 ;   sdjf%nreclast = 1    !  yearly mean
864         ELSEIF( sdjf%nfreqh ==  -1 ) THEN                                          ! monthly mean
865            IF(     sdjf%cltype      == 'monthly' ) THEN   ;   sdjf%nreclast = 1
866            ELSE                                           ;   sdjf%nreclast = 12
867            ENDIF
868         ELSE                                                                       ! higher frequency mean (in hours)
869            IF(     sdjf%cltype      == 'monthly' ) THEN   ;   sdjf%nreclast = NINT( 24 * imonth_len / sdjf%nfreqh )
870            ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   sdjf%nreclast = NINT( 24 * 7          / sdjf%nfreqh )
871            ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   sdjf%nreclast = NINT( 24              / sdjf%nfreqh )
872            ELSE                                           ;   sdjf%nreclast = NINT( 24 * iyear_len  / sdjf%nfreqh )
873            ENDIF
874         ENDIF
875         
876      ENDIF
877      !
878   END SUBROUTINE fld_clopn
879
880
881   SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )
882      !!---------------------------------------------------------------------
883      !!                    ***  ROUTINE fld_fill  ***
884      !!
885      !! ** Purpose :   fill sdf with sdf_n and control print
886      !!----------------------------------------------------------------------
887      TYPE(FLD)  , DIMENSION(:), INTENT(inout) ::   sdf        ! structure of input fields (file informations, fields read)
888      TYPE(FLD_N), DIMENSION(:), INTENT(in   ) ::   sdf_n      ! array of namelist information structures
889      CHARACTER(len=*)         , INTENT(in   ) ::   cdir       ! Root directory for location of flx files
890      CHARACTER(len=*)         , INTENT(in   ) ::   cdcaller   !
891      CHARACTER(len=*)         , INTENT(in   ) ::   cdtitle    !
892      CHARACTER(len=*)         , INTENT(in   ) ::   cdnam      !
893      !
894      INTEGER  ::   jf       ! dummy indices
895      !!---------------------------------------------------------------------
896
897      DO jf = 1, SIZE(sdf)
898         sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
899         sdf(jf)%clname     = "not yet defined"
900         sdf(jf)%nfreqh     = sdf_n(jf)%nfreqh
901         sdf(jf)%clvar      = sdf_n(jf)%clvar
902         sdf(jf)%ln_tint    = sdf_n(jf)%ln_tint
903         sdf(jf)%ln_clim    = sdf_n(jf)%ln_clim
904         sdf(jf)%cltype     = sdf_n(jf)%cltype
905         sdf(jf)%num        = -1
906         sdf(jf)%wgtname    = " "
907         IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
908         sdf(jf)%lsmname = " "
909         IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 )   sdf(jf)%lsmname = TRIM( cdir )//TRIM( sdf_n(jf)%lname )
910         sdf(jf)%vcomp      = sdf_n(jf)%vcomp
911         sdf(jf)%rotn(:)    = .TRUE.   ! pretend to be rotated -> won't try to rotate data before the first call to fld_get
912         IF( sdf(jf)%cltype(1:4) == 'week' .AND. nn_leapy == 0  )   &
913            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs nn_leapy = 1')
914         IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim )   &
915            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.')
916         sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
917      END DO
918
919      IF(lwp) THEN      ! control print
920         WRITE(numout,*)
921         WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
922         WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
923         WRITE(numout,*) '          '//TRIM( cdnam )//' Namelist'
924         WRITE(numout,*) '          list of files and frequency (>0: in hours ; <0 in months)'
925         DO jf = 1, SIZE(sdf)
926            WRITE(numout,*) '               root filename: '  , TRIM( sdf(jf)%clrootname ),   &
927               &                          ' variable name: '  , TRIM( sdf(jf)%clvar      )
928            WRITE(numout,*) '               frequency: '      ,       sdf(jf)%nfreqh      ,   &
929               &                          ' time interp: '    ,       sdf(jf)%ln_tint     ,   &
930               &                          ' climatology: '    ,       sdf(jf)%ln_clim     ,   &
931               &                          ' weights    : '    , TRIM( sdf(jf)%wgtname    ),   &
932               &                          ' pairing    : '    , TRIM( sdf(jf)%vcomp      ),   &
933               &                          ' data type: '      ,       sdf(jf)%cltype      ,   &
934               &                          ' land/sea mask:'   , TRIM( sdf(jf)%lsmname    )
935            call flush(numout)
936         END DO
937      ENDIF
938     
939   END SUBROUTINE fld_fill
940
941
942   SUBROUTINE wgt_list( sd, kwgt )
943      !!---------------------------------------------------------------------
944      !!                    ***  ROUTINE wgt_list  ***
945      !!
946      !! ** Purpose :   search array of WGTs and find a weights file
947      !!                entry, or return a new one adding it to the end
948      !!                if it is a new entry, the weights data is read in and
949      !!                restructured (fld_weight)
950      !!----------------------------------------------------------------------
951      TYPE( FLD ), INTENT(in   ) ::   sd        ! field with name of weights file
952      INTEGER    , INTENT(inout) ::   kwgt      ! index of weights
953      !!
954      INTEGER ::   kw, nestid   ! local integer
955      LOGICAL ::   found        ! local logical
956      !!----------------------------------------------------------------------
957      !
958      !! search down linked list
959      !! weights filename is either present or we hit the end of the list
960      found = .FALSE.
961
962      !! because agrif nest part of filenames are now added in iom_open
963      !! to distinguish between weights files on the different grids, need to track
964      !! nest number explicitly
965      nestid = 0
966#if defined key_agrif
967      nestid = Agrif_Fixed()
968#endif
969      DO kw = 1, nxt_wgt-1
970         IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
971             ref_wgts(kw)%nestid == nestid) THEN
972            kwgt = kw
973            found = .TRUE.
974            EXIT
975         ENDIF
976      END DO
977      IF( .NOT.found ) THEN
978         kwgt = nxt_wgt
979         CALL fld_weight( sd )
980      ENDIF
981      !
982   END SUBROUTINE wgt_list
983
984
985   SUBROUTINE wgt_print( )
986      !!---------------------------------------------------------------------
987      !!                    ***  ROUTINE wgt_print  ***
988      !!
989      !! ** Purpose :   print the list of known weights
990      !!----------------------------------------------------------------------
991      INTEGER ::   kw   !
992      !!----------------------------------------------------------------------
993      !
994      DO kw = 1, nxt_wgt-1
995         WRITE(numout,*) 'weight file:  ',TRIM(ref_wgts(kw)%wgtname)
996         WRITE(numout,*) '      ddims:  ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
997         WRITE(numout,*) '     numwgt:  ',ref_wgts(kw)%numwgt
998         WRITE(numout,*) '     jpiwgt:  ',ref_wgts(kw)%jpiwgt
999         WRITE(numout,*) '     jpjwgt:  ',ref_wgts(kw)%jpjwgt
1000         WRITE(numout,*) '    botleft:  ',ref_wgts(kw)%botleft
1001         WRITE(numout,*) '   topright:  ',ref_wgts(kw)%topright
1002         IF( ref_wgts(kw)%cyclic ) THEN
1003            WRITE(numout,*) '       cyclical'
1004            IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) '              with overlap of ', ref_wgts(kw)%overlap
1005         ELSE
1006            WRITE(numout,*) '       not cyclical'
1007         ENDIF
1008         IF( ASSOCIATED(ref_wgts(kw)%data_wgt) )  WRITE(numout,*) '       allocated'
1009      END DO
1010      !
1011   END SUBROUTINE wgt_print
1012
1013
1014   SUBROUTINE fld_weight( sd )
1015      !!---------------------------------------------------------------------
1016      !!                    ***  ROUTINE fld_weight  ***
1017      !!
1018      !! ** Purpose :   create a new WGT structure and fill in data from 
1019      !!                file, restructuring as required
1020      !!----------------------------------------------------------------------
1021      TYPE( FLD ), INTENT(in) ::   sd   ! field with name of weights file
1022      !!
1023      INTEGER                           ::   jn            ! dummy loop indices
1024      INTEGER                           ::   inum          ! temporary logical unit
1025      INTEGER                           ::   id            ! temporary variable id
1026      INTEGER                           ::   ipk           ! temporary vertical dimension
1027      CHARACTER (len=5)                 ::   aname
1028      INTEGER , DIMENSION(:), ALLOCATABLE ::   ddims
1029      INTEGER , POINTER, DIMENSION(:,:) ::   data_src
1030      REAL(wp), POINTER, DIMENSION(:,:) ::   data_tmp
1031      LOGICAL                           ::   cyclical
1032      INTEGER                           ::   zwrap      ! local integer
1033      !!----------------------------------------------------------------------
1034      !
1035      CALL wrk_alloc( jpi,jpj, data_src )   ! integer
1036      CALL wrk_alloc( jpi,jpj, data_tmp )
1037      !
1038      IF( nxt_wgt > tot_wgts ) THEN
1039        CALL ctl_stop("fld_weight: weights array size exceeded, increase tot_wgts")
1040      ENDIF
1041      !
1042      !! new weights file entry, add in extra information
1043      !! a weights file represents a 2D grid of a certain shape, so we assume that the current
1044      !! input data file is representative of all other files to be opened and processed with the
1045      !! current weights file
1046
1047      !! open input data file (non-model grid)
1048      CALL iom_open( sd%clname, inum, ldiof =  LEN(TRIM(sd%wgtname)) > 0 )
1049
1050      !! get dimensions
1051      IF ( SIZE(sd%fnow, 3) > 1 ) THEN
1052         ALLOCATE( ddims(4) )
1053      ELSE
1054         ALLOCATE( ddims(3) )
1055      ENDIF
1056      id = iom_varid( inum, sd%clvar, ddims )
1057
1058      !! close it
1059      CALL iom_close( inum )
1060
1061      !! now open the weights file
1062
1063      CALL iom_open ( sd%wgtname, inum )   ! interpolation weights
1064      IF ( inum > 0 ) THEN
1065
1066         !! determine whether we have an east-west cyclic grid
1067         !! from global attribute called "ew_wrap" in the weights file
1068         !! note that if not found, iom_getatt returns -999 and cyclic with no overlap is assumed
1069         !! since this is the most common forcing configuration
1070
1071         CALL iom_getatt(inum, 'ew_wrap', zwrap)
1072         IF( zwrap >= 0 ) THEN
1073            cyclical = .TRUE.
1074         ELSE IF( zwrap == -999 ) THEN
1075            cyclical = .TRUE.
1076            zwrap = 0
1077         ELSE
1078            cyclical = .FALSE.
1079         ENDIF
1080
1081         ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
1082         ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
1083         ref_wgts(nxt_wgt)%wgtname = sd%wgtname
1084         ref_wgts(nxt_wgt)%overlap = zwrap
1085         ref_wgts(nxt_wgt)%cyclic = cyclical
1086         ref_wgts(nxt_wgt)%nestid = 0
1087#if defined key_agrif
1088         ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
1089#endif
1090         !! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
1091         !! for each weight wgtNN there is an integer array srcNN which gives the point in
1092         !! the input data grid which is to be multiplied by the weight
1093         !! they are both arrays on the model grid so the result of the multiplication is
1094         !! added into an output array on the model grid as a running sum
1095
1096         !! two possible cases: bilinear (4 weights) or bicubic (16 weights)
1097         id = iom_varid(inum, 'src05', ldstop=.FALSE.)
1098         IF( id <= 0) THEN
1099            ref_wgts(nxt_wgt)%numwgt = 4
1100         ELSE
1101            ref_wgts(nxt_wgt)%numwgt = 16
1102         ENDIF
1103
1104         ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
1105         ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
1106         ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
1107
1108         DO jn = 1,4
1109            aname = ' '
1110            WRITE(aname,'(a3,i2.2)') 'src',jn
1111            data_tmp(:,:) = 0
1112            CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
1113            data_src(:,:) = INT(data_tmp(:,:))
1114            ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
1115            ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
1116         END DO
1117
1118         DO jn = 1, ref_wgts(nxt_wgt)%numwgt
1119            aname = ' '
1120            WRITE(aname,'(a3,i2.2)') 'wgt',jn
1121            ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
1122            CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
1123         END DO
1124         CALL iom_close (inum)
1125 
1126         ! find min and max indices in grid
1127         ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1128         ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1129         ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1130         ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1131
1132         ! and therefore dimensions of the input box
1133         ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
1134         ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
1135
1136         ! shift indexing of source grid
1137         ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
1138         ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
1139
1140         ! create input grid, give it a halo to allow gradient calculations
1141         ! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
1142         ! a more robust solution will be given in next release
1143         ipk =  SIZE(sd%fnow, 3)
1144         ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
1145         IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
1146
1147         nxt_wgt = nxt_wgt + 1
1148
1149      ELSE
1150         CALL ctl_stop( '    fld_weight : unable to read the file ' )
1151      ENDIF
1152
1153      DEALLOCATE (ddims )
1154
1155      CALL wrk_dealloc( jpi,jpj, data_src )   ! integer
1156      CALL wrk_dealloc( jpi,jpj, data_tmp )
1157      !
1158   END SUBROUTINE fld_weight
1159
1160
1161   SUBROUTINE apply_seaoverland(clmaskfile,zfieldo,jpi1_lsm,jpi2_lsm,jpj1_lsm, &
1162                          &      jpj2_lsm,itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
1163      !!---------------------------------------------------------------------
1164      !!                    ***  ROUTINE apply_seaoverland  ***
1165      !!
1166      !! ** Purpose :   avoid spurious fluxes in coastal or near-coastal areas
1167      !!                due to the wrong usage of "land" values from the coarse
1168      !!                atmospheric model when spatial interpolation is required
1169      !!      D. Delrosso INGV         
1170      !!----------------------------------------------------------------------
1171      INTEGER                                   :: inum,jni,jnj,jnz,jc                  ! temporary indices
1172      INTEGER,                   INTENT(in)     :: itmpi,itmpj,itmpz                    ! lengths
1173      INTEGER,                   INTENT(in)     :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm  ! temporary indices
1174      INTEGER, DIMENSION(3),     INTENT(in)     :: rec1_lsm,recn_lsm                    ! temporary arrays for start and length
1175      REAL(wp),DIMENSION (:,:,:),INTENT(inout)  :: zfieldo                              ! input/output array for seaoverland application
1176      REAL(wp),DIMENSION (:,:,:),ALLOCATABLE    :: zslmec1                              ! temporary array for land point detection
1177      REAL(wp),DIMENSION (:,:),  ALLOCATABLE    :: zfieldn                              ! array of forcing field with undeff for land points
1178      REAL(wp),DIMENSION (:,:),  ALLOCATABLE    :: zfield                               ! array of forcing field
1179      CHARACTER (len=100),       INTENT(in)     :: clmaskfile                           ! land/sea mask file name
1180      !!---------------------------------------------------------------------
1181      ALLOCATE ( zslmec1(itmpi,itmpj,itmpz) )
1182      ALLOCATE ( zfieldn(itmpi,itmpj) )
1183      ALLOCATE ( zfield(itmpi,itmpj) )
1184
1185      ! Retrieve the land sea mask data
1186      CALL iom_open( clmaskfile, inum )
1187      SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1188      CASE(1)
1189           CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm)
1190      CASE DEFAULT
1191           CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm)
1192      END SELECT
1193      CALL iom_close( inum )
1194
1195      DO jnz=1,rec1_lsm(3)                            !! Loop over k dimension
1196
1197         DO jni=1,itmpi                               !! copy the original field into a tmp array
1198            DO jnj=1,itmpj                            !! substituting undeff over land points
1199            zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz)
1200               IF ( zslmec1(jni,jnj,jnz) == 1. ) THEN
1201                  zfieldn(jni,jnj) = undeff_lsm
1202               ENDIF
1203            END DO
1204         END DO
1205 
1206      CALL seaoverland(zfieldn,itmpi,itmpj,zfield)
1207      DO jc=1,nn_lsm
1208         CALL seaoverland(zfield,itmpi,itmpj,zfield)
1209      END DO
1210
1211      !   Check for Undeff and substitute original values
1212      IF(ANY(zfield==undeff_lsm)) THEN
1213         DO jni=1,itmpi
1214            DO jnj=1,itmpj
1215               IF (zfield(jni,jnj)==undeff_lsm) THEN
1216                  zfield(jni,jnj) = zfieldo(jni,jnj,jnz)
1217               ENDIF
1218            ENDDO
1219         ENDDO
1220      ENDIF
1221
1222      zfieldo(:,:,jnz)=zfield(:,:)
1223
1224      END DO                          !! End Loop over k dimension
1225
1226      DEALLOCATE ( zslmec1 )
1227      DEALLOCATE ( zfieldn )
1228      DEALLOCATE ( zfield )
1229
1230   END SUBROUTINE apply_seaoverland 
1231
1232
1233   SUBROUTINE seaoverland(zfieldn,ileni,ilenj,zfield)
1234      !!---------------------------------------------------------------------
1235      !!                    ***  ROUTINE seaoverland  ***
1236      !!
1237      !! ** Purpose :   create shifted matrices for seaoverland application 
1238      !!      D. Delrosso INGV
1239      !!----------------------------------------------------------------------
1240      INTEGER,INTENT(in)                       :: ileni,ilenj              ! lengths
1241      REAL,DIMENSION (ileni,ilenj),INTENT(in)  :: zfieldn                  ! array of forcing field with undeff for land points
1242      REAL,DIMENSION (ileni,ilenj),INTENT(out) :: zfield                   ! array of forcing field
1243      REAL,DIMENSION (ileni,ilenj)             :: zmat1,zmat2,zmat3,zmat4  ! temporary arrays for seaoverland application
1244      REAL,DIMENSION (ileni,ilenj)             :: zmat5,zmat6,zmat7,zmat8  ! temporary arrays for seaoverland application
1245      REAL,DIMENSION (ileni,ilenj)             :: zlsm2d                   ! temporary arrays for seaoverland application
1246      REAL,DIMENSION (ileni,ilenj,8)           :: zlsm3d                   ! temporary arrays for seaoverland application
1247      LOGICAL,DIMENSION (ileni,ilenj,8)        :: ll_msknan3d              ! logical mask for undeff detection
1248      LOGICAL,DIMENSION (ileni,ilenj)          :: ll_msknan2d              ! logical mask for undeff detection
1249      !!----------------------------------------------------------------------
1250      zmat8 = eoshift(zfieldn   ,  SHIFT=-1, BOUNDARY = (/zfieldn(:,1)/)    ,DIM=2)
1251      zmat1 = eoshift(zmat8     ,  SHIFT=-1, BOUNDARY = (/zmat8(1,:)/)      ,DIM=1)
1252      zmat2 = eoshift(zfieldn   ,  SHIFT=-1, BOUNDARY = (/zfieldn(1,:)/)    ,DIM=1)
1253      zmat4 = eoshift(zfieldn   ,  SHIFT= 1, BOUNDARY = (/zfieldn(:,ilenj)/),DIM=2)
1254      zmat3 = eoshift(zmat4     ,  SHIFT=-1, BOUNDARY = (/zmat4(1,:)/)      ,DIM=1)
1255      zmat5 = eoshift(zmat4     ,  SHIFT= 1, BOUNDARY = (/zmat4(ileni,:)/)  ,DIM=1)
1256      zmat6 = eoshift(zfieldn   ,  SHIFT= 1, BOUNDARY = (/zfieldn(ileni,:)/),DIM=1)
1257      zmat7 = eoshift(zmat8     ,  SHIFT= 1, BOUNDARY = (/zmat8(ileni,:)/)  ,DIM=1)
1258
1259      zlsm3d  = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /))
1260      ll_msknan3d = .not.(zlsm3d==undeff_lsm)
1261      ll_msknan2d = .not.(zfieldn==undeff_lsm)  ! FALSE where is Undeff (land)
1262      zlsm2d = (SUM ( zlsm3d, 3 , ll_msknan3d ) )/(MAX(1,(COUNT( ll_msknan3d , 3 ))   ))
1263      WHERE ((COUNT( ll_msknan3d , 3 )) == 0.0_wp)  zlsm2d = undeff_lsm
1264      zfield = MERGE (zfieldn,zlsm2d,ll_msknan2d)
1265   END SUBROUTINE seaoverland
1266
1267
1268   SUBROUTINE fld_interp( num, clvar, kw, kk, dta,  &
1269                          &         nrec, lsmfile)     
1270      !!---------------------------------------------------------------------
1271      !!                    ***  ROUTINE fld_interp  ***
1272      !!
1273      !! ** Purpose :   apply weights to input gridded data to create data
1274      !!                on model grid
1275      !!----------------------------------------------------------------------
1276      INTEGER                   , INTENT(in   ) ::   num     ! stream number
1277      CHARACTER(LEN=*)          , INTENT(in   ) ::   clvar   ! variable name
1278      INTEGER                   , INTENT(in   ) ::   kw      ! weights number
1279      INTEGER                   , INTENT(in   ) ::   kk      ! vertical dimension of kk
1280      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   dta     ! output field on model grid
1281      INTEGER                   , INTENT(in   ) ::   nrec    ! record number to read (ie time slice)
1282      CHARACTER(LEN=*)          , INTENT(in   ) ::   lsmfile ! land sea mask file name
1283      !!
1284      REAL(wp),DIMENSION(:,:,:),ALLOCATABLE     ::   ztmp_fly_dta,zfieldo                  ! temporary array of values on input grid
1285      INTEGER, DIMENSION(3)                     ::   rec1,recn                             ! temporary arrays for start and length
1286      INTEGER, DIMENSION(3)                     ::   rec1_lsm,recn_lsm                     ! temporary arrays for start and length in case of seaoverland
1287      INTEGER                                   ::   ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2       ! temporary indices
1288      INTEGER                                   ::   jk, jn, jm, jir, jjr                  ! loop counters
1289      INTEGER                                   ::   ni, nj                                ! lengths
1290      INTEGER                                   ::   jpimin,jpiwid                         ! temporary indices
1291      INTEGER                                   ::   jpimin_lsm,jpiwid_lsm                 ! temporary indices
1292      INTEGER                                   ::   jpjmin,jpjwid                         ! temporary indices
1293      INTEGER                                   ::   jpjmin_lsm,jpjwid_lsm                 ! temporary indices
1294      INTEGER                                   ::   jpi1,jpi2,jpj1,jpj2                   ! temporary indices
1295      INTEGER                                   ::   jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm   ! temporary indices
1296      INTEGER                                   ::   itmpi,itmpj,itmpz                     ! lengths
1297     
1298      !!----------------------------------------------------------------------
1299      !
1300      !! for weighted interpolation we have weights at four corners of a box surrounding
1301      !! a model grid point, each weight is multiplied by a grid value (bilinear case)
1302      !! or by a grid value and gradients at the corner point (bicubic case)
1303      !! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
1304
1305      !! sub grid from non-model input grid which encloses all grid points in this nemo process
1306      jpimin = ref_wgts(kw)%botleft(1)
1307      jpjmin = ref_wgts(kw)%botleft(2)
1308      jpiwid = ref_wgts(kw)%jpiwgt
1309      jpjwid = ref_wgts(kw)%jpjwgt
1310
1311      !! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients
1312      rec1(1) = MAX( jpimin-1, 1 )
1313      rec1(2) = MAX( jpjmin-1, 1 )
1314      rec1(3) = 1
1315      recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
1316      recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1317      recn(3) = kk
1318
1319      !! where we need to put it in the non-nemo grid fly_dta
1320      !! note that jpi1 and jpj1 only differ from 1 when jpimin and jpjmin are 1
1321      !! (ie at the extreme west or south of the whole input grid) and similarly for jpi2 and jpj2
1322      jpi1 = 2 + rec1(1) - jpimin
1323      jpj1 = 2 + rec1(2) - jpjmin
1324      jpi2 = jpi1 + recn(1) - 1
1325      jpj2 = jpj1 + recn(2) - 1
1326
1327
1328      IF( LEN( TRIM(lsmfile) ) > 0 ) THEN
1329      !! indeces for ztmp_fly_dta
1330      ! --------------------------
1331         rec1_lsm(1)=MAX(rec1(1)-nn_lsm,1)  ! starting index for enlarged external data, x direction
1332         rec1_lsm(2)=MAX(rec1(2)-nn_lsm,1)  ! starting index for enlarged external data, y direction
1333         rec1_lsm(3) = 1                    ! vertical dimension
1334         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
1335         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
1336         recn_lsm(3) = kk                   ! number of vertical levels in the input file
1337
1338      !  Avoid out of bound
1339         jpimin_lsm = MAX( rec1_lsm(1)+1, 1 )
1340         jpjmin_lsm = MAX( rec1_lsm(2)+1, 1 )
1341         jpiwid_lsm = MIN( recn_lsm(1)-2,ref_wgts(kw)%ddims(1)-rec1(1)+1)
1342         jpjwid_lsm = MIN( recn_lsm(2)-2,ref_wgts(kw)%ddims(2)-rec1(2)+1)
1343
1344         jpi1_lsm = 2+rec1_lsm(1)-jpimin_lsm
1345         jpj1_lsm = 2+rec1_lsm(2)-jpjmin_lsm
1346         jpi2_lsm = jpi1_lsm + recn_lsm(1) - 1
1347         jpj2_lsm = jpj1_lsm + recn_lsm(2) - 1
1348
1349
1350         itmpi=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),1)
1351         itmpj=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),2)
1352         itmpz=kk
1353         ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz))
1354         ztmp_fly_dta(:,:,:) = 0.0
1355         SELECT CASE( SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1356         CASE(1)
1357              CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1),   &
1358                 &                                                                nrec, rec1_lsm, recn_lsm)
1359         CASE DEFAULT
1360              CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),   &
1361                 &                                                                nrec, rec1_lsm, recn_lsm)
1362         END SELECT
1363         CALL apply_seaoverland(lsmfile,ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),                  &
1364                 &                                      jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm,                  &
1365                 &                                      itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
1366
1367
1368         ! Relative indeces for remapping
1369         ii_lsm1 = (rec1(1)-rec1_lsm(1))+1
1370         ii_lsm2 = (ii_lsm1+recn(1))-1
1371         ij_lsm1 = (rec1(2)-rec1_lsm(2))+1
1372         ij_lsm2 = (ij_lsm1+recn(2))-1
1373
1374         ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1375         ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:) = ztmp_fly_dta(ii_lsm1:ii_lsm2,ij_lsm1:ij_lsm2,:)
1376         DEALLOCATE(ztmp_fly_dta)
1377
1378      ELSE
1379         
1380         ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1381         SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) )
1382         CASE(1)
1383              CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn)
1384         CASE DEFAULT
1385              CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn)
1386         END SELECT
1387      ENDIF
1388     
1389
1390      !! first four weights common to both bilinear and bicubic
1391      !! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft
1392      !! note that we have to offset by 1 into fly_dta array because of halo
1393      dta(:,:,:) = 0.0
1394      DO jk = 1,4
1395        DO jn = 1, jpj
1396          DO jm = 1,jpi
1397            ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1398            nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1399            dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,:)
1400          END DO
1401        END DO
1402      END DO
1403
1404      IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
1405
1406        !! fix up halo points that we couldnt read from file
1407        IF( jpi1 == 2 ) THEN
1408           ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:)
1409        ENDIF
1410        IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1411           ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:)
1412        ENDIF
1413        IF( jpj1 == 2 ) THEN
1414           ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:)
1415        ENDIF
1416        IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
1417           ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:)
1418        ENDIF
1419
1420        !! if data grid is cyclic we can do better on east-west edges
1421        !! but have to allow for whether first and last columns are coincident
1422        IF( ref_wgts(kw)%cyclic ) THEN
1423           rec1(2) = MAX( jpjmin-1, 1 )
1424           recn(1) = 1
1425           recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1426           jpj1 = 2 + rec1(2) - jpjmin
1427           jpj2 = jpj1 + recn(2) - 1
1428           IF( jpi1 == 2 ) THEN
1429              rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap
1430              SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1431              CASE(1)
1432                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1433              CASE DEFAULT
1434                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1435              END SELECT     
1436              ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1437           ENDIF
1438           IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1439              rec1(1) = 1 + ref_wgts(kw)%overlap
1440              SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1441              CASE(1)
1442                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1443              CASE DEFAULT
1444                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1445              END SELECT
1446              ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1447           ENDIF
1448        ENDIF
1449
1450        ! gradient in the i direction
1451        DO jk = 1,4
1452          DO jn = 1, jpj
1453            DO jm = 1,jpi
1454              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1455              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1456              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 *         &
1457                               (ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:))
1458            END DO
1459          END DO
1460        END DO
1461
1462        ! gradient in the j 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+8) * 0.5 *         &
1469                               (ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:))
1470            END DO
1471          END DO
1472        END DO
1473
1474         ! gradient in the ij 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+12) * 0.25 * ( &
1481                               (ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni  ,nj+2,:)) -   &
1482                               (ref_wgts(kw)%fly_dta(ni+2,nj  ,:) - ref_wgts(kw)%fly_dta(ni  ,nj  ,:)))
1483               END DO
1484            END DO
1485         END DO
1486         !
1487      END IF
1488      !
1489   END SUBROUTINE fld_interp
1490
1491
1492   FUNCTION ksec_week( cdday )
1493      !!---------------------------------------------------------------------
1494      !!                    ***  FUNCTION kshift_week ***
1495      !!
1496      !! ** Purpose : 
1497      !!---------------------------------------------------------------------
1498      CHARACTER(len=*), INTENT(in)   ::   cdday   !3 first letters of the first day of the weekly file
1499      !!
1500      INTEGER                        ::   ksec_week  ! output variable
1501      INTEGER                        ::   ijul       !temp variable
1502      INTEGER                        ::   ishift     !temp variable
1503      CHARACTER(len=3),DIMENSION(7)  ::   cl_week 
1504      !!----------------------------------------------------------------------
1505      cl_week = (/"sun","sat","fri","thu","wed","tue","mon"/)
1506      DO ijul = 1, 7
1507         IF( cl_week(ijul) == TRIM(cdday) ) EXIT
1508      END DO
1509      IF( ijul .GT. 7 )   CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): '//TRIM(cdday) )
1510      !
1511      ishift = ijul * NINT(rday)
1512      !
1513      ksec_week = nsec_week + ishift
1514      ksec_week = MOD( ksec_week, 7*NINT(rday) )
1515      !
1516   END FUNCTION ksec_week
1517
1518   !!======================================================================
1519END MODULE fldread
Note: See TracBrowser for help on using the repository browser.