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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/fldread.F90 @ 10922

Last change on this file since 10922 was 10922, checked in by acc, 5 years ago
2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert IOM, LDF, OBS and SBC directories and compatibility changes elsewhere that these changes enforce. Changes pass SETTE and compare with original trunk results. Outstanding issues (currently with work-arounds) in DIU/step_diu.F90 and fld_bdy_interp within SBC/fldread.F90; proper soltions pending
Property svn:keywords set to `Id`
File size: 108.8 KB

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

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

Download in other formats: