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

Context Navigation

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 3152

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

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

Download in other formats: