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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

fldread.F90 in branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 5891

Last change on this file since 5891 was 5891, checked in by jamesharle, 8 years ago
bugfix in bdy_interp
Property svn:keywords set to `Id`
File size: 108.3 KB

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

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

Download in other formats: