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 @ 5620

Last change on this file since 5620 was 5620, checked in by jamesharle, 9 years ago
Merge with r5619 of trunk, update to unstructured BDY interpolation in fldread.F90. Structured BDY interpolation incomplete.
Property svn:keywords set to `Id`
File size: 99.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) ; CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1 ), nrec )
762	CASE DEFAULT ;
763
764	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
765	! Do we include something here to adjust barotropic velocities !
766	! in case of a depth difference between bdy files and !
767	! bathymetry in the case ln_full_vel = .false. and jpk_bdy>0? !
768	! [as the enveloping and parital cells could change H !
769	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
770
771	IF( PRESENT(jpk_bdy) .AND. jpk_bdy>0 ) THEN ! boundary data not on model grid: vertical interpolation
772	CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:jpk_bdy), nrec )
773	SELECT CASE( igrd )
774	CASE(1)
775	CALL iom_get ( num, jpdom_unknown, 'gdept', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
776	CALL iom_get ( num, jpdom_unknown, 'e3t', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
777	CASE(2)
778	CALL iom_get ( num, jpdom_unknown, 'gdepu', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
779	CALL iom_get ( num, jpdom_unknown, 'e3u', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
780	CASE(3)
781	CALL iom_get ( num, jpdom_unknown, 'gdepv', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
782	CALL iom_get ( num, jpdom_unknown, 'e3v', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
783	END SELECT
784	CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar )
785
786	#if defined key_bdy
787	CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl)
788	#endif
789	ELSE ! boundary data assumed to be on model grid
790	CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec )
791	IF (ipj==1) THEN ! we assume that this is an un-structured open boundary file
792	DO ib = 1, ipi
793	DO ik = 1, ipk
794	dta(ib,1,ik) = dta_read(map%ptr(ib),1,ik)
795	END DO
796	END DO
797	ELSE ! we assume that this is a structured open boundary file
798	DO ib = 1, ipi
799	jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
800	ji=map%ptr(ib)-(jj-1)*ilendta
801	DO ik = 1, ipk
802	dta(ib,1,ik) = dta_read(ji,jj,ik)
803	END DO
804	END DO
805	ENDIF
806	ENDIF ! PRESENT(jpk_bdy)
807	END SELECT
808
809	END SUBROUTINE fld_map
810
811	#if defined key_bdy
812	SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl)
813
814	!!---------------------------------------------------------------------
815	!! * ROUTINE fld_bdy_interp *
816	!!
817	!! ** Purpose : on the fly vertical interpolation to allow the use of
818	!! boundary data from non-native vertical grid
819	!!----------------------------------------------------------------------
820	USE bdy_oce, ONLY: idx_bdy ! indexing for map <-> ij transformation
821
822	REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read ! work space for global data
823	REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read_z ! work space for global data
824	REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read_dz ! work space for global data
825	REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional)
826	TYPE(MAP_POINTER) , INTENT(in ) :: map ! global-to-local mapping indices
827	LOGICAL , INTENT(in), OPTIONAL :: fvl ! grid type, set number and number of vertical levels in the bdy data
828	INTEGER , INTENT(in) :: igrd, ibdy, jpk_bdy ! number of levels in bdy data
829	INTEGER :: jpkm1_bdy ! number of levels in bdy data minus 1
830	REAL(wp) , INTENT(in) :: fv ! fillvalue and alternative -ABS(fv)
831	!!
832	INTEGER :: ipi ! length of boundary data on local process
833	INTEGER :: ipj ! length of dummy dimension ( = 1 )
834	INTEGER :: ipk ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk )
835	INTEGER :: ilendta ! length of data in file
836	INTEGER :: ib, ik, ikk! loop counters
837	INTEGER :: ji, jj ! loop counters
838	REAL(wp) :: zl, zi, zh, zz, zdz ! tmp variable for current depth and interpolation factor
839	REAL(wp) :: fv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(fv)
840	!!---------------------------------------------------------------------
841
842	ipi = SIZE( dta, 1 )
843	ipj = SIZE( dta_read, 2 )
844	ipk = SIZE( dta, 3 )
845	ilendta = SIZE( dta_read, 1 )
846	jpkm1_bdy = jpk_bdy-1
847
848	fv_alt = -ABS(fv) ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later
849	!
850	IF ( map%ll_unstruc ) THEN ! unstructured open boundary data file
851	DO ib = 1, ipi
852	DO ik = 1, jpk_bdy
853	IF( ( dta_read(map%ptr(ib),1,ik) == fv ) ) THEN
854	dta_read_z(map%ptr(ib),1,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data
855	dta_read_dz(map%ptr(ib),1,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct
856	ENDIF
857	ENDDO
858	ENDDO
859
860	DO ib = 1, ipi
861	DO ik = 1, ipk
862	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?
863	IF( zl < dta_read_z(map%ptr(ib),1,1) ) THEN ! above the first level of external data
864	dta(ib,1,ik) = dta_read(map%ptr(ib),1,1)
865	ELSEIF( zl > MAXVAL(dta_read_z(map%ptr(ib),1,:),1) ) THEN ! below the last level of external data
866	dta(ib,1,ik) = dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1))
867	ELSE ! inbetween : vertical interpolation between ikk & ikk+1
868	DO ikk = 1, jpkm1_bdy ! when gdept_0(ikk) < zl < gdept_0(ikk+1)
869	IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) &
870	& .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN
871	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))
872	dta(ib,1,ik) = dta_read(map%ptr(ib),1,ikk) + &
873	& ( dta_read(map%ptr(ib),1,ikk+1) - dta_read(map%ptr(ib),1,ikk) ) * zi
874	ENDIF
875	END DO
876	ENDIF
877	END DO
878	END DO
879
880	IF(igrd == 2) THEN ! do we need to adjust the transport term?
881	DO ib = 1, ipi
882	zh = SUM(dta_read_dz(map%ptr(ib),1,:) )
883	ztrans = 0._wp
884	ztrans_new = 0._wp
885	DO ik = 1, jpk_bdy
886	ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik)
887	ENDDO
888	DO ik = 1, ipk
889	zdz = e3u_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)
890	ztrans_new = ztrans_new + dta(ib,1,ik) * zdz
891	ENDDO
892	DO ik = 1, ipk
893	zdz = e3u_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)
894	zz = hur(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd))
895	IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
896	dta(ib,1,ik) = dta(ib,1,ik) + ( ztrans - ztrans_new ) * ( zdz * zz )
897	ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
898	dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * ( zdz * zz )
899	ENDIF
900	ENDDO
901	ENDDO
902	ENDIF
903
904	IF(igrd == 3) THEN ! do we need to adjust the transport term?
905	DO ib = 1, ipi
906	zh = SUM(dta_read_dz(map%ptr(ib),1,:) )
907	ztrans = 0._wp
908	ztrans_new = 0._wp
909	DO ik = 1, jpk_bdy
910	ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik)
911	ENDDO
912	DO ik = 1, ipk
913	zdz = e3v_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)
914	ztrans_new = ztrans_new + dta(ib,1,ik) * zdz
915	ENDDO
916	DO ik = 1, ipk
917	zdz = e3v_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)
918	zz = hvr(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd))
919	IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
920	dta(ib,1,ik) = dta(ib,1,ik) + ( ztrans - ztrans_new ) * ( zdz * zz )
921	ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
922	dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * ( zdz * zz )
923	ENDIF
924	ENDDO
925	ENDDO
926	ENDIF
927
928	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
929	! At this point write out a single velocity profile/dz/H for model and input data !
930	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
931
932	ELSE ! structured open boundary file
933	DO ib = 1, ipi
934	jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
935	ji=map%ptr(ib)-(jj-1)*ilendta
936	DO ik = 1, jpk_bdy
937	IF( ( dta_read(ji,jj,ik) == fv ) ) THEN
938	dta_read_z(ji,jj,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data
939	dta_read_dz(ji,jj,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct
940	ENDIF
941	! dta(ib,1,ik) = fv_alt ! put fillvalue into new field as if all goes well all wet points will be replaced
942	ENDDO
943	ENDDO
944	!
945	DO ib = 1, ipi
946	jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
947	ji=map%ptr(ib)-(jj-1)*ilendta
948	DO ik = 1, ipk
949	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?
950	IF( zl < dta_read_z(ji,jj,1) ) THEN ! above the first level of external data
951	dta(ib,1,ik) = dta_read(ji,jj,1)
952	ELSEIF( zl > MAXVAL(dta_read_z(ji,ji,:),1) ) THEN ! below the last level of external data
953	dta(ib,1,ik) = dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1))
954	ELSE ! inbetween : vertical interpolation between ikk & ikk+1
955	DO ikk = 1, jpkm1_bdy ! when gdept_0(ikk) < zl < gdept_0(ikk+1)
956	IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) &
957	& .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN
958	zi = ( zl - dta_read_z(ji,jj,ikk) ) / (dta_read_z(ji,jj,ikk+1)-dta_read_z(ji,jj,ikk))
959	dta(ib,1,ik) = dta_read(ji,jj,ikk) + &
960	& ( dta_read(ji,jj,ikk+1) - dta_read(ji,jj,ikk) ) * zi
961	ENDIF
962	END DO
963	ENDIF
964	END DO
965	END DO
966	ENDIF
967
968	END SUBROUTINE fld_bdy_interp
969
970	! SUBROUTINE fld_bdy_conserve(dta_read, dta_read_z, map, jpk_bdy, igrd, ibdy, fv, dta)
971
972	! END SUBROUTINE fld_bdy_conserve
973
974	#endif
975
976	SUBROUTINE fld_rot( kt, sd )
977	!!---------------------------------------------------------------------
978	!! * ROUTINE fld_rot *
979	!!
980	!! ** Purpose : Vector fields may need to be rotated onto the local grid direction
981	!!----------------------------------------------------------------------
982	INTEGER , INTENT(in ) :: kt ! ocean time step
983	TYPE(FLD), INTENT(inout), DIMENSION(:) :: sd ! input field related variables
984	!!
985	INTEGER :: ju,jv,jk,jn ! loop indices
986	INTEGER :: imf ! size of the structure sd
987	INTEGER :: ill ! character length
988	INTEGER :: iv ! indice of V component
989	REAL(wp), POINTER, DIMENSION(:,:) :: utmp, vtmp ! temporary arrays for vector rotation
990	CHARACTER (LEN=100) :: clcomp ! dummy weight name
991	!!---------------------------------------------------------------------
992
993	CALL wrk_alloc( jpi,jpj, utmp, vtmp )
994
995	!! (sga: following code should be modified so that pairs arent searched for each time
996	!
997	imf = SIZE( sd )
998	DO ju = 1, imf
999	ill = LEN_TRIM( sd(ju)%vcomp )
1000	DO jn = 2-COUNT((/sd(ju)%ln_tint/)), 2
1001	IF( ill > 0 .AND. .NOT. sd(ju)%rotn(jn) ) THEN ! find vector rotations required
1002	IF( sd(ju)%vcomp(1:1) == 'U' ) THEN ! east-west component has symbolic name starting with 'U'
1003	! look for the north-south component which has same symbolic name but with 'U' replaced with 'V'
1004	clcomp = 'V' // sd(ju)%vcomp(2:ill) ! works even if ill == 1
1005	iv = -1
1006	DO jv = 1, imf
1007	IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) ) iv = jv
1008	END DO
1009	IF( iv > 0 ) THEN ! fields ju and iv are two components which need to be rotated together
1010	DO jk = 1, SIZE( sd(ju)%fnow, 3 )
1011	IF( sd(ju)%ln_tint )THEN
1012	CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->i', utmp(:,:) )
1013	CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->j', vtmp(:,:) )
1014	sd(ju)%fdta(:,:,jk,jn) = utmp(:,:) ; sd(iv)%fdta(:,:,jk,jn) = vtmp(:,:)
1015	ELSE
1016	CALL rot_rep( sd(ju)%fnow(:,:,jk ), sd(iv)%fnow(:,:,jk ), 'T', 'en->i', utmp(:,:) )
1017	CALL rot_rep( sd(ju)%fnow(:,:,jk ), sd(iv)%fnow(:,:,jk ), 'T', 'en->j', vtmp(:,:) )
1018	sd(ju)%fnow(:,:,jk ) = utmp(:,:) ; sd(iv)%fnow(:,:,jk ) = vtmp(:,:)
1019	ENDIF
1020	END DO
1021	sd(ju)%rotn(jn) = .TRUE. ! vector was rotated
1022	IF( lwp .AND. kt == nit000 ) WRITE(numout,*) &
1023	& 'fld_read: vector pair ('//TRIM(sd(ju)%clvar)//', '//TRIM(sd(iv)%clvar)//') rotated on to model grid'
1024	ENDIF
1025	ENDIF
1026	ENDIF
1027	END DO
1028	END DO
1029	!
1030	CALL wrk_dealloc( jpi,jpj, utmp, vtmp )
1031	!
1032	END SUBROUTINE fld_rot
1033
1034
1035	SUBROUTINE fld_clopn( sdjf, kyear, kmonth, kday, ldstop )
1036	!!---------------------------------------------------------------------
1037	!! * ROUTINE fld_clopn *
1038	!!
1039	!! ** Purpose : update the file name and open the file
1040	!!----------------------------------------------------------------------
1041	TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables
1042	INTEGER, OPTIONAL, INTENT(in ) :: kyear ! year value
1043	INTEGER, OPTIONAL, INTENT(in ) :: kmonth ! month value
1044	INTEGER, OPTIONAL, INTENT(in ) :: kday ! day value
1045	LOGICAL, OPTIONAL, INTENT(in ) :: ldstop ! stop if open to read a non-existing file (default = .TRUE.)
1046	!!
1047	LOGICAL :: llprevyr ! are we reading previous year file?
1048	LOGICAL :: llprevmth ! are we reading previous month file?
1049	INTEGER :: iyear, imonth, iday ! first day of the current file in yyyy mm dd
1050	INTEGER :: isec_week ! number of seconds since start of the weekly file
1051	INTEGER :: indexyr ! year undex (O/1/2: previous/current/next)
1052	INTEGER :: iyear_len, imonth_len ! length (days) of iyear and imonth !
1053	CHARACTER(len = 256):: clname ! temporary file name
1054	!!----------------------------------------------------------------------
1055	IF( PRESENT(kyear) ) THEN ! use given values
1056	iyear = kyear
1057	imonth = kmonth
1058	iday = kday
1059	ELSE ! use current day values
1060	IF ( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of the week
1061	isec_week = ksec_week( sdjf%cltype(6:8) ) ! second since the beginning of the week
1062	llprevmth = isec_week > nsec_month ! longer time since beginning of the week than the month
1063	llprevyr = llprevmth .AND. nmonth == 1
1064	ELSE
1065	isec_week = 0
1066	llprevmth = .FALSE.
1067	llprevyr = .FALSE.
1068	ENDIF
1069	iyear = nyear - COUNT((/llprevyr /))
1070	imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
1071	iday = nday + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
1072	ENDIF
1073
1074	! build the new filename if not climatological data
1075	clname=TRIM(sdjf%clrootname)
1076	!
1077	! note that sdjf%ln_clim is is only acting on the presence of the year in the file name
1078	IF( .NOT. sdjf%ln_clim ) THEN
1079	WRITE(clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), iyear ! add year
1080	IF( sdjf%cltype /= 'yearly' ) WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname ), imonth ! add month
1081	ELSE
1082	! build the new filename if climatological data
1083	IF( sdjf%cltype /= 'yearly' ) WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), imonth ! add month
1084	ENDIF
1085	IF( sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) &
1086	& WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname ), iday ! add day
1087	!
1088	IF( TRIM(clname) /= TRIM(sdjf%clname) .OR. sdjf%num == 0 ) THEN ! new file to be open
1089
1090	sdjf%clname = TRIM(clname)
1091	IF( sdjf%num /= 0 ) CALL iom_close( sdjf%num ) ! close file if already open
1092	CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 )
1093
1094	! find the last record to be read -> update sdjf%nreclast
1095	indexyr = iyear - nyear + 1
1096	iyear_len = nyear_len( indexyr )
1097	SELECT CASE ( indexyr )
1098	CASE ( 0 ) ; imonth_len = 31 ! previous year -> imonth = 12
1099	CASE ( 1 ) ; imonth_len = nmonth_len(imonth)
1100	CASE ( 2 ) ; imonth_len = 31 ! next year -> imonth = 1
1101	END SELECT
1102
1103	! last record to be read in the current file
1104	IF ( sdjf%nfreqh == -12 ) THEN ; sdjf%nreclast = 1 ! yearly mean
1105	ELSEIF( sdjf%nfreqh == -1 ) THEN ! monthly mean
1106	IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nreclast = 1
1107	ELSE ; sdjf%nreclast = 12
1108	ENDIF
1109	ELSE ! higher frequency mean (in hours)
1110	IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nreclast = NINT( 24 * imonth_len / sdjf%nfreqh )
1111	ELSEIF( sdjf%cltype(1:4) == 'week' ) THEN ; sdjf%nreclast = NINT( 24 * 7 / sdjf%nfreqh )
1112	ELSEIF( sdjf%cltype == 'daily' ) THEN ; sdjf%nreclast = NINT( 24 / sdjf%nfreqh )
1113	ELSE ; sdjf%nreclast = NINT( 24 * iyear_len / sdjf%nfreqh )
1114	ENDIF
1115	ENDIF
1116
1117	ENDIF
1118	!
1119	END SUBROUTINE fld_clopn
1120
1121
1122	SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )
1123	!!---------------------------------------------------------------------
1124	!! * ROUTINE fld_fill *
1125	!!
1126	!! ** Purpose : fill sdf with sdf_n and control print
1127	!!----------------------------------------------------------------------
1128	TYPE(FLD) , DIMENSION(:), INTENT(inout) :: sdf ! structure of input fields (file informations, fields read)
1129	TYPE(FLD_N), DIMENSION(:), INTENT(in ) :: sdf_n ! array of namelist information structures
1130	CHARACTER(len=*) , INTENT(in ) :: cdir ! Root directory for location of flx files
1131	CHARACTER(len=*) , INTENT(in ) :: cdcaller !
1132	CHARACTER(len=*) , INTENT(in ) :: cdtitle !
1133	CHARACTER(len=*) , INTENT(in ) :: cdnam !
1134	!
1135	INTEGER :: jf ! dummy indices
1136	!!---------------------------------------------------------------------
1137
1138	DO jf = 1, SIZE(sdf)
1139	sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
1140	sdf(jf)%clname = "not yet defined"
1141	sdf(jf)%nfreqh = sdf_n(jf)%nfreqh
1142	sdf(jf)%clvar = sdf_n(jf)%clvar
1143	sdf(jf)%ln_tint = sdf_n(jf)%ln_tint
1144	sdf(jf)%ln_clim = sdf_n(jf)%ln_clim
1145	sdf(jf)%cltype = sdf_n(jf)%cltype
1146	sdf(jf)%num = -1
1147	sdf(jf)%wgtname = " "
1148	IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
1149	sdf(jf)%lsmname = " "
1150	IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 ) sdf(jf)%lsmname = TRIM( cdir )//TRIM( sdf_n(jf)%lname )
1151	sdf(jf)%vcomp = sdf_n(jf)%vcomp
1152	sdf(jf)%rotn(:) = .TRUE. ! pretend to be rotated -> won't try to rotate data before the first call to fld_get
1153	IF( sdf(jf)%cltype(1:4) == 'week' .AND. nn_leapy == 0 ) &
1154	& CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs nn_leapy = 1')
1155	IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim ) &
1156	& CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.')
1157	sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
1158	sdf(jf)%igrd = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during bdy_dta_init
1159	sdf(jf)%ibdy = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during bdy_dta_init
1160	END DO
1161
1162	IF(lwp) THEN ! control print
1163	WRITE(numout,*)
1164	WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
1165	WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
1166	WRITE(numout,*) ' '//TRIM( cdnam )//' Namelist'
1167	WRITE(numout,*) ' list of files and frequency (>0: in hours ; <0 in months)'
1168	DO jf = 1, SIZE(sdf)
1169	WRITE(numout,*) ' root filename: ' , TRIM( sdf(jf)%clrootname ), &
1170	& ' variable name: ' , TRIM( sdf(jf)%clvar )
1171	WRITE(numout,*) ' frequency: ' , sdf(jf)%nfreqh , &
1172	& ' time interp: ' , sdf(jf)%ln_tint , &
1173	& ' climatology: ' , sdf(jf)%ln_clim , &
1174	& ' weights : ' , TRIM( sdf(jf)%wgtname ), &
1175	& ' pairing : ' , TRIM( sdf(jf)%vcomp ), &
1176	& ' data type: ' , sdf(jf)%cltype , &
1177	& ' land/sea mask:' , TRIM( sdf(jf)%lsmname )
1178	call flush(numout)
1179	END DO
1180	ENDIF
1181
1182	END SUBROUTINE fld_fill
1183
1184
1185	SUBROUTINE wgt_list( sd, kwgt )
1186	!!---------------------------------------------------------------------
1187	!! * ROUTINE wgt_list *
1188	!!
1189	!! ** Purpose : search array of WGTs and find a weights file
1190	!! entry, or return a new one adding it to the end
1191	!! if it is a new entry, the weights data is read in and
1192	!! restructured (fld_weight)
1193	!!----------------------------------------------------------------------
1194	TYPE( FLD ), INTENT(in ) :: sd ! field with name of weights file
1195	INTEGER , INTENT(inout) :: kwgt ! index of weights
1196	!!
1197	INTEGER :: kw, nestid ! local integer
1198	LOGICAL :: found ! local logical
1199	!!----------------------------------------------------------------------
1200	!
1201	!! search down linked list
1202	!! weights filename is either present or we hit the end of the list
1203	found = .FALSE.
1204
1205	!! because agrif nest part of filenames are now added in iom_open
1206	!! to distinguish between weights files on the different grids, need to track
1207	!! nest number explicitly
1208	nestid = 0
1209	#if defined key_agrif
1210	nestid = Agrif_Fixed()
1211	#endif
1212	DO kw = 1, nxt_wgt-1
1213	IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
1214	ref_wgts(kw)%nestid == nestid) THEN
1215	kwgt = kw
1216	found = .TRUE.
1217	EXIT
1218	ENDIF
1219	END DO
1220	IF( .NOT.found ) THEN
1221	kwgt = nxt_wgt
1222	CALL fld_weight( sd )
1223	ENDIF
1224	!
1225	END SUBROUTINE wgt_list
1226
1227
1228	SUBROUTINE wgt_print( )
1229	!!---------------------------------------------------------------------
1230	!! * ROUTINE wgt_print *
1231	!!
1232	!! ** Purpose : print the list of known weights
1233	!!----------------------------------------------------------------------
1234	INTEGER :: kw !
1235	!!----------------------------------------------------------------------
1236	!
1237	DO kw = 1, nxt_wgt-1
1238	WRITE(numout,*) 'weight file: ',TRIM(ref_wgts(kw)%wgtname)
1239	WRITE(numout,*) ' ddims: ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
1240	WRITE(numout,*) ' numwgt: ',ref_wgts(kw)%numwgt
1241	WRITE(numout,*) ' jpiwgt: ',ref_wgts(kw)%jpiwgt
1242	WRITE(numout,*) ' jpjwgt: ',ref_wgts(kw)%jpjwgt
1243	WRITE(numout,*) ' botleft: ',ref_wgts(kw)%botleft
1244	WRITE(numout,*) ' topright: ',ref_wgts(kw)%topright
1245	IF( ref_wgts(kw)%cyclic ) THEN
1246	WRITE(numout,*) ' cyclical'
1247	IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) ' with overlap of ', ref_wgts(kw)%overlap
1248	ELSE
1249	WRITE(numout,*) ' not cyclical'
1250	ENDIF
1251	IF( ASSOCIATED(ref_wgts(kw)%data_wgt) ) WRITE(numout,*) ' allocated'
1252	END DO
1253	!
1254	END SUBROUTINE wgt_print
1255
1256
1257	SUBROUTINE fld_weight( sd )
1258	!!---------------------------------------------------------------------
1259	!! * ROUTINE fld_weight *
1260	!!
1261	!! ** Purpose : create a new WGT structure and fill in data from
1262	!! file, restructuring as required
1263	!!----------------------------------------------------------------------
1264	TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file
1265	!!
1266	INTEGER :: jn ! dummy loop indices
1267	INTEGER :: inum ! temporary logical unit
1268	INTEGER :: id ! temporary variable id
1269	INTEGER :: ipk ! temporary vertical dimension
1270	CHARACTER (len=5) :: aname
1271	INTEGER , DIMENSION(:), ALLOCATABLE :: ddims
1272	INTEGER , POINTER, DIMENSION(:,:) :: data_src
1273	REAL(wp), POINTER, DIMENSION(:,:) :: data_tmp
1274	LOGICAL :: cyclical
1275	INTEGER :: zwrap ! local integer
1276	!!----------------------------------------------------------------------
1277	!
1278	CALL wrk_alloc( jpi,jpj, data_src ) ! integer
1279	CALL wrk_alloc( jpi,jpj, data_tmp )
1280	!
1281	IF( nxt_wgt > tot_wgts ) THEN
1282	CALL ctl_stop("fld_weight: weights array size exceeded, increase tot_wgts")
1283	ENDIF
1284	!
1285	!! new weights file entry, add in extra information
1286	!! a weights file represents a 2D grid of a certain shape, so we assume that the current
1287	!! input data file is representative of all other files to be opened and processed with the
1288	!! current weights file
1289
1290	!! open input data file (non-model grid)
1291	CALL iom_open( sd%clname, inum, ldiof = LEN(TRIM(sd%wgtname)) > 0 )
1292
1293	!! get dimensions
1294	IF ( SIZE(sd%fnow, 3) > 1 ) THEN
1295	ALLOCATE( ddims(4) )
1296	ELSE
1297	ALLOCATE( ddims(3) )
1298	ENDIF
1299	id = iom_varid( inum, sd%clvar, ddims )
1300
1301	!! close it
1302	CALL iom_close( inum )
1303
1304	!! now open the weights file
1305
1306	CALL iom_open ( sd%wgtname, inum ) ! interpolation weights
1307	IF ( inum > 0 ) THEN
1308
1309	!! determine whether we have an east-west cyclic grid
1310	!! from global attribute called "ew_wrap" in the weights file
1311	!! note that if not found, iom_getatt returns -999 and cyclic with no overlap is assumed
1312	!! since this is the most common forcing configuration
1313
1314	CALL iom_getatt(inum, 'ew_wrap', zwrap)
1315	IF( zwrap >= 0 ) THEN
1316	cyclical = .TRUE.
1317	ELSE IF( zwrap == -999 ) THEN
1318	cyclical = .TRUE.
1319	zwrap = 0
1320	ELSE
1321	cyclical = .FALSE.
1322	ENDIF
1323
1324	ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
1325	ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
1326	ref_wgts(nxt_wgt)%wgtname = sd%wgtname
1327	ref_wgts(nxt_wgt)%overlap = zwrap
1328	ref_wgts(nxt_wgt)%cyclic = cyclical
1329	ref_wgts(nxt_wgt)%nestid = 0
1330	#if defined key_agrif
1331	ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
1332	#endif
1333	!! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
1334	!! for each weight wgtNN there is an integer array srcNN which gives the point in
1335	!! the input data grid which is to be multiplied by the weight
1336	!! they are both arrays on the model grid so the result of the multiplication is
1337	!! added into an output array on the model grid as a running sum
1338
1339	!! two possible cases: bilinear (4 weights) or bicubic (16 weights)
1340	id = iom_varid(inum, 'src05', ldstop=.FALSE.)
1341	IF( id <= 0) THEN
1342	ref_wgts(nxt_wgt)%numwgt = 4
1343	ELSE
1344	ref_wgts(nxt_wgt)%numwgt = 16
1345	ENDIF
1346
1347	ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
1348	ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
1349	ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
1350
1351	DO jn = 1,4
1352	aname = ' '
1353	WRITE(aname,'(a3,i2.2)') 'src',jn
1354	data_tmp(:,:) = 0
1355	CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
1356	data_src(:,:) = INT(data_tmp(:,:))
1357	ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
1358	ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
1359	END DO
1360
1361	DO jn = 1, ref_wgts(nxt_wgt)%numwgt
1362	aname = ' '
1363	WRITE(aname,'(a3,i2.2)') 'wgt',jn
1364	ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
1365	CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
1366	END DO
1367	CALL iom_close (inum)
1368
1369	! find min and max indices in grid
1370	ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1371	ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1372	ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1373	ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1374
1375	! and therefore dimensions of the input box
1376	ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
1377	ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
1378
1379	! shift indexing of source grid
1380	ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
1381	ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
1382
1383	! create input grid, give it a halo to allow gradient calculations
1384	! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
1385	! a more robust solution will be given in next release
1386	ipk = SIZE(sd%fnow, 3)
1387	ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
1388	IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
1389
1390	nxt_wgt = nxt_wgt + 1
1391
1392	ELSE
1393	CALL ctl_stop( ' fld_weight : unable to read the file ' )
1394	ENDIF
1395
1396	DEALLOCATE (ddims )
1397
1398	CALL wrk_dealloc( jpi,jpj, data_src ) ! integer
1399	CALL wrk_dealloc( jpi,jpj, data_tmp )
1400	!
1401	END SUBROUTINE fld_weight
1402
1403
1404	SUBROUTINE apply_seaoverland(clmaskfile,zfieldo,jpi1_lsm,jpi2_lsm,jpj1_lsm, &
1405	& jpj2_lsm,itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
1406	!!---------------------------------------------------------------------
1407	!! * ROUTINE apply_seaoverland *
1408	!!
1409	!! ** Purpose : avoid spurious fluxes in coastal or near-coastal areas
1410	!! due to the wrong usage of "land" values from the coarse
1411	!! atmospheric model when spatial interpolation is required
1412	!! D. Delrosso INGV
1413	!!----------------------------------------------------------------------
1414	INTEGER :: inum,jni,jnj,jnz,jc ! temporary indices
1415	INTEGER, INTENT(in) :: itmpi,itmpj,itmpz ! lengths
1416	INTEGER, INTENT(in) :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices
1417	INTEGER, DIMENSION(3), INTENT(in) :: rec1_lsm,recn_lsm ! temporary arrays for start and length
1418	REAL(wp),DIMENSION (:,:,:),INTENT(inout) :: zfieldo ! input/output array for seaoverland application
1419	REAL(wp),DIMENSION (:,:,:),ALLOCATABLE :: zslmec1 ! temporary array for land point detection
1420	REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfieldn ! array of forcing field with undeff for land points
1421	REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfield ! array of forcing field
1422	CHARACTER (len=100), INTENT(in) :: clmaskfile ! land/sea mask file name
1423	!!---------------------------------------------------------------------
1424	ALLOCATE ( zslmec1(itmpi,itmpj,itmpz) )
1425	ALLOCATE ( zfieldn(itmpi,itmpj) )
1426	ALLOCATE ( zfield(itmpi,itmpj) )
1427
1428	! Retrieve the land sea mask data
1429	CALL iom_open( clmaskfile, inum )
1430	SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1431	CASE(1)
1432	CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm)
1433	CASE DEFAULT
1434	CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm)
1435	END SELECT
1436	CALL iom_close( inum )
1437
1438	DO jnz=1,rec1_lsm(3) !! Loop over k dimension
1439
1440	DO jni=1,itmpi !! copy the original field into a tmp array
1441	DO jnj=1,itmpj !! substituting undeff over land points
1442	zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz)
1443	IF ( zslmec1(jni,jnj,jnz) == 1. ) THEN
1444	zfieldn(jni,jnj) = undeff_lsm
1445	ENDIF
1446	END DO
1447	END DO
1448
1449	CALL seaoverland(zfieldn,itmpi,itmpj,zfield)
1450	DO jc=1,nn_lsm
1451	CALL seaoverland(zfield,itmpi,itmpj,zfield)
1452	END DO
1453
1454	! Check for Undeff and substitute original values
1455	IF(ANY(zfield==undeff_lsm)) THEN
1456	DO jni=1,itmpi
1457	DO jnj=1,itmpj
1458	IF (zfield(jni,jnj)==undeff_lsm) THEN
1459	zfield(jni,jnj) = zfieldo(jni,jnj,jnz)
1460	ENDIF
1461	ENDDO
1462	ENDDO
1463	ENDIF
1464
1465	zfieldo(:,:,jnz)=zfield(:,:)
1466
1467	END DO !! End Loop over k dimension
1468
1469	DEALLOCATE ( zslmec1 )
1470	DEALLOCATE ( zfieldn )
1471	DEALLOCATE ( zfield )
1472
1473	END SUBROUTINE apply_seaoverland
1474
1475
1476	SUBROUTINE seaoverland(zfieldn,ileni,ilenj,zfield)
1477	!!---------------------------------------------------------------------
1478	!! * ROUTINE seaoverland *
1479	!!
1480	!! ** Purpose : create shifted matrices for seaoverland application
1481	!! D. Delrosso INGV
1482	!!----------------------------------------------------------------------
1483	INTEGER,INTENT(in) :: ileni,ilenj ! lengths
1484	REAL,DIMENSION (ileni,ilenj),INTENT(in) :: zfieldn ! array of forcing field with undeff for land points
1485	REAL,DIMENSION (ileni,ilenj),INTENT(out) :: zfield ! array of forcing field
1486	REAL,DIMENSION (ileni,ilenj) :: zmat1,zmat2,zmat3,zmat4 ! temporary arrays for seaoverland application
1487	REAL,DIMENSION (ileni,ilenj) :: zmat5,zmat6,zmat7,zmat8 ! temporary arrays for seaoverland application
1488	REAL,DIMENSION (ileni,ilenj) :: zlsm2d ! temporary arrays for seaoverland application
1489	REAL,DIMENSION (ileni,ilenj,8) :: zlsm3d ! temporary arrays for seaoverland application
1490	LOGICAL,DIMENSION (ileni,ilenj,8) :: ll_msknan3d ! logical mask for undeff detection
1491	LOGICAL,DIMENSION (ileni,ilenj) :: ll_msknan2d ! logical mask for undeff detection
1492	!!----------------------------------------------------------------------
1493	zmat8 = eoshift(zfieldn , SHIFT=-1, BOUNDARY = (/zfieldn(:,1)/) ,DIM=2)
1494	zmat1 = eoshift(zmat8 , SHIFT=-1, BOUNDARY = (/zmat8(1,:)/) ,DIM=1)
1495	zmat2 = eoshift(zfieldn , SHIFT=-1, BOUNDARY = (/zfieldn(1,:)/) ,DIM=1)
1496	zmat4 = eoshift(zfieldn , SHIFT= 1, BOUNDARY = (/zfieldn(:,ilenj)/),DIM=2)
1497	zmat3 = eoshift(zmat4 , SHIFT=-1, BOUNDARY = (/zmat4(1,:)/) ,DIM=1)
1498	zmat5 = eoshift(zmat4 , SHIFT= 1, BOUNDARY = (/zmat4(ileni,:)/) ,DIM=1)
1499	zmat6 = eoshift(zfieldn , SHIFT= 1, BOUNDARY = (/zfieldn(ileni,:)/),DIM=1)
1500	zmat7 = eoshift(zmat8 , SHIFT= 1, BOUNDARY = (/zmat8(ileni,:)/) ,DIM=1)
1501
1502	zlsm3d = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /))
1503	ll_msknan3d = .not.(zlsm3d==undeff_lsm)
1504	ll_msknan2d = .not.(zfieldn==undeff_lsm) ! FALSE where is Undeff (land)
1505	zlsm2d = (SUM ( zlsm3d, 3 , ll_msknan3d ) )/(MAX(1,(COUNT( ll_msknan3d , 3 )) ))
1506	WHERE ((COUNT( ll_msknan3d , 3 )) == 0.0_wp) zlsm2d = undeff_lsm
1507	zfield = MERGE (zfieldn,zlsm2d,ll_msknan2d)
1508	END SUBROUTINE seaoverland
1509
1510
1511	SUBROUTINE fld_interp( num, clvar, kw, kk, dta, &
1512	& nrec, lsmfile)
1513	!!---------------------------------------------------------------------
1514	!! * ROUTINE fld_interp *
1515	!!
1516	!! ** Purpose : apply weights to input gridded data to create data
1517	!! on model grid
1518	!!----------------------------------------------------------------------
1519	INTEGER , INTENT(in ) :: num ! stream number
1520	CHARACTER(LEN=*) , INTENT(in ) :: clvar ! variable name
1521	INTEGER , INTENT(in ) :: kw ! weights number
1522	INTEGER , INTENT(in ) :: kk ! vertical dimension of kk
1523	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: dta ! output field on model grid
1524	INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice)
1525	CHARACTER(LEN=*) , INTENT(in ) :: lsmfile ! land sea mask file name
1526	!!
1527	REAL(wp),DIMENSION(:,:,:),ALLOCATABLE :: ztmp_fly_dta,zfieldo ! temporary array of values on input grid
1528	INTEGER, DIMENSION(3) :: rec1,recn ! temporary arrays for start and length
1529	INTEGER, DIMENSION(3) :: rec1_lsm,recn_lsm ! temporary arrays for start and length in case of seaoverland
1530	INTEGER :: ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2 ! temporary indices
1531	INTEGER :: jk, jn, jm, jir, jjr ! loop counters
1532	INTEGER :: ni, nj ! lengths
1533	INTEGER :: jpimin,jpiwid ! temporary indices
1534	INTEGER :: jpimin_lsm,jpiwid_lsm ! temporary indices
1535	INTEGER :: jpjmin,jpjwid ! temporary indices
1536	INTEGER :: jpjmin_lsm,jpjwid_lsm ! temporary indices
1537	INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices
1538	INTEGER :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices
1539	INTEGER :: itmpi,itmpj,itmpz ! lengths
1540
1541	!!----------------------------------------------------------------------
1542	!
1543	!! for weighted interpolation we have weights at four corners of a box surrounding
1544	!! a model grid point, each weight is multiplied by a grid value (bilinear case)
1545	!! or by a grid value and gradients at the corner point (bicubic case)
1546	!! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
1547
1548	!! sub grid from non-model input grid which encloses all grid points in this nemo process
1549	jpimin = ref_wgts(kw)%botleft(1)
1550	jpjmin = ref_wgts(kw)%botleft(2)
1551	jpiwid = ref_wgts(kw)%jpiwgt
1552	jpjwid = ref_wgts(kw)%jpjwgt
1553
1554	!! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients
1555	rec1(1) = MAX( jpimin-1, 1 )
1556	rec1(2) = MAX( jpjmin-1, 1 )
1557	rec1(3) = 1
1558	recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
1559	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1560	recn(3) = kk
1561
1562	!! where we need to put it in the non-nemo grid fly_dta
1563	!! note that jpi1 and jpj1 only differ from 1 when jpimin and jpjmin are 1
1564	!! (ie at the extreme west or south of the whole input grid) and similarly for jpi2 and jpj2
1565	jpi1 = 2 + rec1(1) - jpimin
1566	jpj1 = 2 + rec1(2) - jpjmin
1567	jpi2 = jpi1 + recn(1) - 1
1568	jpj2 = jpj1 + recn(2) - 1
1569
1570
1571	IF( LEN( TRIM(lsmfile) ) > 0 ) THEN
1572	!! indeces for ztmp_fly_dta
1573	! --------------------------
1574	rec1_lsm(1)=MAX(rec1(1)-nn_lsm,1) ! starting index for enlarged external data, x direction
1575	rec1_lsm(2)=MAX(rec1(2)-nn_lsm,1) ! starting index for enlarged external data, y direction
1576	rec1_lsm(3) = 1 ! vertical dimension
1577	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
1578	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
1579	recn_lsm(3) = kk ! number of vertical levels in the input file
1580
1581	! Avoid out of bound
1582	jpimin_lsm = MAX( rec1_lsm(1)+1, 1 )
1583	jpjmin_lsm = MAX( rec1_lsm(2)+1, 1 )
1584	jpiwid_lsm = MIN( recn_lsm(1)-2,ref_wgts(kw)%ddims(1)-rec1(1)+1)
1585	jpjwid_lsm = MIN( recn_lsm(2)-2,ref_wgts(kw)%ddims(2)-rec1(2)+1)
1586
1587	jpi1_lsm = 2+rec1_lsm(1)-jpimin_lsm
1588	jpj1_lsm = 2+rec1_lsm(2)-jpjmin_lsm
1589	jpi2_lsm = jpi1_lsm + recn_lsm(1) - 1
1590	jpj2_lsm = jpj1_lsm + recn_lsm(2) - 1
1591
1592
1593	itmpi=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),1)
1594	itmpj=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),2)
1595	itmpz=kk
1596	ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz))
1597	ztmp_fly_dta(:,:,:) = 0.0
1598	SELECT CASE( SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1599	CASE(1)
1600	CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), &
1601	& nrec, rec1_lsm, recn_lsm)
1602	CASE DEFAULT
1603	CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), &
1604	& nrec, rec1_lsm, recn_lsm)
1605	END SELECT
1606	CALL apply_seaoverland(lsmfile,ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), &
1607	& jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm, &
1608	& itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
1609
1610
1611	! Relative indeces for remapping
1612	ii_lsm1 = (rec1(1)-rec1_lsm(1))+1
1613	ii_lsm2 = (ii_lsm1+recn(1))-1
1614	ij_lsm1 = (rec1(2)-rec1_lsm(2))+1
1615	ij_lsm2 = (ij_lsm1+recn(2))-1
1616
1617	ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1618	ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:) = ztmp_fly_dta(ii_lsm1:ii_lsm2,ij_lsm1:ij_lsm2,:)
1619	DEALLOCATE(ztmp_fly_dta)
1620
1621	ELSE
1622
1623	ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1624	SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) )
1625	CASE(1)
1626	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn)
1627	CASE DEFAULT
1628	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn)
1629	END SELECT
1630	ENDIF
1631
1632
1633	!! first four weights common to both bilinear and bicubic
1634	!! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft
1635	!! note that we have to offset by 1 into fly_dta array because of halo
1636	dta(:,:,:) = 0.0
1637	DO jk = 1,4
1638	DO jn = 1, jpj
1639	DO jm = 1,jpi
1640	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1641	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1642	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,:)
1643	END DO
1644	END DO
1645	END DO
1646
1647	IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
1648
1649	!! fix up halo points that we couldnt read from file
1650	IF( jpi1 == 2 ) THEN
1651	ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:)
1652	ENDIF
1653	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1654	ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:)
1655	ENDIF
1656	IF( jpj1 == 2 ) THEN
1657	ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:)
1658	ENDIF
1659	IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
1660	ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:)
1661	ENDIF
1662
1663	!! if data grid is cyclic we can do better on east-west edges
1664	!! but have to allow for whether first and last columns are coincident
1665	IF( ref_wgts(kw)%cyclic ) THEN
1666	rec1(2) = MAX( jpjmin-1, 1 )
1667	recn(1) = 1
1668	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1669	jpj1 = 2 + rec1(2) - jpjmin
1670	jpj2 = jpj1 + recn(2) - 1
1671	IF( jpi1 == 2 ) THEN
1672	rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap
1673	SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1674	CASE(1)
1675	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1676	CASE DEFAULT
1677	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1678	END SELECT
1679	ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1680	ENDIF
1681	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1682	rec1(1) = 1 + ref_wgts(kw)%overlap
1683	SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1684	CASE(1)
1685	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1686	CASE DEFAULT
1687	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1688	END SELECT
1689	ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1690	ENDIF
1691	ENDIF
1692
1693	! gradient in the i direction
1694	DO jk = 1,4
1695	DO jn = 1, jpj
1696	DO jm = 1,jpi
1697	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1698	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1699	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 * &
1700	(ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:))
1701	END DO
1702	END DO
1703	END DO
1704
1705	! gradient in the j direction
1706	DO jk = 1,4
1707	DO jn = 1, jpj
1708	DO jm = 1,jpi
1709	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1710	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1711	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 * &
1712	(ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:))
1713	END DO
1714	END DO
1715	END DO
1716
1717	! gradient in the ij direction
1718	DO jk = 1,4
1719	DO jn = 1, jpj
1720	DO jm = 1,jpi
1721	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1722	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1723	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
1724	(ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni ,nj+2,:)) - &
1725	(ref_wgts(kw)%fly_dta(ni+2,nj ,:) - ref_wgts(kw)%fly_dta(ni ,nj ,:)))
1726	END DO
1727	END DO
1728	END DO
1729	!
1730	END IF
1731	!
1732	END SUBROUTINE fld_interp
1733
1734
1735	FUNCTION ksec_week( cdday )
1736	!!---------------------------------------------------------------------
1737	!! * FUNCTION kshift_week *
1738	!!
1739	!! ** Purpose :
1740	!!---------------------------------------------------------------------
1741	CHARACTER(len=*), INTENT(in) :: cdday !3 first letters of the first day of the weekly file
1742	!!
1743	INTEGER :: ksec_week ! output variable
1744	INTEGER :: ijul !temp variable
1745	INTEGER :: ishift !temp variable
1746	CHARACTER(len=3),DIMENSION(7) :: cl_week
1747	!!----------------------------------------------------------------------
1748	cl_week = (/"sun","sat","fri","thu","wed","tue","mon"/)
1749	DO ijul = 1, 7
1750	IF( cl_week(ijul) == TRIM(cdday) ) EXIT
1751	END DO
1752	IF( ijul .GT. 7 ) CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): '//TRIM(cdday) )
1753	!
1754	ishift = ijul * NINT(rday)
1755	!
1756	ksec_week = nsec_week + ishift
1757	ksec_week = MOD( ksec_week, 7*NINT(rday) )
1758	!
1759	END FUNCTION ksec_week
1760
1761	!!======================================================================
1762	END MODULE fldread

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

Download in other formats: