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

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

fldread.F90 in branches/dev_1784_WEEK/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/dev_1784_WEEK/NEMO/OPA_SRC/SBC/fldread.F90 @ 2059

Last change on this file since 2059 was 2059, checked in by cbricaud, 14 years ago
add modifications to have the possibility to read weekly files with fldread
Property svn:keywords set to `Id`
File size: 55.9 KB

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

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

Download in other formats: