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

Last change on this file since 6079 was 6079, checked in by jamesharle, 8 years ago
merge to trunk@5936
Property svn:keywords set to `Id`
File size: 108.6 KB

Line
1	MODULE fldread
2	!!======================================================================
3	!! * MODULE fldread *
4	!! Ocean forcing: read input field for surface boundary condition
5	!!=====================================================================
6	!! History : 2.0 ! 06-2006 (S. Masson, G. Madec) Original code
7	!! ! 05-2008 (S. Alderson) Modified for Interpolation in memory
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	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	# include "domzgr_substitute.h90"
110
111	!$AGRIF_END_DO_NOT_TREAT
112
113	!!----------------------------------------------------------------------
114	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
115	!! $Id$
116	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
117	!!----------------------------------------------------------------------
118	CONTAINS
119
120	SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset, jpk_bdy, fvl )
121	!!---------------------------------------------------------------------
122	!! * ROUTINE fld_read *
123	!!
124	!! ** Purpose : provide at each time step the surface ocean fluxes
125	!! (momentum, heat, freshwater and runoff)
126	!!
127	!! ** Method : READ each input fields in NetCDF files using IOM
128	!! and intepolate it to the model time-step.
129	!! Several assumptions are made on the input file:
130	!! blahblahblah....
131	!!----------------------------------------------------------------------
132	INTEGER , INTENT(in ) :: kt ! ocean time step
133	INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step)
134	TYPE(FLD), INTENT(inout), DIMENSION(:) :: sd ! input field related variables
135	TYPE(MAP_POINTER),INTENT(in), OPTIONAL, DIMENSION(:) :: map ! global-to-local mapping indices
136	INTEGER , INTENT(in ), OPTIONAL :: kit ! subcycle timestep for timesplitting option
137	INTEGER , INTENT(in ), OPTIONAL :: kt_offset ! provide fields at time other than "now"
138	! kt_offset = -1 => fields at "before" time level
139	! kt_offset = +1 => fields at "after" time level
140	! etc.
141	!!
142	INTEGER , INTENT(in ), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data
143	LOGICAL , INTENT(in ), OPTIONAL :: fvl ! number of vertical levels in the BDY data
144	!!
145	INTEGER :: itmp ! temporary variable
146	INTEGER :: imf ! size of the structure sd
147	INTEGER :: jf ! dummy indices
148	INTEGER :: isecend ! number of second since Jan. 1st 00h of nit000 year at nitend
149	INTEGER :: isecsbc ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step
150	INTEGER :: it_offset ! local time offset variable
151	LOGICAL :: llnxtyr ! open next year file?
152	LOGICAL :: llnxtmth ! open next month file?
153	LOGICAL :: llstop ! stop is the file does not exist
154	LOGICAL :: ll_firstcall ! true if this is the first call to fld_read for this set of fields
155	REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation
156	REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation
157	CHARACTER(LEN=1000) :: clfmt ! write format
158	TYPE(MAP_POINTER) :: imap ! global-to-local mapping indices
159	!!---------------------------------------------------------------------
160	ll_firstcall = kt == nit000
161	IF( PRESENT(kit) ) ll_firstcall = ll_firstcall .and. kit == 1
162
163	IF ( nn_components == jp_iam_sas ) THEN ; it_offset = nn_fsbc
164	ELSE ; it_offset = 0
165	ENDIF
166	IF( PRESENT(kt_offset) ) it_offset = kt_offset
167
168	imap%ptr => NULL()
169
170	! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
171	IF( present(kit) ) THEN ! ignore kn_fsbc in this case
172	isecsbc = nsec_year + nsec1jan000 + (kit+it_offset)*NINT( rdt/REAL(nn_baro,wp) )
173	ELSE ! middle of sbc time step
174	isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdttra(1)) + it_offset * NINT(rdttra(1))
175	ENDIF
176	imf = SIZE( sd )
177	!
178	IF( ll_firstcall ) THEN ! initialization
179	DO jf = 1, imf
180	IF( PRESENT(map) ) imap = map(jf)
181	IF( PRESENT(jpk_bdy) ) THEN
182	CALL fld_init( kn_fsbc, sd(jf), imap, jpk_bdy, fvl ) ! read each before field (put them in after as they will be swapped)
183	ELSE
184	CALL fld_init( kn_fsbc, sd(jf), imap ) ! read each before field (put them in after as they will be swapped)
185	ENDIF
186	END DO
187	IF( lwp ) CALL wgt_print() ! control print
188	ENDIF
189	! ! ====================================== !
190	IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN ! update field at each kn_fsbc time-step !
191	! ! ====================================== !
192	!
193	DO jf = 1, imf ! --- loop over field --- !
194
195	IF( isecsbc > sd(jf)%nrec_a(2) .OR. ll_firstcall ) THEN ! read/update the after data?
196
197	IF( PRESENT(map) ) imap = map(jf) ! temporary definition of map
198
199	sd(jf)%nrec_b(:) = sd(jf)%nrec_a(:) ! swap before record informations
200	sd(jf)%rotn(1) = sd(jf)%rotn(2) ! swap before rotate informations
201	IF( sd(jf)%ln_tint ) sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2) ! swap before record field
202
203	CALL fld_rec( kn_fsbc, sd(jf), kt_offset = it_offset, kit = kit ) ! update after record informations
204
205	! if kn_fsbc*rdttra is larger than nfreqh (which is kind of odd),
206	! it is possible that the before value is no more the good one... we have to re-read it
207	! if before is not the last record of the file currently opened and after is the first record to be read
208	! in a new file which means after = 1 (the file to be opened corresponds to the current time)
209	! or after = nreclast + 1 (the file to be opened corresponds to a future time step)
210	IF( .NOT. ll_firstcall .AND. sd(jf)%ln_tint .AND. sd(jf)%nrec_b(1) /= sd(jf)%nreclast &
211	& .AND. MOD( sd(jf)%nrec_a(1), sd(jf)%nreclast ) == 1 ) THEN
212	itmp = sd(jf)%nrec_a(1) ! temporary storage
213	sd(jf)%nrec_a(1) = sd(jf)%nreclast ! read the last record of the file currently opened
214	CALL fld_get( sd(jf), imap ) ! read after data
215	sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2) ! re-swap before record field
216	sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1) ! update before record informations
217	sd(jf)%nrec_b(2) = sd(jf)%nrec_a(2) - NINT( sd(jf)%nfreqh * 3600 ) ! assume freq to be in hours in this case
218	sd(jf)%rotn(1) = sd(jf)%rotn(2) ! update before rotate informations
219	sd(jf)%nrec_a(1) = itmp ! move back to after record
220	ENDIF
221
222	CALL fld_clopn( sd(jf) ) ! Do we need to open a new year/month/week/day file?
223
224	IF( sd(jf)%ln_tint ) THEN
225
226	! if kn_fsbc*rdttra is larger than nfreqh (which is kind of odd),
227	! it is possible that the before value is no more the good one... we have to re-read it
228	! if before record is not just just before the after record...
229	IF( .NOT. ll_firstcall .AND. MOD( sd(jf)%nrec_a(1), sd(jf)%nreclast ) /= 1 &
230	& .AND. sd(jf)%nrec_b(1) /= sd(jf)%nrec_a(1) - 1 ) THEN
231	sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) - 1 ! move back to before record
232	CALL fld_get( sd(jf), imap ) ! read after data
233	sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2) ! re-swap before record field
234	sd(jf)%nrec_b(1) = sd(jf)%nrec_a(1) ! update before record informations
235	sd(jf)%nrec_b(2) = sd(jf)%nrec_a(2) - NINT( sd(jf)%nfreqh * 3600 ) ! assume freq to be in hours in this case
236	sd(jf)%rotn(1) = sd(jf)%rotn(2) ! update before rotate informations
237	sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) + 1 ! move back to after record
238	ENDIF
239
240	! do we have to change the year/month/week/day of the forcing field??
241	! if we do time interpolation we will need to open next year/month/week/day file before the end of the current
242	! 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)
243	! will be larger than the record number that should be read for current year/month/week/day
244	! do we need next file data?
245	IF( sd(jf)%nrec_a(1) > sd(jf)%nreclast ) THEN
246
247	sd(jf)%nrec_a(1) = sd(jf)%nrec_a(1) - sd(jf)%nreclast !
248
249	IF( .NOT. ( sd(jf)%ln_clim .AND. sd(jf)%cltype == 'yearly' ) ) THEN ! close/open the current/new file
250
251	llnxtmth = sd(jf)%cltype == 'monthly' .OR. nday == nmonth_len(nmonth) ! open next month file?
252	llnxtyr = sd(jf)%cltype == 'yearly' .OR. (nmonth == 12 .AND. llnxtmth) ! open next year file?
253
254	! if the run finishes at the end of the current year/month/week/day, we will allow next
255	! year/month/week/day file to be not present. If the run continue further than the current
256	! year/month/week/day, next year/month/week/day file must exist
257	isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdttra(1)) ! second at the end of the run
258	llstop = isecend > sd(jf)%nrec_a(2) ! read more than 1 record of next year
259	! we suppose that the date of next file is next day (should be ok even for weekly files...)
260	CALL fld_clopn( sd(jf), nyear + COUNT((/llnxtyr /)) , &
261	& nmonth + COUNT((/llnxtmth/)) - 12 * COUNT((/llnxtyr /)), &
262	& nday + 1 - nmonth_len(nmonth) * COUNT((/llnxtmth/)), llstop )
263
264	IF( sd(jf)%num <= 0 .AND. .NOT. llstop ) THEN ! next year file does not exist
265	CALL ctl_warn('next year/month/week/day file: '//TRIM(sd(jf)%clname)// &
266	& ' not present -> back to current year/month/day')
267	CALL fld_clopn( sd(jf) ) ! back to the current year/month/day
268	sd(jf)%nrec_a(1) = sd(jf)%nreclast ! force to read the last record in the current year file
269	ENDIF
270
271	ENDIF
272	ENDIF ! open need next file?
273
274	ENDIF ! temporal interpolation?
275
276	! read after data
277	IF( PRESENT(jpk_bdy) ) THEN
278	CALL fld_get( sd(jf), imap, jpk_bdy, fvl)
279	ELSE
280	CALL fld_get( sd(jf), imap )
281	ENDIF
282	ENDIF ! read new data?
283	END DO ! --- end loop over field --- !
284
285	CALL fld_rot( kt, sd ) ! rotate vector before/now/after fields if needed
286
287	DO jf = 1, imf ! --- loop over field --- !
288	!
289	IF( sd(jf)%ln_tint ) THEN ! temporal interpolation
290	IF(lwp .AND. kt - nit000 <= 100 ) THEN
291	clfmt = "('fld_read: var ', a, ' kt = ', i8, ' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," // &
292	& "', records b/a: ', i6.4, '/', i6.4, ' (days ', f9.4,'/', f9.4, ')')"
293	WRITE(numout, clfmt) TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday, &
294	& 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
295	WRITE(numout, *) 'it_offset is : ',it_offset
296	ENDIF
297	! temporal interpolation weights
298	ztinta = REAL( isecsbc - sd(jf)%nrec_b(2), wp ) / REAL( sd(jf)%nrec_a(2) - sd(jf)%nrec_b(2), wp )
299	ztintb = 1. - ztinta
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	IF ( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of the week
1217	isec_week = ksec_week( sdjf%cltype(6:8) )- (86400 * 8 )
1218	llprevmth = isec_week > nsec_month ! longer time since beginning of the week than the month
1219	llprevyr = llprevmth .AND. nmonth == 1
1220	iyear = nyear - COUNT((/llprevyr /))
1221	imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
1222	iday = nday + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
1223	ENDIF
1224	ELSE ! use current day values
1225	IF ( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of the week
1226	isec_week = ksec_week( sdjf%cltype(6:8) ) ! second since the beginning of the week
1227	llprevmth = isec_week > nsec_month ! longer time since beginning of the week than the month
1228	llprevyr = llprevmth .AND. nmonth == 1
1229	ELSE
1230	isec_week = 0
1231	llprevmth = .FALSE.
1232	llprevyr = .FALSE.
1233	ENDIF
1234	iyear = nyear - COUNT((/llprevyr /))
1235	imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
1236	iday = nday + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
1237	ENDIF
1238
1239	! build the new filename if not climatological data
1240	clname=TRIM(sdjf%clrootname)
1241	!
1242	! note that sdjf%ln_clim is is only acting on the presence of the year in the file name
1243	IF( .NOT. sdjf%ln_clim ) THEN
1244	WRITE(clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), iyear ! add year
1245	IF( sdjf%cltype /= 'yearly' ) WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname ), imonth ! add month
1246	ELSE
1247	! build the new filename if climatological data
1248	IF( sdjf%cltype /= 'yearly' ) WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), imonth ! add month
1249	ENDIF
1250	IF( sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) &
1251	& WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname ), iday ! add day
1252	!
1253	IF( TRIM(clname) /= TRIM(sdjf%clname) .OR. sdjf%num == 0 ) THEN ! new file to be open
1254
1255	sdjf%clname = TRIM(clname)
1256	IF( sdjf%num /= 0 ) CALL iom_close( sdjf%num ) ! close file if already open
1257	CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 )
1258
1259	! find the last record to be read -> update sdjf%nreclast
1260	indexyr = iyear - nyear + 1
1261	iyear_len = nyear_len( indexyr )
1262	SELECT CASE ( indexyr )
1263	CASE ( 0 ) ; imonth_len = 31 ! previous year -> imonth = 12
1264	CASE ( 1 ) ; imonth_len = nmonth_len(imonth)
1265	CASE ( 2 ) ; imonth_len = 31 ! next year -> imonth = 1
1266	END SELECT
1267
1268	! last record to be read in the current file
1269	IF ( sdjf%nfreqh == -12 ) THEN ; sdjf%nreclast = 1 ! yearly mean
1270	ELSEIF( sdjf%nfreqh == -1 ) THEN ! monthly mean
1271	IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nreclast = 1
1272	ELSE ; sdjf%nreclast = 12
1273	ENDIF
1274	ELSE ! higher frequency mean (in hours)
1275	IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nreclast = NINT( 24 * imonth_len / sdjf%nfreqh )
1276	ELSEIF( sdjf%cltype(1:4) == 'week' ) THEN ; sdjf%nreclast = NINT( 24 * 7 / sdjf%nfreqh )
1277	ELSEIF( sdjf%cltype == 'daily' ) THEN ; sdjf%nreclast = NINT( 24 / sdjf%nfreqh )
1278	ELSE ; sdjf%nreclast = NINT( 24 * iyear_len / sdjf%nfreqh )
1279	ENDIF
1280	ENDIF
1281
1282	ENDIF
1283	!
1284	END SUBROUTINE fld_clopn
1285
1286
1287	SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )
1288	!!---------------------------------------------------------------------
1289	!! * ROUTINE fld_fill *
1290	!!
1291	!! ** Purpose : fill sdf with sdf_n and control print
1292	!!----------------------------------------------------------------------
1293	TYPE(FLD) , DIMENSION(:), INTENT(inout) :: sdf ! structure of input fields (file informations, fields read)
1294	TYPE(FLD_N), DIMENSION(:), INTENT(in ) :: sdf_n ! array of namelist information structures
1295	CHARACTER(len=*) , INTENT(in ) :: cdir ! Root directory for location of flx files
1296	CHARACTER(len=*) , INTENT(in ) :: cdcaller !
1297	CHARACTER(len=*) , INTENT(in ) :: cdtitle !
1298	CHARACTER(len=*) , INTENT(in ) :: cdnam !
1299	!
1300	INTEGER :: jf ! dummy indices
1301	!!---------------------------------------------------------------------
1302
1303	DO jf = 1, SIZE(sdf)
1304	sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
1305	sdf(jf)%clname = "not yet defined"
1306	sdf(jf)%nfreqh = sdf_n(jf)%nfreqh
1307	sdf(jf)%clvar = sdf_n(jf)%clvar
1308	sdf(jf)%ln_tint = sdf_n(jf)%ln_tint
1309	sdf(jf)%ln_clim = sdf_n(jf)%ln_clim
1310	sdf(jf)%cltype = sdf_n(jf)%cltype
1311	sdf(jf)%num = -1
1312	sdf(jf)%wgtname = " "
1313	IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
1314	sdf(jf)%lsmname = " "
1315	IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 ) sdf(jf)%lsmname = TRIM( cdir )//TRIM( sdf_n(jf)%lname )
1316	sdf(jf)%vcomp = sdf_n(jf)%vcomp
1317	sdf(jf)%rotn(:) = .TRUE. ! pretend to be rotated -> won't try to rotate data before the first call to fld_get
1318	IF( sdf(jf)%cltype(1:4) == 'week' .AND. nn_leapy == 0 ) &
1319	& CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs nn_leapy = 1')
1320	IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim ) &
1321	& CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.')
1322	sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
1323	END DO
1324
1325	IF(lwp) THEN ! control print
1326	WRITE(numout,*)
1327	WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
1328	WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
1329	WRITE(numout,*) ' '//TRIM( cdnam )//' Namelist'
1330	WRITE(numout,*) ' list of files and frequency (>0: in hours ; <0 in months)'
1331	DO jf = 1, SIZE(sdf)
1332	WRITE(numout,*) ' root filename: ' , TRIM( sdf(jf)%clrootname ), &
1333	& ' variable name: ' , TRIM( sdf(jf)%clvar )
1334	WRITE(numout,*) ' frequency: ' , sdf(jf)%nfreqh , &
1335	& ' time interp: ' , sdf(jf)%ln_tint , &
1336	& ' climatology: ' , sdf(jf)%ln_clim , &
1337	& ' weights : ' , TRIM( sdf(jf)%wgtname ), &
1338	& ' pairing : ' , TRIM( sdf(jf)%vcomp ), &
1339	& ' data type: ' , sdf(jf)%cltype , &
1340	& ' land/sea mask:' , TRIM( sdf(jf)%lsmname )
1341	call flush(numout)
1342	END DO
1343	ENDIF
1344
1345	END SUBROUTINE fld_fill
1346
1347
1348	SUBROUTINE wgt_list( sd, kwgt )
1349	!!---------------------------------------------------------------------
1350	!! * ROUTINE wgt_list *
1351	!!
1352	!! ** Purpose : search array of WGTs and find a weights file
1353	!! entry, or return a new one adding it to the end
1354	!! if it is a new entry, the weights data is read in and
1355	!! restructured (fld_weight)
1356	!!----------------------------------------------------------------------
1357	TYPE( FLD ), INTENT(in ) :: sd ! field with name of weights file
1358	INTEGER , INTENT(inout) :: kwgt ! index of weights
1359	!!
1360	INTEGER :: kw, nestid ! local integer
1361	LOGICAL :: found ! local logical
1362	!!----------------------------------------------------------------------
1363	!
1364	!! search down linked list
1365	!! weights filename is either present or we hit the end of the list
1366	found = .FALSE.
1367
1368	!! because agrif nest part of filenames are now added in iom_open
1369	!! to distinguish between weights files on the different grids, need to track
1370	!! nest number explicitly
1371	nestid = 0
1372	#if defined key_agrif
1373	nestid = Agrif_Fixed()
1374	#endif
1375	DO kw = 1, nxt_wgt-1
1376	IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
1377	ref_wgts(kw)%nestid == nestid) THEN
1378	kwgt = kw
1379	found = .TRUE.
1380	EXIT
1381	ENDIF
1382	END DO
1383	IF( .NOT.found ) THEN
1384	kwgt = nxt_wgt
1385	CALL fld_weight( sd )
1386	ENDIF
1387	!
1388	END SUBROUTINE wgt_list
1389
1390
1391	SUBROUTINE wgt_print( )
1392	!!---------------------------------------------------------------------
1393	!! * ROUTINE wgt_print *
1394	!!
1395	!! ** Purpose : print the list of known weights
1396	!!----------------------------------------------------------------------
1397	INTEGER :: kw !
1398	!!----------------------------------------------------------------------
1399	!
1400	DO kw = 1, nxt_wgt-1
1401	WRITE(numout,*) 'weight file: ',TRIM(ref_wgts(kw)%wgtname)
1402	WRITE(numout,*) ' ddims: ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
1403	WRITE(numout,*) ' numwgt: ',ref_wgts(kw)%numwgt
1404	WRITE(numout,*) ' jpiwgt: ',ref_wgts(kw)%jpiwgt
1405	WRITE(numout,*) ' jpjwgt: ',ref_wgts(kw)%jpjwgt
1406	WRITE(numout,*) ' botleft: ',ref_wgts(kw)%botleft
1407	WRITE(numout,*) ' topright: ',ref_wgts(kw)%topright
1408	IF( ref_wgts(kw)%cyclic ) THEN
1409	WRITE(numout,*) ' cyclical'
1410	IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) ' with overlap of ', ref_wgts(kw)%overlap
1411	ELSE
1412	WRITE(numout,*) ' not cyclical'
1413	ENDIF
1414	IF( ASSOCIATED(ref_wgts(kw)%data_wgt) ) WRITE(numout,*) ' allocated'
1415	END DO
1416	!
1417	END SUBROUTINE wgt_print
1418
1419
1420	SUBROUTINE fld_weight( sd )
1421	!!---------------------------------------------------------------------
1422	!! * ROUTINE fld_weight *
1423	!!
1424	!! ** Purpose : create a new WGT structure and fill in data from
1425	!! file, restructuring as required
1426	!!----------------------------------------------------------------------
1427	TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file
1428	!!
1429	INTEGER :: jn ! dummy loop indices
1430	INTEGER :: inum ! temporary logical unit
1431	INTEGER :: id ! temporary variable id
1432	INTEGER :: ipk ! temporary vertical dimension
1433	CHARACTER (len=5) :: aname
1434	INTEGER , DIMENSION(:), ALLOCATABLE :: ddims
1435	INTEGER , POINTER, DIMENSION(:,:) :: data_src
1436	REAL(wp), POINTER, DIMENSION(:,:) :: data_tmp
1437	LOGICAL :: cyclical
1438	INTEGER :: zwrap ! local integer
1439	!!----------------------------------------------------------------------
1440	!
1441	CALL wrk_alloc( jpi,jpj, data_src ) ! integer
1442	CALL wrk_alloc( jpi,jpj, data_tmp )
1443	!
1444	IF( nxt_wgt > tot_wgts ) THEN
1445	CALL ctl_stop("fld_weight: weights array size exceeded, increase tot_wgts")
1446	ENDIF
1447	!
1448	!! new weights file entry, add in extra information
1449	!! a weights file represents a 2D grid of a certain shape, so we assume that the current
1450	!! input data file is representative of all other files to be opened and processed with the
1451	!! current weights file
1452
1453	!! open input data file (non-model grid)
1454	CALL iom_open( sd%clname, inum, ldiof = LEN(TRIM(sd%wgtname)) > 0 )
1455
1456	!! get dimensions
1457	IF ( SIZE(sd%fnow, 3) > 1 ) THEN
1458	ALLOCATE( ddims(4) )
1459	ELSE
1460	ALLOCATE( ddims(3) )
1461	ENDIF
1462	id = iom_varid( inum, sd%clvar, ddims )
1463
1464	!! close it
1465	CALL iom_close( inum )
1466
1467	!! now open the weights file
1468
1469	CALL iom_open ( sd%wgtname, inum ) ! interpolation weights
1470	IF ( inum > 0 ) THEN
1471
1472	!! determine whether we have an east-west cyclic grid
1473	!! from global attribute called "ew_wrap" in the weights file
1474	!! note that if not found, iom_getatt returns -999 and cyclic with no overlap is assumed
1475	!! since this is the most common forcing configuration
1476
1477	CALL iom_getatt(inum, 'ew_wrap', zwrap)
1478	IF( zwrap >= 0 ) THEN
1479	cyclical = .TRUE.
1480	ELSE IF( zwrap == -999 ) THEN
1481	cyclical = .TRUE.
1482	zwrap = 0
1483	ELSE
1484	cyclical = .FALSE.
1485	ENDIF
1486
1487	ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
1488	ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
1489	ref_wgts(nxt_wgt)%wgtname = sd%wgtname
1490	ref_wgts(nxt_wgt)%overlap = zwrap
1491	ref_wgts(nxt_wgt)%cyclic = cyclical
1492	ref_wgts(nxt_wgt)%nestid = 0
1493	#if defined key_agrif
1494	ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
1495	#endif
1496	!! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
1497	!! for each weight wgtNN there is an integer array srcNN which gives the point in
1498	!! the input data grid which is to be multiplied by the weight
1499	!! they are both arrays on the model grid so the result of the multiplication is
1500	!! added into an output array on the model grid as a running sum
1501
1502	!! two possible cases: bilinear (4 weights) or bicubic (16 weights)
1503	id = iom_varid(inum, 'src05', ldstop=.FALSE.)
1504	IF( id <= 0) THEN
1505	ref_wgts(nxt_wgt)%numwgt = 4
1506	ELSE
1507	ref_wgts(nxt_wgt)%numwgt = 16
1508	ENDIF
1509
1510	ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
1511	ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
1512	ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
1513
1514	DO jn = 1,4
1515	aname = ' '
1516	WRITE(aname,'(a3,i2.2)') 'src',jn
1517	data_tmp(:,:) = 0
1518	CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
1519	data_src(:,:) = INT(data_tmp(:,:))
1520	ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
1521	ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
1522	END DO
1523
1524	DO jn = 1, ref_wgts(nxt_wgt)%numwgt
1525	aname = ' '
1526	WRITE(aname,'(a3,i2.2)') 'wgt',jn
1527	ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
1528	CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
1529	END DO
1530	CALL iom_close (inum)
1531
1532	! find min and max indices in grid
1533	ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1534	ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1535	ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1536	ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1537
1538	! and therefore dimensions of the input box
1539	ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
1540	ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
1541
1542	! shift indexing of source grid
1543	ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
1544	ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
1545
1546	! create input grid, give it a halo to allow gradient calculations
1547	! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
1548	! a more robust solution will be given in next release
1549	ipk = SIZE(sd%fnow, 3)
1550	ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
1551	IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
1552
1553	nxt_wgt = nxt_wgt + 1
1554
1555	ELSE
1556	CALL ctl_stop( ' fld_weight : unable to read the file ' )
1557	ENDIF
1558
1559	DEALLOCATE (ddims )
1560
1561	CALL wrk_dealloc( jpi,jpj, data_src ) ! integer
1562	CALL wrk_dealloc( jpi,jpj, data_tmp )
1563	!
1564	END SUBROUTINE fld_weight
1565
1566
1567	SUBROUTINE apply_seaoverland(clmaskfile,zfieldo,jpi1_lsm,jpi2_lsm,jpj1_lsm, &
1568	& jpj2_lsm,itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
1569	!!---------------------------------------------------------------------
1570	!! * ROUTINE apply_seaoverland *
1571	!!
1572	!! ** Purpose : avoid spurious fluxes in coastal or near-coastal areas
1573	!! due to the wrong usage of "land" values from the coarse
1574	!! atmospheric model when spatial interpolation is required
1575	!! D. Delrosso INGV
1576	!!----------------------------------------------------------------------
1577	INTEGER :: inum,jni,jnj,jnz,jc ! temporary indices
1578	INTEGER, INTENT(in) :: itmpi,itmpj,itmpz ! lengths
1579	INTEGER, INTENT(in) :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices
1580	INTEGER, DIMENSION(3), INTENT(in) :: rec1_lsm,recn_lsm ! temporary arrays for start and length
1581	REAL(wp),DIMENSION (:,:,:),INTENT(inout) :: zfieldo ! input/output array for seaoverland application
1582	REAL(wp),DIMENSION (:,:,:),ALLOCATABLE :: zslmec1 ! temporary array for land point detection
1583	REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfieldn ! array of forcing field with undeff for land points
1584	REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfield ! array of forcing field
1585	CHARACTER (len=100), INTENT(in) :: clmaskfile ! land/sea mask file name
1586	!!---------------------------------------------------------------------
1587	ALLOCATE ( zslmec1(itmpi,itmpj,itmpz) )
1588	ALLOCATE ( zfieldn(itmpi,itmpj) )
1589	ALLOCATE ( zfield(itmpi,itmpj) )
1590
1591	! Retrieve the land sea mask data
1592	CALL iom_open( clmaskfile, inum )
1593	SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1594	CASE(1)
1595	CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm)
1596	CASE DEFAULT
1597	CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm)
1598	END SELECT
1599	CALL iom_close( inum )
1600
1601	DO jnz=1,rec1_lsm(3) !! Loop over k dimension
1602
1603	DO jni=1,itmpi !! copy the original field into a tmp array
1604	DO jnj=1,itmpj !! substituting undeff over land points
1605	zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz)
1606	IF ( zslmec1(jni,jnj,jnz) == 1. ) THEN
1607	zfieldn(jni,jnj) = undeff_lsm
1608	ENDIF
1609	END DO
1610	END DO
1611
1612	CALL seaoverland(zfieldn,itmpi,itmpj,zfield)
1613	DO jc=1,nn_lsm
1614	CALL seaoverland(zfield,itmpi,itmpj,zfield)
1615	END DO
1616
1617	! Check for Undeff and substitute original values
1618	IF(ANY(zfield==undeff_lsm)) THEN
1619	DO jni=1,itmpi
1620	DO jnj=1,itmpj
1621	IF (zfield(jni,jnj)==undeff_lsm) THEN
1622	zfield(jni,jnj) = zfieldo(jni,jnj,jnz)
1623	ENDIF
1624	ENDDO
1625	ENDDO
1626	ENDIF
1627
1628	zfieldo(:,:,jnz)=zfield(:,:)
1629
1630	END DO !! End Loop over k dimension
1631
1632	DEALLOCATE ( zslmec1 )
1633	DEALLOCATE ( zfieldn )
1634	DEALLOCATE ( zfield )
1635
1636	END SUBROUTINE apply_seaoverland
1637
1638
1639	SUBROUTINE seaoverland(zfieldn,ileni,ilenj,zfield)
1640	!!---------------------------------------------------------------------
1641	!! * ROUTINE seaoverland *
1642	!!
1643	!! ** Purpose : create shifted matrices for seaoverland application
1644	!! D. Delrosso INGV
1645	!!----------------------------------------------------------------------
1646	INTEGER,INTENT(in) :: ileni,ilenj ! lengths
1647	REAL,DIMENSION (ileni,ilenj),INTENT(in) :: zfieldn ! array of forcing field with undeff for land points
1648	REAL,DIMENSION (ileni,ilenj),INTENT(out) :: zfield ! array of forcing field
1649	REAL,DIMENSION (ileni,ilenj) :: zmat1,zmat2,zmat3,zmat4 ! temporary arrays for seaoverland application
1650	REAL,DIMENSION (ileni,ilenj) :: zmat5,zmat6,zmat7,zmat8 ! temporary arrays for seaoverland application
1651	REAL,DIMENSION (ileni,ilenj) :: zlsm2d ! temporary arrays for seaoverland application
1652	REAL,DIMENSION (ileni,ilenj,8) :: zlsm3d ! temporary arrays for seaoverland application
1653	LOGICAL,DIMENSION (ileni,ilenj,8) :: ll_msknan3d ! logical mask for undeff detection
1654	LOGICAL,DIMENSION (ileni,ilenj) :: ll_msknan2d ! logical mask for undeff detection
1655	!!----------------------------------------------------------------------
1656	zmat8 = eoshift(zfieldn , SHIFT=-1, BOUNDARY = (/zfieldn(:,1)/) ,DIM=2)
1657	zmat1 = eoshift(zmat8 , SHIFT=-1, BOUNDARY = (/zmat8(1,:)/) ,DIM=1)
1658	zmat2 = eoshift(zfieldn , SHIFT=-1, BOUNDARY = (/zfieldn(1,:)/) ,DIM=1)
1659	zmat4 = eoshift(zfieldn , SHIFT= 1, BOUNDARY = (/zfieldn(:,ilenj)/),DIM=2)
1660	zmat3 = eoshift(zmat4 , SHIFT=-1, BOUNDARY = (/zmat4(1,:)/) ,DIM=1)
1661	zmat5 = eoshift(zmat4 , SHIFT= 1, BOUNDARY = (/zmat4(ileni,:)/) ,DIM=1)
1662	zmat6 = eoshift(zfieldn , SHIFT= 1, BOUNDARY = (/zfieldn(ileni,:)/),DIM=1)
1663	zmat7 = eoshift(zmat8 , SHIFT= 1, BOUNDARY = (/zmat8(ileni,:)/) ,DIM=1)
1664
1665	zlsm3d = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /))
1666	ll_msknan3d = .not.(zlsm3d==undeff_lsm)
1667	ll_msknan2d = .not.(zfieldn==undeff_lsm) ! FALSE where is Undeff (land)
1668	zlsm2d = (SUM ( zlsm3d, 3 , ll_msknan3d ) )/(MAX(1,(COUNT( ll_msknan3d , 3 )) ))
1669	WHERE ((COUNT( ll_msknan3d , 3 )) == 0.0_wp) zlsm2d = undeff_lsm
1670	zfield = MERGE (zfieldn,zlsm2d,ll_msknan2d)
1671	END SUBROUTINE seaoverland
1672
1673
1674	SUBROUTINE fld_interp( num, clvar, kw, kk, dta, &
1675	& nrec, lsmfile)
1676	!!---------------------------------------------------------------------
1677	!! * ROUTINE fld_interp *
1678	!!
1679	!! ** Purpose : apply weights to input gridded data to create data
1680	!! on model grid
1681	!!----------------------------------------------------------------------
1682	INTEGER , INTENT(in ) :: num ! stream number
1683	CHARACTER(LEN=*) , INTENT(in ) :: clvar ! variable name
1684	INTEGER , INTENT(in ) :: kw ! weights number
1685	INTEGER , INTENT(in ) :: kk ! vertical dimension of kk
1686	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: dta ! output field on model grid
1687	INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice)
1688	CHARACTER(LEN=*) , INTENT(in ) :: lsmfile ! land sea mask file name
1689	!!
1690	REAL(wp),DIMENSION(:,:,:),ALLOCATABLE :: ztmp_fly_dta,zfieldo ! temporary array of values on input grid
1691	INTEGER, DIMENSION(3) :: rec1,recn ! temporary arrays for start and length
1692	INTEGER, DIMENSION(3) :: rec1_lsm,recn_lsm ! temporary arrays for start and length in case of seaoverland
1693	INTEGER :: ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2 ! temporary indices
1694	INTEGER :: jk, jn, jm, jir, jjr ! loop counters
1695	INTEGER :: ni, nj ! lengths
1696	INTEGER :: jpimin,jpiwid ! temporary indices
1697	INTEGER :: jpimin_lsm,jpiwid_lsm ! temporary indices
1698	INTEGER :: jpjmin,jpjwid ! temporary indices
1699	INTEGER :: jpjmin_lsm,jpjwid_lsm ! temporary indices
1700	INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices
1701	INTEGER :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices
1702	INTEGER :: itmpi,itmpj,itmpz ! lengths
1703
1704	!!----------------------------------------------------------------------
1705	!
1706	!! for weighted interpolation we have weights at four corners of a box surrounding
1707	!! a model grid point, each weight is multiplied by a grid value (bilinear case)
1708	!! or by a grid value and gradients at the corner point (bicubic case)
1709	!! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
1710
1711	!! sub grid from non-model input grid which encloses all grid points in this nemo process
1712	jpimin = ref_wgts(kw)%botleft(1)
1713	jpjmin = ref_wgts(kw)%botleft(2)
1714	jpiwid = ref_wgts(kw)%jpiwgt
1715	jpjwid = ref_wgts(kw)%jpjwgt
1716
1717	!! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients
1718	rec1(1) = MAX( jpimin-1, 1 )
1719	rec1(2) = MAX( jpjmin-1, 1 )
1720	rec1(3) = 1
1721	recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
1722	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1723	recn(3) = kk
1724
1725	!! where we need to put it in the non-nemo grid fly_dta
1726	!! note that jpi1 and jpj1 only differ from 1 when jpimin and jpjmin are 1
1727	!! (ie at the extreme west or south of the whole input grid) and similarly for jpi2 and jpj2
1728	jpi1 = 2 + rec1(1) - jpimin
1729	jpj1 = 2 + rec1(2) - jpjmin
1730	jpi2 = jpi1 + recn(1) - 1
1731	jpj2 = jpj1 + recn(2) - 1
1732
1733
1734	IF( LEN( TRIM(lsmfile) ) > 0 ) THEN
1735	!! indeces for ztmp_fly_dta
1736	! --------------------------
1737	rec1_lsm(1)=MAX(rec1(1)-nn_lsm,1) ! starting index for enlarged external data, x direction
1738	rec1_lsm(2)=MAX(rec1(2)-nn_lsm,1) ! starting index for enlarged external data, y direction
1739	rec1_lsm(3) = 1 ! vertical dimension
1740	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
1741	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
1742	recn_lsm(3) = kk ! number of vertical levels in the input file
1743
1744	! Avoid out of bound
1745	jpimin_lsm = MAX( rec1_lsm(1)+1, 1 )
1746	jpjmin_lsm = MAX( rec1_lsm(2)+1, 1 )
1747	jpiwid_lsm = MIN( recn_lsm(1)-2,ref_wgts(kw)%ddims(1)-rec1(1)+1)
1748	jpjwid_lsm = MIN( recn_lsm(2)-2,ref_wgts(kw)%ddims(2)-rec1(2)+1)
1749
1750	jpi1_lsm = 2+rec1_lsm(1)-jpimin_lsm
1751	jpj1_lsm = 2+rec1_lsm(2)-jpjmin_lsm
1752	jpi2_lsm = jpi1_lsm + recn_lsm(1) - 1
1753	jpj2_lsm = jpj1_lsm + recn_lsm(2) - 1
1754
1755
1756	itmpi=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),1)
1757	itmpj=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),2)
1758	itmpz=kk
1759	ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz))
1760	ztmp_fly_dta(:,:,:) = 0.0
1761	SELECT CASE( SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1762	CASE(1)
1763	CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), &
1764	& nrec, rec1_lsm, recn_lsm)
1765	CASE DEFAULT
1766	CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), &
1767	& nrec, rec1_lsm, recn_lsm)
1768	END SELECT
1769	CALL apply_seaoverland(lsmfile,ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), &
1770	& jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm, &
1771	& itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
1772
1773
1774	! Relative indeces for remapping
1775	ii_lsm1 = (rec1(1)-rec1_lsm(1))+1
1776	ii_lsm2 = (ii_lsm1+recn(1))-1
1777	ij_lsm1 = (rec1(2)-rec1_lsm(2))+1
1778	ij_lsm2 = (ij_lsm1+recn(2))-1
1779
1780	ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1781	ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:) = ztmp_fly_dta(ii_lsm1:ii_lsm2,ij_lsm1:ij_lsm2,:)
1782	DEALLOCATE(ztmp_fly_dta)
1783
1784	ELSE
1785
1786	ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1787	SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) )
1788	CASE(1)
1789	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn)
1790	CASE DEFAULT
1791	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn)
1792	END SELECT
1793	ENDIF
1794
1795
1796	!! first four weights common to both bilinear and bicubic
1797	!! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft
1798	!! note that we have to offset by 1 into fly_dta array because of halo
1799	dta(:,:,:) = 0.0
1800	DO jk = 1,4
1801	DO jn = 1, jpj
1802	DO jm = 1,jpi
1803	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1804	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1805	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,:)
1806	END DO
1807	END DO
1808	END DO
1809
1810	IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
1811
1812	!! fix up halo points that we couldnt read from file
1813	IF( jpi1 == 2 ) THEN
1814	ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:)
1815	ENDIF
1816	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1817	ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:)
1818	ENDIF
1819	IF( jpj1 == 2 ) THEN
1820	ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:)
1821	ENDIF
1822	IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
1823	ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:)
1824	ENDIF
1825
1826	!! if data grid is cyclic we can do better on east-west edges
1827	!! but have to allow for whether first and last columns are coincident
1828	IF( ref_wgts(kw)%cyclic ) THEN
1829	rec1(2) = MAX( jpjmin-1, 1 )
1830	recn(1) = 1
1831	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1832	jpj1 = 2 + rec1(2) - jpjmin
1833	jpj2 = jpj1 + recn(2) - 1
1834	IF( jpi1 == 2 ) THEN
1835	rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap
1836	SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1837	CASE(1)
1838	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1839	CASE DEFAULT
1840	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1841	END SELECT
1842	ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1843	ENDIF
1844	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1845	rec1(1) = 1 + ref_wgts(kw)%overlap
1846	SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1847	CASE(1)
1848	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1849	CASE DEFAULT
1850	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1851	END SELECT
1852	ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1853	ENDIF
1854	ENDIF
1855
1856	! gradient in the i direction
1857	DO jk = 1,4
1858	DO jn = 1, jpj
1859	DO jm = 1,jpi
1860	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1861	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1862	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 * &
1863	(ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:))
1864	END DO
1865	END DO
1866	END DO
1867
1868	! gradient in the j direction
1869	DO jk = 1,4
1870	DO jn = 1, jpj
1871	DO jm = 1,jpi
1872	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1873	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1874	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 * &
1875	(ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:))
1876	END DO
1877	END DO
1878	END DO
1879
1880	! gradient in the ij direction
1881	DO jk = 1,4
1882	DO jn = 1, jpj
1883	DO jm = 1,jpi
1884	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1885	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1886	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
1887	(ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni ,nj+2,:)) - &
1888	(ref_wgts(kw)%fly_dta(ni+2,nj ,:) - ref_wgts(kw)%fly_dta(ni ,nj ,:)))
1889	END DO
1890	END DO
1891	END DO
1892	!
1893	END IF
1894	!
1895	END SUBROUTINE fld_interp
1896
1897
1898	FUNCTION ksec_week( cdday )
1899	!!---------------------------------------------------------------------
1900	!! * FUNCTION kshift_week *
1901	!!
1902	!! ** Purpose :
1903	!!---------------------------------------------------------------------
1904	CHARACTER(len=*), INTENT(in) :: cdday !3 first letters of the first day of the weekly file
1905	!!
1906	INTEGER :: ksec_week ! output variable
1907	INTEGER :: ijul !temp variable
1908	INTEGER :: ishift !temp variable
1909	CHARACTER(len=3),DIMENSION(7) :: cl_week
1910	!!----------------------------------------------------------------------
1911	cl_week = (/"sun","sat","fri","thu","wed","tue","mon"/)
1912	DO ijul = 1, 7
1913	IF( cl_week(ijul) == TRIM(cdday) ) EXIT
1914	END DO
1915	IF( ijul .GT. 7 ) CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): '//TRIM(cdday) )
1916	!
1917	ishift = ijul * NINT(rday)
1918	!
1919	ksec_week = nsec_week + ishift
1920	ksec_week = MOD( ksec_week, 7*NINT(rday) )
1921	!
1922	END FUNCTION ksec_week
1923
1924	!!======================================================================
1925	END MODULE fldread

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

Download in other formats: