source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 4371

Last change on this file since 4371 was 4371, checked in by rfurner, 7 years ago

Fix to unitialised value in fld read, see ticket 1222

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