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

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

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

Context Navigation

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

Last change on this file since 5700 was 5700, checked in by jamesharle, 9 years ago
Bug fixes for structured boundary conditions using BDY on-the-fly interpolation. Structured boundary conditions example now matches the unstructured equivalent.
Property svn:keywords set to `Id`
File size: 102.4 KB

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

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

Download in other formats: