Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

fldread.F90

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

Download in other formats: