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 NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/SBC/fldread.F90 @ 12558

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

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

Download in other formats: