source: branches/UKMO/dev_r7573_xios_write/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 8310

Last change on this file since 8310 was 8310, checked in by andmirek, 4 years ago

#1882 works with AGRIF and few small fixes/changes

File size: 85.2 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   PUBLIC   fld_clopn
35
36   TYPE, PUBLIC ::   FLD_N      !: Namelist field informations
37      CHARACTER(len = 256) ::   clname      ! generic name of the NetCDF flux file
38      REAL(wp)             ::   nfreqh      ! frequency of each flux file
39      CHARACTER(len = 34)  ::   clvar       ! generic name of the variable in the NetCDF flux file
40      LOGICAL              ::   ln_tint     ! time interpolation or not (T/F)
41      LOGICAL              ::   ln_clim     ! climatology or not (T/F)
42      CHARACTER(len = 8)   ::   cltype      ! type of data file 'daily', 'monthly' or yearly'
43      CHARACTER(len = 256) ::   wname       ! generic name of a NetCDF weights file to be used, blank if not
44      CHARACTER(len = 34)  ::   vcomp       ! symbolic component name if a vector that needs rotation
45      !                                     ! a string starting with "U" or "V" for each component   
46      !                                     ! chars 2 onwards identify which components go together 
47      CHARACTER(len = 34)  ::   lname       ! generic name of a NetCDF land/sea mask file to be used, blank if not
48      !                                     ! 0=sea 1=land
49   END TYPE FLD_N
50
51   TYPE, PUBLIC ::   FLD        !: Input field related variables
52      CHARACTER(len = 256)            ::   clrootname   ! generic name of the NetCDF file
53      CHARACTER(len = 256)            ::   clname       ! current name of the NetCDF file
54      REAL(wp)                        ::   nfreqh       ! frequency of each flux file
55      CHARACTER(len = 34)             ::   clvar        ! generic name of the variable in the NetCDF flux file
56      LOGICAL                         ::   ln_tint      ! time interpolation or not (T/F)
57      LOGICAL                         ::   ln_clim      ! climatology or not (T/F)
58      CHARACTER(len = 8)              ::   cltype       ! type of data file 'daily', 'monthly' or yearly'
59      INTEGER                         ::   num          ! iom id of the jpfld files to be read
60      INTEGER , DIMENSION(2)          ::   nrec_b       ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year)
61      INTEGER , DIMENSION(2)          ::   nrec_a       ! after  record (1: index, 2: second since Jan. 1st 00h of nit000 year)
62      REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:  ) ::   fnow       ! input fields interpolated to now time step
63      REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) ::   fdta       ! 2 consecutive record of input fields
64      CHARACTER(len = 256)            ::   wgtname      ! current name of the NetCDF weight file acting as a key
65      !                                                 ! into the WGTLIST structure
66      CHARACTER(len = 34)             ::   vcomp        ! symbolic name for a vector component that needs rotation
67      LOGICAL, DIMENSION(2)           ::   rotn         ! flag to indicate whether before/after field has been rotated
68      INTEGER                         ::   nreclast     ! last record to be read in the current file
69      CHARACTER(len = 256)            ::   lsmname      ! current name of the NetCDF mask file acting as a key
70   END TYPE FLD
71
72   TYPE, PUBLIC ::   MAP_POINTER      !: Map from input data file to local domain
73      INTEGER, POINTER, DIMENSION(:)  ::  ptr           ! Array of integer pointers to 1D arrays
74      LOGICAL                         ::  ll_unstruc    ! Unstructured (T) or structured (F) boundary data file
75   END TYPE MAP_POINTER
76
77!$AGRIF_DO_NOT_TREAT
78
79   !! keep list of all weights variables so they're only read in once
80   !! need to add AGRIF directives not to process this structure
81   !! also need to force wgtname to include AGRIF nest number
82   TYPE         ::   WGT        !: Input weights related variables
83      CHARACTER(len = 256)                    ::   wgtname      ! current name of the NetCDF weight file
84      INTEGER , DIMENSION(2)                  ::   ddims        ! shape of input grid
85      INTEGER , DIMENSION(2)                  ::   botleft      ! top left corner of box in input grid containing
86      !                                                         ! current processor grid
87      INTEGER , DIMENSION(2)                  ::   topright     ! top right corner of box
88      INTEGER                                 ::   jpiwgt       ! width of box on input grid
89      INTEGER                                 ::   jpjwgt       ! height of box on input grid
90      INTEGER                                 ::   numwgt       ! number of weights (4=bilinear, 16=bicubic)
91      INTEGER                                 ::   nestid       ! for agrif, keep track of nest we're in
92      INTEGER                                 ::   overlap      ! =0 when cyclic grid has no overlapping EW columns
93      !                                                         ! =>1 when they have one or more overlapping columns     
94      !                                                         ! =-1 not cyclic
95      LOGICAL                                 ::   cyclic       ! east-west cyclic or not
96      INTEGER,  DIMENSION(:,:,:), POINTER     ::   data_jpi     ! array of source integers
97      INTEGER,  DIMENSION(:,:,:), POINTER     ::   data_jpj     ! array of source integers
98      REAL(wp), DIMENSION(:,:,:), POINTER     ::   data_wgt     ! array of weights on model grid
99      REAL(wp), DIMENSION(:,:,:), POINTER     ::   fly_dta      ! array of values on input grid
100      REAL(wp), DIMENSION(:,:,:), POINTER     ::   col          ! temporary array for reading in columns
101   END TYPE WGT
102
103   INTEGER,     PARAMETER             ::   tot_wgts = 10
104   TYPE( WGT ), DIMENSION(tot_wgts)   ::   ref_wgts     ! array of wgts
105   INTEGER                            ::   nxt_wgt = 1  ! point to next available space in ref_wgts array
106   REAL(wp), PARAMETER                ::   undeff_lsm = -999.00_wp
107
108!$AGRIF_END_DO_NOT_TREAT
109
110   !!----------------------------------------------------------------------
111   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
112   !! $Id$
113   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
114   !!----------------------------------------------------------------------
115CONTAINS
116
117   SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset )
118      !!---------------------------------------------------------------------
119      !!                    ***  ROUTINE fld_read  ***
120      !!                   
121      !! ** Purpose :   provide at each time step the surface ocean fluxes
122      !!                (momentum, heat, freshwater and runoff)
123      !!
124      !! ** Method  :   READ each input fields in NetCDF files using IOM
125      !!      and intepolate it to the model time-step.
126      !!         Several assumptions are made on the input file:
127      !!      blahblahblah....
128      !!----------------------------------------------------------------------
129      INTEGER  , INTENT(in   )               ::   kt        ! ocean time step
130      INTEGER  , INTENT(in   )               ::   kn_fsbc   ! sbc computation period (in time step)
131      TYPE(FLD), INTENT(inout), DIMENSION(:) ::   sd        ! input field related variables
132      TYPE(MAP_POINTER),INTENT(in), OPTIONAL, DIMENSION(:) ::   map   ! global-to-local mapping indices
133      INTEGER  , INTENT(in   ), OPTIONAL     ::   kit       ! subcycle timestep for timesplitting option
134      INTEGER  , INTENT(in   ), OPTIONAL     ::   kt_offset ! provide fields at time other than "now"
135                                                            !   kt_offset = -1 => fields at "before" time level
136                                                            !   kt_offset = +1 => fields at "after"  time level
137                                                            !   etc.
138      !!
139      INTEGER  ::   itmp       ! temporary variable
140      INTEGER  ::   imf        ! size of the structure sd
141      INTEGER  ::   jf         ! dummy indices
142      INTEGER  ::   isecend    ! number of second since Jan. 1st 00h of nit000 year at nitend
143      INTEGER  ::   isecsbc    ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step
144      INTEGER  ::   it_offset  ! local time offset variable
145      LOGICAL  ::   llnxtyr    ! open next year  file?
146      LOGICAL  ::   llnxtmth   ! open next month file?
147      LOGICAL  ::   llstop     ! stop is the file does not exist
148      LOGICAL  ::   ll_firstcall ! true if this is the first call to fld_read for this set of fields
149      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation
150      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation
151      CHARACTER(LEN=1000) ::   clfmt   ! write format
152      TYPE(MAP_POINTER) ::   imap   ! global-to-local mapping indices
153      !!---------------------------------------------------------------------
154      ll_firstcall = kt == nit000
155      IF( PRESENT(kit) )   ll_firstcall = ll_firstcall .and. kit == 1
156
157      IF ( nn_components == jp_iam_sas ) THEN   ;   it_offset = nn_fsbc
158      ELSE                                      ;   it_offset = 0
159      ENDIF
160      IF( PRESENT(kt_offset) )   it_offset = kt_offset
161
162      imap%ptr => NULL()
163
164      ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
165      IF( present(kit) ) THEN   ! ignore kn_fsbc in this case
166         isecsbc = nsec_year + nsec1jan000 + (kit+it_offset)*NINT( rdt/REAL(nn_baro,wp) )
167      ELSE                      ! middle of sbc time step
168         isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdttra(1)) + it_offset * NINT(rdttra(1))
169      ENDIF
170      imf = SIZE( sd )
171      !
172      IF( ll_firstcall ) THEN                      ! initialization
173         DO jf = 1, imf 
174            IF( PRESENT(map) ) imap = map(jf)
175            CALL fld_init( kn_fsbc, sd(jf), imap )  ! read each before field (put them in after as they will be swapped)
176         END DO
177         IF( lwp ) CALL wgt_print()                ! control print
178      ENDIF
179      !                                            ! ====================================== !
180      IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! update field at each kn_fsbc time-step !
181         !                                         ! ====================================== !
182         !
183         DO jf = 1, imf                            ! ---   loop over field   --- !
184           
185            IF( isecsbc > sd(jf)%nrec_a(2) .OR. ll_firstcall ) THEN    ! read/update the after data?
186
187               IF( PRESENT(map) )   imap = map(jf)   ! temporary definition of map
188
189               sd(jf)%nrec_b(:) = sd(jf)%nrec_a(:)                                  ! swap before record informations
190               sd(jf)%rotn(1) = sd(jf)%rotn(2)                                      ! swap before rotate informations
191               IF( sd(jf)%ln_tint )   sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)   ! swap before record field
192
193               CALL fld_rec( kn_fsbc, sd(jf), kt_offset = it_offset, kit = kit )    ! update after record informations
194
195               ! if kn_fsbc*rdttra is larger than nfreqh (which is kind of odd),
196               ! it is possible that the before value is no more the good one... we have to re-read it
197               ! if before is not the last record of the file currently opened and after is the first record to be read
198               ! in a new file which means after = 1 (the file to be opened corresponds to the current time)
199               ! or after = nreclast + 1 (the file to be opened corresponds to a future time step)
200               IF( .NOT. ll_firstcall .AND. sd(jf)%ln_tint .AND. sd(jf)%nrec_b(1) /= sd(jf)%nreclast &
201                  &                   .AND. MOD( sd(jf)%nrec_a(1), sd(jf)%nreclast ) == 1 ) THEN
202                  itmp = sd(jf)%nrec_a(1)                       ! temporary storage
203                  sd(jf)%nrec_a(1) = sd(jf)%nreclast            ! read the last record of the file currently opened
204                  CALL fld_get( sd(jf), imap )                  ! read after data
205                  sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)   ! re-swap before record field
206                  sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1)           ! update before record informations
207                  sd(jf)%nrec_b(2) = sd(jf)%nrec_a(2) - NINT( sd(jf)%nfreqh * 3600 )  ! assume freq to be in hours in this case
208                  sd(jf)%rotn(1)   = sd(jf)%rotn(2)             ! update before rotate informations
209                  sd(jf)%nrec_a(1) = itmp                       ! move back to after record
210               ENDIF
211
212               CALL fld_clopn( sd(jf) )   ! Do we need to open a new year/month/week/day file?
213               
214               IF( sd(jf)%ln_tint ) THEN
215                 
216                  ! if kn_fsbc*rdttra is larger than nfreqh (which is kind of odd),
217                  ! it is possible that the before value is no more the good one... we have to re-read it
218                  ! if before record is not just just before the after record...
219                  IF( .NOT. ll_firstcall .AND. MOD( sd(jf)%nrec_a(1), sd(jf)%nreclast ) /= 1 &
220                     &                   .AND. sd(jf)%nrec_b(1) /= sd(jf)%nrec_a(1) - 1 ) THEN   
221                     sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) - 1       ! move back to before record
222                     CALL fld_get( sd(jf), imap )                  ! read after data
223                     sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)   ! re-swap before record field
224                     sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1)           ! update before record informations
225                     sd(jf)%nrec_b(2) = sd(jf)%nrec_a(2) - NINT( sd(jf)%nfreqh * 3600 )  ! assume freq to be in hours in this case
226                     sd(jf)%rotn(1)   = sd(jf)%rotn(2)             ! update before rotate informations
227                     sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) + 1       ! move back to after record
228                  ENDIF
229
230                  ! do we have to change the year/month/week/day of the forcing field??
231                  ! if we do time interpolation we will need to open next year/month/week/day file before the end of the current
232                  ! one. If so, we are still before the end of the year/month/week/day when calling fld_rec so sd(jf)%nrec_a(1)
233                  ! will be larger than the record number that should be read for current year/month/week/day
234                  ! do we need next file data?
235                  IF( sd(jf)%nrec_a(1) > sd(jf)%nreclast ) THEN
236                     
237                     sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) - sd(jf)%nreclast   !
238                     
239                     IF( .NOT. ( sd(jf)%ln_clim .AND. sd(jf)%cltype == 'yearly' ) ) THEN   ! close/open the current/new file
240                       
241                        llnxtmth = sd(jf)%cltype == 'monthly' .OR. nday == nmonth_len(nmonth)      ! open next month file?
242                        llnxtyr  = sd(jf)%cltype == 'yearly'  .OR. (nmonth == 12 .AND. llnxtmth)   ! open next year  file?
243
244                        ! if the run finishes at the end of the current year/month/week/day, we will allow next
245                        ! year/month/week/day file to be not present. If the run continue further than the current
246                        ! year/month/week/day, next year/month/week/day file must exist
247                        isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdttra(1))   ! second at the end of the run
248                        llstop = isecend > sd(jf)%nrec_a(2)                                   ! read more than 1 record of next year
249                        ! we suppose that the date of next file is next day (should be ok even for weekly files...)
250                        CALL fld_clopn( sd(jf), nyear  + COUNT((/llnxtyr /))                                           ,         &
251                           &                    nmonth + COUNT((/llnxtmth/)) - 12                 * COUNT((/llnxtyr /)),         &
252                           &                    nday   + 1                   - nmonth_len(nmonth) * COUNT((/llnxtmth/)), llstop )
253
254                        IF( sd(jf)%num <= 0 .AND. .NOT. llstop ) THEN    ! next year file does not exist
255                           CALL ctl_warn('next year/month/week/day file: '//TRIM(sd(jf)%clname)//     &
256                              &     ' not present -> back to current year/month/day')
257                           CALL fld_clopn( sd(jf) )       ! back to the current year/month/day
258                           sd(jf)%nrec_a(1) = sd(jf)%nreclast     ! force to read the last record in the current year file
259                        ENDIF
260                       
261                     ENDIF
262                  ENDIF   ! open need next file?
263                 
264               ENDIF   ! temporal interpolation?
265
266               ! read after data
267               CALL fld_get( sd(jf), imap )
268
269            ENDIF   ! read new data?
270         END DO                                    ! --- end loop over field --- !
271
272         CALL fld_rot( kt, sd )                    ! rotate vector before/now/after fields if needed
273
274         DO jf = 1, imf                            ! ---   loop over field   --- !
275            !
276            IF( sd(jf)%ln_tint ) THEN              ! temporal interpolation
277               IF(lwp .AND. kt - nit000 <= 100 ) THEN
278                  clfmt = "('fld_read: var ', a, ' kt = ', i8, ' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
279                     &    "', records b/a: ', i6.4, '/', i6.4, ' (days ', f9.4,'/', f9.4, ')')"
280                  WRITE(numout, clfmt)  TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,   &           
281                     & sd(jf)%nrec_b(1), sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday
282                  WRITE(numout, *) 'it_offset is : ',it_offset
283               ENDIF
284               ! temporal interpolation weights
285               ztinta =  REAL( isecsbc - sd(jf)%nrec_b(2), wp ) / REAL( sd(jf)%nrec_a(2) - sd(jf)%nrec_b(2), wp )
286               ztintb =  1. - ztinta
287!CDIR COLLAPSE
288               sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:,1) + ztinta * sd(jf)%fdta(:,:,:,2)
289            ELSE   ! nothing to do...
290               IF(lwp .AND. kt - nit000 <= 100 ) THEN
291                  clfmt = "('fld_read: var ', a, ' kt = ', i8,' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
292                     &    "', record: ', i6.4, ' (days ', f9.4, ' <-> ', f9.4, ')')"
293                  WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,    &
294                     &                 sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday
295               ENDIF
296            ENDIF
297            !
298            IF( kt == nitend - kn_fsbc + 1 )   CALL iom_close( sd(jf)%num )   ! Close the input files
299
300         END DO                                    ! --- end loop over field --- !
301         !
302         !                                         ! ====================================== !
303      ENDIF                                        ! update field at each kn_fsbc time-step !
304      !                                            ! ====================================== !
305      !
306   END SUBROUTINE fld_read
307
308
309   SUBROUTINE fld_init( kn_fsbc, sdjf, map )
310      !!---------------------------------------------------------------------
311      !!                    ***  ROUTINE fld_init  ***
312      !!
313      !! ** Purpose :  - first call to fld_rec to define before values
314      !!               - if time interpolation, read before data
315      !!----------------------------------------------------------------------
316      INTEGER  , INTENT(in   ) ::   kn_fsbc   ! sbc computation period (in time step)
317      TYPE(FLD), INTENT(inout) ::   sdjf      ! input field related variables
318      TYPE(MAP_POINTER),INTENT(in) ::   map   ! global-to-local mapping indices
319      !!
320      LOGICAL :: llprevyr              ! are we reading previous year  file?
321      LOGICAL :: llprevmth             ! are we reading previous month file?
322      LOGICAL :: llprevweek            ! are we reading previous week  file?
323      LOGICAL :: llprevday             ! are we reading previous day   file?
324      LOGICAL :: llprev                ! llprevyr .OR. llprevmth .OR. llprevweek .OR. llprevday
325      INTEGER :: idvar                 ! variable id
326      INTEGER :: inrec                 ! number of record existing for this variable
327      INTEGER :: iyear, imonth, iday   ! first day of the current file in yyyy mm dd
328      INTEGER :: isec_week             ! number of seconds since start of the weekly file
329      CHARACTER(LEN=1000) ::   clfmt   ! write format
330      !!---------------------------------------------------------------------
331      llprevyr   = .FALSE.
332      llprevmth  = .FALSE.
333      llprevweek = .FALSE.
334      llprevday  = .FALSE.
335      isec_week  = 0
336           
337      ! define record informations
338      CALL fld_rec( kn_fsbc, sdjf, ldbefore = .TRUE. )  ! return before values in sdjf%nrec_a (as we will swap it later)
339
340      ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
341
342      IF( sdjf%ln_tint ) THEN ! we need to read the previous record and we will put it in the current record structure
343
344         IF( sdjf%nrec_a(1) == 0  ) THEN   ! we redefine record sdjf%nrec_a(1) with the last record of previous year file
345            IF    ( sdjf%nfreqh == -12 ) THEN   ! yearly mean
346               IF( sdjf%cltype == 'yearly' ) THEN             ! yearly file
347                  sdjf%nrec_a(1) = 1                                                       ! force to read the unique record
348                  llprevyr  = .NOT. sdjf%ln_clim                                           ! use previous year  file?
349               ELSE
350                  CALL ctl_stop( "fld_init: yearly mean file must be in a yearly type of file: "//TRIM(sdjf%clrootname) )
351               ENDIF
352            ELSEIF( sdjf%nfreqh ==  -1 ) THEN   ! monthly mean
353               IF( sdjf%cltype == 'monthly' ) THEN            ! monthly file
354                  sdjf%nrec_a(1) = 1                                                       ! force to read the unique record
355                  llprevmth = .TRUE.                                                       ! use previous month file?
356                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file?
357               ELSE                                           ! yearly file
358                  sdjf%nrec_a(1) = 12                                                      ! force to read december mean
359                  llprevyr = .NOT. sdjf%ln_clim                                            ! use previous year  file?
360               ENDIF
361            ELSE                                ! higher frequency mean (in hours)
362               IF    ( sdjf%cltype      == 'monthly' ) THEN   ! monthly file
363                  sdjf%nrec_a(1) = NINT( 24 * nmonth_len(nmonth-1) / sdjf%nfreqh )         ! last record of previous month
364                  llprevmth = .TRUE.                                                       ! use previous month file?
365                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file?
366               ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ! weekly file
367                  llprevweek = .TRUE.                                                      ! use previous week  file?
368                  sdjf%nrec_a(1) = NINT( 24 * 7 / sdjf%nfreqh )                            ! last record of previous week
369                  isec_week = NINT(rday) * 7                                               ! add a shift toward previous week
370               ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ! daily file
371                  sdjf%nrec_a(1) = NINT( 24 / sdjf%nfreqh )                                ! last record of previous day
372                  llprevday = .TRUE.                                                       ! use previous day   file?
373                  llprevmth = llprevday .AND. nday   == 1                                  ! use previous month file?
374                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file?
375               ELSE                                           ! yearly file
376                  sdjf%nrec_a(1) = NINT( 24 * nyear_len(0) / sdjf%nfreqh )                 ! last record of previous year
377                  llprevyr = .NOT. sdjf%ln_clim                                            ! use previous year  file?
378               ENDIF
379            ENDIF
380         ENDIF
381         !
382         IF ( sdjf%cltype(1:4) == 'week' ) THEN
383            isec_week = isec_week + ksec_week( sdjf%cltype(6:8) )   ! second since the beginning of the week
384            llprevmth = isec_week > nsec_month                      ! longer time since the beginning of the week than the month
385            llprevyr  = llprevmth .AND. nmonth == 1
386         ENDIF
387         llprev = llprevyr .OR. llprevmth .OR. llprevweek .OR. llprevday
388         !
389         iyear  = nyear  - COUNT((/llprevyr /))
390         imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
391         iday   = nday   - COUNT((/llprevday/)) + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
392         !
393         CALL fld_clopn( sdjf, iyear, imonth, iday, .NOT. llprev )
394
395         ! if previous year/month/day file does not exist, we switch to the current year/month/day
396         IF( llprev .AND. sdjf%num <= 0 ) THEN
397            CALL ctl_warn( 'previous year/month/week/day file: '//TRIM(sdjf%clrootname)//   &
398               &           ' not present -> back to current year/month/week/day' )
399            ! we force to read the first record of the current year/month/day instead of last record of previous year/month/day
400            llprev = .FALSE.
401            sdjf%nrec_a(1) = 1
402            CALL fld_clopn( sdjf )
403         ENDIF
404         
405         IF( llprev ) THEN   ! check if the record sdjf%nrec_a(1) exists in the file
406            idvar = iom_varid( sdjf%num, sdjf%clvar )                                        ! id of the variable sdjf%clvar
407            IF( idvar <= 0 )   RETURN
408            inrec = iom_file( sdjf%num )%dimsz( iom_file( sdjf%num )%ndims(idvar), idvar )   ! size of the last dim of idvar
409            sdjf%nrec_a(1) = MIN( sdjf%nrec_a(1), inrec )   ! make sure we select an existing record
410         ENDIF
411
412         ! read before data in after arrays(as we will swap it later)
413         CALL fld_get( sdjf, map )
414
415         clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')"
416         IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_a(1), REAL(sdjf%nrec_a(2),wp)/rday
417
418      ENDIF
419      !
420   END SUBROUTINE fld_init
421
422
423   SUBROUTINE fld_rec( kn_fsbc, sdjf, ldbefore, kit, kt_offset )
424      !!---------------------------------------------------------------------
425      !!                    ***  ROUTINE fld_rec  ***
426      !!
427      !! ** Purpose : Compute
428      !!              if sdjf%ln_tint = .TRUE.
429      !!                  nrec_a: record number and its time (nrec_b is obtained from nrec_a when swapping)
430      !!              if sdjf%ln_tint = .FALSE.
431      !!                  nrec_a(1): record number
432      !!                  nrec_b(2) and nrec_a(2): time of the beginning and end of the record (for print only)
433      !!----------------------------------------------------------------------
434      INTEGER  , INTENT(in   )           ::   kn_fsbc   ! sbc computation period (in time step)
435      TYPE(FLD), INTENT(inout)           ::   sdjf      ! input field related variables
436      LOGICAL  , INTENT(in   ), OPTIONAL ::   ldbefore  ! sent back before record values (default = .FALSE.)
437      INTEGER  , INTENT(in   ), OPTIONAL ::   kit       ! index of barotropic subcycle
438                                                        ! used only if sdjf%ln_tint = .TRUE.
439      INTEGER  , INTENT(in   ), OPTIONAL ::   kt_offset ! Offset of required time level compared to "now"
440                                                        !   time level in units of time steps.
441      !!
442      LOGICAL  ::   llbefore    ! local definition of ldbefore
443      INTEGER  ::   iendrec     ! end of this record (in seconds)
444      INTEGER  ::   imth        ! month number
445      INTEGER  ::   ifreq_sec   ! frequency mean (in seconds)
446      INTEGER  ::   isec_week   ! number of seconds since the start of the weekly file
447      INTEGER  ::   it_offset   ! local time offset variable
448      REAL(wp) ::   ztmp        ! temporary variable
449      !!----------------------------------------------------------------------
450      !
451      ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
452      !
453      IF( PRESENT(ldbefore) ) THEN   ;   llbefore = ldbefore .AND. sdjf%ln_tint   ! needed only if sdjf%ln_tint = .TRUE.
454      ELSE                           ;   llbefore = .FALSE.
455      ENDIF
456      !
457      IF ( nn_components == jp_iam_sas ) THEN   ;   it_offset = nn_fsbc
458      ELSE                                      ;   it_offset = 0
459      ENDIF
460      IF( PRESENT(kt_offset) )   it_offset = kt_offset
461      IF( PRESENT(kit) ) THEN   ;   it_offset = ( kit + it_offset ) * NINT( rdt/REAL(nn_baro,wp) )
462      ELSE                      ;   it_offset =         it_offset   * NINT(       rdttra(1)      )
463      ENDIF
464      !
465      !                                      ! =========== !
466      IF    ( sdjf%nfreqh == -12 ) THEN      ! yearly mean
467         !                                   ! =========== !
468         !
469         IF( sdjf%ln_tint ) THEN                 ! time interpolation, shift by 1/2 record
470            !
471            !                  INT( ztmp )
472            !                     /|\
473            !                    1 |    *----
474            !                    0 |----(             
475            !                      |----+----|--> time
476            !                      0   /|\   1   (nday/nyear_len(1))
477            !                           |   
478            !                           |   
479            !       forcing record :    1
480            !                           
481            ztmp = REAL( nsec_year, wp ) / ( REAL( nyear_len(1), wp ) * rday ) + 0.5 &
482           &       + REAL( it_offset, wp ) / ( REAL( nyear_len(1), wp ) * rday )
483            sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
484            ! swap at the middle of the year
485            IF( llbefore ) THEN   ;   sdjf%nrec_a(2) = nsec1jan000 - (1 - INT(ztmp)) * NINT(0.5 * rday) * nyear_len(0) + &
486                                    & INT(ztmp) * NINT( 0.5 * rday) * nyear_len(1) 
487            ELSE                  ;   sdjf%nrec_a(2) = nsec1jan000 + (1 - INT(ztmp)) * NINT(0.5 * rday) * nyear_len(1) + &
488                                    & INT(ztmp) * INT(rday) * nyear_len(1) + INT(ztmp) * NINT( 0.5 * rday) * nyear_len(2) 
489            ENDIF
490         ELSE                                    ! no time interpolation
491            sdjf%nrec_a(1) = 1
492            sdjf%nrec_a(2) = NINT(rday) * nyear_len(1) + nsec1jan000   ! swap at the end    of the year
493            sdjf%nrec_b(2) = nsec1jan000                               ! beginning of the year (only for print)
494         ENDIF
495         !
496         !                                   ! ============ !
497      ELSEIF( sdjf%nfreqh ==  -1 ) THEN      ! monthly mean !
498         !                                   ! ============ !
499         !
500         IF( sdjf%ln_tint ) THEN                 ! time interpolation, shift by 1/2 record
501            !
502            !                  INT( ztmp )
503            !                     /|\
504            !                    1 |    *----
505            !                    0 |----(             
506            !                      |----+----|--> time
507            !                      0   /|\   1   (nday/nmonth_len(nmonth))
508            !                           |   
509            !                           |   
510            !       forcing record :  nmonth
511            !                           
512            ztmp = REAL( nsec_month, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) + 0.5 &
513           &       + REAL( it_offset, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday )
514            imth = nmonth + INT( ztmp ) - COUNT((/llbefore/))
515            IF( sdjf%cltype == 'monthly' ) THEN   ;   sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
516            ELSE                                  ;   sdjf%nrec_a(1) = imth
517            ENDIF
518            sdjf%nrec_a(2) = nmonth_half(   imth ) + nsec1jan000   ! swap at the middle of the month
519         ELSE                                    ! no time interpolation
520            IF( sdjf%cltype == 'monthly' ) THEN   ;   sdjf%nrec_a(1) = 1
521            ELSE                                  ;   sdjf%nrec_a(1) = nmonth
522            ENDIF
523            sdjf%nrec_a(2) =  nmonth_end(nmonth  ) + nsec1jan000   ! swap at the end    of the month
524            sdjf%nrec_b(2) =  nmonth_end(nmonth-1) + nsec1jan000   ! beginning of the month (only for print)
525         ENDIF
526         !
527         !                                   ! ================================ !
528      ELSE                                   ! higher frequency mean (in hours)
529         !                                   ! ================================ !
530         !
531         ifreq_sec = NINT( sdjf%nfreqh * 3600 )                                         ! frequency mean (in seconds)
532         IF( sdjf%cltype(1:4) == 'week' )   isec_week = ksec_week( sdjf%cltype(6:8) )   ! since the first day of the current week
533         ! number of second since the beginning of the file
534         IF(     sdjf%cltype      == 'monthly' ) THEN   ;   ztmp = REAL(nsec_month,wp)  ! since the first day of the current month
535         ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   ztmp = REAL(isec_week ,wp)  ! since the first day of the current week
536         ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   ztmp = REAL(nsec_day  ,wp)  ! since 00h of the current day
537         ELSE                                           ;   ztmp = REAL(nsec_year ,wp)  ! since 00h on Jan 1 of the current year
538         ENDIF
539         ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdttra(1) + REAL( it_offset, wp )  ! centrered in the middle of sbc time step
540         ztmp = ztmp + 0.01 * rdttra(1)                                                 ! avoid truncation error
541         IF( sdjf%ln_tint ) THEN                ! time interpolation, shift by 1/2 record
542            !
543            !          INT( ztmp/ifreq_sec + 0.5 )
544            !                     /|\
545            !                    2 |        *-----(
546            !                    1 |  *-----(
547            !                    0 |--(             
548            !                      |--+--|--+--|--+--|--> time
549            !                      0 /|\ 1 /|\ 2 /|\ 3    (ztmp/ifreq_sec)
550            !                         |     |     |
551            !                         |     |     |
552            !       forcing record :  1     2     3
553            !                   
554            ztmp= ztmp / REAL(ifreq_sec, wp) + 0.5
555         ELSE                                   ! no time interpolation
556            !
557            !           INT( ztmp/ifreq_sec )
558            !                     /|\
559            !                    2 |           *-----(
560            !                    1 |     *-----(
561            !                    0 |-----(             
562            !                      |--+--|--+--|--+--|--> time
563            !                      0 /|\ 1 /|\ 2 /|\ 3    (ztmp/ifreq_sec)
564            !                         |     |     |
565            !                         |     |     |
566            !       forcing record :  1     2     3
567            !                           
568            ztmp= ztmp / REAL(ifreq_sec, wp)
569         ENDIF
570         sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))   ! record number to be read
571
572         iendrec = ifreq_sec * sdjf%nrec_a(1) + nsec1jan000       ! end of this record (in second)
573         ! add the number of seconds between 00h Jan 1 and the end of previous month/week/day (ok if nmonth=1)
574         IF( sdjf%cltype      == 'monthly' )   iendrec = iendrec + NINT(rday) * SUM(nmonth_len(1:nmonth -1))
575         IF( sdjf%cltype(1:4) == 'week'    )   iendrec = iendrec + ( nsec_year - isec_week )
576         IF( sdjf%cltype      == 'daily'   )   iendrec = iendrec + NINT(rday) * ( nday_year - 1 )
577         IF( sdjf%ln_tint ) THEN
578             sdjf%nrec_a(2) = iendrec - ifreq_sec / 2        ! swap at the middle of the record
579         ELSE
580             sdjf%nrec_a(2) = iendrec                        ! swap at the end    of the record
581             sdjf%nrec_b(2) = iendrec - ifreq_sec            ! beginning of the record (only for print)
582         ENDIF
583         !
584      ENDIF
585      !
586   END SUBROUTINE fld_rec
587
588
589   SUBROUTINE fld_get( sdjf, map )
590      !!---------------------------------------------------------------------
591      !!                    ***  ROUTINE fld_get  ***
592      !!
593      !! ** Purpose :   read the data
594      !!----------------------------------------------------------------------
595      TYPE(FLD), INTENT(inout) ::   sdjf   ! input field related variables
596      TYPE(MAP_POINTER),INTENT(in) ::   map   ! global-to-local mapping indices
597      !!
598      INTEGER                  ::   ipk    ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
599      INTEGER                  ::   iw     ! index into wgts array
600      INTEGER                  ::   ipdom  ! index of the domain
601      INTEGER                  ::   idvar  ! variable ID
602      INTEGER                  ::   idmspc ! number of spatial dimensions
603      LOGICAL                  ::   lmoor  ! C1D case: point data
604      !!---------------------------------------------------------------------
605      !
606      ipk = SIZE( sdjf%fnow, 3 )
607      !
608      IF( ASSOCIATED(map%ptr) ) THEN
609         IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map )
610         ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1), map )
611         ENDIF
612      ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
613         CALL wgt_list( sdjf, iw )
614         IF( sdjf%ln_tint ) THEN   ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk  , sdjf%fdta(:,:,:,2), & 
615              & sdjf%nrec_a(1), sdjf%lsmname )
616         ELSE                      ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk  , sdjf%fnow(:,:,:  ), &
617              & sdjf%nrec_a(1), sdjf%lsmname )
618         ENDIF
619      ELSE
620         IF( SIZE(sdjf%fnow, 1) == jpi ) THEN  ;  ipdom = jpdom_data
621         ELSE                                  ;  ipdom = jpdom_unknown
622         ENDIF
623         ! C1D case: If product of spatial dimensions == ipk, then x,y are of
624         ! size 1 (point/mooring data): this must be read onto the central grid point
625         idvar  = iom_varid( sdjf%num, sdjf%clvar )
626         idmspc = iom_file( sdjf%num )%ndims( idvar )
627         IF( iom_file( sdjf%num )%luld( idvar ) )   idmspc = idmspc - 1
628         lmoor  = (idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk)
629         !
630         SELECT CASE( ipk )
631         CASE(1)
632            IF( lk_c1d .AND. lmoor ) THEN
633               IF( sdjf%ln_tint ) THEN
634                  CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fdta(2,2,1,2), sdjf%nrec_a(1) )
635                  CALL lbc_lnk( sdjf%fdta(:,:,1,2),'Z',1. )
636               ELSE
637                  CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fnow(2,2,1  ), sdjf%nrec_a(1) )
638                  CALL lbc_lnk( sdjf%fnow(:,:,1  ),'Z',1. )
639               ENDIF
640            ELSE
641               IF( sdjf%ln_tint ) THEN   ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,1,2), sdjf%nrec_a(1) )
642               ELSE                      ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,1  ), sdjf%nrec_a(1) )
643               ENDIF
644            ENDIF
645         CASE DEFAULT
646            IF (lk_c1d .AND. lmoor ) THEN
647               IF( sdjf%ln_tint ) THEN
648                  CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fdta(2,2,:,2), sdjf%nrec_a(1) )
649                  CALL lbc_lnk( sdjf%fdta(:,:,:,2),'Z',1. )
650               ELSE
651                  CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fnow(2,2,:  ), sdjf%nrec_a(1) )
652                  CALL lbc_lnk( sdjf%fnow(:,:,:  ),'Z',1. )
653               ENDIF
654            ELSE
655               IF( sdjf%ln_tint ) THEN   ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) )
656               ELSE                      ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1) )
657               ENDIF
658            ENDIF
659         END SELECT
660      ENDIF
661      !
662      sdjf%rotn(2) = .false.   ! vector not yet rotated
663
664   END SUBROUTINE fld_get
665
666   SUBROUTINE fld_map( num, clvar, dta, nrec, map )
667      !!---------------------------------------------------------------------
668      !!                    ***  ROUTINE fld_map  ***
669      !!
670      !! ** Purpose :   read global data from file and map onto local data
671      !!                using a general mapping (for open boundaries)
672      !!----------------------------------------------------------------------
673#if defined key_bdy
674      USE bdy_oce, ONLY:  dta_global, dta_global2         ! workspace to read in global data arrays
675#endif
676      INTEGER                   , INTENT(in ) ::   num     ! stream number
677      CHARACTER(LEN=*)          , INTENT(in ) ::   clvar   ! variable name
678      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta   ! output field on model grid (2 dimensional)
679      INTEGER                   , INTENT(in ) ::   nrec    ! record number to read (ie time slice)
680      TYPE(MAP_POINTER)         , INTENT(in ) ::   map     ! global-to-local mapping indices
681      !!
682      INTEGER                                 ::   ipi      ! length of boundary data on local process
683      INTEGER                                 ::   ipj      ! length of dummy dimension ( = 1 )
684      INTEGER                                 ::   ipk      ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk )
685      INTEGER                                 ::   ilendta  ! length of data in file
686      INTEGER                                 ::   idvar    ! variable ID
687      INTEGER                                 ::   ib, ik, ji, jj   ! loop counters
688      INTEGER                                 ::   ierr
689      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read  ! work space for global data
690      !!---------------------------------------------------------------------
691           
692      ipi = SIZE( dta, 1 )
693      ipj = 1
694      ipk = SIZE( dta, 3 )
695
696      idvar   = iom_varid( num, clvar )
697      ilendta = iom_file(num)%dimsz(1,idvar)
698
699#if defined key_bdy
700      ipj = iom_file(num)%dimsz(2,idvar)
701      IF ( map%ll_unstruc) THEN ! unstructured open boundary data file
702         dta_read => dta_global
703      ELSE                      ! structured open boundary data file
704         dta_read => dta_global2
705      ENDIF
706#endif
707
708      IF(lwp) WRITE(numout,*) 'Dim size for ',TRIM(clvar),' is ', ilendta
709      IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk
710
711      SELECT CASE( ipk )
712      CASE(1)        ;   CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1    ), nrec )
713      CASE DEFAULT   ;   CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec )
714      END SELECT
715      !
716      IF ( map%ll_unstruc ) THEN ! unstructured open boundary data file
717         DO ib = 1, ipi
718            DO ik = 1, ipk
719               dta(ib,1,ik) =  dta_read(map%ptr(ib),1,ik)
720            END DO
721         END DO
722      ELSE                       ! structured open boundary data file
723         DO ib = 1, ipi
724            jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
725            ji=map%ptr(ib)-(jj-1)*ilendta
726            DO ik = 1, ipk
727               dta(ib,1,ik) =  dta_read(ji,jj,ik)
728            END DO
729         END DO
730      ENDIF
731
732   END SUBROUTINE fld_map
733
734
735   SUBROUTINE fld_rot( kt, sd )
736      !!---------------------------------------------------------------------
737      !!                    ***  ROUTINE fld_rot  ***
738      !!
739      !! ** Purpose :   Vector fields may need to be rotated onto the local grid direction
740      !!----------------------------------------------------------------------
741      INTEGER  , INTENT(in   )               ::   kt        ! ocean time step
742      TYPE(FLD), INTENT(inout), DIMENSION(:) ::   sd        ! input field related variables
743      !!
744      INTEGER                           ::   ju,jv,jk,jn  ! loop indices
745      INTEGER                           ::   imf          ! size of the structure sd
746      INTEGER                           ::   ill          ! character length
747      INTEGER                           ::   iv           ! indice of V component
748      REAL(wp), POINTER, DIMENSION(:,:) ::   utmp, vtmp   ! temporary arrays for vector rotation
749      CHARACTER (LEN=100)               ::   clcomp       ! dummy weight name
750      !!---------------------------------------------------------------------
751
752      CALL wrk_alloc( jpi,jpj, utmp, vtmp )
753
754      !! (sga: following code should be modified so that pairs arent searched for each time
755      !
756      imf = SIZE( sd )
757      DO ju = 1, imf
758         ill = LEN_TRIM( sd(ju)%vcomp )
759         DO jn = 2-COUNT((/sd(ju)%ln_tint/)), 2
760            IF( ill > 0 .AND. .NOT. sd(ju)%rotn(jn) ) THEN   ! find vector rotations required             
761               IF( sd(ju)%vcomp(1:1) == 'U' ) THEN      ! east-west component has symbolic name starting with 'U'
762                  ! look for the north-south component which has same symbolic name but with 'U' replaced with 'V'
763                  clcomp = 'V' // sd(ju)%vcomp(2:ill)   ! works even if ill == 1
764                  iv = -1
765                  DO jv = 1, imf
766                     IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) )   iv = jv
767                  END DO
768                  IF( iv > 0 ) THEN   ! fields ju and iv are two components which need to be rotated together
769                     DO jk = 1, SIZE( sd(ju)%fnow, 3 )
770                        IF( sd(ju)%ln_tint )THEN
771                           CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->i', utmp(:,:) )
772                           CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->j', vtmp(:,:) )
773                           sd(ju)%fdta(:,:,jk,jn) = utmp(:,:)   ;   sd(iv)%fdta(:,:,jk,jn) = vtmp(:,:)
774                        ELSE
775                           CALL rot_rep( sd(ju)%fnow(:,:,jk  ), sd(iv)%fnow(:,:,jk  ), 'T', 'en->i', utmp(:,:) )
776                           CALL rot_rep( sd(ju)%fnow(:,:,jk  ), sd(iv)%fnow(:,:,jk  ), 'T', 'en->j', vtmp(:,:) )
777                           sd(ju)%fnow(:,:,jk   ) = utmp(:,:)   ;   sd(iv)%fnow(:,:,jk   ) = vtmp(:,:)
778                        ENDIF
779                     END DO
780                     sd(ju)%rotn(jn) = .TRUE.               ! vector was rotated
781                     IF( lwp .AND. kt == nit000 )   WRITE(numout,*)   &
782                        &   'fld_read: vector pair ('//TRIM(sd(ju)%clvar)//', '//TRIM(sd(iv)%clvar)//') rotated on to model grid'
783                  ENDIF
784               ENDIF
785            ENDIF
786         END DO
787       END DO
788      !
789      CALL wrk_dealloc( jpi,jpj, utmp, vtmp )
790      !
791   END SUBROUTINE fld_rot
792
793
794   SUBROUTINE fld_clopn( sdjf, kyear, kmonth, kday, ldstop )
795      !!---------------------------------------------------------------------
796      !!                    ***  ROUTINE fld_clopn  ***
797      !!
798      !! ** Purpose :   update the file name and open the file
799      !!----------------------------------------------------------------------
800      TYPE(FLD)        , INTENT(inout) ::   sdjf     ! input field related variables
801      INTEGER, OPTIONAL, INTENT(in   ) ::   kyear    ! year value
802      INTEGER, OPTIONAL, INTENT(in   ) ::   kmonth   ! month value
803      INTEGER, OPTIONAL, INTENT(in   ) ::   kday     ! day value
804      LOGICAL, OPTIONAL, INTENT(in   ) ::   ldstop   ! stop if open to read a non-existing file (default = .TRUE.)
805      !!
806      LOGICAL :: llprevyr              ! are we reading previous year  file?
807      LOGICAL :: llprevmth             ! are we reading previous month file?
808      INTEGER :: iyear, imonth, iday   ! first day of the current file in yyyy mm dd
809      INTEGER :: isec_week             ! number of seconds since start of the weekly file
810      INTEGER :: indexyr               ! year undex (O/1/2: previous/current/next)
811      INTEGER :: iyear_len, imonth_len ! length (days) of iyear and imonth             !
812      CHARACTER(len = 256)::   clname  ! temporary file name
813      !!----------------------------------------------------------------------
814      IF( PRESENT(kyear) ) THEN                             ! use given values
815         iyear = kyear
816         imonth = kmonth
817         iday = kday
818         IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week
819            isec_week = ksec_week( sdjf%cltype(6:8) )- (86400 * 8 ) 
820            llprevmth  = isec_week > nsec_month             ! longer time since beginning of the week than the month
821            llprevyr   = llprevmth .AND. nmonth == 1
822            iyear  = nyear  - COUNT((/llprevyr /))
823            imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
824            iday   = nday   + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
825         ENDIF
826      ELSE                                                  ! use current day values
827         IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week
828            isec_week  = ksec_week( sdjf%cltype(6:8) )      ! second since the beginning of the week
829            llprevmth  = isec_week > nsec_month             ! longer time since beginning of the week than the month
830            llprevyr   = llprevmth .AND. nmonth == 1
831         ELSE
832            isec_week  = 0
833            llprevmth  = .FALSE.
834            llprevyr   = .FALSE.
835         ENDIF
836         iyear  = nyear  - COUNT((/llprevyr /))
837         imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
838         iday   = nday   + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
839      ENDIF
840
841      ! build the new filename if not climatological data
842      clname=TRIM(sdjf%clrootname)
843      !
844      ! note that sdjf%ln_clim is is only acting on the presence of the year in the file name
845      IF( .NOT. sdjf%ln_clim ) THEN   
846                                         WRITE(clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), iyear    ! add year
847         IF( sdjf%cltype /= 'yearly' )   WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname          ), imonth   ! add month
848      ELSE
849         ! build the new filename if climatological data
850         IF( sdjf%cltype /= 'yearly' )   WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), imonth   ! add month
851      ENDIF
852      IF( sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) &
853            &                            WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname          ), iday     ! add day
854      !
855      IF( TRIM(clname) /= TRIM(sdjf%clname) .OR. sdjf%num == 0 ) THEN   ! new file to be open
856
857         sdjf%clname = TRIM(clname)
858         IF( sdjf%num /= 0 )   CALL iom_close( sdjf%num )   ! close file if already open
859         CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof =  LEN(TRIM(sdjf%wgtname)) > 0 )
860
861         ! find the last record to be read -> update sdjf%nreclast
862         indexyr = iyear - nyear + 1
863         iyear_len = nyear_len( indexyr )
864         SELECT CASE ( indexyr )
865         CASE ( 0 )   ;   imonth_len = 31   ! previous year -> imonth = 12
866         CASE ( 1 )   ;   imonth_len = nmonth_len(imonth) 
867         CASE ( 2 )   ;   imonth_len = 31   ! next     year -> imonth = 1
868         END SELECT
869         
870         ! last record to be read in the current file
871         IF    ( sdjf%nfreqh == -12 ) THEN                 ;   sdjf%nreclast = 1    !  yearly mean
872         ELSEIF( sdjf%nfreqh ==  -1 ) THEN                                          ! monthly mean
873            IF(     sdjf%cltype      == 'monthly' ) THEN   ;   sdjf%nreclast = 1
874            ELSE                                           ;   sdjf%nreclast = 12
875            ENDIF
876         ELSE                                                                       ! higher frequency mean (in hours)
877            IF(     sdjf%cltype      == 'monthly' ) THEN   ;   sdjf%nreclast = NINT( 24 * imonth_len / sdjf%nfreqh )
878            ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   sdjf%nreclast = NINT( 24 * 7          / sdjf%nfreqh )
879            ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   sdjf%nreclast = NINT( 24              / sdjf%nfreqh )
880            ELSE                                           ;   sdjf%nreclast = NINT( 24 * iyear_len  / sdjf%nfreqh )
881            ENDIF
882         ENDIF
883         
884      ENDIF
885      !
886   END SUBROUTINE fld_clopn
887
888
889   SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )
890      !!---------------------------------------------------------------------
891      !!                    ***  ROUTINE fld_fill  ***
892      !!
893      !! ** Purpose :   fill sdf with sdf_n and control print
894      !!----------------------------------------------------------------------
895      TYPE(FLD)  , DIMENSION(:), INTENT(inout) ::   sdf        ! structure of input fields (file informations, fields read)
896      TYPE(FLD_N), DIMENSION(:), INTENT(in   ) ::   sdf_n      ! array of namelist information structures
897      CHARACTER(len=*)         , INTENT(in   ) ::   cdir       ! Root directory for location of flx files
898      CHARACTER(len=*)         , INTENT(in   ) ::   cdcaller   !
899      CHARACTER(len=*)         , INTENT(in   ) ::   cdtitle    !
900      CHARACTER(len=*)         , INTENT(in   ) ::   cdnam      !
901      !
902      INTEGER  ::   jf       ! dummy indices
903      !!---------------------------------------------------------------------
904
905      DO jf = 1, SIZE(sdf)
906         sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
907         sdf(jf)%clname     = "not yet defined"
908         sdf(jf)%nfreqh     = sdf_n(jf)%nfreqh
909         sdf(jf)%clvar      = sdf_n(jf)%clvar
910         sdf(jf)%ln_tint    = sdf_n(jf)%ln_tint
911         sdf(jf)%ln_clim    = sdf_n(jf)%ln_clim
912         sdf(jf)%cltype     = sdf_n(jf)%cltype
913         sdf(jf)%num        = -1
914         sdf(jf)%wgtname    = " "
915         IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
916         sdf(jf)%lsmname = " "
917         IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 )   sdf(jf)%lsmname = TRIM( cdir )//TRIM( sdf_n(jf)%lname )
918         sdf(jf)%vcomp      = sdf_n(jf)%vcomp
919         sdf(jf)%rotn(:)    = .TRUE.   ! pretend to be rotated -> won't try to rotate data before the first call to fld_get
920         IF( sdf(jf)%cltype(1:4) == 'week' .AND. nn_leapy == 0  )   &
921            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs nn_leapy = 1')
922         IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim )   &
923            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.')
924         sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
925      END DO
926
927      IF(lwp) THEN      ! control print
928         WRITE(numout,*)
929         WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
930         WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
931         WRITE(numout,*) '          '//TRIM( cdnam )//' Namelist'
932         WRITE(numout,*) '          list of files and frequency (>0: in hours ; <0 in months)'
933         DO jf = 1, SIZE(sdf)
934            WRITE(numout,*) '               root filename: '  , TRIM( sdf(jf)%clrootname ),   &
935               &                          ' variable name: '  , TRIM( sdf(jf)%clvar      )
936            WRITE(numout,*) '               frequency: '      ,       sdf(jf)%nfreqh      ,   &
937               &                          ' time interp: '    ,       sdf(jf)%ln_tint     ,   &
938               &                          ' climatology: '    ,       sdf(jf)%ln_clim     ,   &
939               &                          ' weights    : '    , TRIM( sdf(jf)%wgtname    ),   &
940               &                          ' pairing    : '    , TRIM( sdf(jf)%vcomp      ),   &
941               &                          ' data type: '      ,       sdf(jf)%cltype      ,   &
942               &                          ' land/sea mask:'   , TRIM( sdf(jf)%lsmname    )
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                          ! 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=jpi2_lsm-jpi1_lsm+1
1358         itmpj=jpj2_lsm-jpj1_lsm+1
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.