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

Context Navigation

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 2633

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

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

Download in other formats: