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/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 7412

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

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

Download in other formats: