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 @ 2818

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

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

Download in other formats: