Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

fldread.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
Property svn:keywords set to `Id`
File size: 108.6 KB

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

Download in other formats: