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

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

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

Context Navigation

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

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

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

Download in other formats: