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

Context Navigation

source: trunk/NEMO/OPA_SRC/SBC/fldread.F90 @ 2468

Last change on this file since 2468 was 2323, checked in by smasson, 14 years ago
bugfix: temporal shift in fldread in the trunk, see ticket:745
Property svn:keywords set to `Id`
File size: 52.5 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 phycst ! ???
18	USE in_out_manager ! I/O manager
19	USE iom ! I/O manager library
20	USE geo2ocean ! for vector rotation on to model grid
21
22	IMPLICIT NONE
23	PRIVATE
24
25	TYPE, PUBLIC :: FLD_N !: Namelist field informations
26	CHARACTER(len = 256) :: clname ! generic name of the NetCDF flux file
27	INTEGER :: nfreqh ! frequency of each flux file
28	CHARACTER(len = 34) :: clvar ! generic name of the variable in the NetCDF flux file
29	LOGICAL :: ln_tint ! time interpolation or not (T/F)
30	LOGICAL :: ln_clim ! climatology or not (T/F)
31	CHARACTER(len = 7) :: cltype ! type of data file 'daily', 'monthly' or yearly'
32	CHARACTER(len = 34) :: wname ! generic name of a NetCDF weights file to be used, blank if not
33	CHARACTER(len = 34) :: vcomp ! symbolic component name if a vector that needs rotation
34	! a string starting with "U" or "V" for each component
35	! chars 2 onwards identify which components go together
36	END TYPE FLD_N
37
38	TYPE, PUBLIC :: FLD !: Input field related variables
39	CHARACTER(len = 256) :: clrootname ! generic name of the NetCDF file
40	CHARACTER(len = 256) :: clname ! current name of the NetCDF file
41	INTEGER :: nfreqh ! frequency of each flux file
42	CHARACTER(len = 34) :: clvar ! generic name of the variable in the NetCDF flux file
43	LOGICAL :: ln_tint ! time interpolation or not (T/F)
44	LOGICAL :: ln_clim ! climatology or not (T/F)
45	CHARACTER(len = 7) :: cltype ! type of data file 'daily', 'monthly' or yearly'
46	INTEGER :: num ! iom id of the jpfld files to be read
47	INTEGER :: nswap_sec ! swapping time in second since Jan. 1st 00h of nit000 year
48	INTEGER , DIMENSION(2) :: nrec_b ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year)
49	INTEGER , DIMENSION(2) :: nrec_a ! after record (1: index, 2: second since Jan. 1st 00h of nit000 year)
50	REAL(wp) , ALLOCATABLE, DIMENSION(:,:) :: fnow ! input fields interpolated to now time step
51	REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:) :: fdta ! 2 consecutive record of input fields
52	CHARACTER(len = 256) :: wgtname ! current name of the NetCDF weight file acting as a key
53	! into the WGTLIST structure
54	CHARACTER(len = 34) :: vcomp ! symbolic name for a vector component that needs rotation
55	LOGICAL , DIMENSION(2) :: rotn ! flag to indicate whether field has been rotated
56	END TYPE FLD
57
58	!$AGRIF_DO_NOT_TREAT
59
60	!! keep list of all weights variables so they're only read in once
61	!! need to add AGRIF directives not to process this structure
62	!! also need to force wgtname to include AGRIF nest number
63	TYPE :: WGT !: Input weights related variables
64	CHARACTER(len = 256) :: wgtname ! current name of the NetCDF weight file
65	INTEGER , DIMENSION(2) :: ddims ! shape of input grid
66	INTEGER , DIMENSION(2) :: botleft ! top left corner of box in input grid containing
67	! current processor grid
68	INTEGER , DIMENSION(2) :: topright ! top right corner of box
69	INTEGER :: jpiwgt ! width of box on input grid
70	INTEGER :: jpjwgt ! height of box on input grid
71	INTEGER :: numwgt ! number of weights (4=bilinear, 16=bicubic)
72	INTEGER :: nestid ! for agrif, keep track of nest we're in
73	INTEGER :: offset ! =0 when cyclic grid has coincident first/last columns,
74	! =1 when they assumed to be one grid spacing apart
75	! =-1 otherwise
76	LOGICAL :: cyclic ! east-west cyclic or not
77	INTEGER, DIMENSION(:,:,:), POINTER :: data_jpi ! array of source integers
78	INTEGER, DIMENSION(:,:,:), POINTER :: data_jpj ! array of source integers
79	REAL(wp), DIMENSION(:,:,:), POINTER :: data_wgt ! array of weights on model grid
80	REAL(wp), DIMENSION(:,:), POINTER :: fly_dta ! array of values on input grid
81	REAL(wp), DIMENSION(:,:), POINTER :: col2 ! temporary array for reading in columns
82	END TYPE WGT
83
84	INTEGER, PARAMETER :: tot_wgts = 10
85	TYPE( WGT ), DIMENSION(tot_wgts) :: ref_wgts ! array of wgts
86	INTEGER :: nxt_wgt = 1 ! point to next available space in ref_wgts array
87
88	!$AGRIF_END_DO_NOT_TREAT
89
90	PUBLIC fld_read, fld_fill ! called by sbc... modules
91
92	!!----------------------------------------------------------------------
93	!! OPA 9.0 , LOCEAN-IPSL (2006)
94	!! $Id$
95	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
96	!!----------------------------------------------------------------------
97
98	CONTAINS
99
100	SUBROUTINE fld_read( kt, kn_fsbc, sd )
101	!!---------------------------------------------------------------------
102	!! * ROUTINE fld_read *
103	!!
104	!! ** Purpose : provide at each time step the surface ocean fluxes
105	!! (momentum, heat, freshwater and runoff)
106	!!
107	!! ** Method : READ each input fields in NetCDF files using IOM
108	!! and intepolate it to the model time-step.
109	!! Several assumptions are made on the input file:
110	!! blahblahblah....
111	!!----------------------------------------------------------------------
112	INTEGER , INTENT(in ) :: kt ! ocean time step
113	INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step)
114	TYPE(FLD), INTENT(inout), DIMENSION(:) :: sd ! input field related variables
115	!!
116	CHARACTER (LEN=34) :: acomp ! dummy weight name
117	INTEGER :: kf, nf ! dummy indices
118	INTEGER :: imf ! size of the structure sd
119	REAL(wp), DIMENSION(jpi,jpj) :: utmp, vtmp! temporary arrays for vector rotation
120
121	INTEGER :: jf ! dummy indices
122	INTEGER :: kw ! index into wgts array
123	INTEGER :: ireclast ! last record to be read in the current year file
124	INTEGER :: isecend ! number of second since Jan. 1st 00h of nit000 year at nitend
125	INTEGER :: isecsbc ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step
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	isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdttra(1)) ! centrered in the middle of sbc time step
136	!
137	! ! ===================== !
138	DO jf = 1, imf ! LOOP OVER FIELD !
139	! ! ===================== !
140	!
141	IF( kt == nit000 ) CALL fld_init( kn_fsbc, sd(jf) )
142	!
143	! read/update the after data?
144	IF( isecsbc > sd(jf)%nswap_sec ) THEN
145
146	IF( sd(jf)%ln_tint ) THEN ! time interpolation: swap before record field
147	!CDIR COLLAPSE
148	sd(jf)%fdta(:,:,1) = sd(jf)%fdta(:,:,2)
149	sd(jf)%rotn(1) = sd(jf)%rotn(2)
150	ENDIF
151
152	! update record informations
153	CALL fld_rec( kn_fsbc, sd(jf) )
154
155	! do we have to change the year/month/day of the forcing field??
156	IF( sd(jf)%ln_tint ) THEN
157	! if we do time interpolation we will need to open next year/month/day file before the end of the current one
158	! 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
159	! larger than the record number that should be read for current year/month/day (for ex. 13 for monthly mean file)
160
161	! last record to be read in the current file
162	IF( sd(jf)%nfreqh == -1 ) THEN ; ireclast = 12
163	ELSE
164	IF( sd(jf)%cltype == 'monthly' ) THEN ; ireclast = 24 * nmonth_len(nmonth) / 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( isecsbc - 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( kn_fsbc, 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	INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step)
317	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
318	!!
319	LOGICAL :: llprevyr ! are we reading previous year file?
320	LOGICAL :: llprevmth ! are we reading previous month 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	llprevday = .FALSE.
335
336	! define record informations
337	CALL fld_rec( kn_fsbc, sdjf )
338
339	! Note: shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
340
341	IF( sdjf%ln_tint ) THEN ! we need to read the previous record and we will put it in the current record structure
342
343	IF( sdjf%nrec_b(1) == 0 ) THEN ! we redefine record sdjf%nrec_b(1) with the last record of previous year file
344	IF( sdjf%nfreqh == -1 ) THEN ! monthly mean
345	IF( sdjf%cltype == 'monthly' ) THEN ! monthly file
346	sdjf%nrec_b(1) = 1 ! force to read the unique record
347	llprevmth = .NOT. sdjf%ln_clim ! use previous month file?
348	llprevyr = llprevmth .AND. nmonth == 1 ! use previous year file?
349	ELSE ! yearly file
350	sdjf%nrec_b(1) = 12 ! force to read december mean
351	llprevyr = .NOT. sdjf%ln_clim ! use previous year file?
352	ENDIF
353	ELSE
354	IF( sdjf%cltype == 'monthly' ) THEN ! monthly file
355	sdjf%nrec_b(1) = 24 * nmonth_len(nmonth-1) / sdjf%nfreqh ! last record of previous month
356	llprevmth = .NOT. sdjf%ln_clim ! use previous month file?
357	llprevyr = llprevmth .AND. nmonth == 1 ! use previous year file?
358	ELSEIF( sdjf%cltype == 'daily' ) THEN ! daily file
359	sdjf%nrec_b(1) = 24 / sdjf%nfreqh ! last record of previous day
360	llprevday = .NOT. sdjf%ln_clim ! use previous day file?
361	llprevmth = llprevday .AND. nday == 1 ! use previous month file?
362	llprevyr = llprevmth .AND. nmonth == 1 ! use previous year file?
363	ELSE ! yearly file
364	sdjf%nrec_b(1) = 24 * nyear_len(0) / sdjf%nfreqh ! last record of year month
365	llprevyr = .NOT. sdjf%ln_clim ! use previous year file?
366	ENDIF
367	ENDIF
368	ENDIF
369	llprev = llprevyr .OR. llprevmth .OR. llprevday
370
371	CALL fld_clopn( sdjf, nyear - COUNT((/llprevyr /)) , &
372	& nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /)), &
373	& nday - COUNT((/llprevday/)) + nmonth_len(nmonth-1) * COUNT((/llprevmth/)), .NOT. llprev )
374
375	! if previous year/month/day file does not exist, we switch to the current year/month/day
376	IF( llprev .AND. sdjf%num <= 0 ) THEN
377	CALL ctl_warn( 'previous year/month/day file: '//TRIM(sdjf%clname)//' not present -> back to current year/month/day')
378	! we force to read the first record of the current year/month/day instead of last record of previous year/month/day
379	llprev = .false.
380	sdjf%nrec_b(1) = 1
381	CALL fld_clopn( sdjf, nyear, nmonth, nday )
382	ENDIF
383
384	IF( llprev ) THEN ! check if the last record sdjf%nrec_n(1) exists in the file
385	idvar = iom_varid( sdjf%num, sdjf%clvar ) ! id of the variable sdjf%clvar
386	IF( idvar <= 0 ) RETURN
387	inrec = iom_file( sdjf%num )%dimsz( iom_file( sdjf%num )%ndims(idvar), idvar ) ! size of the last dim of idvar
388	sdjf%nrec_b(1) = MIN( sdjf%nrec_b(1), inrec ) ! make sure we select an existing record
389	ENDIF
390
391	! read before data into sdjf%fdta(:,:,2) because we will swap data in the following part of fld_read
392	IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
393	CALL wgt_list( sdjf, kwgt )
394	CALL fld_interp( sdjf%num, sdjf%clvar, kwgt, sdjf%fdta(:,:,2), sdjf%nrec_b(1) )
395	ELSE
396	CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,2), sdjf%nrec_b(1) )
397	ENDIF
398	sdjf%rotn(2) = .FALSE.
399
400	clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i4, ' at time = ', f7.2, ' days')"
401	IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_b(1), REAL(sdjf%nrec_b(2),wp)/rday
402
403	IF( llprev ) CALL iom_close( sdjf%num ) ! close previous year file (-> redefine sdjf%num to 0)
404
405	ENDIF
406
407	IF( sdjf%num <= 0 ) CALL fld_clopn( sdjf, nyear, nmonth, nday ) ! make sure current year/month/day file is opened
408
409	sdjf%nswap_sec = nsec_year + nsec1jan000 - 1 ! force read/update the after data in the following part of fld_read
410
411	END SUBROUTINE fld_init
412
413
414	SUBROUTINE fld_rec( kn_fsbc, sdjf )
415	!!---------------------------------------------------------------------
416	!! * ROUTINE fld_rec *
417	!!
418	!! ** Purpose : compute nrec_a, nrec_b and nswap_sec
419	!!
420	!! ** Method :
421	!!----------------------------------------------------------------------
422	INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step)
423	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
424	!!
425	INTEGER :: irec ! record number
426	INTEGER :: isecd ! rday
427	REAL(wp) :: ztmp ! temporary variable
428	INTEGER :: ifreq_sec ! frequency mean (in seconds)
429	!!----------------------------------------------------------------------
430	!
431	! Note: shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
432	!
433	IF( sdjf%nfreqh == -1 ) THEN ! monthly mean
434	!
435	IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record
436	!
437	! INT( ztmp )
438	! /\|\
439	! 1 \| *----
440	! 0 \|----(
441	! \|----+----\|--> time
442	! 0 /\|\ 1 (nday/nmonth_len(nmonth))
443	! \|
444	! \|
445	! forcing record : nmonth
446	!
447	ztmp = REAL( nday, wp ) / REAL( nmonth_len(nmonth), wp ) + 0.5
448	ELSE
449	ztmp = 0.e0
450	ENDIF
451	irec = nmonth + INT( ztmp )
452
453	IF( sdjf%ln_tint ) THEN ; sdjf%nswap_sec = nmonth_half(irec) + nsec1jan000 ! swap at the middle of the month
454	ELSE ; sdjf%nswap_sec = nmonth_end (irec) + nsec1jan000 ! swap at the end of the month
455	ENDIF
456
457	sdjf%nrec_a(:) = (/ irec, nmonth_half(irec) + nsec1jan000 /) ! define after record number and time
458	irec = irec - 1 ! move back to previous record
459	sdjf%nrec_b(:) = (/ irec, nmonth_half(irec) + nsec1jan000 /) ! define before record number and time
460	!
461	ELSE ! higher frequency mean (in hours)
462	!
463	ifreq_sec = sdjf%nfreqh * 3600 ! frequency mean (in seconds)
464	! number of second since the beginning of the file
465	IF( sdjf%cltype == 'monthly' ) THEN ; ztmp = REAL(nsec_month,wp) ! since 00h on the 1st day of the current month
466	ELSEIF( sdjf%cltype == 'daily' ) THEN ; ztmp = REAL(nsec_day ,wp) ! since 00h of the current day
467	ELSE ; ztmp = REAL(nsec_year ,wp) ! since 00h on Jan 1 of the current year
468	ENDIF
469	ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdttra(1) ! shift time to be centrered in the middle of sbc time step
470	ztmp = ztmp + 0.01 * rdttra(1) ! add 0.01 time step to avoid truncation error
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 / REAL(ifreq_sec, wp) + 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 / REAL(ifreq_sec, wp)
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 == 'daily' ) & ! add the number of seconds between 00h Jan 1 and the end of previous day
508	sdjf%nrec_a(2) = sdjf%nrec_a(2) + isecd * ( nday_year - 1 )
509
510	! before record index and second since Jan. 1st 00h of nit000 year
511	irec = irec - 1. ! move back to previous record
512	sdjf%nrec_b(:) = (/ irec, ifreq_sec * irec - ifreq_sec / 2 + nsec1jan000 /)
513	IF( sdjf%cltype == 'monthly' ) & ! add the number of seconds between 00h Jan 1 and the end of previous month
514	sdjf%nrec_b(2) = sdjf%nrec_b(2) + isecd * SUM(nmonth_len(1:nmonth -1)) ! ok if nmonth=1
515	IF( sdjf%cltype == 'daily' ) & ! add the number of seconds between 00h Jan 1 and the end of previous day
516	sdjf%nrec_b(2) = sdjf%nrec_b(2) + isecd * ( nday_year - 1 )
517
518	! swapping time in second since Jan. 1st 00h of nit000 year
519	IF( sdjf%ln_tint ) THEN ; sdjf%nswap_sec = sdjf%nrec_a(2) ! swap at the middle of the record
520	ELSE ; sdjf%nswap_sec = sdjf%nrec_a(2) + ifreq_sec / 2 ! swap at the end of the record
521	ENDIF
522	!
523	ENDIF
524	!
525	END SUBROUTINE fld_rec
526
527
528	SUBROUTINE fld_clopn( sdjf, kyear, kmonth, kday, ldstop )
529	!!---------------------------------------------------------------------
530	!! * ROUTINE fld_clopn *
531	!!
532	!! ** Purpose : update the file name and open the file
533	!!
534	!! ** Method :
535	!!----------------------------------------------------------------------
536	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
537	INTEGER , INTENT(in ) :: kyear ! year value
538	INTEGER , INTENT(in ) :: kmonth ! month value
539	INTEGER , INTENT(in ) :: kday ! day value
540	LOGICAL , INTENT(in ), OPTIONAL :: ldstop ! stop if open to read a non-existing file (default = .TRUE.)
541
542	IF( sdjf%num /= 0 ) CALL iom_close( sdjf%num ) ! close file if already open
543	! build the new filename if not climatological data
544	IF( .NOT. sdjf%ln_clim ) THEN ; WRITE(sdjf%clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), kyear ! add year
545	IF( sdjf%cltype /= 'yearly' ) WRITE(sdjf%clname, '(a,"m" ,i2.2)' ) TRIM( sdjf%clname ), kmonth ! add month
546	IF( sdjf%cltype == 'daily' ) WRITE(sdjf%clname, '(a,"d" ,i2.2)' ) TRIM( sdjf%clname ), kday ! add day
547	ENDIF
548	CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 )
549	!
550	END SUBROUTINE fld_clopn
551
552
553	SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )
554	!!---------------------------------------------------------------------
555	!! * ROUTINE fld_fill *
556	!!
557	!! ** Purpose : fill sdf with sdf_n and control print
558	!!
559	!! ** Method :
560	!!----------------------------------------------------------------------
561	TYPE(FLD) , DIMENSION(:), INTENT(inout) :: sdf ! structure of input fields (file informations, fields read)
562	TYPE(FLD_N), DIMENSION(:), INTENT(in ) :: sdf_n ! array of namelist information structures
563	CHARACTER(len=*) , INTENT(in ) :: cdir ! Root directory for location of flx files
564	CHARACTER(len=*) , INTENT(in ) :: cdcaller !
565	CHARACTER(len=*) , INTENT(in ) :: cdtitle !
566	CHARACTER(len=*) , INTENT(in ) :: cdnam !
567	!
568	INTEGER :: jf ! dummy indices
569	!!---------------------------------------------------------------------
570
571	DO jf = 1, SIZE(sdf)
572	sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
573	sdf(jf)%nfreqh = sdf_n(jf)%nfreqh
574	sdf(jf)%clvar = sdf_n(jf)%clvar
575	sdf(jf)%ln_tint = sdf_n(jf)%ln_tint
576	sdf(jf)%ln_clim = sdf_n(jf)%ln_clim
577	IF( sdf(jf)%nfreqh == -1. ) THEN ; sdf(jf)%cltype = 'yearly'
578	ELSE ; sdf(jf)%cltype = sdf_n(jf)%cltype
579	ENDIF
580	sdf(jf)%wgtname = " "
581	IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
582	sdf(jf)%vcomp = sdf_n(jf)%vcomp
583	END DO
584
585	IF(lwp) THEN ! control print
586	WRITE(numout,*)
587	WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
588	WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
589	WRITE(numout,*) ' '//TRIM( cdnam )//' Namelist'
590	WRITE(numout,*) ' list of files and frequency (>0: in hours ; <0 in months)'
591	DO jf = 1, SIZE(sdf)
592	WRITE(numout,*) ' root filename: ' , TRIM( sdf(jf)%clrootname ), &
593	& ' variable name: ' , TRIM( sdf(jf)%clvar )
594	WRITE(numout,*) ' frequency: ' , sdf(jf)%nfreqh , &
595	& ' time interp: ' , sdf(jf)%ln_tint , &
596	& ' climatology: ' , sdf(jf)%ln_clim , &
597	& ' weights : ' , TRIM( sdf(jf)%wgtname ), &
598	& ' pairing : ' , TRIM( sdf(jf)%vcomp ), &
599	& ' data type: ' , sdf(jf)%cltype
600	END DO
601	ENDIF
602
603	END SUBROUTINE fld_fill
604
605
606	SUBROUTINE wgt_list( sd, kwgt )
607	!!---------------------------------------------------------------------
608	!! * ROUTINE wgt_list *
609	!!
610	!! ** Purpose : search array of WGTs and find a weights file
611	!! entry, or return a new one adding it to the end
612	!! if it is a new entry, the weights data is read in and
613	!! restructured (fld_weight)
614	!!
615	!! ** Method :
616	!!----------------------------------------------------------------------
617	TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file
618	INTEGER, INTENT(inout) :: kwgt ! index of weights
619	!!
620	INTEGER :: kw
621	INTEGER :: nestid
622	LOGICAL :: found
623	!!----------------------------------------------------------------------
624	!
625	!! search down linked list
626	!! weights filename is either present or we hit the end of the list
627	found = .FALSE.
628
629	!! because agrif nest part of filenames are now added in iom_open
630	!! to distinguish between weights files on the different grids, need to track
631	!! nest number explicitly
632	nestid = 0
633	#if defined key_agrif
634	nestid = Agrif_Fixed()
635	#endif
636	DO kw = 1, nxt_wgt-1
637	IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
638	ref_wgts(kw)%nestid == nestid) THEN
639	kwgt = kw
640	found = .TRUE.
641	EXIT
642	ENDIF
643	END DO
644	IF( .NOT.found ) THEN
645	kwgt = nxt_wgt
646	CALL fld_weight( sd )
647	ENDIF
648
649	END SUBROUTINE wgt_list
650
651	SUBROUTINE wgt_print( )
652	!!---------------------------------------------------------------------
653	!! * ROUTINE wgt_print *
654	!!
655	!! ** Purpose : print the list of known weights
656	!!
657	!! ** Method :
658	!!----------------------------------------------------------------------
659	!!
660	INTEGER :: kw
661	!!----------------------------------------------------------------------
662	!
663
664	DO kw = 1, nxt_wgt-1
665	WRITE(numout,*) 'weight file: ',TRIM(ref_wgts(kw)%wgtname)
666	WRITE(numout,*) ' ddims: ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
667	WRITE(numout,*) ' numwgt: ',ref_wgts(kw)%numwgt
668	WRITE(numout,*) ' jpiwgt: ',ref_wgts(kw)%jpiwgt
669	WRITE(numout,*) ' jpjwgt: ',ref_wgts(kw)%jpjwgt
670	WRITE(numout,*) ' botleft: ',ref_wgts(kw)%botleft
671	WRITE(numout,*) ' topright: ',ref_wgts(kw)%topright
672	IF( ref_wgts(kw)%cyclic ) THEN
673	WRITE(numout,*) ' cyclical'
674	IF( ref_wgts(kw)%offset > 0 ) WRITE(numout,*) ' with offset'
675	ELSE
676	WRITE(numout,*) ' not cyclical'
677	ENDIF
678	IF( ASSOCIATED(ref_wgts(kw)%data_wgt) ) WRITE(numout,*) ' allocated'
679	END DO
680
681	END SUBROUTINE wgt_print
682
683	SUBROUTINE fld_weight( sd )
684	!!---------------------------------------------------------------------
685	!! * ROUTINE fld_weight *
686	!!
687	!! ** Purpose : create a new WGT structure and fill in data from
688	!! file, restructuring as required
689	!!
690	!! ** Method :
691	!!----------------------------------------------------------------------
692	TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file
693	!!
694	INTEGER :: jn ! dummy loop indices
695	INTEGER :: inum ! temporary logical unit
696	INTEGER :: id ! temporary variable id
697	CHARACTER (len=5) :: aname
698	INTEGER , DIMENSION(3) :: ddims
699	INTEGER , DIMENSION(jpi, jpj) :: data_src
700	REAL(wp), DIMENSION(jpi, jpj) :: data_tmp
701	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: line2, lines ! temporary array to read 2 lineumns
702	CHARACTER (len=34) :: lonvar
703	LOGICAL :: cyclical
704	REAL(wp) :: resid, dlon ! temporary array to read 2 lineumns
705	INTEGER :: offset ! temporary integer
706	!!----------------------------------------------------------------------
707	!
708	IF( nxt_wgt > tot_wgts ) THEN
709	CALL ctl_stop("fld_weights: weights array size exceeded, increase tot_wgts")
710	ENDIF
711	!
712	!! new weights file entry, add in extra information
713	!! a weights file represents a 2D grid of a certain shape, so we assume that the current
714	!! input data file is representative of all other files to be opened and processed with the
715	!! current weights file
716
717	!! open input data file (non-model grid)
718	CALL iom_open( sd%clname, inum, ldiof = LEN(TRIM(sd%wgtname)) > 0 )
719
720	!! get dimensions
721	id = iom_varid( inum, sd%clvar, ddims )
722
723	!! check for an east-west cyclic grid
724	!! try to guess name of longitude variable
725
726	lonvar = 'nav_lon'
727	id = iom_varid(inum, TRIM(lonvar), ldstop=.FALSE.)
728	IF( id <= 0 ) THEN
729	lonvar = 'lon'
730	id = iom_varid(inum, TRIM(lonvar), ldstop=.FALSE.)
731	ENDIF
732
733	offset = -1
734	cyclical = .FALSE.
735	IF( id > 0 ) THEN
736	!! found a longitude variable
737	!! now going to assume that grid is regular so we can read a single row
738
739	!! because input array is 2d, have to present iom with 2d array even though we only need 1d slice
740	!! worse, we cant pass line2(:,1) to iom_get since this is treated as a 1d array which doesnt match input file
741	ALLOCATE( lines(ddims(1),2) )
742	CALL iom_get(inum, jpdom_unknown, lonvar, lines(:,:), 1, kstart=(/1,1/), kcount=(/ddims(1),2/) )
743
744	!! find largest grid spacing
745	lines(1:ddims(1)-1,2) = lines(2:ddims(1),1) - lines(1:ddims(1)-1,1)
746	dlon = MAXVAL( lines(1:ddims(1)-1,2) )
747
748	resid = ABS(ABS(lines(ddims(1),1)-lines(1,1))-360.0)
749	IF( resid < rsmall ) THEN
750	!! end rows overlap in longitude
751	offset = 0
752	cyclical = .TRUE.
753	ELSEIF( resid < 2.0*dlon ) THEN
754	!! also call it cyclic if difference between end points is less than twice dlon from 360
755	offset = 1
756	cyclical = .TRUE.
757	ENDIF
758
759	DEALLOCATE( lines )
760
761	ELSE
762	!! guessing failed
763	!! read in first and last columns of data variable
764	!! since we dont know the name of the longitude variable (or even if there is one)
765	!! we assume that if these two columns are equal, file is cyclic east-west
766
767	!! because input array is 2d, have to present iom with 2d array even though we only need 1d slice
768	!! worse, we cant pass line2(1,:) to iom_get since this is treated as a 1d array which doesnt match input file
769	ALLOCATE( lines(2,ddims(2)), line2(2,ddims(2)) )
770	CALL iom_get(inum, jpdom_unknown, sd%clvar, line2(:,:), 1, kstart=(/1,1/), kcount=(/2,ddims(2)/) )
771	lines(2,:) = line2(1,:)
772
773	CALL iom_get(inum, jpdom_unknown, sd%clvar, line2(:,:), 1, kstart=(/ddims(1)-1,1/), kcount=(/2,ddims(2)/) )
774	lines(1,:) = line2(2,:)
775
776	resid = SUM( ABS(lines(1,:) - lines(2,:)) )
777	IF( resid < ddims(2)*rsmall ) THEN
778	offset = 0
779	cyclical = .TRUE.
780	ENDIF
781
782	DEALLOCATE( lines, line2 )
783	ENDIF
784
785	!! close it
786	CALL iom_close( inum )
787
788	!! now open the weights file
789
790	CALL iom_open ( sd%wgtname, inum ) ! interpolation weights
791	IF ( inum > 0 ) THEN
792
793	ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
794	ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
795	ref_wgts(nxt_wgt)%wgtname = sd%wgtname
796	ref_wgts(nxt_wgt)%offset = -1
797	ref_wgts(nxt_wgt)%cyclic = .FALSE.
798	IF( cyclical ) THEN
799	ref_wgts(nxt_wgt)%offset = offset
800	ref_wgts(nxt_wgt)%cyclic = .TRUE.
801	ENDIF
802	ref_wgts(nxt_wgt)%nestid = 0
803	#if defined key_agrif
804	ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
805	#endif
806	!! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
807	!! for each weight wgtNN there is an integer array srcNN which gives the point in
808	!! the input data grid which is to be multiplied by the weight
809	!! they are both arrays on the model grid so the result of the multiplication is
810	!! added into an output array on the model grid as a running sum
811
812	!! two possible cases: bilinear (4 weights) or bicubic (16 weights)
813	id = iom_varid(inum, 'src05', ldstop=.FALSE.)
814	IF( id <= 0) THEN
815	ref_wgts(nxt_wgt)%numwgt = 4
816	ELSE
817	ref_wgts(nxt_wgt)%numwgt = 16
818	ENDIF
819
820	ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
821	ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
822	ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
823
824	DO jn = 1,4
825	aname = ' '
826	WRITE(aname,'(a3,i2.2)') 'src',jn
827	data_tmp(:,:) = 0
828	CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
829	data_src(:,:) = INT(data_tmp(:,:))
830	ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
831	ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
832	END DO
833
834	DO jn = 1, ref_wgts(nxt_wgt)%numwgt
835	aname = ' '
836	WRITE(aname,'(a3,i2.2)') 'wgt',jn
837	ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
838	CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
839	END DO
840	CALL iom_close (inum)
841
842	! find min and max indices in grid
843	ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
844	ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
845	ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
846	ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
847
848	! and therefore dimensions of the input box
849	ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
850	ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
851
852	! shift indexing of source grid
853	ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
854	ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
855
856	! create input grid, give it a halo to allow gradient calculations
857	! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
858	! a more robust solution will be given in next release
859	ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3) )
860	IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col2(2,ref_wgts(nxt_wgt)%jpjwgt+3) )
861
862	nxt_wgt = nxt_wgt + 1
863
864	ELSE
865	CALL ctl_stop( ' fld_weight : unable to read the file ' )
866	ENDIF
867
868	END SUBROUTINE fld_weight
869
870	SUBROUTINE fld_interp(num, clvar, kw, dta, nrec)
871	!!---------------------------------------------------------------------
872	!! * ROUTINE fld_interp *
873	!!
874	!! ** Purpose : apply weights to input gridded data to create data
875	!! on model grid
876	!!
877	!! ** Method :
878	!!----------------------------------------------------------------------
879	INTEGER, INTENT(in) :: num ! stream number
880	CHARACTER(LEN=*), INTENT(in) :: clvar ! variable name
881	INTEGER, INTENT(in) :: kw ! weights number
882	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: dta ! output field on model grid
883	INTEGER, INTENT(in) :: nrec ! record number to read (ie time slice)
884	!!
885	INTEGER, DIMENSION(2) :: rec1,recn ! temporary arrays for start and length
886	INTEGER :: jk, jn, jm ! loop counters
887	INTEGER :: ni, nj ! lengths
888	INTEGER :: jpimin,jpiwid ! temporary indices
889	INTEGER :: jpjmin,jpjwid ! temporary indices
890	INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices
891	!!----------------------------------------------------------------------
892	!
893
894	!! for weighted interpolation we have weights at four corners of a box surrounding
895	!! a model grid point, each weight is multiplied by a grid value (bilinear case)
896	!! or by a grid value and gradients at the corner point (bicubic case)
897	!! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
898
899	!! sub grid where we already have weights
900	jpimin = ref_wgts(kw)%botleft(1)
901	jpjmin = ref_wgts(kw)%botleft(2)
902	jpiwid = ref_wgts(kw)%jpiwgt
903	jpjwid = ref_wgts(kw)%jpjwgt
904
905	!! what we need to read into sub grid in order to calculate gradients
906	rec1(1) = MAX( jpimin-1, 1 )
907	rec1(2) = MAX( jpjmin-1, 1 )
908	recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
909	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
910
911	!! where we need to read it to
912	jpi1 = 2 + rec1(1) - jpimin
913	jpj1 = 2 + rec1(2) - jpjmin
914	jpi2 = jpi1 + recn(1) - 1
915	jpj2 = jpj1 + recn(2) - 1
916
917	ref_wgts(kw)%fly_dta(:,:) = 0.0
918	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2), nrec, rec1, recn)
919
920	!! first four weights common to both bilinear and bicubic
921	!! note that we have to offset by 1 into fly_dta array because of halo
922	dta(:,:) = 0.0
923	DO jk = 1,4
924	DO jn = 1, jpj
925	DO jm = 1,jpi
926	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
927	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
928	dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1)
929	END DO
930	END DO
931	END DO
932
933	IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
934
935	!! fix up halo points that we couldnt read from file
936	IF( jpi1 == 2 ) THEN
937	ref_wgts(kw)%fly_dta(jpi1-1,:) = ref_wgts(kw)%fly_dta(jpi1,:)
938	ENDIF
939	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
940	ref_wgts(kw)%fly_dta(jpi2+1,:) = ref_wgts(kw)%fly_dta(jpi2,:)
941	ENDIF
942	IF( jpj1 == 2 ) THEN
943	ref_wgts(kw)%fly_dta(:,jpj1-1) = ref_wgts(kw)%fly_dta(:,jpj1)
944	ENDIF
945	IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
946	ref_wgts(kw)%fly_dta(:,jpj2+1) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2) - ref_wgts(kw)%fly_dta(:,jpj2-1)
947	ENDIF
948
949	!! if data grid is cyclic we can do better on east-west edges
950	!! but have to allow for whether first and last columns are coincident
951	IF( ref_wgts(kw)%cyclic ) THEN
952	rec1(2) = MAX( jpjmin-1, 1 )
953	recn(1) = 2
954	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
955	jpj1 = 2 + rec1(2) - jpjmin
956	jpj2 = jpj1 + recn(2) - 1
957	IF( jpi1 == 2 ) THEN
958	rec1(1) = ref_wgts(kw)%ddims(1) - 1
959	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2), nrec, rec1, recn)
960	ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2) = ref_wgts(kw)%col2(ref_wgts(kw)%offset+1,jpj1:jpj2)
961	ENDIF
962	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
963	rec1(1) = 1
964	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2), nrec, rec1, recn)
965	ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2) = ref_wgts(kw)%col2(2-ref_wgts(kw)%offset,jpj1:jpj2)
966	ENDIF
967	ENDIF
968
969	! gradient in the i direction
970	DO jk = 1,4
971	DO jn = 1, jpj
972	DO jm = 1,jpi
973	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
974	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
975	dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 * &
976	(ref_wgts(kw)%fly_dta(ni+2,nj+1) - ref_wgts(kw)%fly_dta(ni,nj+1))
977	END DO
978	END DO
979	END DO
980
981	! gradient in the j direction
982	DO jk = 1,4
983	DO jn = 1, jpj
984	DO jm = 1,jpi
985	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
986	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
987	dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 * &
988	(ref_wgts(kw)%fly_dta(ni+1,nj+2) - ref_wgts(kw)%fly_dta(ni+1,nj))
989	END DO
990	END DO
991	END DO
992
993	! gradient in the ij direction
994	DO jk = 1,4
995	DO jn = 1, jpj
996	DO jm = 1,jpi
997	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
998	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
999	dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
1000	(ref_wgts(kw)%fly_dta(ni+2,nj+2) - ref_wgts(kw)%fly_dta(ni ,nj+2)) - &
1001	(ref_wgts(kw)%fly_dta(ni+2,nj ) - ref_wgts(kw)%fly_dta(ni ,nj )))
1002	END DO
1003	END DO
1004	END DO
1005
1006	END IF
1007
1008	END SUBROUTINE fld_interp
1009
1010	END MODULE fldread

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

Download in other formats: