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

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

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

Download in other formats: