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

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

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

Download in other formats: