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

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

source: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 5220

Last change on this file since 5220 was 5220, checked in by smasson, 9 years ago

dev_r5218_CNRS17_coupling: first update

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