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

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

fldread.F90 in branches/2013/dev_LOCEAN_CMCC_INGV_2013/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2013/dev_LOCEAN_CMCC_INGV_2013/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 4230

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

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

Download in other formats: