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 trunk/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 3851

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

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

Download in other formats: