Context Navigation

Back to Ticket #632

Ticket #632: fldread.F90_new

File fldread.F90_new, 54.7 KB (added by cbricaud, 14 years ago)

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

Download in other formats:

Original Format

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

Context Navigation

Ticket #632: fldread.F90_new

Download in other formats: