Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

fldread.F90 @ 7029

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

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

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

Context Navigation

source: branches/NERC/dev_r6998_ORCHESTRA/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 7029

Download in other formats: