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 trunk/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago
Update NEMOGCM from branch nemo_v3_3_beta
Property svn:keywords set to `Id`
File size: 61.5 KB

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

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

Download in other formats: