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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

fldread.F90 in branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 2865

Last change on this file since 2865 was 2865, checked in by davestorkey, 13 years ago
Updates for dynspg_exp option. Implement time_offset functionality in obc_dta. Add option to specify boundaries in the namelist. Re-activate obc_vol option. Update to namelist control of tidal harmonics.
Property svn:keywords set to `Id`
File size: 67.0 KB

Line
1	MODULE fldread
2	!!======================================================================
3	!! * MODULE fldread *
4	!! Ocean forcing: read input field for surface boundary condition
5	!!=====================================================================
6	!! History : 2.0 ! 06-2006 (S. Masson, G. Madec) Original code
7	!! ! 05-2008 (S. Alderson) Modified for Interpolation in memory
8	!! ! from input grid to model grid
9	!!----------------------------------------------------------------------
10
11	!!----------------------------------------------------------------------
12	!! fld_read : read input fields used for the computation of the
13	!! surface boundary condition
14	!!----------------------------------------------------------------------
15	USE oce ! ocean dynamics and tracers
16	USE dom_oce ! ocean space and time domain
17	USE phycst ! ???
18	USE in_out_manager ! I/O manager
19	USE iom ! I/O manager library
20	USE geo2ocean ! for vector rotation on to model grid
21	USE lib_mpp ! MPP library
22	USE ioipsl, ONLY : ymds2ju, ju2ymds ! for calendar
23
24	IMPLICIT NONE
25	PRIVATE
26
27	PUBLIC fld_map ! routine called by tides_init
28
29	TYPE, PUBLIC :: FLD_N !: Namelist field informations
30	CHARACTER(len = 256) :: clname ! generic name of the NetCDF flux file
31	INTEGER :: nfreqh ! frequency of each flux file
32	CHARACTER(len = 34) :: clvar ! generic name of the variable in the NetCDF flux file
33	LOGICAL :: ln_tint ! time interpolation or not (T/F)
34	LOGICAL :: ln_clim ! climatology or not (T/F)
35	CHARACTER(len = 8) :: cltype ! type of data file 'daily', 'monthly' or yearly'
36	CHARACTER(len = 34) :: wname ! generic name of a NetCDF weights file to be used, blank if not
37	CHARACTER(len = 34) :: vcomp ! symbolic component name if a vector that needs rotation
38	! ! a string starting with "U" or "V" for each component
39	! ! chars 2 onwards identify which components go together
40	END TYPE FLD_N
41
42	TYPE, PUBLIC :: FLD !: Input field related variables
43	CHARACTER(len = 256) :: clrootname ! generic name of the NetCDF file
44	CHARACTER(len = 256) :: clname ! current name of the NetCDF file
45	INTEGER :: nfreqh ! frequency of each flux file
46	CHARACTER(len = 34) :: clvar ! generic name of the variable in the NetCDF flux file
47	LOGICAL :: ln_tint ! time interpolation or not (T/F)
48	LOGICAL :: ln_clim ! climatology or not (T/F)
49	CHARACTER(len = 8) :: cltype ! type of data file 'daily', 'monthly' or yearly'
50	INTEGER :: num ! iom id of the jpfld files to be read
51	INTEGER , DIMENSION(2) :: nrec_b ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year)
52	INTEGER , DIMENSION(2) :: nrec_a ! after record (1: index, 2: second since Jan. 1st 00h of nit000 year)
53	REAL(wp) , ALLOCATABLE, DIMENSION(:,:,: ) :: fnow ! input fields interpolated to now time step
54	REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) :: fdta ! 2 consecutive record of input fields
55	CHARACTER(len = 256) :: wgtname ! current name of the NetCDF weight file acting as a key
56	! ! into the WGTLIST structure
57	CHARACTER(len = 34) :: vcomp ! symbolic name for a vector component that needs rotation
58	LOGICAL :: rotn ! flag to indicate whether field has been rotated
59	END TYPE FLD
60
61	TYPE, PUBLIC :: MAP_POINTER !: Array of integer pointers to 1D arrays
62	INTEGER, POINTER :: ptr(:)
63	END TYPE MAP_POINTER
64
65	!$AGRIF_DO_NOT_TREAT
66
67	!! keep list of all weights variables so they're only read in once
68	!! need to add AGRIF directives not to process this structure
69	!! also need to force wgtname to include AGRIF nest number
70	TYPE :: WGT !: Input weights related variables
71	CHARACTER(len = 256) :: wgtname ! current name of the NetCDF weight file
72	INTEGER , DIMENSION(2) :: ddims ! shape of input grid
73	INTEGER , DIMENSION(2) :: botleft ! top left corner of box in input grid containing
74	! ! current processor grid
75	INTEGER , DIMENSION(2) :: topright ! top right corner of box
76	INTEGER :: jpiwgt ! width of box on input grid
77	INTEGER :: jpjwgt ! height of box on input grid
78	INTEGER :: numwgt ! number of weights (4=bilinear, 16=bicubic)
79	INTEGER :: nestid ! for agrif, keep track of nest we're in
80	INTEGER :: overlap ! =0 when cyclic grid has no overlapping EW columns
81	! ! =>1 when they have one or more overlapping columns
82	! ! =-1 not cyclic
83	LOGICAL :: cyclic ! east-west cyclic or not
84	INTEGER, DIMENSION(:,:,:), POINTER :: data_jpi ! array of source integers
85	INTEGER, DIMENSION(:,:,:), POINTER :: data_jpj ! array of source integers
86	REAL(wp), DIMENSION(:,:,:), POINTER :: data_wgt ! array of weights on model grid
87	REAL(wp), DIMENSION(:,:,:), POINTER :: fly_dta ! array of values on input grid
88	REAL(wp), DIMENSION(:,:,:), POINTER :: col ! temporary array for reading in columns
89	END TYPE WGT
90
91	INTEGER, PARAMETER :: tot_wgts = 10
92	TYPE( WGT ), DIMENSION(tot_wgts) :: ref_wgts ! array of wgts
93	INTEGER :: nxt_wgt = 1 ! point to next available space in ref_wgts array
94
95	!$AGRIF_END_DO_NOT_TREAT
96
97	PUBLIC fld_read, fld_fill ! called by sbc... modules
98
99	!!----------------------------------------------------------------------
100	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
101	!! $Id$
102	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
103	!!----------------------------------------------------------------------
104	CONTAINS
105
106	SUBROUTINE fld_read( kt, kn_fsbc, sd, map, jit, time_offset )
107	!!---------------------------------------------------------------------
108	!! * ROUTINE fld_read *
109	!!
110	!! ** Purpose : provide at each time step the surface ocean fluxes
111	!! (momentum, heat, freshwater and runoff)
112	!!
113	!! ** Method : READ each input fields in NetCDF files using IOM
114	!! and intepolate it to the model time-step.
115	!! Several assumptions are made on the input file:
116	!! blahblahblah....
117	!!----------------------------------------------------------------------
118	INTEGER , INTENT(in ) :: kt ! ocean time step
119	INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step)
120	TYPE(FLD), INTENT(inout), DIMENSION(:) :: sd ! input field related variables
121	TYPE(MAP_POINTER),INTENT(in), OPTIONAL, DIMENSION(:) :: map ! global-to-local mapping index
122	INTEGER , INTENT(in ), OPTIONAL :: jit ! subcycle timestep for timesplitting option
123	INTEGER , INTENT(in ), OPTIONAL :: time_offset ! provide fields at time other than "now"
124	! time_offset = -1 => fields at "before" time level
125	! time_offset = +1 => fields at "after" time levels
126	! etc.
127	!!
128	INTEGER :: imf ! size of the structure sd
129	INTEGER :: jf ! dummy indices
130	INTEGER :: ireclast ! last record to be read in the current year file
131	INTEGER :: isecend ! number of second since Jan. 1st 00h of nit000 year at nitend
132	INTEGER :: isecsbc ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step
133	INTEGER :: time_add ! local time_offset variable
134	LOGICAL :: llnxtyr ! open next year file?
135	LOGICAL :: llnxtmth ! open next month file?
136	LOGICAL :: llstop ! stop is the file does not exist
137	LOGICAL :: ll_firstcall ! true if this is the first call to fld_read for this set of fields
138	REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation
139	REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation
140	CHARACTER(LEN=1000) :: clfmt ! write format
141	!!---------------------------------------------------------------------
142	ll_firstcall = .false.
143	IF( PRESENT(jit) ) THEN
144	IF(kt == nit000 .and. jit == 1) ll_firstcall = .true.
145	ELSE
146	IF(kt == nit000) ll_firstcall = .true.
147	ENDIF
148
149	time_add = 0
150	IF( PRESENT(time_offset) ) THEN
151	time_add = time_offset
152	ENDIF
153
154	! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
155	IF( present(jit) ) THEN
156	! ignore kn_fsbc in this case
157	isecsbc = nsec_year + nsec1jan000 + (jit+time_add)*rdt/REAL(nn_baro,wp)
158	ELSE
159	isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdttra(1)) + time_add * rdttra(1) ! middle of sbc time step
160	ENDIF
161	imf = SIZE( sd )
162	!
163	IF( ll_firstcall ) THEN ! initialization
164	IF( PRESENT(map) ) THEN
165	DO jf = 1, imf
166	CALL fld_init( kn_fsbc, sd(jf), map(jf)%ptr ) ! read each before field (put them in after as they will be swapped)
167	END DO
168	ELSE
169	DO jf = 1, imf
170	CALL fld_init( kn_fsbc, sd(jf) ) ! read each before field (put them in after as they will be swapped)
171	END DO
172	ENDIF
173	IF( lwp ) CALL wgt_print() ! control print
174	CALL fld_rot( kt, sd ) ! rotate vector fiels if needed
175	ENDIF
176	! ! ====================================== !
177	IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN ! update field at each kn_fsbc time-step !
178	! ! ====================================== !
179	!
180	DO jf = 1, imf ! --- loop over field --- !
181
182	IF( isecsbc > sd(jf)%nrec_a(2) .OR. ll_firstcall ) THEN ! read/update the after data?
183
184	IF( sd(jf)%ln_tint ) THEN ! swap before record field and informations
185	sd(jf)%nrec_b(:) = sd(jf)%nrec_a(:)
186	!CDIR COLLAPSE
187	sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)
188	ENDIF
189
190	IF( PRESENT(jit) ) THEN
191	CALL fld_rec( kn_fsbc, sd(jf), jit=jit ) ! update record informations
192	ELSE
193	CALL fld_rec( kn_fsbc, sd(jf) ) ! update record informations
194	ENDIF
195
196	! do we have to change the year/month/week/day of the forcing field??
197	IF( sd(jf)%ln_tint ) THEN
198	! if we do time interpolation we will need to open next year/month/week/day file before the end of the current
199	! 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)
200	! will be larger than the record number that should be read for current year/month/week/day
201
202	! last record to be read in the current file
203	IF ( sd(jf)%nfreqh == -12 ) THEN ; ireclast = 1 ! yearly mean
204	ELSEIF( sd(jf)%nfreqh == -1 ) THEN ! monthly mean
205	IF( sd(jf)%cltype == 'monthly' ) THEN ; ireclast = 1
206	ELSE ; ireclast = 12
207	ENDIF
208	ELSE ! higher frequency mean (in hours)
209	IF( sd(jf)%cltype == 'monthly' ) THEN ; ireclast = 24 * nmonth_len(nmonth) / sd(jf)%nfreqh
210	ELSEIF( sd(jf)%cltype(1:4) == 'week' ) THEN ; ireclast = 24 * 7 / sd(jf)%nfreqh
211	ELSEIF( sd(jf)%cltype == 'daily' ) THEN ; ireclast = 24 / sd(jf)%nfreqh
212	ELSE ; ireclast = 24 * nyear_len( 1 ) / sd(jf)%nfreqh
213	ENDIF
214	ENDIF
215
216	! do we need next file data?
217	IF( sd(jf)%nrec_a(1) > ireclast ) THEN
218
219	sd(jf)%nrec_a(1) = 1 ! force to read the first record of the next file
220
221	IF( .NOT. sd(jf)%ln_clim ) THEN ! close the current file and open a new one.
222
223	llnxtmth = sd(jf)%cltype == 'monthly' .OR. nday == nmonth_len(nmonth) ! open next month file?
224	llnxtyr = sd(jf)%cltype == 'yearly' .OR. (nmonth == 12 .AND. llnxtmth) ! open next year file?
225
226	! if the run finishes at the end of the current year/month/week/day, we will allow next
227	! year/month/week/day file to be not present. If the run continue further than the current
228	! year/month/week/day, next year/month/week/day file must exist
229	isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdttra(1)) ! second at the end of the run
230	llstop = isecend > sd(jf)%nrec_a(2) ! read more than 1 record of next year
231
232	CALL fld_clopn( sd(jf), nyear + COUNT((/llnxtyr /)) , &
233	& nmonth + COUNT((/llnxtmth/)) - 12 * COUNT((/llnxtyr /)), &
234	& nday + 1 - nmonth_len(nmonth) * COUNT((/llnxtmth/)), llstop )
235
236	IF( sd(jf)%num <= 0 .AND. .NOT. llstop ) THEN ! next year file does not exist
237	CALL ctl_warn('next year/month/week/day file: '//TRIM(sd(jf)%clname)// &
238	& ' not present -> back to current year/month/day')
239	CALL fld_clopn( sd(jf), nyear, nmonth, nday ) ! back to the current year/month/day
240	sd(jf)%nrec_a(1) = ireclast ! force to read the last record to be read in the current year file
241	ENDIF
242
243	ENDIF
244	ENDIF
245
246	ELSE
247	! if we are not doing time interpolation, we must change the year/month/week/day of the file just after
248	! switching to the NEW year/month/week/day. If it is the case, we are at the beginning of the
249	! year/month/week/day when calling fld_rec so sd(jf)%nrec_a(1) = 1
250	IF( sd(jf)%nrec_a(1) == 1 .AND. .NOT. ( sd(jf)%ln_clim .AND. sd(jf)%cltype == 'yearly' ) ) &
251	& CALL fld_clopn( sd(jf), nyear, nmonth, nday )
252	ENDIF
253
254	! read after data
255	IF( PRESENT(map) ) THEN
256	CALL fld_get( sd(jf), map(jf)%ptr )
257	ELSE
258	CALL fld_get( sd(jf) )
259	ENDIF
260
261	ENDIF
262	END DO ! --- end loop over field --- !
263
264	CALL fld_rot( kt, sd ) ! rotate vector fiels if needed
265
266	DO jf = 1, imf ! --- loop over field --- !
267	!
268	IF( sd(jf)%ln_tint ) THEN ! temporal interpolation
269	IF(lwp .AND. kt - nit000 <= 100 ) THEN
270	clfmt = "('fld_read: var ', a, ' kt = ', i8, ' (', f7.2,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," // &
271	& "', records b/a: ', i4.4, '/', i4.4, ' (days ', f7.2,'/', f7.2, ')')"
272	WRITE(numout, clfmt) TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday, &
273	& 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
274	WRITE(numout, *) 'time_add is : ',time_add
275	ENDIF
276	! temporal interpolation weights
277	ztinta = REAL( isecsbc - sd(jf)%nrec_b(2), wp ) / REAL( sd(jf)%nrec_a(2) - sd(jf)%nrec_b(2), wp )
278	ztintb = 1. - ztinta
279	!CDIR COLLAPSE
280	sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:,1) + ztinta * sd(jf)%fdta(:,:,:,2)
281	ELSE ! nothing to do...
282	IF(lwp .AND. kt - nit000 <= 100 ) THEN
283	clfmt = "('fld_read: var ', a, ' kt = ', i8,' (', f7.2,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," // &
284	& "', record: ', i4.4, ' (days ', f7.2, ' <-> ', f7.2, ')')"
285	WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday, &
286	& sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday
287	ENDIF
288	ENDIF
289	!
290	IF( kt == nitend - kn_fsbc + 1 ) CALL iom_close( sd(jf)%num ) ! Close the input files
291
292	END DO ! --- end loop over field --- !
293	!
294	! ! ====================================== !
295	ENDIF ! update field at each kn_fsbc time-step !
296	! ! ====================================== !
297	!
298	END SUBROUTINE fld_read
299
300
301	SUBROUTINE fld_init( kn_fsbc, sdjf, map )
302	!!---------------------------------------------------------------------
303	!! * ROUTINE fld_init *
304	!!
305	!! ** Purpose : - if time interpolation, read before data
306	!! - open current year file
307	!!----------------------------------------------------------------------
308	INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step)
309	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
310	INTEGER , INTENT(in), OPTIONAL, DIMENSION(:) :: map ! global-to-local mapping indices
311	!!
312	LOGICAL :: llprevyr ! are we reading previous year file?
313	LOGICAL :: llprevmth ! are we reading previous month file?
314	LOGICAL :: llprevweek ! are we reading previous week file?
315	LOGICAL :: llprevday ! are we reading previous day file?
316	LOGICAL :: llprev ! llprevyr .OR. llprevmth .OR. llprevweek .OR. llprevday
317	INTEGER :: idvar ! variable id
318	INTEGER :: inrec ! number of record existing for this variable
319	INTEGER :: iyear, imonth, iday ! first day of the current file in yyyy mm dd
320	INTEGER :: isec_week ! number of seconds since start of the weekly file
321	CHARACTER(LEN=1000) :: clfmt ! write format
322	!!---------------------------------------------------------------------
323
324	! some default definitions...
325	sdjf%num = 0 ! default definition for non-opened file
326	IF( sdjf%ln_clim ) sdjf%clname = TRIM( sdjf%clrootname ) ! file name defaut definition, never change in this case
327	llprevyr = .FALSE.
328	llprevmth = .FALSE.
329	llprevweek = .FALSE.
330	llprevday = .FALSE.
331	isec_week = 0
332
333	IF( sdjf%cltype(1:4) == 'week' .AND. nn_leapy == 0 ) &
334	& CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdjf%clrootname)//') needs nn_leapy = 1')
335	IF( sdjf%cltype(1:4) == 'week' .AND. sdjf%ln_clim ) &
336	& CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdjf%clrootname)//') needs ln_clim = .FALSE.')
337
338	! define record informations
339	CALL fld_rec( kn_fsbc, sdjf, ldbefore = .TRUE. ) ! return before values in sdjf%nrec_a (as we will swap it later)
340
341	! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
342
343	IF( sdjf%ln_tint ) THEN ! we need to read the previous record and we will put it in the current record structure
344
345	IF( sdjf%nrec_a(1) == 0 ) THEN ! we redefine record sdjf%nrec_a(1) with the last record of previous year file
346	IF ( sdjf%nfreqh == -12 ) THEN ! yearly mean
347	IF( sdjf%cltype == 'yearly' ) THEN ! yearly file
348	sdjf%nrec_a(1) = 1 ! force to read the unique record
349	llprevyr = .NOT. sdjf%ln_clim ! use previous year file?
350	ELSE
351	CALL ctl_stop( "fld_init: yearly mean file must be in a yearly type of file: "//TRIM(sdjf%clname) )
352	ENDIF
353	ELSEIF( sdjf%nfreqh == -1 ) THEN ! monthly mean
354	IF( sdjf%cltype == 'monthly' ) THEN ! monthly file
355	sdjf%nrec_a(1) = 1 ! force to read the unique record
356	llprevmth = .TRUE. ! use previous month file?
357	llprevyr = llprevmth .AND. nmonth == 1 ! use previous year file?
358	ELSE ! yearly file
359	sdjf%nrec_a(1) = 12 ! force to read december mean
360	llprevyr = .NOT. sdjf%ln_clim ! use previous year file?
361	ENDIF
362	ELSE ! higher frequency mean (in hours)
363	IF ( sdjf%cltype == 'monthly' ) THEN ! monthly file
364	sdjf%nrec_a(1) = 24 * nmonth_len(nmonth-1) / sdjf%nfreqh ! last record of previous month
365	llprevmth = .TRUE. ! use previous month file?
366	llprevyr = llprevmth .AND. nmonth == 1 ! use previous year file?
367	ELSEIF( sdjf%cltype(1:4) == 'week' ) THEN ! weekly file
368	llprevweek = .TRUE. ! use previous week file?
369	sdjf%nrec_a(1) = 24 * 7 / sdjf%nfreqh ! last record of previous week
370	isec_week = NINT(rday) * 7 ! add a shift toward previous week
371	ELSEIF( sdjf%cltype == 'daily' ) THEN ! daily file
372	sdjf%nrec_a(1) = 24 / sdjf%nfreqh ! last record of previous day
373	llprevday = .TRUE. ! use previous day file?
374	llprevmth = llprevday .AND. nday == 1 ! use previous month file?
375	llprevyr = llprevmth .AND. nmonth == 1 ! use previous year file?
376	ELSE ! yearly file
377	sdjf%nrec_a(1) = 24 * nyear_len(0) / sdjf%nfreqh ! last record of previous year
378	llprevyr = .NOT. sdjf%ln_clim ! use previous year file?
379	ENDIF
380	ENDIF
381	ENDIF
382	IF ( sdjf%cltype(1:4) == 'week' ) THEN
383	isec_week = isec_week + ksec_week( sdjf%cltype(6:8) ) ! second since the beginning of the week
384	llprevmth = isec_week > nsec_month ! longer time since the beginning of the week than the month
385	llprevyr = llprevmth .AND. nmonth == 1
386	ENDIF
387	llprev = llprevyr .OR. llprevmth .OR. llprevweek .OR. llprevday
388	!
389	iyear = nyear - COUNT((/llprevyr /))
390	imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
391	iday = nday - COUNT((/llprevday/)) + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
392	!
393	CALL fld_clopn( sdjf, iyear, imonth, iday, .NOT. llprev )
394
395	! if previous year/month/day file does not exist, we switch to the current year/month/day
396	IF( llprev .AND. sdjf%num <= 0 ) THEN
397	CALL ctl_warn( 'previous year/month/week/day file: '//TRIM(sdjf%clname)// &
398	& ' not present -> back to current year/month/week/day' )
399	! we force to read the first record of the current year/month/day instead of last record of previous year/month/day
400	llprev = .FALSE.
401	sdjf%nrec_a(1) = 1
402	CALL fld_clopn( sdjf, nyear, nmonth, nday )
403	ENDIF
404
405	IF( llprev ) THEN ! check if the last record sdjf%nrec_n(1) exists in the file
406	idvar = iom_varid( sdjf%num, sdjf%clvar ) ! id of the variable sdjf%clvar
407	IF( idvar <= 0 ) RETURN
408	inrec = iom_file( sdjf%num )%dimsz( iom_file( sdjf%num )%ndims(idvar), idvar ) ! size of the last dim of idvar
409	sdjf%nrec_a(1) = MIN( sdjf%nrec_a(1), inrec ) ! make sure we select an existing record
410	ENDIF
411
412	! read before data
413	IF( PRESENT(map) ) THEN
414	CALL fld_get( sdjf, map ) ! read before values in after arrays(as we will swap it later)
415	ELSE
416	CALL fld_get( sdjf ) ! read before values in after arrays(as we will swap it later)
417	ENDIF
418
419	clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i4, ' at time = ', f7.2, ' days')"
420	IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_a(1), REAL(sdjf%nrec_a(2),wp)/rday
421
422	IF( llprev ) CALL iom_close( sdjf%num ) ! force to close previous year file (-> redefine sdjf%num to 0)
423
424	ENDIF
425
426	! make sure current year/month/day file is opened
427	IF( sdjf%num <= 0 ) THEN
428	!
429	IF ( sdjf%cltype(1:4) == 'week' ) THEN
430	isec_week = ksec_week( sdjf%cltype(6:8) ) ! second since the beginning of the week
431	llprevmth = isec_week > nsec_month ! longer time since beginning of the week than the month
432	llprevyr = llprevmth .AND. nmonth == 1
433	ELSE
434	isec_week = 0
435	llprevmth = .FALSE.
436	llprevyr = .FALSE.
437	ENDIF
438	!
439	iyear = nyear - COUNT((/llprevyr /))
440	imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
441	iday = nday + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
442	!
443	CALL fld_clopn( sdjf, iyear, imonth, iday )
444	ENDIF
445	!
446	END SUBROUTINE fld_init
447
448
449	SUBROUTINE fld_rec( kn_fsbc, sdjf, ldbefore, jit )
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 :: jit ! index of barotropic subcycle
464	! used only if sdjf%ln_tint = .TRUE.
465	!!
466	LOGICAL :: llbefore ! local definition of ldbefore
467	INTEGER :: iendrec ! end of this record (in seconds)
468	INTEGER :: imth ! month number
469	INTEGER :: ifreq_sec ! frequency mean (in seconds)
470	INTEGER :: isec_week ! number of seconds since the start of the weekly file
471	REAL(wp) :: ztmp ! temporary variable
472	!!----------------------------------------------------------------------
473	!
474	! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
475	!
476	IF( PRESENT(ldbefore) ) THEN ; llbefore = ldbefore .AND. sdjf%ln_tint ! needed only if sdjf%ln_tint = .TRUE.
477	ELSE ; llbefore = .FALSE.
478	ENDIF
479	!
480	! ! =========== !
481	IF ( sdjf%nfreqh == -12 ) THEN ! yearly mean
482	! ! =========== !
483	!
484	IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record
485	!
486	! INT( ztmp )
487	! /\|\
488	! 1 \| *----
489	! 0 \|----(
490	! \|----+----\|--> time
491	! 0 /\|\ 1 (nday/nyear_len(1))
492	! \|
493	! \|
494	! forcing record : 1
495	!
496	ztmp = REAL( nday, wp ) / REAL( nyear_len(1), wp ) + 0.5
497	IF( PRESENT(jit) ) ztmp = ztmp + jit*rdt/REAL(nn_baro,wp)
498	sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
499	! swap at the middle of the year
500	IF( llbefore ) THEN ; sdjf%nrec_a(2) = nsec1jan000 - NINT(0.5 * rday) * nyear_len(0)
501	ELSE ; sdjf%nrec_a(2) = nsec1jan000 + NINT(0.5 * rday) * nyear_len(1)
502	ENDIF
503	ELSE ! no time interpolation
504	sdjf%nrec_a(1) = 1
505	sdjf%nrec_a(2) = NINT(rday) * nyear_len(1) + nsec1jan000 ! swap at the end of the year
506	sdjf%nrec_b(2) = nsec1jan000 ! beginning of the year (only for print)
507	ENDIF
508	!
509	! ! ============ !
510	ELSEIF( sdjf%nfreqh == -1 ) THEN ! monthly mean !
511	! ! ============ !
512	!
513	IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record
514	!
515	! INT( ztmp )
516	! /\|\
517	! 1 \| *----
518	! 0 \|----(
519	! \|----+----\|--> time
520	! 0 /\|\ 1 (nday/nmonth_len(nmonth))
521	! \|
522	! \|
523	! forcing record : nmonth
524	!
525	ztmp = REAL( nday, wp ) / REAL( nmonth_len(nmonth), wp ) + 0.5
526	IF( PRESENT(jit) ) ztmp = ztmp + jit*rdt/REAL(nn_baro,wp)
527	imth = nmonth + INT( ztmp ) - COUNT((/llbefore/))
528	IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
529	ELSE ; sdjf%nrec_a(1) = imth
530	ENDIF
531	sdjf%nrec_a(2) = nmonth_half( imth ) + nsec1jan000 ! swap at the middle of the month
532	ELSE ! no time interpolation
533	IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nrec_a(1) = 1
534	ELSE ; sdjf%nrec_a(1) = nmonth
535	ENDIF
536	sdjf%nrec_a(2) = nmonth_end(nmonth ) + nsec1jan000 ! swap at the end of the month
537	sdjf%nrec_b(2) = nmonth_end(nmonth-1) + nsec1jan000 ! beginning of the month (only for print)
538	ENDIF
539	!
540	! ! ================================ !
541	ELSE ! higher frequency mean (in hours)
542	! ! ================================ !
543	!
544	ifreq_sec = sdjf%nfreqh * 3600 ! frequency mean (in seconds)
545	IF( sdjf%cltype(1:4) == 'week' ) isec_week = ksec_week( sdjf%cltype(6:8) ) ! since the first day of the current week
546	! number of second since the beginning of the file
547	IF( sdjf%cltype == 'monthly' ) THEN ; ztmp = REAL(nsec_month,wp) ! since the first day of the current month
548	ELSEIF( sdjf%cltype(1:4) == 'week' ) THEN ; ztmp = REAL(isec_week ,wp) ! since the first day of the current week
549	ELSEIF( sdjf%cltype == 'daily' ) THEN ; ztmp = REAL(nsec_day ,wp) ! since 00h of the current day
550	ELSE ; ztmp = REAL(nsec_year ,wp) ! since 00h on Jan 1 of the current year
551	ENDIF
552	ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdttra(1) ! shift time to be centrered in the middle of sbc time step
553	ztmp = ztmp + 0.01 * rdttra(1) ! add 0.01 time step to avoid truncation error
554	IF( PRESENT(jit) ) ztmp = ztmp + jit*rdt/REAL(nn_baro,wp)
555	IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record
556	!
557	! INT( ztmp )
558	! /\|\
559	! 2 \| *-----(
560	! 1 \| *-----(
561	! 0 \|--(
562	! \|--+--\|--+--\|--+--\|--> time
563	! 0 /\|\ 1 /\|\ 2 /\|\ 3 (nsec_year/ifreq_sec) or (nsec_month/ifreq_sec)
564	! \| \| \|
565	! \| \| \|
566	! forcing record : 1 2 3
567	!
568	ztmp= ztmp / REAL(ifreq_sec, wp) + 0.5
569	ELSE ! no time interpolation
570	!
571	! INT( ztmp )
572	! /\|\
573	! 2 \| *-----(
574	! 1 \| *-----(
575	! 0 \|-----(
576	! \|--+--\|--+--\|--+--\|--> time
577	! 0 /\|\ 1 /\|\ 2 /\|\ 3 (nsec_year/ifreq_sec) or (nsec_month/ifreq_sec)
578	! \| \| \|
579	! \| \| \|
580	! forcing record : 1 2 3
581	!
582	ztmp= ztmp / REAL(ifreq_sec, wp)
583	ENDIF
584	sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/)) ! record nomber to be read
585
586	iendrec = ifreq_sec * sdjf%nrec_a(1) + nsec1jan000 ! end of this record (in second)
587	! add the number of seconds between 00h Jan 1 and the end of previous month/week/day (ok if nmonth=1)
588	IF( sdjf%cltype == 'monthly' ) iendrec = iendrec + NINT(rday) * SUM(nmonth_len(1:nmonth -1))
589	IF( sdjf%cltype(1:4) == 'week' ) iendrec = iendrec + ( nsec_year - isec_week )
590	IF( sdjf%cltype == 'daily' ) iendrec = iendrec + NINT(rday) * ( nday_year - 1 )
591	IF( sdjf%ln_tint ) THEN
592	sdjf%nrec_a(2) = iendrec - ifreq_sec / 2 ! swap at the middle of the record
593	ELSE
594	sdjf%nrec_a(2) = iendrec ! swap at the end of the record
595	sdjf%nrec_b(2) = iendrec - ifreq_sec ! beginning of the record (only for print)
596	ENDIF
597	!
598	ENDIF
599	!
600	END SUBROUTINE fld_rec
601
602
603	SUBROUTINE fld_get( sdjf, map )
604	!!---------------------------------------------------------------------
605	!! * ROUTINE fld_get *
606	!!
607	!! ** Purpose : read the data
608	!!----------------------------------------------------------------------
609	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
610	INTEGER , INTENT(in), OPTIONAL, DIMENSION(:) :: map ! global-to-local mapping indices
611	!!
612	INTEGER :: ipk ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
613	INTEGER :: iw ! index into wgts array
614	!!---------------------------------------------------------------------
615
616	ipk = SIZE( sdjf%fnow, 3 )
617
618	IF( PRESENT(map) ) THEN
619	IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map )
620	ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1), map )
621	ENDIF
622	ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
623	CALL wgt_list( sdjf, iw )
624	IF( sdjf%ln_tint ) THEN ; CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) )
625	ELSE ; CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fnow(:,:,: ), sdjf%nrec_a(1) )
626	ENDIF
627	ELSE
628	SELECT CASE( ipk )
629	CASE(1)
630	IF( sdjf%ln_tint ) THEN ; CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,1,2), sdjf%nrec_a(1) )
631	ELSE ; CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fnow(:,:,1 ), sdjf%nrec_a(1) )
632	ENDIF
633	CASE DEFAULT
634	IF( sdjf%ln_tint ) THEN ; CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) )
635	ELSE ; CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1) )
636	ENDIF
637	END SELECT
638	ENDIF
639	!
640	sdjf%rotn = .false. ! vector not yet rotated
641
642	END SUBROUTINE fld_get
643
644	SUBROUTINE fld_map( num, clvar, dta, nrec, map )
645	!!---------------------------------------------------------------------
646	!! * ROUTINE fld_get *
647	!!
648	!! ** Purpose : read global data from file and map onto local data
649	!! using a general mapping (for open boundaries)
650	!!----------------------------------------------------------------------
651	#if defined key_obc
652	USE obc_oce, ONLY: dta_global ! workspace to read in global data arrays
653	#endif
654
655	INTEGER , INTENT(in ) :: num ! stream number
656	CHARACTER(LEN=*) , INTENT(in ) :: clvar ! variable name
657	REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional)
658	INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice)
659	INTEGER, DIMENSION(:) , INTENT(in ) :: map ! global-to-local mapping indices
660	!!
661	INTEGER :: ipi ! length of boundary data on local process
662	INTEGER :: ipj ! length of dummy dimension ( = 1 )
663	INTEGER :: ipk ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk )
664	INTEGER :: ilendta ! length of data in file
665	INTEGER :: idvar ! variable ID
666	INTEGER :: ib, ik ! loop counters
667	INTEGER :: ierr
668	REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read ! work space for global data
669	!!---------------------------------------------------------------------
670
671	#if defined key_obc
672	dta_read => dta_global
673	#endif
674
675	ipi = SIZE( dta, 1 )
676	ipj = 1
677	ipk = SIZE( dta, 3 )
678
679	idvar = iom_varid( num, clvar )
680	ilendta = iom_file(num)%dimsz(1,idvar)
681	IF(lwp) WRITE(numout,*) 'Dim size for ',TRIM(clvar),' is ', ilendta
682	IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk
683
684	SELECT CASE( ipk )
685	CASE(1)
686	CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1 ), nrec )
687	CASE DEFAULT
688	CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec )
689	END SELECT
690	!
691	DO ib = 1, ipi
692	DO ik = 1, ipk
693	dta(ib,1,ik) = dta_read(map(ib),1,ik)
694	END DO
695	END DO
696
697	END SUBROUTINE fld_map
698
699
700	SUBROUTINE fld_rot( kt, sd )
701	!!---------------------------------------------------------------------
702	!! * ROUTINE fld_rot *
703	!!
704	!! ** Purpose : Vector fields may need to be rotated onto the local grid direction
705	!!----------------------------------------------------------------------
706	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
707	USE wrk_nemo, ONLY: utmp => wrk_2d_24, vtmp => wrk_2d_25 ! 2D workspace
708	!!
709	INTEGER , INTENT(in ) :: kt ! ocean time step
710	TYPE(FLD), INTENT(inout), DIMENSION(:) :: sd ! input field related variables
711	!!
712	INTEGER :: ju, jv, jk ! loop indices
713	INTEGER :: imf ! size of the structure sd
714	INTEGER :: ill ! character length
715	INTEGER :: iv ! indice of V component
716	CHARACTER (LEN=100) :: clcomp ! dummy weight name
717	!!---------------------------------------------------------------------
718
719	IF(wrk_in_use(2, 24,25) ) THEN
720	CALL ctl_stop('fld_rot: ERROR: requested workspace arrays are unavailable.') ; RETURN
721	END IF
722
723	!! (sga: following code should be modified so that pairs arent searched for each time
724	!
725	imf = SIZE( sd )
726	DO ju = 1, imf
727	ill = LEN_TRIM( sd(ju)%vcomp )
728	IF( ill > 0 .AND. .NOT. sd(ju)%rotn ) THEN ! find vector rotations required
729	IF( sd(ju)%vcomp(1:1) == 'U' ) THEN ! east-west component has symbolic name starting with 'U'
730	! look for the north-south component which has same symbolic name but with 'U' replaced with 'V'
731	clcomp = 'V' // sd(ju)%vcomp(2:ill) ! works even if ill == 1
732	iv = -1
733	DO jv = 1, imf
734	IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) ) iv = jv
735	END DO
736	IF( iv > 0 ) THEN ! fields ju and iv are two components which need to be rotated together
737	DO jk = 1, SIZE( sd(ju)%fnow, 3 )
738	IF( sd(ju)%ln_tint )THEN
739	CALL rot_rep( sd(ju)%fdta(:,:,jk,2), sd(iv)%fdta(:,:,jk,2), 'T', 'en->i', utmp(:,:) )
740	CALL rot_rep( sd(ju)%fdta(:,:,jk,2), sd(iv)%fdta(:,:,jk,2), 'T', 'en->j', vtmp(:,:) )
741	sd(ju)%fdta(:,:,jk,2) = utmp(:,:) ; sd(iv)%fdta(:,:,jk,2) = vtmp(:,:)
742	ELSE
743	CALL rot_rep( sd(ju)%fnow(:,:,jk ), sd(iv)%fnow(:,:,jk ), 'T', 'en->i', utmp(:,:) )
744	CALL rot_rep( sd(ju)%fnow(:,:,jk ), sd(iv)%fnow(:,:,jk ), 'T', 'en->j', vtmp(:,:) )
745	sd(ju)%fnow(:,:,jk ) = utmp(:,:) ; sd(iv)%fnow(:,:,jk ) = vtmp(:,:)
746	ENDIF
747	END DO
748	sd(ju)%rotn = .TRUE. ! vector was rotated
749	IF( lwp .AND. kt == nit000 ) WRITE(numout,*) &
750	& 'fld_read: vector pair ('//TRIM(sd(ju)%clvar)//', '//TRIM(sd(iv)%clvar)//') rotated on to model grid'
751	ENDIF
752	ENDIF
753	ENDIF
754	END DO
755	!
756	IF(wrk_not_released(2, 24,25) ) CALL ctl_stop('fld_rot: ERROR: failed to release workspace arrays.')
757	!
758	END SUBROUTINE fld_rot
759
760
761	SUBROUTINE fld_clopn( sdjf, kyear, kmonth, kday, ldstop )
762	!!---------------------------------------------------------------------
763	!! * ROUTINE fld_clopn *
764	!!
765	!! ** Purpose : update the file name and open the file
766	!!----------------------------------------------------------------------
767	TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables
768	INTEGER , INTENT(in ) :: kyear ! year value
769	INTEGER , INTENT(in ) :: kmonth ! month value
770	INTEGER , INTENT(in ) :: kday ! day value
771	LOGICAL, OPTIONAL, INTENT(in ) :: ldstop ! stop if open to read a non-existing file (default = .TRUE.)
772	!!----------------------------------------------------------------------
773
774	IF( sdjf%num /= 0 ) CALL iom_close( sdjf%num ) ! close file if already open
775	! build the new filename if not climatological data
776	sdjf%clname=TRIM(sdjf%clrootname)
777	!
778	! note that sdjf%ln_clim is is only acting on presence of the year in the file
779	IF( .NOT. sdjf%ln_clim ) THEN
780	WRITE(sdjf%clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), kyear ! add year
781	IF( sdjf%cltype /= 'yearly' ) WRITE(sdjf%clname, '(a,"m" ,i2.2)' ) TRIM( sdjf%clname ), kmonth ! add month
782	ELSE
783	! build the new filename if climatological data
784	IF( sdjf%cltype /= 'yearly' ) WRITE(sdjf%clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), kmonth ! add month
785	ENDIF
786	IF( sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) &
787	& WRITE(sdjf%clname, '(a,"d" ,i2.2)' ) TRIM( sdjf%clname ), kday ! add day
788	!
789	CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 )
790	!
791	END SUBROUTINE fld_clopn
792
793
794	SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )
795	!!---------------------------------------------------------------------
796	!! * ROUTINE fld_fill *
797	!!
798	!! ** Purpose : fill sdf with sdf_n and control print
799	!!----------------------------------------------------------------------
800	TYPE(FLD) , DIMENSION(:), INTENT(inout) :: sdf ! structure of input fields (file informations, fields read)
801	TYPE(FLD_N), DIMENSION(:), INTENT(in ) :: sdf_n ! array of namelist information structures
802	CHARACTER(len=*) , INTENT(in ) :: cdir ! Root directory for location of flx files
803	CHARACTER(len=*) , INTENT(in ) :: cdcaller !
804	CHARACTER(len=*) , INTENT(in ) :: cdtitle !
805	CHARACTER(len=*) , INTENT(in ) :: cdnam !
806	!
807	INTEGER :: jf ! dummy indices
808	!!---------------------------------------------------------------------
809
810	DO jf = 1, SIZE(sdf)
811	sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
812	sdf(jf)%nfreqh = sdf_n(jf)%nfreqh
813	sdf(jf)%clvar = sdf_n(jf)%clvar
814	sdf(jf)%ln_tint = sdf_n(jf)%ln_tint
815	sdf(jf)%ln_clim = sdf_n(jf)%ln_clim
816	sdf(jf)%cltype = sdf_n(jf)%cltype
817	sdf(jf)%wgtname = " "
818	IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
819	sdf(jf)%vcomp = sdf_n(jf)%vcomp
820	END DO
821
822	IF(lwp) THEN ! control print
823	WRITE(numout,*)
824	WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
825	WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
826	WRITE(numout,*) ' '//TRIM( cdnam )//' Namelist'
827	WRITE(numout,*) ' list of files and frequency (>0: in hours ; <0 in months)'
828	DO jf = 1, SIZE(sdf)
829	WRITE(numout,*) ' root filename: ' , TRIM( sdf(jf)%clrootname ), &
830	& ' variable name: ' , TRIM( sdf(jf)%clvar )
831	WRITE(numout,*) ' frequency: ' , sdf(jf)%nfreqh , &
832	& ' time interp: ' , sdf(jf)%ln_tint , &
833	& ' climatology: ' , sdf(jf)%ln_clim , &
834	& ' weights : ' , TRIM( sdf(jf)%wgtname ), &
835	& ' pairing : ' , TRIM( sdf(jf)%vcomp ), &
836	& ' data type: ' , sdf(jf)%cltype
837	call flush(numout)
838	END DO
839	ENDIF
840
841	END SUBROUTINE fld_fill
842
843
844	SUBROUTINE wgt_list( sd, kwgt )
845	!!---------------------------------------------------------------------
846	!! * ROUTINE wgt_list *
847	!!
848	!! ** Purpose : search array of WGTs and find a weights file
849	!! entry, or return a new one adding it to the end
850	!! if it is a new entry, the weights data is read in and
851	!! restructured (fld_weight)
852	!!----------------------------------------------------------------------
853	TYPE( FLD ), INTENT(in ) :: sd ! field with name of weights file
854	INTEGER , INTENT(inout) :: kwgt ! index of weights
855	!!
856	INTEGER :: kw, nestid ! local integer
857	LOGICAL :: found ! local logical
858	!!----------------------------------------------------------------------
859	!
860	!! search down linked list
861	!! weights filename is either present or we hit the end of the list
862	found = .FALSE.
863
864	!! because agrif nest part of filenames are now added in iom_open
865	!! to distinguish between weights files on the different grids, need to track
866	!! nest number explicitly
867	nestid = 0
868	#if defined key_agrif
869	nestid = Agrif_Fixed()
870	#endif
871	DO kw = 1, nxt_wgt-1
872	IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
873	ref_wgts(kw)%nestid == nestid) THEN
874	kwgt = kw
875	found = .TRUE.
876	EXIT
877	ENDIF
878	END DO
879	IF( .NOT.found ) THEN
880	kwgt = nxt_wgt
881	CALL fld_weight( sd )
882	ENDIF
883	!
884	END SUBROUTINE wgt_list
885
886
887	SUBROUTINE wgt_print( )
888	!!---------------------------------------------------------------------
889	!! * ROUTINE wgt_print *
890	!!
891	!! ** Purpose : print the list of known weights
892	!!----------------------------------------------------------------------
893	INTEGER :: kw !
894	!!----------------------------------------------------------------------
895	!
896	DO kw = 1, nxt_wgt-1
897	WRITE(numout,*) 'weight file: ',TRIM(ref_wgts(kw)%wgtname)
898	WRITE(numout,*) ' ddims: ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
899	WRITE(numout,*) ' numwgt: ',ref_wgts(kw)%numwgt
900	WRITE(numout,*) ' jpiwgt: ',ref_wgts(kw)%jpiwgt
901	WRITE(numout,*) ' jpjwgt: ',ref_wgts(kw)%jpjwgt
902	WRITE(numout,*) ' botleft: ',ref_wgts(kw)%botleft
903	WRITE(numout,*) ' topright: ',ref_wgts(kw)%topright
904	IF( ref_wgts(kw)%cyclic ) THEN
905	WRITE(numout,*) ' cyclical'
906	IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) ' with overlap of ', ref_wgts(kw)%overlap
907	ELSE
908	WRITE(numout,*) ' not cyclical'
909	ENDIF
910	IF( ASSOCIATED(ref_wgts(kw)%data_wgt) ) WRITE(numout,*) ' allocated'
911	END DO
912	!
913	END SUBROUTINE wgt_print
914
915
916	SUBROUTINE fld_weight( sd )
917	!!---------------------------------------------------------------------
918	!! * ROUTINE fld_weight *
919	!!
920	!! ** Purpose : create a new WGT structure and fill in data from
921	!! file, restructuring as required
922	!!----------------------------------------------------------------------
923	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released, iwrk_in_use, iwrk_not_released
924	USE wrk_nemo, ONLY: data_tmp => wrk_2d_1 ! 2D real workspace
925	USE wrk_nemo, ONLY: data_src => iwrk_2d_1 ! 2D integer workspace
926	!!
927	TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file
928	!!
929	INTEGER :: jn ! dummy loop indices
930	INTEGER :: inum ! temporary logical unit
931	INTEGER :: id ! temporary variable id
932	INTEGER :: ipk ! temporary vertical dimension
933	CHARACTER (len=5) :: aname
934	INTEGER , DIMENSION(3) :: ddims
935	LOGICAL :: cyclical
936	INTEGER :: zwrap ! local integer
937	!!----------------------------------------------------------------------
938	!
939	IF( wrk_in_use(2, 1) .OR. iwrk_in_use(2,1) ) THEN
940	CALL ctl_stop('fld_weight: requested workspace arrays are unavailable') ; RETURN
941	ENDIF
942	!
943	IF( nxt_wgt > tot_wgts ) THEN
944	CALL ctl_stop("fld_weight: weights array size exceeded, increase tot_wgts")
945	ENDIF
946	!
947	!! new weights file entry, add in extra information
948	!! a weights file represents a 2D grid of a certain shape, so we assume that the current
949	!! input data file is representative of all other files to be opened and processed with the
950	!! current weights file
951
952	!! open input data file (non-model grid)
953	CALL iom_open( sd%clname, inum, ldiof = LEN(TRIM(sd%wgtname)) > 0 )
954
955	!! get dimensions
956	id = iom_varid( inum, sd%clvar, ddims )
957
958	!! close it
959	CALL iom_close( inum )
960
961	!! now open the weights file
962
963	CALL iom_open ( sd%wgtname, inum ) ! interpolation weights
964	IF ( inum > 0 ) THEN
965
966	!! determine whether we have an east-west cyclic grid
967	!! from global attribute called "ew_wrap" in the weights file
968	!! note that if not found, iom_getatt returns -999 and cyclic with no overlap is assumed
969	!! since this is the most common forcing configuration
970
971	CALL iom_getatt(inum, 'ew_wrap', zwrap)
972	IF( zwrap >= 0 ) THEN
973	cyclical = .TRUE.
974	ELSE IF( zwrap == -999 ) THEN
975	cyclical = .TRUE.
976	zwrap = 0
977	ELSE
978	cyclical = .FALSE.
979	ENDIF
980
981	ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
982	ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
983	ref_wgts(nxt_wgt)%wgtname = sd%wgtname
984	ref_wgts(nxt_wgt)%overlap = zwrap
985	ref_wgts(nxt_wgt)%cyclic = cyclical
986	ref_wgts(nxt_wgt)%nestid = 0
987	#if defined key_agrif
988	ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
989	#endif
990	!! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
991	!! for each weight wgtNN there is an integer array srcNN which gives the point in
992	!! the input data grid which is to be multiplied by the weight
993	!! they are both arrays on the model grid so the result of the multiplication is
994	!! added into an output array on the model grid as a running sum
995
996	!! two possible cases: bilinear (4 weights) or bicubic (16 weights)
997	id = iom_varid(inum, 'src05', ldstop=.FALSE.)
998	IF( id <= 0) THEN
999	ref_wgts(nxt_wgt)%numwgt = 4
1000	ELSE
1001	ref_wgts(nxt_wgt)%numwgt = 16
1002	ENDIF
1003
1004	ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
1005	ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
1006	ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
1007
1008	DO jn = 1,4
1009	aname = ' '
1010	WRITE(aname,'(a3,i2.2)') 'src',jn
1011	data_tmp(:,:) = 0
1012	CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
1013	data_src(:,:) = INT(data_tmp(:,:))
1014	ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
1015	ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
1016	END DO
1017
1018	DO jn = 1, ref_wgts(nxt_wgt)%numwgt
1019	aname = ' '
1020	WRITE(aname,'(a3,i2.2)') 'wgt',jn
1021	ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
1022	CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
1023	END DO
1024	CALL iom_close (inum)
1025
1026	! find min and max indices in grid
1027	ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1028	ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1029	ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1030	ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1031
1032	! and therefore dimensions of the input box
1033	ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
1034	ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
1035
1036	! shift indexing of source grid
1037	ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
1038	ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
1039
1040	! create input grid, give it a halo to allow gradient calculations
1041	! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
1042	! a more robust solution will be given in next release
1043	ipk = SIZE(sd%fnow, 3)
1044	ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
1045	IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
1046
1047	nxt_wgt = nxt_wgt + 1
1048
1049	ELSE
1050	CALL ctl_stop( ' fld_weight : unable to read the file ' )
1051	ENDIF
1052
1053	IF( wrk_not_released(2, 1) .OR. &
1054	iwrk_not_released(2, 1) ) CALL ctl_stop('fld_weight: failed to release workspace arrays')
1055	!
1056	END SUBROUTINE fld_weight
1057
1058
1059	SUBROUTINE fld_interp( num, clvar, kw, kk, dta, nrec )
1060	!!---------------------------------------------------------------------
1061	!! * ROUTINE fld_interp *
1062	!!
1063	!! ** Purpose : apply weights to input gridded data to create data
1064	!! on model grid
1065	!!----------------------------------------------------------------------
1066	INTEGER , INTENT(in ) :: num ! stream number
1067	CHARACTER(LEN=*) , INTENT(in ) :: clvar ! variable name
1068	INTEGER , INTENT(in ) :: kw ! weights number
1069	INTEGER , INTENT(in ) :: kk ! vertical dimension of kk
1070	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: dta ! output field on model grid
1071	INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice)
1072	!!
1073	INTEGER, DIMENSION(3) :: rec1,recn ! temporary arrays for start and length
1074	INTEGER :: jk, jn, jm ! loop counters
1075	INTEGER :: ni, nj ! lengths
1076	INTEGER :: jpimin,jpiwid ! temporary indices
1077	INTEGER :: jpjmin,jpjwid ! temporary indices
1078	INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices
1079	!!----------------------------------------------------------------------
1080	!
1081	!! for weighted interpolation we have weights at four corners of a box surrounding
1082	!! a model grid point, each weight is multiplied by a grid value (bilinear case)
1083	!! or by a grid value and gradients at the corner point (bicubic case)
1084	!! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
1085
1086	!! sub grid from non-model input grid which encloses all grid points in this nemo process
1087	jpimin = ref_wgts(kw)%botleft(1)
1088	jpjmin = ref_wgts(kw)%botleft(2)
1089	jpiwid = ref_wgts(kw)%jpiwgt
1090	jpjwid = ref_wgts(kw)%jpjwgt
1091
1092	!! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients
1093	rec1(1) = MAX( jpimin-1, 1 )
1094	rec1(2) = MAX( jpjmin-1, 1 )
1095	rec1(3) = 1
1096	recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
1097	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1098	recn(3) = kk
1099
1100	!! where we need to put it in the non-nemo grid fly_dta
1101	!! note that jpi1 and jpj1 only differ from 1 when jpimin and jpjmin are 1
1102	!! (ie at the extreme west or south of the whole input grid) and similarly for jpi2 and jpj2
1103	jpi1 = 2 + rec1(1) - jpimin
1104	jpj1 = 2 + rec1(2) - jpjmin
1105	jpi2 = jpi1 + recn(1) - 1
1106	jpj2 = jpj1 + recn(2) - 1
1107
1108	ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1109	SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) )
1110	CASE(1)
1111	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn)
1112	CASE DEFAULT
1113	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn)
1114	END SELECT
1115
1116	!! first four weights common to both bilinear and bicubic
1117	!! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft
1118	!! note that we have to offset by 1 into fly_dta array because of halo
1119	dta(:,:,:) = 0.0
1120	DO jk = 1,4
1121	DO jn = 1, jpj
1122	DO jm = 1,jpi
1123	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1124	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1125	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,:)
1126	END DO
1127	END DO
1128	END DO
1129
1130	IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
1131
1132	!! fix up halo points that we couldnt read from file
1133	IF( jpi1 == 2 ) THEN
1134	ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:)
1135	ENDIF
1136	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1137	ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:)
1138	ENDIF
1139	IF( jpj1 == 2 ) THEN
1140	ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:)
1141	ENDIF
1142	IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
1143	ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:)
1144	ENDIF
1145
1146	!! if data grid is cyclic we can do better on east-west edges
1147	!! but have to allow for whether first and last columns are coincident
1148	IF( ref_wgts(kw)%cyclic ) THEN
1149	rec1(2) = MAX( jpjmin-1, 1 )
1150	recn(1) = 1
1151	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1152	jpj1 = 2 + rec1(2) - jpjmin
1153	jpj2 = jpj1 + recn(2) - 1
1154	IF( jpi1 == 2 ) THEN
1155	rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap
1156	SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1157	CASE(1)
1158	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1159	CASE DEFAULT
1160	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1161	END SELECT
1162	ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1163	ENDIF
1164	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1165	rec1(1) = 1 + ref_wgts(kw)%overlap
1166	SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
1167	CASE(1)
1168	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
1169	CASE DEFAULT
1170	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1171	END SELECT
1172	ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1173	ENDIF
1174	ENDIF
1175
1176	! gradient in the i direction
1177	DO jk = 1,4
1178	DO jn = 1, jpj
1179	DO jm = 1,jpi
1180	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1181	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1182	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 * &
1183	(ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:))
1184	END DO
1185	END DO
1186	END DO
1187
1188	! gradient in the j direction
1189	DO jk = 1,4
1190	DO jn = 1, jpj
1191	DO jm = 1,jpi
1192	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1193	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1194	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 * &
1195	(ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:))
1196	END DO
1197	END DO
1198	END DO
1199
1200	! gradient in the ij direction
1201	DO jk = 1,4
1202	DO jn = 1, jpj
1203	DO jm = 1,jpi
1204	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1205	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1206	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
1207	(ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni ,nj+2,:)) - &
1208	(ref_wgts(kw)%fly_dta(ni+2,nj ,:) - ref_wgts(kw)%fly_dta(ni ,nj ,:)))
1209	END DO
1210	END DO
1211	END DO
1212	!
1213	END IF
1214	!
1215	END SUBROUTINE fld_interp
1216
1217
1218	FUNCTION ksec_week( cdday )
1219	!!---------------------------------------------------------------------
1220	!! * FUNCTION kshift_week *
1221	!!
1222	!! ** Purpose :
1223	!!---------------------------------------------------------------------
1224	CHARACTER(len=*), INTENT(in) :: cdday !3 first letters of the first day of the weekly file
1225	!!
1226	INTEGER :: ksec_week ! output variable
1227	INTEGER :: ijul !temp variable
1228	INTEGER :: ishift !temp variable
1229	CHARACTER(len=3),DIMENSION(7) :: cl_week
1230	!!----------------------------------------------------------------------
1231	cl_week = (/"sun","sat","fri","thu","wed","tue","mon"/)
1232	DO ijul = 1, 7
1233	IF( cl_week(ijul) == TRIM(cdday) ) EXIT
1234	END DO
1235	IF( ijul .GT. 7 ) CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): '//TRIM(cdday) )
1236	!
1237	ishift = ijul * NINT(rday)
1238	!
1239	ksec_week = nsec_week + ishift
1240	ksec_week = MOD( ksec_week, 7*NINT(rday) )
1241	!
1242	END FUNCTION ksec_week
1243
1244	!!======================================================================
1245	END MODULE fldread

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

Download in other formats: