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

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

source: NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/SBC/fldread.F90 @ 12750

Last change on this file since 12750 was 12750, checked in by jcastill, 2 years ago

Comment some excessive warnings

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