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 NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/fldread.F90 @ 12256

Last change on this file since 12256 was 12256, checked in by smasson, 4 years ago
rev12232_dev_r12072_MERGE_OPTION2_2019: bugfix in fldread offset in bdy, see #2263. all sette tests ok
Property svn:keywords set to `Id`
File size: 87.3 KB

Line
1	MODULE fldread
2	!!======================================================================
3	!! * MODULE fldread *
4	!! Ocean forcing: read input field for surface boundary condition
5	!!=====================================================================
6	!! History : 2.0 ! 2006-06 (S. Masson, G. Madec) Original code
7	!! 3.0 ! 2008-05 (S. Alderson) Modified for Interpolation in memory from input grid to model grid
8	!! 3.4 ! 2013-10 (D. Delrosso, P. Oddo) suppression of land point prior to interpolation
9	!! ! 12-2015 (J. Harle) Adding BDY on-the-fly interpolation
10	!!----------------------------------------------------------------------
11
12	!!----------------------------------------------------------------------
13	!! fld_read : read input fields used for the computation of the surface boundary condition
14	!! fld_init : initialization of field read
15	!! fld_def : define the record(s) of the file and its name
16	!! fld_get : read the data
17	!! fld_map : read global data from file and map onto local data using a general mapping (use for open boundaries)
18	!! fld_rot : rotate the vector fields onto the local grid direction
19	!! fld_clopn : close/open the files
20	!! fld_fill : fill the data structure with the associated information read in namelist
21	!! wgt_list : manage the weights used for interpolation
22	!! wgt_print : print the list of known weights
23	!! fld_weight : create a WGT structure and fill in data from file, restructuring as required
24	!! apply_seaoverland : fill land with ocean values
25	!! seaoverland : create shifted matrices for seaoverland application
26	!! fld_interp : apply weights to input gridded data to create data on model grid
27	!! fld_filename : define the filename according to a given date
28	!! ksec_week : function returning seconds between 00h of the beginning of the week and half of the current time step
29	!!----------------------------------------------------------------------
30	USE oce ! ocean dynamics and tracers
31	USE dom_oce ! ocean space and time domain
32	USE phycst ! physical constant
33	USE sbc_oce ! surface boundary conditions : fields
34	USE geo2ocean ! for vector rotation on to model grid
35	!
36	USE in_out_manager ! I/O manager
37	USE iom ! I/O manager library
38	USE ioipsl , ONLY : ymds2ju, ju2ymds ! for calendar
39	USE lib_mpp ! MPP library
40	USE lbclnk ! ocean lateral boundary conditions (C1D case)
41
42	IMPLICIT NONE
43	PRIVATE
44
45	PUBLIC fld_map ! routine called by tides_init
46	PUBLIC fld_read, fld_fill ! called by sbc... modules
47	PUBLIC fld_def
48
49	TYPE, PUBLIC :: FLD_N !: Namelist field informations
50	CHARACTER(len = 256) :: clname ! generic name of the NetCDF flux file
51	REAL(wp) :: freqh ! frequency of each flux file
52	CHARACTER(len = 34) :: clvar ! generic name of the variable in the NetCDF flux file
53	LOGICAL :: ln_tint ! time interpolation or not (T/F)
54	LOGICAL :: ln_clim ! climatology or not (T/F)
55	CHARACTER(len = 8) :: cltype ! type of data file 'daily', 'monthly' or yearly'
56	CHARACTER(len = 256) :: wname ! generic name of a NetCDF weights file to be used, blank if not
57	CHARACTER(len = 34) :: vcomp ! symbolic component name if a vector that needs rotation
58	! ! a string starting with "U" or "V" for each component
59	! ! chars 2 onwards identify which components go together
60	CHARACTER(len = 34) :: lname ! generic name of a NetCDF land/sea mask file to be used, blank if not
61	! ! 0=sea 1=land
62	END TYPE FLD_N
63
64	TYPE, PUBLIC :: FLD !: Input field related variables
65	CHARACTER(len = 256) :: clrootname ! generic name of the NetCDF file
66	CHARACTER(len = 256) :: clname ! current name of the NetCDF file
67	REAL(wp) :: freqh ! frequency of each flux file
68	CHARACTER(len = 34) :: clvar ! generic name of the variable in the NetCDF flux file
69	LOGICAL :: ln_tint ! time interpolation or not (T/F)
70	LOGICAL :: ln_clim ! climatology or not (T/F)
71	CHARACTER(len = 8) :: cltype ! type of data file 'daily', 'monthly' or yearly'
72	INTEGER :: num ! iom id of the jpfld files to be read
73	INTEGER , DIMENSION(2) :: nrec_b ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year)
74	INTEGER , DIMENSION(2) :: nrec_a ! after record (1: index, 2: second since Jan. 1st 00h of nit000 year)
75	INTEGER , ALLOCATABLE, DIMENSION(: ) :: nrecsec !
76	REAL(wp), ALLOCATABLE, DIMENSION(:,:,: ) :: fnow ! input fields interpolated to now time step
77	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: fdta ! 2 consecutive record of input fields
78	CHARACTER(len = 256) :: wgtname ! current name of the NetCDF weight file acting as a key
79	! ! into the WGTLIST structure
80	CHARACTER(len = 34) :: vcomp ! symbolic name for a vector component that needs rotation
81	LOGICAL, DIMENSION(2) :: rotn ! flag to indicate whether before/after field has been rotated
82	INTEGER :: nreclast ! last record to be read in the current file
83	CHARACTER(len = 256) :: lsmname ! current name of the NetCDF mask file acting as a key
84	! !
85	! ! Variables related to BDY
86	INTEGER :: igrd ! grid type for bdy data
87	INTEGER :: ibdy ! bdy set id number
88	INTEGER, POINTER, DIMENSION(:) :: imap ! Array of integer pointers to 1D arrays
89	LOGICAL :: ltotvel ! total velocity or not (T/F)
90	LOGICAL :: lzint ! T if it requires a vertical interpolation
91	END TYPE FLD
92
93	!$AGRIF_DO_NOT_TREAT
94
95	!! keep list of all weights variables so they're only read in once
96	!! need to add AGRIF directives not to process this structure
97	!! also need to force wgtname to include AGRIF nest number
98	TYPE :: WGT !: Input weights related variables
99	CHARACTER(len = 256) :: wgtname ! current name of the NetCDF weight file
100	INTEGER , DIMENSION(2) :: ddims ! shape of input grid
101	INTEGER , DIMENSION(2) :: botleft ! top left corner of box in input grid containing
102	! ! current processor grid
103	INTEGER , DIMENSION(2) :: topright ! top right corner of box
104	INTEGER :: jpiwgt ! width of box on input grid
105	INTEGER :: jpjwgt ! height of box on input grid
106	INTEGER :: numwgt ! number of weights (4=bilinear, 16=bicubic)
107	INTEGER :: nestid ! for agrif, keep track of nest we're in
108	INTEGER :: overlap ! =0 when cyclic grid has no overlapping EW columns
109	! ! =>1 when they have one or more overlapping columns
110	! ! =-1 not cyclic
111	LOGICAL :: cyclic ! east-west cyclic or not
112	INTEGER, DIMENSION(:,:,:), POINTER :: data_jpi ! array of source integers
113	INTEGER, DIMENSION(:,:,:), POINTER :: data_jpj ! array of source integers
114	REAL(wp), DIMENSION(:,:,:), POINTER :: data_wgt ! array of weights on model grid
115	REAL(wp), DIMENSION(:,:,:), POINTER :: fly_dta ! array of values on input grid
116	REAL(wp), DIMENSION(:,:,:), POINTER :: col ! temporary array for reading in columns
117	END TYPE WGT
118
119	INTEGER, PARAMETER :: tot_wgts = 20
120	TYPE( WGT ), DIMENSION(tot_wgts) :: ref_wgts ! array of wgts
121	INTEGER :: nxt_wgt = 1 ! point to next available space in ref_wgts array
122	INTEGER :: nflag = 0
123	REAL(wp), PARAMETER :: undeff_lsm = -999.00_wp
124
125	!$AGRIF_END_DO_NOT_TREAT
126
127	!!----------------------------------------------------------------------
128	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
129	!! $Id$
130	!! Software governed by the CeCILL license (see ./LICENSE)
131	!!----------------------------------------------------------------------
132	CONTAINS
133
134	SUBROUTINE fld_read( kt, kn_fsbc, sd, kit, pt_offset )
135	!!---------------------------------------------------------------------
136	!! * ROUTINE fld_read *
137	!!
138	!! ** Purpose : provide at each time step the surface ocean fluxes
139	!! (momentum, heat, freshwater and runoff)
140	!!
141	!! ** Method : READ each input fields in NetCDF files using IOM
142	!! and intepolate it to the model time-step.
143	!! Several assumptions are made on the input file:
144	!! blahblahblah....
145	!!----------------------------------------------------------------------
146	INTEGER , INTENT(in ) :: kt ! ocean time step
147	INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step)
148	TYPE(FLD), INTENT(inout), DIMENSION(:) :: sd ! input field related variables
149	INTEGER , INTENT(in ), OPTIONAL :: kit ! subcycle timestep for timesplitting option
150	REAL(wp) , INTENT(in ), OPTIONAL :: pt_offset ! provide fields at time other than "now"
151	!!
152	INTEGER :: imf ! size of the structure sd
153	INTEGER :: jf ! dummy indices
154	INTEGER :: isecsbc ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step
155	LOGICAL :: ll_firstcall ! true if this is the first call to fld_read for this set of fields
156	REAL(wp) :: zt_offset ! local time offset variable
157	REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation
158	REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation
159	CHARACTER(LEN=1000) :: clfmt ! write format
160	!!---------------------------------------------------------------------
161	ll_firstcall = kt == nit000
162	IF( PRESENT(kit) ) ll_firstcall = ll_firstcall .and. kit == 1
163
164	IF( nn_components == jp_iam_sas ) THEN ; zt_offset = REAL( nn_fsbc, wp )
165	ELSE ; zt_offset = 0.
166	ENDIF
167	IF( PRESENT(pt_offset) ) zt_offset = pt_offset
168
169	! Note that all varibles starting by nsec_* are shifted time by +1/2 time step to be centrered
170	IF( PRESENT(kit) ) THEN ! ignore kn_fsbc in this case
171	isecsbc = nsec_year + nsec1jan000 + NINT( ( REAL( kit,wp) + zt_offset ) * rdt / REAL(nn_baro,wp) )
172	ELSE ! middle of sbc time step
173	! note: we use kn_fsbc-1 because nsec_year is defined at the middle of the current time step
174	isecsbc = nsec_year + nsec1jan000 + NINT( ( 0.5REAL(kn_fsbc-1,wp) + zt_offset ) rdt )
175	ENDIF
176	imf = SIZE( sd )
177	!
178	IF( ll_firstcall ) THEN ! initialization
179	DO jf = 1, imf
180	IF( TRIM(sd(jf)%clrootname) == 'NOT USED' ) CYCLE
181	CALL fld_init( isecsbc, sd(jf) ) ! read each before field (put them in after as they will be swapped)
182	END DO
183	IF( lwp ) CALL wgt_print() ! control print
184	ENDIF
185	! ! ====================================== !
186	IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN ! update field at each kn_fsbc time-step !
187	! ! ====================================== !
188	!
189	DO jf = 1, imf ! --- loop over field --- !
190	!
191	IF( TRIM(sd(jf)%clrootname) == 'NOT USED' ) CYCLE
192	CALL fld_update( isecsbc, sd(jf) )
193	!
194	END DO ! --- end loop over field --- !
195
196	CALL fld_rot( kt, sd ) ! rotate vector before/now/after fields if needed
197
198	DO jf = 1, imf ! --- loop over field --- !
199	!
200	IF( TRIM(sd(jf)%clrootname) == 'NOT USED' ) CYCLE
201	!
202	IF( sd(jf)%ln_tint ) THEN ! temporal interpolation
203	IF(lwp .AND. kt - nit000 <= 100 ) THEN
204	clfmt = "(' fld_read: var ', a, ' kt = ', i8, ' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," // &
205	& "', records b/a: ', i6.4, '/', i6.4, ' (days ', f9.4,'/', f9.4, ')')"
206	WRITE(numout, clfmt) TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday, &
207	& 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
208	WRITE(numout, *) ' zt_offset is : ',zt_offset
209	ENDIF
210	! temporal interpolation weights
211	ztinta = REAL( isecsbc - sd(jf)%nrec_b(2), wp ) / REAL( sd(jf)%nrec_a(2) - sd(jf)%nrec_b(2), wp )
212	ztintb = 1. - ztinta
213	sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:,1) + ztinta * sd(jf)%fdta(:,:,:,2)
214	ELSE ! nothing to do...
215	IF(lwp .AND. kt - nit000 <= 100 ) THEN
216	clfmt = "(' fld_read: var ', a, ' kt = ', i8,' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," // &
217	& "', record: ', i6.4, ' (days ', f9.4, ' <-> ', f9.4, ')')"
218	WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday, &
219	& sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday
220	ENDIF
221	ENDIF
222	!
223	IF( kt == nitend - kn_fsbc + 1 ) CALL iom_close( sd(jf)%num ) ! Close the input files
224
225	END DO ! --- end loop over field --- !
226	!
227	ENDIF
228	!
229	END SUBROUTINE fld_read
230
231
232	SUBROUTINE fld_init( ksecsbc, sdjf )
233	!!---------------------------------------------------------------------
234	!! * ROUTINE fld_init *
235	!!
236	!! ** Purpose : - first call(s) to fld_def to define before values
237	!! - open file
238	!!----------------------------------------------------------------------
239	INTEGER , INTENT(in ) :: ksecsbc !
240	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
241	!!---------------------------------------------------------------------
242	!
243	IF( nflag == 0 ) nflag = -( HUGE(0) - 10 )
244	!
245	CALL fld_def( sdjf )
246	IF( sdjf%ln_tint .AND. ksecsbc < sdjf%nrecsec(1) ) CALL fld_def( sdjf, ldprev = .TRUE. )
247	!
248	CALL fld_clopn( sdjf )
249	sdjf%nrec_a(:) = (/ 1, nflag /) ! default definition to force flp_update to read the file.
250	!
251	END SUBROUTINE fld_init
252
253
254	SUBROUTINE fld_update( ksecsbc, sdjf )
255	!!---------------------------------------------------------------------
256	!! * ROUTINE fld_update *
257	!!
258	!! ** Purpose : Compute
259	!! if sdjf%ln_tint = .TRUE.
260	!! nrec_a: record number and its time (nrec_b is obtained from nrec_a when swapping)
261	!! if sdjf%ln_tint = .FALSE.
262	!! nrec_a(1): record number
263	!! nrec_b(2) and nrec_a(2): time of the beginning and end of the record
264	!!----------------------------------------------------------------------
265	INTEGER , INTENT(in ) :: ksecsbc !
266	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
267	!
268	INTEGER :: ja ! end of this record (in seconds)
269	!!----------------------------------------------------------------------
270	!
271	IF( ksecsbc > sdjf%nrec_a(2) ) THEN ! --> we need to update after data
272
273	! find where is the new after record... (it is not necessary sdjf%nrec_a(1)+1 )
274	ja = sdjf%nrec_a(1)
275	DO WHILE ( ksecsbc >= sdjf%nrecsec(ja) .AND. ja < sdjf%nreclast ) ! Warning: make sure ja <= sdjf%nreclast in this test
276	ja = ja + 1
277	END DO
278	IF( ksecsbc > sdjf%nrecsec(ja) ) ja = ja + 1 ! in case ksecsbc > sdjf%nrecsec(sdjf%nreclast)
279
280	! if ln_tint and if the new after is not ja+1, we need also to update after data before the swap
281	! so, after the swap, sdjf%nrec_b(2) will still be the closest value located just before ksecsbc
282	IF( sdjf%ln_tint .AND. ( ja > sdjf%nrec_a(1) + 1 .OR. sdjf%nrec_a(2) == nflag ) ) THEN
283	sdjf%nrec_a(:) = (/ ja-1, sdjf%nrecsec(ja-1) /) ! update nrec_a with before information
284	CALL fld_get( sdjf ) ! read after data that will be used as before data
285	ENDIF
286
287	! if after is in the next file...
288	IF( ja > sdjf%nreclast ) THEN
289
290	CALL fld_def( sdjf )
291	IF( ksecsbc > sdjf%nrecsec(sdjf%nreclast) ) CALL fld_def( sdjf, ldnext = .TRUE. )
292	CALL fld_clopn( sdjf ) ! open next file
293
294	! find where is after in this new file
295	ja = 1
296	DO WHILE ( ksecsbc > sdjf%nrecsec(ja) .AND. ja < sdjf%nreclast )
297	ja = ja + 1
298	END DO
299	IF( ksecsbc > sdjf%nrecsec(ja) ) ja = ja + 1 ! in case ksecsbc > sdjf%nrecsec(sdjf%nreclast)
300
301	IF( ja > sdjf%nreclast ) THEN
302	CALL ctl_stop( "STOP", "fld_def: need next-next file? we should not be there... file: "//TRIM(sdjf%clrootname) )
303	ENDIF
304
305	! if ln_tint and if after is not the first record, we must (potentially again) update after data before the swap
306	IF( sdjf%ln_tint .AND. ja > 1 ) THEN
307	IF( sdjf%nrecsec(0) /= nflag ) THEN ! no trick used: after file is not the current file
308	sdjf%nrec_a(:) = (/ ja-1, sdjf%nrecsec(ja-1) /) ! update nrec_a with before information
309	CALL fld_get( sdjf ) ! read after data that will be used as before data
310	ENDIF
311	ENDIF
312
313	ENDIF
314
315	IF( sdjf%ln_tint ) THEN
316	! Swap data
317	sdjf%nrec_b(:) = sdjf%nrec_a(:) ! swap before record informations
318	sdjf%rotn(1) = sdjf%rotn(2) ! swap before rotate informations
319	sdjf%fdta(:,:,:,1) = sdjf%fdta(:,:,:,2) ! swap before record field
320	ELSE
321	sdjf%nrec_b(:) = (/ ja-1, sdjf%nrecsec(ja-1) /) ! only for print
322	ENDIF
323
324	! read new after data
325	sdjf%nrec_a(:) = (/ ja, sdjf%nrecsec(ja) /) ! update nrec_a as it is used by fld_get
326	CALL fld_get( sdjf ) ! read after data (with nrec_a informations)
327
328	ENDIF
329	!
330	END SUBROUTINE fld_update
331
332
333	SUBROUTINE fld_get( sdjf )
334	!!---------------------------------------------------------------------
335	!! * ROUTINE fld_get *
336	!!
337	!! ** Purpose : read the data
338	!!----------------------------------------------------------------------
339	TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables
340	!
341	INTEGER :: ipk ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
342	INTEGER :: iw ! index into wgts array
343	INTEGER :: ipdom ! index of the domain
344	INTEGER :: idvar ! variable ID
345	INTEGER :: idmspc ! number of spatial dimensions
346	LOGICAL :: lmoor ! C1D case: point data
347	!!---------------------------------------------------------------------
348	!
349	ipk = SIZE( sdjf%fnow, 3 )
350	!
351	IF( ASSOCIATED(sdjf%imap) ) THEN
352	IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), &
353	& sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel, sdjf%lzint )
354	ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1), &
355	& sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel, sdjf%lzint )
356	ENDIF
357	ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
358	CALL wgt_list( sdjf, iw )
359	IF( sdjf%ln_tint ) THEN ; CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fdta(:,:,:,2), &
360	& sdjf%nrec_a(1), sdjf%lsmname )
361	ELSE ; CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fnow(:,:,: ), &
362	& sdjf%nrec_a(1), sdjf%lsmname )
363	ENDIF
364	ELSE
365	IF( SIZE(sdjf%fnow, 1) == jpi ) THEN ; ipdom = jpdom_data
366	ELSE ; ipdom = jpdom_unknown
367	ENDIF
368	! C1D case: If product of spatial dimensions == ipk, then x,y are of
369	! size 1 (point/mooring data): this must be read onto the central grid point
370	idvar = iom_varid( sdjf%num, sdjf%clvar )
371	idmspc = iom_file ( sdjf%num )%ndims( idvar )
372	IF( iom_file( sdjf%num )%luld( idvar ) ) idmspc = idmspc - 1
373	lmoor = ( idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk )
374	!
375	SELECT CASE( ipk )
376	CASE(1)
377	IF( lk_c1d .AND. lmoor ) THEN
378	IF( sdjf%ln_tint ) THEN
379	CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fdta(2,2,1,2), sdjf%nrec_a(1) )
380	CALL lbc_lnk( 'fldread', sdjf%fdta(:,:,1,2),'Z',1. )
381	ELSE
382	CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fnow(2,2,1 ), sdjf%nrec_a(1) )
383	CALL lbc_lnk( 'fldread', sdjf%fnow(:,:,1 ),'Z',1. )
384	ENDIF
385	ELSE
386	IF( sdjf%ln_tint ) THEN ; CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,1,2), sdjf%nrec_a(1) )
387	ELSE ; CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,1 ), sdjf%nrec_a(1) )
388	ENDIF
389	ENDIF
390	CASE DEFAULT
391	IF(lk_c1d .AND. lmoor ) THEN
392	IF( sdjf%ln_tint ) THEN
393	CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fdta(2,2,:,2), sdjf%nrec_a(1) )
394	CALL lbc_lnk( 'fldread', sdjf%fdta(:,:,:,2),'Z',1. )
395	ELSE
396	CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fnow(2,2,: ), sdjf%nrec_a(1) )
397	CALL lbc_lnk( 'fldread', sdjf%fnow(:,:,: ),'Z',1. )
398	ENDIF
399	ELSE
400	IF( sdjf%ln_tint ) THEN ; CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) )
401	ELSE ; CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1) )
402	ENDIF
403	ENDIF
404	END SELECT
405	ENDIF
406	!
407	sdjf%rotn(2) = .false. ! vector not yet rotated
408	!
409	END SUBROUTINE fld_get
410
411
412	SUBROUTINE fld_map( knum, cdvar, pdta, krec, kmap, kgrd, kbdy, ldtotvel, ldzint )
413	!!---------------------------------------------------------------------
414	!! * ROUTINE fld_map *
415	!!
416	!! ** Purpose : read global data from file and map onto local data
417	!! using a general mapping (for open boundaries)
418	!!----------------------------------------------------------------------
419	INTEGER , INTENT(in ) :: knum ! stream number
420	CHARACTER(LEN=*) , INTENT(in ) :: cdvar ! variable name
421	REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pdta ! bdy output field on model grid
422	INTEGER , INTENT(in ) :: krec ! record number to read (ie time slice)
423	INTEGER , DIMENSION(:) , INTENT(in ) :: kmap ! global-to-local bdy mapping indices
424	! optional variables used for vertical interpolation:
425	INTEGER, OPTIONAL , INTENT(in ) :: kgrd ! grid type (t, u, v)
426	INTEGER, OPTIONAL , INTENT(in ) :: kbdy ! bdy number
427	LOGICAL, OPTIONAL , INTENT(in ) :: ldtotvel ! true if total ( = barotrop + barocline) velocity
428	LOGICAL, OPTIONAL , INTENT(in ) :: ldzint ! true if 3D variable requires a vertical interpolation
429	!!
430	INTEGER :: ipi ! length of boundary data on local process
431	INTEGER :: ipj ! length of dummy dimension ( = 1 )
432	INTEGER :: ipk ! number of vertical levels of pdta ( 2D: ipk=1 ; 3D: ipk=jpk )
433	INTEGER :: ipkb ! number of vertical levels in boundary data file
434	INTEGER :: idvar ! variable ID
435	INTEGER :: indims ! number of dimensions of the variable
436	INTEGER, DIMENSION(4) :: idimsz ! size of variable dimensions
437	REAL(wp) :: zfv ! fillvalue
438	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zz_read ! work space for global boundary data
439	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zdta_read ! work space local data requiring vertical interpolation
440	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zdta_read_z ! work space local data requiring vertical interpolation
441	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zdta_read_dz ! work space local data requiring vertical interpolation
442	CHARACTER(LEN=1),DIMENSION(3) :: clgrid
443	LOGICAL :: lluld ! is the variable using the unlimited dimension
444	LOGICAL :: llzint ! local value of ldzint
445	!!---------------------------------------------------------------------
446	!
447	clgrid = (/'t','u','v'/)
448	!
449	ipi = SIZE( pdta, 1 )
450	ipj = SIZE( pdta, 2 ) ! must be equal to 1
451	ipk = SIZE( pdta, 3 )
452	!
453	llzint = .FALSE.
454	IF( PRESENT(ldzint) ) llzint = ldzint
455	!
456	idvar = iom_varid( knum, cdvar, kndims = indims, kdimsz = idimsz, lduld = lluld )
457	IF( indims == 4 .OR. ( indims == 3 .AND. .NOT. lluld ) ) THEN ; ipkb = idimsz(3) ! xy(zl)t or xy(zl)
458	ELSE ; ipkb = 1 ! xy or xyt
459	ENDIF
460	!
461	ALLOCATE( zz_read( idimsz(1), idimsz(2), ipkb ) ) ! ++++++++ !!! this can be very big...
462	!
463	IF( ipk == 1 ) THEN
464
465	IF( ipkb /= 1 ) CALL ctl_stop( 'fld_map : we must have ipkb = 1 to read surface data' )
466	CALL iom_get ( knum, jpdom_unknown, cdvar, zz_read(:,:,1), krec ) ! call iom_get with a 2D file
467	CALL fld_map_core( zz_read, kmap, pdta )
468
469	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
470	! Do we include something here to adjust barotropic velocities !
471	! in case of a depth difference between bdy files and !
472	! bathymetry in the case ln_totvel = .false. and ipkb>0? !
473	! [as the enveloping and parital cells could change H] !
474	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
475
476	ELSE
477	!
478	CALL iom_get ( knum, jpdom_unknown, cdvar, zz_read(:,:,:), krec ) ! call iom_get with a 3D file
479	!
480	IF( ipkb /= ipk .OR. llzint ) THEN ! boundary data not on model vertical grid : vertical interpolation
481	!
482	IF( ipk == jpk .AND. iom_varid(knum,'gdep'//clgrid(kgrd)) /= -1 .AND. iom_varid(knum,'e3'//clgrid(kgrd)) /= -1 ) THEN
483
484	ALLOCATE( zdta_read(ipi,ipj,ipkb), zdta_read_z(ipi,ipj,ipkb), zdta_read_dz(ipi,ipj,ipkb) )
485
486	CALL fld_map_core( zz_read, kmap, zdta_read )
487	CALL iom_get ( knum, jpdom_unknown, 'gdep'//clgrid(kgrd), zz_read ) ! read only once? Potential temporal evolution?
488	CALL fld_map_core( zz_read, kmap, zdta_read_z )
489	CALL iom_get ( knum, jpdom_unknown, 'e3'//clgrid(kgrd), zz_read ) ! read only once? Potential temporal evolution?
490	CALL fld_map_core( zz_read, kmap, zdta_read_dz )
491
492	CALL iom_getatt(knum, '_FillValue', zfv, cdvar=cdvar )
493	CALL fld_bdy_interp(zdta_read, zdta_read_z, zdta_read_dz, pdta, kgrd, kbdy, zfv, ldtotvel)
494	DEALLOCATE( zdta_read, zdta_read_z, zdta_read_dz )
495
496	ELSE
497	IF( ipk /= jpk ) CALL ctl_stop( 'fld_map : this should be an impossible case...' )
498	WRITE(ctmp1,*) 'fld_map : vertical interpolation for bdy variable '//TRIM(cdvar)//' requires '
499	IF( iom_varid(knum, 'gdep'//clgrid(kgrd)) == -1 ) CALL ctl_stop( ctmp1//'gdep'//clgrid(kgrd)//' variable' )
500	IF( iom_varid(knum, 'e3'//clgrid(kgrd)) == -1 ) CALL ctl_stop( ctmp1// 'e3'//clgrid(kgrd)//' variable' )
501
502	ENDIF
503	!
504	ELSE ! bdy data assumed to be the same levels as bdy variables
505	!
506	CALL fld_map_core( zz_read, kmap, pdta )
507	!
508	ENDIF ! ipkb /= ipk
509	ENDIF ! ipk == 1
510
511	DEALLOCATE( zz_read )
512
513	END SUBROUTINE fld_map
514
515
516	SUBROUTINE fld_map_core( pdta_read, kmap, pdta_bdy )
517	!!---------------------------------------------------------------------
518	!! * ROUTINE fld_map_core *
519	!!
520	!! ** Purpose : inner core of fld_map
521	!!----------------------------------------------------------------------
522	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdta_read ! global boundary data
523	INTEGER, DIMENSION(: ), INTENT(in ) :: kmap ! global-to-local bdy mapping indices
524	REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pdta_bdy ! bdy output field on model grid
525	!!
526	INTEGER, DIMENSION(3) :: idim_read, idim_bdy ! arrays dimensions
527	INTEGER :: ji, jj, jk, jb ! loop counters
528	INTEGER :: im1
529	!!---------------------------------------------------------------------
530	!
531	idim_read = SHAPE( pdta_read )
532	idim_bdy = SHAPE( pdta_bdy )
533	!
534	! in all cases: idim_bdy(2) == 1 .AND. idim_read(1) * idim_read(2) == idim_bdy(1)
535	! structured BDY with rimwidth > 1 : idim_read(2) == rimwidth /= 1
536	! structured BDY with rimwidth == 1 or unstructured BDY: idim_read(2) == 1
537	!
538	IF( idim_read(2) > 1 ) THEN ! structured BDY with rimwidth > 1
539	DO jk = 1, idim_bdy(3)
540	DO jb = 1, idim_bdy(1)
541	im1 = kmap(jb) - 1
542	jj = im1 / idim_read(1) + 1
543	ji = MOD( im1, idim_read(1) ) + 1
544	pdta_bdy(jb,1,jk) = pdta_read(ji,jj,jk)
545	END DO
546	END DO
547	ELSE
548	DO jk = 1, idim_bdy(3)
549	DO jb = 1, idim_bdy(1) ! horizontal remap of bdy data on the local bdy
550	pdta_bdy(jb,1,jk) = pdta_read(kmap(jb),1,jk)
551	END DO
552	END DO
553	ENDIF
554
555	END SUBROUTINE fld_map_core
556
557
558	SUBROUTINE fld_bdy_interp(pdta_read, pdta_read_z, pdta_read_dz, pdta, kgrd, kbdy, pfv, ldtotvel)
559	!!---------------------------------------------------------------------
560	!! * ROUTINE fld_bdy_interp *
561	!!
562	!! ** Purpose : on the fly vertical interpolation to allow the use of
563	!! boundary data from non-native vertical grid
564	!!----------------------------------------------------------------------
565	USE bdy_oce, ONLY: idx_bdy ! indexing for map <-> ij transformation
566
567	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdta_read ! data read in bdy file
568	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdta_read_z ! depth of the data read in bdy file
569	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdta_read_dz ! thickness of the levels in bdy file
570	REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pdta ! output field on model grid (2 dimensional)
571	REAL(wp) , INTENT(in ) :: pfv ! fillvalue of the data read in bdy file
572	LOGICAL , INTENT(in ) :: ldtotvel ! true if toal ( = barotrop + barocline) velocity
573	INTEGER , INTENT(in ) :: kgrd ! grid type (t, u, v)
574	INTEGER , INTENT(in ) :: kbdy ! bdy number
575	!!
576	INTEGER :: ipi ! length of boundary data on local process
577	INTEGER :: ipkb ! number of vertical levels in boundary data file
578	INTEGER :: jb, ji, jj, jk, jkb ! loop counters
579	REAL(wp) :: zcoef
580	REAL(wp) :: zl, zi, zh ! tmp variable for current depth and interpolation factor
581	REAL(wp) :: zfv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(pfv)
582	REAL(wp), DIMENSION(jpk) :: zdepth, zdhalf ! level and half-level depth
583	!!---------------------------------------------------------------------
584
585	ipi = SIZE( pdta, 1 )
586	ipkb = SIZE( pdta_read, 3 )
587
588	zfv_alt = -ABS(pfv) ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later
589	!
590	WHERE( pdta_read == pfv )
591	pdta_read_z = zfv_alt ! safety: put fillvalue into external depth field so consistent with data
592	pdta_read_dz = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct
593	ENDWHERE
594
595	DO jb = 1, ipi
596	ji = idx_bdy(kbdy)%nbi(jb,kgrd)
597	jj = idx_bdy(kbdy)%nbj(jb,kgrd)
598	zh = SUM(pdta_read_dz(jb,1,:) )
599	!
600	SELECT CASE( kgrd )
601	CASE(1)
602	! depth of T points:
603	zdepth(:) = gdept_n(ji,jj,:)
604	CASE(2)
605	! depth of U points: we must not use gdept_n as we don't want to do a communication
606	! --> copy what is done for gdept_n in domvvl...
607	zdhalf(1) = 0.0_wp
608	zdepth(1) = 0.5_wp * e3uw_n(ji,jj,1)
609	DO jk = 2, jpk ! vertical sum
610	! zcoef = umask - wumask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
611	! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
612	! ! 0.5 where jk = mikt
613	!!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ??
614	zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) )
615	zdhalf(jk) = zdhalf(jk-1) + e3u_n(ji,jj,jk-1)
616	zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5 * e3uw_n(ji,jj,jk)) &
617	& + (1-zcoef) * ( zdepth(jk-1) + e3uw_n(ji,jj,jk))
618	END DO
619	CASE(3)
620	! depth of V points: we must not use gdept_n as we don't want to do a communication
621	! --> copy what is done for gdept_n in domvvl...
622	zdhalf(1) = 0.0_wp
623	zdepth(1) = 0.5_wp * e3vw_n(ji,jj,1)
624	DO jk = 2, jpk ! vertical sum
625	! zcoef = vmask - wvmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
626	! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
627	! ! 0.5 where jk = mikt
628	!!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ??
629	zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) )
630	zdhalf(jk) = zdhalf(jk-1) + e3v_n(ji,jj,jk-1)
631	zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5 * e3vw_n(ji,jj,jk)) &
632	& + (1-zcoef) * ( zdepth(jk-1) + e3vw_n(ji,jj,jk))
633	END DO
634	END SELECT
635	!
636	DO jk = 1, jpk
637	IF( zdepth(jk) < pdta_read_z(jb,1, 1) ) THEN ! above the first level of external data
638	pdta(jb,1,jk) = pdta_read(jb,1,1)
639	ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkb) ) THEN ! below the last level of external data
640	pdta(jb,1,jk) = pdta_read(jb,1,MAXLOC(pdta_read_z(jb,1,:),1))
641	ELSE ! inbetween: vertical interpolation between jkb & jkb+1
642	DO jkb = 1, ipkb-1 ! when gdept_n(jkb) < zdepth(jk) < gdept_n(jkb+1)
643	IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) <= 0._wp ) &
644	& .AND. ( pdta_read_z(jb,1,jkb+1) /= zfv_alt) ) THEN ! linear interpolation between 2 levels
645	zi = ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) / ( pdta_read_z(jb,1,jkb+1) - pdta_read_z(jb,1,jkb) )
646	pdta(jb,1,jk) = pdta_read(jb,1,jkb) + ( pdta_read (jb,1,jkb+1) - pdta_read (jb,1,jkb) ) * zi
647	ENDIF
648	END DO
649	ENDIF
650	END DO ! jpk
651	!
652	END DO ! ipi
653
654	IF(kgrd == 2) THEN ! do we need to adjust the transport term?
655	DO jb = 1, ipi
656	ji = idx_bdy(kbdy)%nbi(jb,kgrd)
657	jj = idx_bdy(kbdy)%nbj(jb,kgrd)
658	zh = SUM(pdta_read_dz(jb,1,:) )
659	ztrans = 0._wp
660	ztrans_new = 0._wp
661	DO jkb = 1, ipkb ! calculate transport on input grid
662	ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)
663	ENDDO
664	DO jk = 1, jpk ! calculate transport on model grid
665	ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3u_n(ji,jj,jk ) * umask(ji,jj,jk)
666	ENDDO
667	DO jk = 1, jpk ! make transport correction
668	IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data
669	pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu_n(ji,jj) ) * umask(ji,jj,jk)
670	ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
671	pdta(jb,1,jk) = pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hu_n(ji,jj) * umask(ji,jj,jk)
672	ENDIF
673	ENDDO
674	ENDDO
675	ENDIF
676
677	IF(kgrd == 3) THEN ! do we need to adjust the transport term?
678	DO jb = 1, ipi
679	ji = idx_bdy(kbdy)%nbi(jb,kgrd)
680	jj = idx_bdy(kbdy)%nbj(jb,kgrd)
681	zh = SUM(pdta_read_dz(jb,1,:) )
682	ztrans = 0._wp
683	ztrans_new = 0._wp
684	DO jkb = 1, ipkb ! calculate transport on input grid
685	ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)
686	ENDDO
687	DO jk = 1, jpk ! calculate transport on model grid
688	ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3v_n(ji,jj,jk ) * vmask(ji,jj,jk)
689	ENDDO
690	DO jk = 1, jpk ! make transport correction
691	IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data
692	pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv_n(ji,jj) ) * vmask(ji,jj,jk)
693	ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
694	pdta(jb,1,jk) = pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hv_n(ji,jj) * vmask(ji,jj,jk)
695	ENDIF
696	ENDDO
697	ENDDO
698	ENDIF
699
700	END SUBROUTINE fld_bdy_interp
701
702
703	SUBROUTINE fld_rot( kt, sd )
704	!!---------------------------------------------------------------------
705	!! * ROUTINE fld_rot *
706	!!
707	!! ** Purpose : Vector fields may need to be rotated onto the local grid direction
708	!!----------------------------------------------------------------------
709	INTEGER , INTENT(in ) :: kt ! ocean time step
710	TYPE(FLD), DIMENSION(:), INTENT(inout) :: sd ! input field related variables
711	!
712	INTEGER :: ju, jv, jk, jn ! loop indices
713	INTEGER :: imf ! size of the structure sd
714	INTEGER :: ill ! character length
715	INTEGER :: iv ! indice of V component
716	CHARACTER (LEN=100) :: clcomp ! dummy weight name
717	REAL(wp), DIMENSION(jpi,jpj) :: utmp, vtmp ! temporary arrays for vector rotation
718	!!---------------------------------------------------------------------
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	IF( TRIM(sd(ju)%clrootname) == 'NOT USED' ) CYCLE
725	ill = LEN_TRIM( sd(ju)%vcomp )
726	DO jn = 2-COUNT((/sd(ju)%ln_tint/)), 2
727	IF( ill > 0 .AND. .NOT. sd(ju)%rotn(jn) ) THEN ! find vector rotations required
728	IF( sd(ju)%vcomp(1:1) == 'U' ) THEN ! east-west component has symbolic name starting with 'U'
729	! look for the north-south component which has same symbolic name but with 'U' replaced with 'V'
730	clcomp = 'V' // sd(ju)%vcomp(2:ill) ! works even if ill == 1
731	iv = -1
732	DO jv = 1, imf
733	IF( TRIM(sd(jv)%clrootname) == 'NOT USED' ) CYCLE
734	IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) ) iv = jv
735	END DO
736	IF( iv > 0 ) THEN ! fields ju and iv are two components which need to be rotated together
737	DO jk = 1, SIZE( sd(ju)%fnow, 3 )
738	IF( sd(ju)%ln_tint )THEN
739	CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->i', utmp(:,:) )
740	CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->j', vtmp(:,:) )
741	sd(ju)%fdta(:,:,jk,jn) = utmp(:,:) ; sd(iv)%fdta(:,:,jk,jn) = vtmp(:,:)
742	ELSE
743	CALL rot_rep( sd(ju)%fnow(:,:,jk ), sd(iv)%fnow(:,:,jk ), 'T', 'en->i', utmp(:,:) )
744	CALL rot_rep( sd(ju)%fnow(:,:,jk ), sd(iv)%fnow(:,:,jk ), 'T', 'en->j', vtmp(:,:) )
745	sd(ju)%fnow(:,:,jk ) = utmp(:,:) ; sd(iv)%fnow(:,:,jk ) = vtmp(:,:)
746	ENDIF
747	END DO
748	sd(ju)%rotn(jn) = .TRUE. ! vector was rotated
749	IF( lwp .AND. kt == nit000 ) WRITE(numout,*) &
750	& 'fld_read: vector pair ('//TRIM(sd(ju)%clvar)//', '//TRIM(sd(iv)%clvar)//') rotated on to model grid'
751	ENDIF
752	ENDIF
753	ENDIF
754	END DO
755	END DO
756	!
757	END SUBROUTINE fld_rot
758
759
760	SUBROUTINE fld_def( sdjf, ldprev, ldnext )
761	!!---------------------------------------------------------------------
762	!! * ROUTINE fld_def *
763	!!
764	!! ** Purpose : define the record(s) of the file and its name
765	!!----------------------------------------------------------------------
766	TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables
767	LOGICAL, OPTIONAL, INTENT(in ) :: ldprev !
768	LOGICAL, OPTIONAL, INTENT(in ) :: ldnext !
769	!
770	INTEGER :: jt
771	INTEGER :: idaysec ! number of seconds in 1 day = NINT(rday)
772	INTEGER :: iyr, imt, idy, isecwk
773	INTEGER :: indexyr, indexmt
774	INTEGER :: ireclast
775	INTEGER :: ishift, istart
776	INTEGER, DIMENSION(2) :: isave
777	REAL(wp) :: zfreqs
778	LOGICAL :: llprev, llnext, llstop
779	LOGICAL :: llprevmt, llprevyr
780	LOGICAL :: llnextmt, llnextyr
781	!!----------------------------------------------------------------------
782	idaysec = NINT(rday)
783	!
784	IF( PRESENT(ldprev) ) THEN ; llprev = ldprev
785	ELSE ; llprev = .FALSE.
786	ENDIF
787	IF( PRESENT(ldnext) ) THEN ; llnext = ldnext
788	ELSE ; llnext = .FALSE.
789	ENDIF
790
791	! current file parameters
792	IF( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of the current week
793	isecwk = ksec_week( sdjf%cltype(6:8) ) ! seconds between the beginning of the week and half of current time step
794	llprevmt = isecwk > nsec_month ! longer time since beginning of the current week than the current month
795	llprevyr = llprevmt .AND. nmonth == 1
796	iyr = nyear - COUNT((/llprevyr/))
797	imt = nmonth - COUNT((/llprevmt/)) + 12 * COUNT((/llprevyr/))
798	idy = nday + nmonth_len(nmonth-1) * COUNT((/llprevmt/)) - isecwk / idaysec
799	isecwk = nsec_year - isecwk ! seconds between 00h jan 1st of current year and current week beginning
800	ELSE
801	iyr = nyear
802	imt = nmonth
803	idy = nday
804	isecwk = 0
805	ENDIF
806
807	! previous file parameters
808	IF( llprev ) THEN
809	IF( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of previous week
810	isecwk = isecwk + 7 * idaysec ! seconds between the beginning of previous week and half of the time step
811	llprevmt = isecwk > nsec_month ! longer time since beginning of the previous week than the current month
812	llprevyr = llprevmt .AND. nmonth == 1
813	iyr = nyear - COUNT((/llprevyr/))
814	imt = nmonth - COUNT((/llprevmt/)) + 12 * COUNT((/llprevyr/))
815	idy = nday + nmonth_len(nmonth-1) * COUNT((/llprevmt/)) - isecwk / idaysec
816	isecwk = nsec_year - isecwk ! seconds between 00h jan 1st of current year and previous week beginning
817	ELSE
818	idy = nday - COUNT((/ sdjf%cltype == 'daily' /))
819	imt = nmonth - COUNT((/ sdjf%cltype == 'monthly' .OR. idy == 0 /))
820	iyr = nyear - COUNT((/ sdjf%cltype == 'yearly' .OR. imt == 0 /))
821	IF( idy == 0 ) idy = nmonth_len(imt)
822	IF( imt == 0 ) imt = 12
823	isecwk = 0
824	ENDIF
825	ENDIF
826
827	! next file parameters
828	IF( llnext ) THEN
829	IF( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of next week
830	isecwk = 7 * idaysec - isecwk ! seconds between half of the time step and the beginning of next week
831	llnextmt = isecwk > ( nmonth_len(nmonth)*idaysec - nsec_month ) ! larger than the seconds to the end of the month
832	llnextyr = llnextmt .AND. nmonth == 12
833	iyr = nyear + COUNT((/llnextyr/))
834	imt = nmonth + COUNT((/llnextmt/)) - 12 * COUNT((/llnextyr/))
835	idy = nday - nmonth_len(nmonth) * COUNT((/llnextmt/)) + isecwk / idaysec + 1
836	isecwk = nsec_year + isecwk ! seconds between 00h jan 1st of current year and next week beginning
837	ELSE
838	idy = nday + COUNT((/ sdjf%cltype == 'daily' /))
839	imt = nmonth + COUNT((/ sdjf%cltype == 'monthly' .OR. idy > nmonth_len(nmonth) /))
840	iyr = nyear + COUNT((/ sdjf%cltype == 'yearly' .OR. imt == 13 /))
841	IF( idy > nmonth_len(nmonth) ) idy = 1
842	IF( imt == 13 ) imt = 1
843	isecwk = 0
844	ENDIF
845	ENDIF
846	!
847	! find the last record to be read -> update sdjf%nreclast
848	indexyr = iyr - nyear + 1 ! which year are we looking for? previous(0), current(1) or next(2)?
849	indexmt = imt + 12 * ( indexyr - 1 ) ! which month are we looking for (relatively to current year)?
850	!
851	! Last record to be read in the current file
852	! Predefine the number of record in the file according of its type.
853	! We could compare this number with the number of records in the file and make a stop if the 2 numbers do not match...
854	! However this would be much less fexible (e.g. for tests) and will force to rewite input files according to nleapy...
855	IF ( NINT(sdjf%freqh) == -12 ) THEN ; ireclast = 1 ! yearly mean: consider only 1 record
856	ELSEIF( NINT(sdjf%freqh) == -1 ) THEN ! monthly mean:
857	IF( sdjf%cltype == 'monthly' ) THEN ; ireclast = 1 ! consider that the file has 1 record
858	ELSE ; ireclast = 12 ! consider that the file has 12 record
859	ENDIF
860	ELSE ! higher frequency mean (in hours)
861	IF( sdjf%cltype == 'monthly' ) THEN ; ireclast = NINT( 24. * REAL(nmonth_len(indexmt), wp) / sdjf%freqh )
862	ELSEIF( sdjf%cltype(1:4) == 'week' ) THEN ; ireclast = NINT( 24. * 7. / sdjf%freqh )
863	ELSEIF( sdjf%cltype == 'daily' ) THEN ; ireclast = NINT( 24. / sdjf%freqh )
864	ELSE ; ireclast = NINT( 24. * REAL( nyear_len(indexyr), wp) / sdjf%freqh )
865	ENDIF
866	ENDIF
867
868	sdjf%nreclast = ireclast
869	! Allocate arrays for beginning/middle/end of each record (seconds since Jan. 1st 00h of nit000 year)
870	IF( ALLOCATED(sdjf%nrecsec) ) DEALLOCATE( sdjf%nrecsec )
871	ALLOCATE( sdjf%nrecsec( 0:ireclast ) )
872	!
873	IF ( NINT(sdjf%freqh) == -12 ) THEN ! yearly mean and yearly file
874	SELECT CASE( indexyr )
875	CASE(0) ; sdjf%nrecsec(0) = nsec1jan000 - nyear_len( 0 ) * idaysec
876	CASE(1) ; sdjf%nrecsec(0) = nsec1jan000
877	CASE(2) ; sdjf%nrecsec(0) = nsec1jan000 + nyear_len( 1 ) * idaysec
878	ENDSELECT
879	sdjf%nrecsec(1) = sdjf%nrecsec(0) + nyear_len( indexyr ) * idaysec
880	ELSEIF( NINT(sdjf%freqh) == -1 ) THEN ! monthly mean:
881	IF( sdjf%cltype == 'monthly' ) THEN ! monthly file
882	sdjf%nrecsec(0 ) = nsec1jan000 + nmonth_beg(indexmt )
883	sdjf%nrecsec(1 ) = nsec1jan000 + nmonth_beg(indexmt+1)
884	ELSE ! yearly file
885	ishift = 12 * ( indexyr - 1 )
886	sdjf%nrecsec(0:12) = nsec1jan000 + nmonth_beg(1+ishift:13+ishift)
887	ENDIF
888	ELSE ! higher frequency mean (in hours)
889	IF( sdjf%cltype == 'monthly' ) THEN ; istart = nsec1jan000 + nmonth_beg(indexmt)
890	ELSEIF( sdjf%cltype(1:4) == 'week' ) THEN ; istart = nsec1jan000 + isecwk
891	ELSEIF( sdjf%cltype == 'daily' ) THEN ; istart = nsec1jan000 + nmonth_beg(indexmt) + ( idy - 1 ) * idaysec
892	ELSEIF( indexyr == 0 ) THEN ; istart = nsec1jan000 - nyear_len( 0 ) * idaysec
893	ELSEIF( indexyr == 2 ) THEN ; istart = nsec1jan000 + nyear_len( 1 ) * idaysec
894	ELSE ; istart = nsec1jan000
895	ENDIF
896	zfreqs = sdjf%freqh * rhhmm * rmmss
897	DO jt = 0, sdjf%nreclast
898	sdjf%nrecsec(jt) = istart + NINT( zfreqs * REAL(jt,wp) )
899	END DO
900	ENDIF
901	!
902	IF( sdjf%ln_tint ) THEN ! record time defined in the middle of the record
903	sdjf%nrecsec(1:sdjf%nreclast) = 0.5 * ( sdjf%nrecsec(0:sdjf%nreclast-1) + sdjf%nrecsec(1:sdjf%nreclast) )
904	END IF
905	!
906	sdjf%clname = fld_filename( sdjf, idy, imt, iyr )
907	!
908	END SUBROUTINE fld_def
909
910
911	SUBROUTINE fld_clopn( sdjf )
912	!!---------------------------------------------------------------------
913	!! * ROUTINE fld_clopn *
914	!!
915	!! ** Purpose : close/open the files
916	!!----------------------------------------------------------------------
917	TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables
918	!
919	INTEGER, DIMENSION(2) :: isave
920	LOGICAL :: llprev, llnext, llstop
921	!!----------------------------------------------------------------------
922	!
923	llprev = sdjf%nrecsec(sdjf%nreclast) < nsec000_1jan000 ! file ends before the beginning of the job -> file may not exist
924	llnext = sdjf%nrecsec( 0 ) > nsecend_1jan000 ! file begins after the end of the job -> file may not exist
925
926	llstop = sdjf%ln_clim .OR. .NOT. ( llprev .OR. llnext )
927
928	IF( sdjf%num <= 0 .OR. .NOT. sdjf%ln_clim ) THEN
929	IF( sdjf%num > 0 ) CALL iom_close( sdjf%num ) ! close file if already open
930	CALL iom_open( sdjf%clname, sdjf%num, ldstop = llstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 )
931	ENDIF
932	!
933	IF( sdjf%num <= 0 .AND. .NOT. llstop ) THEN ! file not found but we do accept this...
934	!
935	IF( llprev ) THEN ! previous file does not exist : go back to current and accept to read only the first record
936	CALL ctl_warn('previous file: '//TRIM(sdjf%clname)//' not found -> go back to current year/month/week/day file')
937	isave(1:2) = sdjf%nrecsec(sdjf%nreclast-1:sdjf%nreclast) ! save previous file info
938	CALL fld_def( sdjf ) ! go back to current file
939	sdjf%nreclast = 1 ! force to use only the first record (do as if other were not existing...)
940	sdjf%nrecsec(0:1) = isave(1:2)
941	ENDIF
942	!
943	IF( llnext ) THEN ! next file does not exist : go back to current and accept to read only the last record
944	CALL ctl_warn('next file: '//TRIM(sdjf%clname)//' not found -> go back to current year/month/week/day file')
945	isave(1:2) = sdjf%nrecsec(0:1) ! save next file info
946	CALL fld_def( sdjf ) ! go back to current file
947	! -> read last record but keep record info from the first record of next file
948	sdjf%nrecsec(sdjf%nreclast-1:sdjf%nreclast) = isave(1:2)
949	sdjf%nrecsec(0:sdjf%nreclast-2) = nflag
950	ENDIF
951	!
952	CALL iom_open( sdjf%clname, sdjf%num, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 )
953	!
954	ENDIF
955	!
956	END SUBROUTINE fld_clopn
957
958
959	SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam, knoprint )
960	!!---------------------------------------------------------------------
961	!! * ROUTINE fld_fill *
962	!!
963	!! ** Purpose : fill the data structure (sdf) with the associated information
964	!! read in namelist (sdf_n) and control print
965	!!----------------------------------------------------------------------
966	TYPE(FLD) , DIMENSION(:) , INTENT(inout) :: sdf ! structure of input fields (file informations, fields read)
967	TYPE(FLD_N), DIMENSION(:) , INTENT(in ) :: sdf_n ! array of namelist information structures
968	CHARACTER(len=*) , INTENT(in ) :: cdir ! Root directory for location of flx files
969	CHARACTER(len=*) , INTENT(in ) :: cdcaller ! name of the calling routine
970	CHARACTER(len=*) , INTENT(in ) :: cdtitle ! description of the calling routine
971	CHARACTER(len=*) , INTENT(in ) :: cdnam ! name of the namelist from which sdf_n comes
972	INTEGER , OPTIONAL, INTENT(in ) :: knoprint ! no calling routine information printed
973	!
974	INTEGER :: jf ! dummy indices
975	!!---------------------------------------------------------------------
976	!
977	DO jf = 1, SIZE(sdf)
978	sdf(jf)%clrootname = sdf_n(jf)%clname
979	IF( TRIM(sdf_n(jf)%clname) /= 'NOT USED' ) sdf(jf)%clrootname = TRIM( cdir )//sdf(jf)%clrootname
980	sdf(jf)%clname = "not yet defined"
981	sdf(jf)%freqh = sdf_n(jf)%freqh
982	sdf(jf)%clvar = sdf_n(jf)%clvar
983	sdf(jf)%ln_tint = sdf_n(jf)%ln_tint
984	sdf(jf)%ln_clim = sdf_n(jf)%ln_clim
985	sdf(jf)%cltype = sdf_n(jf)%cltype
986	sdf(jf)%num = -1
987	sdf(jf)%wgtname = " "
988	IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) sdf(jf)%wgtname = TRIM( cdir )//sdf_n(jf)%wname
989	sdf(jf)%lsmname = " "
990	IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 ) sdf(jf)%lsmname = TRIM( cdir )//sdf_n(jf)%lname
991	sdf(jf)%vcomp = sdf_n(jf)%vcomp
992	sdf(jf)%rotn(:) = .TRUE. ! pretend to be rotated -> won't try to rotate data before the first call to fld_get
993	IF( sdf(jf)%cltype(1:4) == 'week' .AND. nn_leapy == 0 ) &
994	& CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs nn_leapy = 1')
995	IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim ) &
996	& CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.')
997	sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
998	sdf(jf)%igrd = 0
999	sdf(jf)%ibdy = 0
1000	sdf(jf)%imap => NULL()
1001	sdf(jf)%ltotvel = .FALSE.
1002	sdf(jf)%lzint = .FALSE.
1003	END DO
1004	!
1005	IF(lwp) THEN ! control print
1006	WRITE(numout,*)
1007	IF( .NOT.PRESENT( knoprint) ) THEN
1008	WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
1009	WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
1010	ENDIF
1011	WRITE(numout,*) ' fld_fill : fill data structure with information from namelist '//TRIM( cdnam )
1012	WRITE(numout,*) ' ~~~~~~~~'
1013	WRITE(numout,*) ' list of files and frequency (>0: in hours ; <0 in months)'
1014	DO jf = 1, SIZE(sdf)
1015	WRITE(numout,*) ' root filename: ' , TRIM( sdf(jf)%clrootname ), ' variable name: ', TRIM( sdf(jf)%clvar )
1016	WRITE(numout,*) ' frequency: ' , sdf(jf)%freqh , &
1017	& ' time interp: ' , sdf(jf)%ln_tint , &
1018	& ' climatology: ' , sdf(jf)%ln_clim
1019	WRITE(numout,*) ' weights: ' , TRIM( sdf(jf)%wgtname ), &
1020	& ' pairing: ' , TRIM( sdf(jf)%vcomp ), &
1021	& ' data type: ' , sdf(jf)%cltype , &
1022	& ' land/sea mask:' , TRIM( sdf(jf)%lsmname )
1023	call flush(numout)
1024	END DO
1025	ENDIF
1026	!
1027	END SUBROUTINE fld_fill
1028
1029
1030	SUBROUTINE wgt_list( sd, kwgt )
1031	!!---------------------------------------------------------------------
1032	!! * ROUTINE wgt_list *
1033	!!
1034	!! ** Purpose : search array of WGTs and find a weights file entry,
1035	!! or return a new one adding it to the end if new entry.
1036	!! the weights data is read in and restructured (fld_weight)
1037	!!----------------------------------------------------------------------
1038	TYPE( FLD ), INTENT(in ) :: sd ! field with name of weights file
1039	INTEGER , INTENT(inout) :: kwgt ! index of weights
1040	!
1041	INTEGER :: kw, nestid ! local integer
1042	LOGICAL :: found ! local logical
1043	!!----------------------------------------------------------------------
1044	!
1045	!! search down linked list
1046	!! weights filename is either present or we hit the end of the list
1047	found = .FALSE.
1048	!
1049	!! because agrif nest part of filenames are now added in iom_open
1050	!! to distinguish between weights files on the different grids, need to track
1051	!! nest number explicitly
1052	nestid = 0
1053	#if defined key_agrif
1054	nestid = Agrif_Fixed()
1055	#endif
1056	DO kw = 1, nxt_wgt-1
1057	IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
1058	ref_wgts(kw)%nestid == nestid) THEN
1059	kwgt = kw
1060	found = .TRUE.
1061	EXIT
1062	ENDIF
1063	END DO
1064	IF( .NOT.found ) THEN
1065	kwgt = nxt_wgt
1066	CALL fld_weight( sd )
1067	ENDIF
1068	!
1069	END SUBROUTINE wgt_list
1070
1071
1072	SUBROUTINE wgt_print( )
1073	!!---------------------------------------------------------------------
1074	!! * ROUTINE wgt_print *
1075	!!
1076	!! ** Purpose : print the list of known weights
1077	!!----------------------------------------------------------------------
1078	INTEGER :: kw !
1079	!!----------------------------------------------------------------------
1080	!
1081	DO kw = 1, nxt_wgt-1
1082	WRITE(numout,*) 'weight file: ',TRIM(ref_wgts(kw)%wgtname)
1083	WRITE(numout,*) ' ddims: ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
1084	WRITE(numout,*) ' numwgt: ',ref_wgts(kw)%numwgt
1085	WRITE(numout,*) ' jpiwgt: ',ref_wgts(kw)%jpiwgt
1086	WRITE(numout,*) ' jpjwgt: ',ref_wgts(kw)%jpjwgt
1087	WRITE(numout,*) ' botleft: ',ref_wgts(kw)%botleft
1088	WRITE(numout,*) ' topright: ',ref_wgts(kw)%topright
1089	IF( ref_wgts(kw)%cyclic ) THEN
1090	WRITE(numout,*) ' cyclical'
1091	IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) ' with overlap of ', ref_wgts(kw)%overlap
1092	ELSE
1093	WRITE(numout,*) ' not cyclical'
1094	ENDIF
1095	IF( ASSOCIATED(ref_wgts(kw)%data_wgt) ) WRITE(numout,*) ' allocated'
1096	END DO
1097	!
1098	END SUBROUTINE wgt_print
1099
1100
1101	SUBROUTINE fld_weight( sd )
1102	!!---------------------------------------------------------------------
1103	!! * ROUTINE fld_weight *
1104	!!
1105	!! ** Purpose : create a new WGT structure and fill in data from file,
1106	!! restructuring as required
1107	!!----------------------------------------------------------------------
1108	TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file
1109	!!
1110	INTEGER :: jn ! dummy loop indices
1111	INTEGER :: inum ! local logical unit
1112	INTEGER :: id ! local variable id
1113	INTEGER :: ipk ! local vertical dimension
1114	INTEGER :: zwrap ! local integer
1115	LOGICAL :: cyclical !
1116	CHARACTER (len=5) :: aname !
1117	INTEGER , DIMENSION(:), ALLOCATABLE :: ddims
1118	INTEGER, DIMENSION(jpi,jpj) :: data_src
1119	REAL(wp), DIMENSION(jpi,jpj) :: data_tmp
1120	!!----------------------------------------------------------------------
1121	!
1122	IF( nxt_wgt > tot_wgts ) THEN
1123	CALL ctl_stop("fld_weight: weights array size exceeded, increase tot_wgts")
1124	ENDIF
1125	!
1126	!! new weights file entry, add in extra information
1127	!! a weights file represents a 2D grid of a certain shape, so we assume that the current
1128	!! input data file is representative of all other files to be opened and processed with the
1129	!! current weights file
1130
1131	!! open input data file (non-model grid)
1132	CALL iom_open( sd%clname, inum, ldiof = LEN(TRIM(sd%wgtname)) > 0 )
1133
1134	!! get dimensions: we consider 2D data as 3D data with vertical dim size = 1
1135	IF( SIZE(sd%fnow, 3) > 0 ) THEN
1136	ALLOCATE( ddims(4) )
1137	ELSE
1138	ALLOCATE( ddims(3) )
1139	ENDIF
1140	id = iom_varid( inum, sd%clvar, ddims )
1141
1142	!! close it
1143	CALL iom_close( inum )
1144
1145	!! now open the weights file
1146
1147	CALL iom_open ( sd%wgtname, inum ) ! interpolation weights
1148	IF( inum > 0 ) THEN
1149
1150	!! determine whether we have an east-west cyclic grid
1151	!! from global attribute called "ew_wrap" in the weights file
1152	!! note that if not found, iom_getatt returns -999 and cyclic with no overlap is assumed
1153	!! since this is the most common forcing configuration
1154
1155	CALL iom_getatt(inum, 'ew_wrap', zwrap)
1156	IF( zwrap >= 0 ) THEN
1157	cyclical = .TRUE.
1158	ELSE IF( zwrap == -999 ) THEN
1159	cyclical = .TRUE.
1160	zwrap = 0
1161	ELSE
1162	cyclical = .FALSE.
1163	ENDIF
1164
1165	ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
1166	ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
1167	ref_wgts(nxt_wgt)%wgtname = sd%wgtname
1168	ref_wgts(nxt_wgt)%overlap = zwrap
1169	ref_wgts(nxt_wgt)%cyclic = cyclical
1170	ref_wgts(nxt_wgt)%nestid = 0
1171	#if defined key_agrif
1172	ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
1173	#endif
1174	!! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
1175	!! for each weight wgtNN there is an integer array srcNN which gives the point in
1176	!! the input data grid which is to be multiplied by the weight
1177	!! they are both arrays on the model grid so the result of the multiplication is
1178	!! added into an output array on the model grid as a running sum
1179
1180	!! two possible cases: bilinear (4 weights) or bicubic (16 weights)
1181	id = iom_varid(inum, 'src05', ldstop=.FALSE.)
1182	IF( id <= 0) THEN
1183	ref_wgts(nxt_wgt)%numwgt = 4
1184	ELSE
1185	ref_wgts(nxt_wgt)%numwgt = 16
1186	ENDIF
1187
1188	ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
1189	ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
1190	ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
1191
1192	DO jn = 1,4
1193	aname = ' '
1194	WRITE(aname,'(a3,i2.2)') 'src',jn
1195	data_tmp(:,:) = 0
1196	CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
1197	data_src(:,:) = INT(data_tmp(:,:))
1198	ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
1199	ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
1200	END DO
1201
1202	DO jn = 1, ref_wgts(nxt_wgt)%numwgt
1203	aname = ' '
1204	WRITE(aname,'(a3,i2.2)') 'wgt',jn
1205	ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
1206	CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
1207	END DO
1208	CALL iom_close (inum)
1209
1210	! find min and max indices in grid
1211	ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1212	ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1213	ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
1214	ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
1215
1216	! and therefore dimensions of the input box
1217	ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
1218	ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
1219
1220	! shift indexing of source grid
1221	ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
1222	ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
1223
1224	! create input grid, give it a halo to allow gradient calculations
1225	! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
1226	! a more robust solution will be given in next release
1227	ipk = SIZE(sd%fnow, 3)
1228	ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
1229	IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
1230	!
1231	nxt_wgt = nxt_wgt + 1
1232	!
1233	ELSE
1234	CALL ctl_stop( ' fld_weight : unable to read the file ' )
1235	ENDIF
1236
1237	DEALLOCATE (ddims )
1238	!
1239	END SUBROUTINE fld_weight
1240
1241
1242	SUBROUTINE apply_seaoverland( clmaskfile, zfieldo, jpi1_lsm, jpi2_lsm, jpj1_lsm, &
1243	& jpj2_lsm, itmpi, itmpj, itmpz, rec1_lsm, recn_lsm )
1244	!!---------------------------------------------------------------------
1245	!! * ROUTINE apply_seaoverland *
1246	!!
1247	!! ** Purpose : avoid spurious fluxes in coastal or near-coastal areas
1248	!! due to the wrong usage of "land" values from the coarse
1249	!! atmospheric model when spatial interpolation is required
1250	!! D. Delrosso INGV
1251	!!----------------------------------------------------------------------
1252	INTEGER, INTENT(in ) :: itmpi,itmpj,itmpz ! lengths
1253	INTEGER, INTENT(in ) :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices
1254	INTEGER, DIMENSION(3), INTENT(in ) :: rec1_lsm,recn_lsm ! temporary arrays for start and length
1255	REAL(wp),DIMENSION (:,:,:),INTENT(inout) :: zfieldo ! input/output array for seaoverland application
1256	CHARACTER (len=100), INTENT(in ) :: clmaskfile ! land/sea mask file name
1257	!
1258	INTEGER :: inum,jni,jnj,jnz,jc ! local indices
1259	REAL(wp),DIMENSION (:,:,:),ALLOCATABLE :: zslmec1 ! local array for land point detection
1260	REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfieldn ! array of forcing field with undeff for land points
1261	REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfield ! array of forcing field
1262	!!---------------------------------------------------------------------
1263	!
1264	ALLOCATE ( zslmec1(itmpi,itmpj,itmpz), zfieldn(itmpi,itmpj), zfield(itmpi,itmpj) )
1265	!
1266	! Retrieve the land sea mask data
1267	CALL iom_open( clmaskfile, inum )
1268	SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1269	CASE(1)
1270	CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm)
1271	CASE DEFAULT
1272	CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm)
1273	END SELECT
1274	CALL iom_close( inum )
1275	!
1276	DO jnz=1,rec1_lsm(3) !! Loop over k dimension
1277	!
1278	DO jni = 1, itmpi !! copy the original field into a tmp array
1279	DO jnj = 1, itmpj !! substituting undeff over land points
1280	zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz)
1281	IF( zslmec1(jni,jnj,jnz) == 1. ) zfieldn(jni,jnj) = undeff_lsm
1282	END DO
1283	END DO
1284	!
1285	CALL seaoverland( zfieldn, itmpi, itmpj, zfield )
1286	DO jc = 1, nn_lsm
1287	CALL seaoverland( zfield, itmpi, itmpj, zfield )
1288	END DO
1289	!
1290	! Check for Undeff and substitute original values
1291	IF( ANY(zfield==undeff_lsm) ) THEN
1292	DO jni = 1, itmpi
1293	DO jnj = 1, itmpj
1294	IF( zfield(jni,jnj)==undeff_lsm ) zfield(jni,jnj) = zfieldo(jni,jnj,jnz)
1295	END DO
1296	END DO
1297	ENDIF
1298	!
1299	zfieldo(:,:,jnz) = zfield(:,:)
1300	!
1301	END DO !! End Loop over k dimension
1302	!
1303	DEALLOCATE ( zslmec1, zfieldn, zfield )
1304	!
1305	END SUBROUTINE apply_seaoverland
1306
1307
1308	SUBROUTINE seaoverland( zfieldn, ileni, ilenj, zfield )
1309	!!---------------------------------------------------------------------
1310	!! * ROUTINE seaoverland *
1311	!!
1312	!! ** Purpose : create shifted matrices for seaoverland application
1313	!! D. Delrosso INGV
1314	!!----------------------------------------------------------------------
1315	INTEGER , INTENT(in ) :: ileni,ilenj ! lengths
1316	REAL, DIMENSION (ileni,ilenj), INTENT(in ) :: zfieldn ! array of forcing field with undeff for land points
1317	REAL, DIMENSION (ileni,ilenj), INTENT( out) :: zfield ! array of forcing field
1318	!
1319	REAL , DIMENSION (ileni,ilenj) :: zmat1, zmat2, zmat3, zmat4 ! local arrays
1320	REAL , DIMENSION (ileni,ilenj) :: zmat5, zmat6, zmat7, zmat8 ! - -
1321	REAL , DIMENSION (ileni,ilenj) :: zlsm2d ! - -
1322	REAL , DIMENSION (ileni,ilenj,8) :: zlsm3d ! - -
1323	LOGICAL, DIMENSION (ileni,ilenj,8) :: ll_msknan3d ! logical mask for undeff detection
1324	LOGICAL, DIMENSION (ileni,ilenj) :: ll_msknan2d ! logical mask for undeff detection
1325	!!----------------------------------------------------------------------
1326	zmat8 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(:,1)/) , DIM=2 )
1327	zmat1 = eoshift( zmat8 , SHIFT=-1 , BOUNDARY = (/zmat8(1,:)/) , DIM=1 )
1328	zmat2 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(1,:)/) , DIM=1 )
1329	zmat4 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(:,ilenj)/) , DIM=2 )
1330	zmat3 = eoshift( zmat4 , SHIFT=-1 , BOUNDARY = (/zmat4(1,:)/) , DIM=1 )
1331	zmat5 = eoshift( zmat4 , SHIFT= 1 , BOUNDARY = (/zmat4(ileni,:)/) , DIM=1 )
1332	zmat6 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(ileni,:)/) , DIM=1 )
1333	zmat7 = eoshift( zmat8 , SHIFT= 1 , BOUNDARY = (/zmat8(ileni,:)/) , DIM=1 )
1334	!
1335	zlsm3d = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /))
1336	ll_msknan3d = .NOT.( zlsm3d == undeff_lsm )
1337	ll_msknan2d = .NOT.( zfieldn == undeff_lsm ) ! FALSE where is Undeff (land)
1338	zlsm2d = SUM( zlsm3d, 3 , ll_msknan3d ) / MAX( 1 , COUNT( ll_msknan3d , 3 ) )
1339	WHERE( COUNT( ll_msknan3d , 3 ) == 0._wp ) zlsm2d = undeff_lsm
1340	zfield = MERGE( zfieldn, zlsm2d, ll_msknan2d )
1341	!
1342	END SUBROUTINE seaoverland
1343
1344
1345	SUBROUTINE fld_interp( num, clvar, kw, kk, dta, &
1346	& nrec, lsmfile)
1347	!!---------------------------------------------------------------------
1348	!! * ROUTINE fld_interp *
1349	!!
1350	!! ** Purpose : apply weights to input gridded data to create data
1351	!! on model grid
1352	!!----------------------------------------------------------------------
1353	INTEGER , INTENT(in ) :: num ! stream number
1354	CHARACTER(LEN=*) , INTENT(in ) :: clvar ! variable name
1355	INTEGER , INTENT(in ) :: kw ! weights number
1356	INTEGER , INTENT(in ) :: kk ! vertical dimension of kk
1357	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: dta ! output field on model grid
1358	INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice)
1359	CHARACTER(LEN=*) , INTENT(in ) :: lsmfile ! land sea mask file name
1360	!
1361	INTEGER, DIMENSION(3) :: rec1, recn ! temporary arrays for start and length
1362	INTEGER, DIMENSION(3) :: rec1_lsm, recn_lsm ! temporary arrays for start and length in case of seaoverland
1363	INTEGER :: ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2 ! temporary indices
1364	INTEGER :: jk, jn, jm, jir, jjr ! loop counters
1365	INTEGER :: ni, nj ! lengths
1366	INTEGER :: jpimin,jpiwid ! temporary indices
1367	INTEGER :: jpimin_lsm,jpiwid_lsm ! temporary indices
1368	INTEGER :: jpjmin,jpjwid ! temporary indices
1369	INTEGER :: jpjmin_lsm,jpjwid_lsm ! temporary indices
1370	INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices
1371	INTEGER :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices
1372	INTEGER :: itmpi,itmpj,itmpz ! lengths
1373	REAL(wp),DIMENSION(:,:,:), ALLOCATABLE :: ztmp_fly_dta ! local array of values on input grid
1374	!!----------------------------------------------------------------------
1375	!
1376	!! for weighted interpolation we have weights at four corners of a box surrounding
1377	!! a model grid point, each weight is multiplied by a grid value (bilinear case)
1378	!! or by a grid value and gradients at the corner point (bicubic case)
1379	!! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
1380
1381	!! sub grid from non-model input grid which encloses all grid points in this nemo process
1382	jpimin = ref_wgts(kw)%botleft(1)
1383	jpjmin = ref_wgts(kw)%botleft(2)
1384	jpiwid = ref_wgts(kw)%jpiwgt
1385	jpjwid = ref_wgts(kw)%jpjwgt
1386
1387	!! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients
1388	rec1(1) = MAX( jpimin-1, 1 )
1389	rec1(2) = MAX( jpjmin-1, 1 )
1390	rec1(3) = 1
1391	recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
1392	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1393	recn(3) = kk
1394
1395	!! where we need to put it in the non-nemo grid fly_dta
1396	!! note that jpi1 and jpj1 only differ from 1 when jpimin and jpjmin are 1
1397	!! (ie at the extreme west or south of the whole input grid) and similarly for jpi2 and jpj2
1398	jpi1 = 2 + rec1(1) - jpimin
1399	jpj1 = 2 + rec1(2) - jpjmin
1400	jpi2 = jpi1 + recn(1) - 1
1401	jpj2 = jpj1 + recn(2) - 1
1402
1403
1404	IF( LEN( TRIM(lsmfile) ) > 0 ) THEN
1405	!! indeces for ztmp_fly_dta
1406	! --------------------------
1407	rec1_lsm(1)=MAX(rec1(1)-nn_lsm,1) ! starting index for enlarged external data, x direction
1408	rec1_lsm(2)=MAX(rec1(2)-nn_lsm,1) ! starting index for enlarged external data, y direction
1409	rec1_lsm(3) = 1 ! vertical dimension
1410	recn_lsm(1)=MIN(rec1(1)-rec1_lsm(1)+recn(1)+nn_lsm,ref_wgts(kw)%ddims(1)-rec1_lsm(1)) ! n points in x direction
1411	recn_lsm(2)=MIN(rec1(2)-rec1_lsm(2)+recn(2)+nn_lsm,ref_wgts(kw)%ddims(2)-rec1_lsm(2)) ! n points in y direction
1412	recn_lsm(3) = kk ! number of vertical levels in the input file
1413
1414	! Avoid out of bound
1415	jpimin_lsm = MAX( rec1_lsm(1)+1, 1 )
1416	jpjmin_lsm = MAX( rec1_lsm(2)+1, 1 )
1417	jpiwid_lsm = MIN( recn_lsm(1)-2,ref_wgts(kw)%ddims(1)-rec1(1)+1)
1418	jpjwid_lsm = MIN( recn_lsm(2)-2,ref_wgts(kw)%ddims(2)-rec1(2)+1)
1419
1420	jpi1_lsm = 2+rec1_lsm(1)-jpimin_lsm
1421	jpj1_lsm = 2+rec1_lsm(2)-jpjmin_lsm
1422	jpi2_lsm = jpi1_lsm + recn_lsm(1) - 1
1423	jpj2_lsm = jpj1_lsm + recn_lsm(2) - 1
1424
1425
1426	itmpi=jpi2_lsm-jpi1_lsm+1
1427	itmpj=jpj2_lsm-jpj1_lsm+1
1428	itmpz=kk
1429	ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz))
1430	ztmp_fly_dta(:,:,:) = 0.0
1431	SELECT CASE( SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
1432	CASE(1)
1433	CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), &
1434	& nrec, rec1_lsm, recn_lsm)
1435	CASE DEFAULT
1436	CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), &
1437	& nrec, rec1_lsm, recn_lsm)
1438	END SELECT
1439	CALL apply_seaoverland(lsmfile,ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), &
1440	& jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm, &
1441	& itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
1442
1443
1444	! Relative indeces for remapping
1445	ii_lsm1 = (rec1(1)-rec1_lsm(1))+1
1446	ii_lsm2 = (ii_lsm1+recn(1))-1
1447	ij_lsm1 = (rec1(2)-rec1_lsm(2))+1
1448	ij_lsm2 = (ij_lsm1+recn(2))-1
1449
1450	ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1451	ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:) = ztmp_fly_dta(ii_lsm1:ii_lsm2,ij_lsm1:ij_lsm2,:)
1452	DEALLOCATE(ztmp_fly_dta)
1453
1454	ELSE
1455
1456	ref_wgts(kw)%fly_dta(:,:,:) = 0.0
1457	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn)
1458	ENDIF
1459
1460
1461	!! first four weights common to both bilinear and bicubic
1462	!! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft
1463	!! note that we have to offset by 1 into fly_dta array because of halo
1464	dta(:,:,:) = 0.0
1465	DO jk = 1,4
1466	DO jn = 1, jpj
1467	DO jm = 1,jpi
1468	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1469	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1470	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,:)
1471	END DO
1472	END DO
1473	END DO
1474
1475	IF(ref_wgts(kw)%numwgt .EQ. 16) THEN
1476
1477	!! fix up halo points that we couldnt read from file
1478	IF( jpi1 == 2 ) THEN
1479	ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:)
1480	ENDIF
1481	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1482	ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:)
1483	ENDIF
1484	IF( jpj1 == 2 ) THEN
1485	ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:)
1486	ENDIF
1487	IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
1488	ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:)
1489	ENDIF
1490
1491	!! if data grid is cyclic we can do better on east-west edges
1492	!! but have to allow for whether first and last columns are coincident
1493	IF( ref_wgts(kw)%cyclic ) THEN
1494	rec1(2) = MAX( jpjmin-1, 1 )
1495	recn(1) = 1
1496	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
1497	jpj1 = 2 + rec1(2) - jpjmin
1498	jpj2 = jpj1 + recn(2) - 1
1499	IF( jpi1 == 2 ) THEN
1500	rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap
1501	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1502	ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1503	ENDIF
1504	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
1505	rec1(1) = 1 + ref_wgts(kw)%overlap
1506	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
1507	ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
1508	ENDIF
1509	ENDIF
1510
1511	! gradient in the i direction
1512	DO jk = 1,4
1513	DO jn = 1, jpj
1514	DO jm = 1,jpi
1515	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1516	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1517	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 * &
1518	(ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:))
1519	END DO
1520	END DO
1521	END DO
1522
1523	! gradient in the j direction
1524	DO jk = 1,4
1525	DO jn = 1, jpj
1526	DO jm = 1,jpi
1527	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1528	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1529	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 * &
1530	(ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:))
1531	END DO
1532	END DO
1533	END DO
1534
1535	! gradient in the ij direction
1536	DO jk = 1,4
1537	DO jn = 1, jpj
1538	DO jm = 1,jpi
1539	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1540	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1541	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
1542	(ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni ,nj+2,:)) - &
1543	(ref_wgts(kw)%fly_dta(ni+2,nj ,:) - ref_wgts(kw)%fly_dta(ni ,nj ,:)))
1544	END DO
1545	END DO
1546	END DO
1547	!
1548	ENDIF
1549	!
1550	END SUBROUTINE fld_interp
1551
1552
1553	FUNCTION fld_filename( sdjf, kday, kmonth, kyear )
1554	!!---------------------------------------------------------------------
1555	!! * FUNCTION fld_filename *
1556	!!
1557	!! ** Purpose : define the filename according to a given date
1558	!!---------------------------------------------------------------------
1559	TYPE(FLD), INTENT(in) :: sdjf ! input field related variables
1560	INTEGER , INTENT(in) :: kday, kmonth, kyear
1561	!
1562	CHARACTER(len = 256) :: clname, fld_filename
1563	!!---------------------------------------------------------------------
1564
1565
1566	! build the new filename if not climatological data
1567	clname=TRIM(sdjf%clrootname)
1568	!
1569	! note that sdjf%ln_clim is is only acting on the presence of the year in the file name
1570	IF( .NOT. sdjf%ln_clim ) THEN
1571	WRITE(clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), kyear ! add year
1572	IF( sdjf%cltype /= 'yearly' ) WRITE(clname, '(a, "m",i2.2)' ) TRIM( clname ), kmonth ! add month
1573	ELSE
1574	! build the new filename if climatological data
1575	IF( sdjf%cltype /= 'yearly' ) WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), kmonth ! add month
1576	ENDIF
1577	IF( sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) &
1578	& WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname ), kday ! add day
1579
1580	fld_filename = clname
1581
1582	END FUNCTION fld_filename
1583
1584
1585	FUNCTION ksec_week( cdday )
1586	!!---------------------------------------------------------------------
1587	!! * FUNCTION ksec_week *
1588	!!
1589	!! ** Purpose : seconds between 00h of the beginning of the week and half of the current time step
1590	!!---------------------------------------------------------------------
1591	CHARACTER(len=*), INTENT(in) :: cdday ! first 3 letters of the first day of the weekly file
1592	!!
1593	INTEGER :: ksec_week ! output variable
1594	INTEGER :: ijul, ishift ! local integer
1595	CHARACTER(len=3),DIMENSION(7) :: cl_week
1596	!!----------------------------------------------------------------------
1597	cl_week = (/"sun","sat","fri","thu","wed","tue","mon"/)
1598	DO ijul = 1, 7
1599	IF( cl_week(ijul) == TRIM(cdday) ) EXIT
1600	END DO
1601	IF( ijul .GT. 7 ) CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): '//TRIM(cdday) )
1602	!
1603	ishift = ijul * NINT(rday)
1604	!
1605	ksec_week = nsec_monday + ishift
1606	ksec_week = MOD( ksec_week, 7*NINT(rday) )
1607	!
1608	END FUNCTION ksec_week
1609
1610	!!======================================================================
1611	END MODULE fldread

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

Download in other formats: