Context Navigation

fldread.F90 @ 11080

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

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

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

Context Navigation

source: branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 11080

Download in other formats: