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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

fldread.F90 in branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 4792

Last change on this file since 4792 was 4792, checked in by jamesharle, 10 years ago
Updates to code after first successful test + merge with HEAD of trunk
Property svn:keywords set to `Id`
File size: 94.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: