source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 5752

Last change on this file since 5752 was 5752, checked in by jpaul, 5 years ago

fix tipying error, see ticket #1598

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