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

source: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 5777

Last change on this file since 5777 was 5777, checked in by gm, 9 years ago

#1593: LDF-ADV, III. Phasing of the improvements/simplifications of ADV & LDF momentum trends (see wiki page)

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