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

Context Navigation

source: tags/nemo_v3_2_2/NEMO/OPA_SRC/SBC/fldread.F90 @ 2700

Last change on this file since 2700 was 2700, checked in by smasson, 13 years ago
bugfix for fldread in v3_2_2, see ticket#810
Property svn:keywords set to `Id`
File size: 52.8 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	!
136	! ! ===================== !
137	DO jf = 1, imf ! LOOP OVER FIELD !
138	! ! ===================== !
139	!
140	! centrered in the middle of sbc time step if doing time interpolation
141	isecsbc = nsec_year + nsec1jan000 + COUNT((/ sd(jf)%ln_tint /)) * NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdttra(1))
142	!
143	IF( kt == nit000 ) CALL fld_init( kn_fsbc, sd(jf) )
144	!
145	! read/update the after data?
146	IF( isecsbc > sd(jf)%nswap_sec ) THEN
147
148	IF( sd(jf)%ln_tint ) THEN ! time interpolation: swap before record field
149	!CDIR COLLAPSE
150	sd(jf)%fdta(:,:,1) = sd(jf)%fdta(:,:,2)
151	sd(jf)%rotn(1) = sd(jf)%rotn(2)
152	ENDIF
153
154	! update record informations
155	CALL fld_rec( kn_fsbc, sd(jf) )
156
157	! do we have to change the year/month/day of the forcing field??
158	IF( sd(jf)%ln_tint ) THEN
159	! if we do time interpolation we will need to open next year/month/day file before the end of the current one
160	! 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
161	! larger than the record number that should be read for current year/month/day (for ex. 13 for monthly mean file)
162
163	! last record to be read in the current file
164	IF( sd(jf)%nfreqh == -1 ) THEN ; ireclast = 12
165	ELSE
166	IF( sd(jf)%cltype == 'monthly' ) THEN ; ireclast = 24 * nmonth_len(nmonth) / sd(jf)%nfreqh
167	ELSEIF( sd(jf)%cltype == 'daily' ) THEN ; ireclast = 24 / sd(jf)%nfreqh
168	ELSE ; ireclast = 24 * nyear_len( 1 ) / sd(jf)%nfreqh
169	ENDIF
170	ENDIF
171
172	! do we need next file data?
173	IF( sd(jf)%nrec_a(1) > ireclast ) THEN
174
175	sd(jf)%nrec_a(1) = 1 ! force to read the first record of the next file
176
177	IF( .NOT. sd(jf)%ln_clim ) THEN ! close the current file and open a new one.
178
179	llnxtmth = sd(jf)%cltype == 'monthly' .OR. nday == nmonth_len(nmonth) ! open next month file?
180	llnxtyr = sd(jf)%cltype == 'yearly' .OR. (nmonth == 12 .AND. llnxtmth) ! open next year file?
181
182	! if the run finishes at the end of the current year/month/day, we will allow next year/month/day file to be
183	! not present. If the run continue further than the current year/month/day, next year/month/day file must exist
184	isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdttra(1)) ! second at the end of the run
185	llstop = isecend > sd(jf)%nswap_sec ! read more than 1 record of next year
186
187	CALL fld_clopn( sd(jf), nyear + COUNT((/llnxtyr /)) , &
188	& nmonth + COUNT((/llnxtmth/)) - 12 * COUNT((/llnxtyr /)), &
189	& nday + 1 - nmonth_len(nmonth) * COUNT((/llnxtmth/)), llstop )
190
191	IF( sd(jf)%num <= 0 .AND. .NOT. llstop ) THEN ! next year file does not exist
192	CALL ctl_warn('next year/month/day file: '//TRIM(sd(jf)%clname)// &
193	& ' not present -> back to current year/month/day')
194	CALL fld_clopn( sd(jf), nyear, nmonth, nday ) ! back to the current year/month/day
195	sd(jf)%nrec_a(1) = ireclast ! force to read the last record to be read in the current year file
196	ENDIF
197
198	ENDIF
199	ENDIF
200
201	ELSE
202	! if we are not doing time interpolation, we must change the year/month/day of the file just after switching
203	! to the NEW year/month/day. If it is the case, we are at the beginning of the year/month/day when calling
204	! fld_rec so sd(jf)%nrec_a(1) = 1
205	IF( sd(jf)%nrec_a(1) == 1 .AND. .NOT. sd(jf)%ln_clim ) CALL fld_clopn( sd(jf), nyear, nmonth, nday )
206	ENDIF
207
208	! read after data
209	IF( LEN(TRIM(sd(jf)%wgtname)) > 0 ) THEN
210	CALL wgt_list( sd(jf), kw )
211	CALL fld_interp( sd(jf)%num, sd(jf)%clvar, kw, sd(jf)%fdta(:,:,2), sd(jf)%nrec_a(1) )
212	ELSE
213	CALL iom_get( sd(jf)%num, jpdom_data, sd(jf)%clvar, sd(jf)%fdta(:,:,2), sd(jf)%nrec_a(1) )
214	ENDIF
215	sd(jf)%rotn(2) = .FALSE.
216
217	ENDIF
218	! ! ===================== !
219	END DO ! END LOOP OVER FIELD !
220	! ! ===================== !
221
222	IF( kt == nit000 .AND. lwp ) CALL wgt_print()
223
224	!! Vector fields may need to be rotated onto the local grid direction
225	!! This has to happen before the time interpolations
226	!! (sga: following code should be modified so that pairs arent searched for each time
227
228	DO jf = 1, imf
229	!! find vector rotations required
230	IF( LEN(TRIM(sd(jf)%vcomp)) > 0 ) THEN
231	!! east-west component has symbolic name starting with 'U'
232	IF( sd(jf)%vcomp(1:1) == 'U' ) THEN
233	!! found an east-west component, look for the north-south component
234	!! which has same symbolic name but with 'U' replaced with 'V'
235	nf = LEN_TRIM( sd(jf)%vcomp )
236	IF( nf == 1) THEN
237	acomp = 'V'
238	ELSE
239	acomp = 'V' // sd(jf)%vcomp(2:nf)
240	ENDIF
241	kf = -1
242	DO nf = 1, imf
243	IF( TRIM(sd(nf)%vcomp) == TRIM(acomp) ) kf = nf
244	END DO
245	IF( kf > 0 ) THEN
246	!! fields jf,kf are two components which need to be rotated together
247	DO nf = 1,2
248	!! check each time level of this pair
249	IF( .NOT. sd(jf)%rotn(nf) .AND. .NOT. sd(kf)%rotn(nf) ) THEN
250	utmp(:,:) = 0.0
251	vtmp(:,:) = 0.0
252	CALL rot_rep( sd(jf)%fdta(:,:,nf), sd(kf)%fdta(:,:,nf), 'T', 'en->i', utmp(:,:) )
253	CALL rot_rep( sd(jf)%fdta(:,:,nf), sd(kf)%fdta(:,:,nf), 'T', 'en->j', vtmp(:,:) )
254	sd(jf)%fdta(:,:,nf) = utmp(:,:)
255	sd(kf)%fdta(:,:,nf) = vtmp(:,:)
256	sd(jf)%rotn(nf) = .TRUE.
257	sd(kf)%rotn(nf) = .TRUE.
258	IF( lwp .AND. kt == nit000 ) &
259	WRITE(numout,*) 'fld_read: vector pair (', &
260	TRIM(sd(jf)%clvar),',',TRIM(sd(kf)%clvar), &
261	') rotated on to model grid'
262	ENDIF
263	END DO
264	ENDIF
265	ENDIF
266	ENDIF
267	END DO
268
269	isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdttra(1)) ! used only for temporal interpolation
270
271	! ! ===================== !
272	DO jf = 1, imf ! LOOP OVER FIELD !
273	! ! ===================== !
274	!
275	! update field at each kn_fsbc time-step
276	IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
277	!
278	IF( sd(jf)%ln_tint ) THEN
279	IF(lwp .AND. kt - nit000 <= 100 ) THEN
280	clfmt = "('fld_read: var ', a, ' kt = ', i8,' Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," // &
281	& "' records b/a: ', i4.4, '/', i4.4, ' (', f7.2,'/', f7.2, ' days)')"
282	WRITE(numout, clfmt) TRIM( sd(jf)%clvar ), kt, nyear, nmonth, nday, &
283	& 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
284	ENDIF
285	!
286	ztinta = REAL( isecsbc - sd(jf)%nrec_b(2), wp ) / REAL( sd(jf)%nrec_a(2) - sd(jf)%nrec_b(2), wp )
287	ztintb = 1. - ztinta
288	!CDIR COLLAPSE
289	sd(jf)%fnow(:,:) = ztintb * sd(jf)%fdta(:,:,1) + ztinta * sd(jf)%fdta(:,:,2)
290	ELSE
291	IF(lwp .AND. kt - nit000 <= 100 ) THEN
292	clfmt = "('fld_read: var ', a, ' kt = ', i8,' Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," // &
293	& "' record: ', i4.4, ' at ', f7.2, ' day')"
294	WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, nyear, nmonth, nday, sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_a(2),wp)/rday
295	ENDIF
296	!CDIR COLLAPSE
297	sd(jf)%fnow(:,:) = sd(jf)%fdta(:,:,2) ! piecewise constant field
298
299	ENDIF
300	!
301	ENDIF
302
303	IF( kt == nitend ) CALL iom_close( sd(jf)%num ) ! Close the input files
304
305	! ! ===================== !
306	END DO ! END LOOP OVER FIELD !
307	! ! ===================== !
308	END SUBROUTINE fld_read
309
310
311	SUBROUTINE fld_init( kn_fsbc, sdjf )
312	!!---------------------------------------------------------------------
313	!! * ROUTINE fld_init *
314	!!
315	!! ** Purpose : - if time interpolation, read before data
316	!! - open current year file
317	!!
318	!! ** Method :
319	!!----------------------------------------------------------------------
320	INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step)
321	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
322	!!
323	LOGICAL :: llprevyr ! are we reading previous year file?
324	LOGICAL :: llprevmth ! are we reading previous month file?
325	LOGICAL :: llprevday ! are we reading previous day file?
326	LOGICAL :: llprev ! llprevyr .OR. llprevmth .OR. llprevday
327	INTEGER :: idvar ! variable id
328	INTEGER :: inrec ! number of record existing for this variable
329	INTEGER :: kwgt
330	CHARACTER(LEN=1000) :: clfmt ! write format
331	!!---------------------------------------------------------------------
332
333	! some default definitions...
334	sdjf%num = 0 ! default definition for non-opened file
335	IF( sdjf%ln_clim ) sdjf%clname = TRIM( sdjf%clrootname ) ! file name defaut definition, never change in this case
336	llprevyr = .FALSE.
337	llprevmth = .FALSE.
338	llprevday = .FALSE.
339
340	! define record informations
341	CALL fld_rec( kn_fsbc, sdjf )
342
343	! Note: shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
344
345	IF( sdjf%ln_tint ) THEN ! we need to read the previous record and we will put it in the current record structure
346
347	IF( sdjf%nrec_b(1) == 0 ) THEN ! we redefine record sdjf%nrec_b(1) with the last record of previous year file
348	IF( sdjf%nfreqh == -1 ) THEN ! monthly mean
349	IF( sdjf%cltype == 'monthly' ) THEN ! monthly file
350	sdjf%nrec_b(1) = 1 ! force to read the unique record
351	llprevmth = .NOT. sdjf%ln_clim ! use previous month file?
352	llprevyr = llprevmth .AND. nmonth == 1 ! use previous year file?
353	ELSE ! yearly file
354	sdjf%nrec_b(1) = 12 ! force to read december mean
355	llprevyr = .NOT. sdjf%ln_clim ! use previous year file?
356	ENDIF
357	ELSE
358	IF( sdjf%cltype == 'monthly' ) THEN ! monthly file
359	sdjf%nrec_b(1) = 24 * nmonth_len(nmonth-1) / sdjf%nfreqh ! last record of previous month
360	llprevmth = .NOT. sdjf%ln_clim ! use previous month file?
361	llprevyr = llprevmth .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. 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( kn_fsbc, sdjf )
419	!!---------------------------------------------------------------------
420	!! * ROUTINE fld_rec *
421	!!
422	!! ** Purpose : compute nrec_a, nrec_b and nswap_sec
423	!!
424	!! ** Method :
425	!!----------------------------------------------------------------------
426	INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step)
427	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
428	!!
429	INTEGER :: irec ! record number
430	INTEGER :: isecd ! rday
431	REAL(wp) :: ztmp ! temporary variable
432	INTEGER :: ifreq_sec ! frequency mean (in seconds)
433	!!----------------------------------------------------------------------
434	!
435	! Note: shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
436	!
437	IF( sdjf%nfreqh == -1 ) THEN ! monthly mean
438	!
439	IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record
440	!
441	! INT( ztmp )
442	! /\|\
443	! 1 \| *----
444	! 0 \|----(
445	! \|----+----\|--> time
446	! 0 /\|\ 1 (nday/nmonth_len(nmonth))
447	! \|
448	! \|
449	! forcing record : nmonth
450	!
451	ztmp = REAL( nday, wp ) / REAL( nmonth_len(nmonth), wp ) + 0.5
452	ELSE
453	ztmp = 0.e0
454	ENDIF
455	irec = nmonth + INT( ztmp )
456
457	IF( sdjf%ln_tint ) THEN ; sdjf%nswap_sec = nmonth_half(irec) + nsec1jan000 ! swap at the middle of the month
458	ELSE ; sdjf%nswap_sec = nmonth_end (irec) + nsec1jan000 ! swap at the end of the month
459	ENDIF
460
461	sdjf%nrec_a(:) = (/ irec, nmonth_half(irec) + nsec1jan000 /) ! define after record number and time
462	irec = irec - 1 ! move back to previous record
463	sdjf%nrec_b(:) = (/ irec, nmonth_half(irec) + nsec1jan000 /) ! define before record number and time
464	!
465	ELSE ! higher frequency mean (in hours)
466	!
467	ifreq_sec = sdjf%nfreqh * 3600 ! frequency mean (in seconds)
468	! number of second since the beginning of the file
469	IF( sdjf%cltype == 'monthly' ) THEN ; ztmp = REAL(nsec_month,wp) ! since 00h on the 1st day of the current month
470	ELSEIF( sdjf%cltype == 'daily' ) THEN ; ztmp = REAL(nsec_day ,wp) ! since 00h of the current day
471	ELSE ; ztmp = REAL(nsec_year ,wp) ! since 00h on Jan 1 of the current year
472	ENDIF
473	ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdttra(1) ! shift time to be centrered in the middle of sbc time step
474	ztmp = ztmp + 0.01 * rdttra(1) ! add 0.01 time step to avoid truncation error
475	IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record
476	!
477	! INT( ztmp )
478	! /\|\
479	! 2 \| *-----(
480	! 1 \| *-----(
481	! 0 \|--(
482	! \|--+--\|--+--\|--+--\|--> time
483	! 0 /\|\ 1 /\|\ 2 /\|\ 3 (nsec_year/ifreq_sec) or (nsec_month/ifreq_sec)
484	! \| \| \|
485	! \| \| \|
486	! forcing record : 1 2 3
487	!
488	ztmp = ztmp / REAL(ifreq_sec, wp) + 0.5
489	ELSE
490	!
491	! INT( ztmp )
492	! /\|\
493	! 2 \| *-----(
494	! 1 \| *-----(
495	! 0 \|-----(
496	! \|--+--\|--+--\|--+--\|--> time
497	! 0 /\|\ 1 /\|\ 2 /\|\ 3 (nsec_year/ifreq_sec) or (nsec_month/ifreq_sec)
498	! \| \| \|
499	! \| \| \|
500	! forcing record : 1 2 3
501	!
502	ztmp = ztmp / REAL(ifreq_sec, wp)
503	ENDIF
504	irec = 1 + INT( ztmp )
505
506	isecd = NINT(rday)
507	! after record index and second since Jan. 1st 00h of nit000 year
508	sdjf%nrec_a(:) = (/ irec, ifreq_sec * irec - ifreq_sec / 2 + nsec1jan000 /)
509	IF( sdjf%cltype == 'monthly' ) & ! add the number of seconds between 00h Jan 1 and the end of previous month
510	sdjf%nrec_a(2) = sdjf%nrec_a(2) + isecd * SUM(nmonth_len(1:nmonth -1)) ! ok if nmonth=1
511	IF( sdjf%cltype == 'daily' ) & ! add the number of seconds between 00h Jan 1 and the end of previous day
512	sdjf%nrec_a(2) = sdjf%nrec_a(2) + isecd * ( nday_year - 1 )
513
514	! before record index and second since Jan. 1st 00h of nit000 year
515	irec = irec - 1. ! move back to previous record
516	sdjf%nrec_b(:) = (/ irec, ifreq_sec * irec - ifreq_sec / 2 + nsec1jan000 /)
517	IF( sdjf%cltype == 'monthly' ) & ! add the number of seconds between 00h Jan 1 and the end of previous month
518	sdjf%nrec_b(2) = sdjf%nrec_b(2) + isecd * SUM(nmonth_len(1:nmonth -1)) ! ok if nmonth=1
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
546	IF( sdjf%num /= 0 ) CALL iom_close( sdjf%num ) ! close file if already open
547	! build the new filename if not climatological data
548	IF( .NOT. sdjf%ln_clim ) THEN ; WRITE(sdjf%clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), kyear ! add year
549	IF( sdjf%cltype /= 'yearly' ) WRITE(sdjf%clname, '(a,"m" ,i2.2)' ) TRIM( sdjf%clname ), kmonth ! add month
550	IF( sdjf%cltype == 'daily' ) WRITE(sdjf%clname, '(a,"d" ,i2.2)' ) TRIM( sdjf%clname ), kday ! add day
551	ENDIF
552	CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 )
553	!
554	END SUBROUTINE fld_clopn
555
556
557	SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )
558	!!---------------------------------------------------------------------
559	!! * ROUTINE fld_fill *
560	!!
561	!! ** Purpose : fill sdf with sdf_n and control print
562	!!
563	!! ** Method :
564	!!----------------------------------------------------------------------
565	TYPE(FLD) , DIMENSION(:), INTENT(inout) :: sdf ! structure of input fields (file informations, fields read)
566	TYPE(FLD_N), DIMENSION(:), INTENT(in ) :: sdf_n ! array of namelist information structures
567	CHARACTER(len=*) , INTENT(in ) :: cdir ! Root directory for location of flx files
568	CHARACTER(len=*) , INTENT(in ) :: cdcaller !
569	CHARACTER(len=*) , INTENT(in ) :: cdtitle !
570	CHARACTER(len=*) , INTENT(in ) :: cdnam !
571	!
572	INTEGER :: jf ! dummy indices
573	!!---------------------------------------------------------------------
574
575	DO jf = 1, SIZE(sdf)
576	sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
577	sdf(jf)%nfreqh = sdf_n(jf)%nfreqh
578	sdf(jf)%clvar = sdf_n(jf)%clvar
579	sdf(jf)%ln_tint = sdf_n(jf)%ln_tint
580	sdf(jf)%ln_clim = sdf_n(jf)%ln_clim
581	IF( sdf(jf)%nfreqh == -1. ) THEN ; sdf(jf)%cltype = 'yearly'
582	ELSE ; sdf(jf)%cltype = sdf_n(jf)%cltype
583	ENDIF
584	sdf(jf)%wgtname = " "
585	IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
586	sdf(jf)%vcomp = sdf_n(jf)%vcomp
587	END DO
588
589	IF(lwp) THEN ! control print
590	WRITE(numout,*)
591	WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
592	WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
593	WRITE(numout,*) ' '//TRIM( cdnam )//' Namelist'
594	WRITE(numout,*) ' list of files and frequency (>0: in hours ; <0 in months)'
595	DO jf = 1, SIZE(sdf)
596	WRITE(numout,*) ' root filename: ' , TRIM( sdf(jf)%clrootname ), &
597	& ' variable name: ' , TRIM( sdf(jf)%clvar )
598	WRITE(numout,*) ' frequency: ' , sdf(jf)%nfreqh , &
599	& ' time interp: ' , sdf(jf)%ln_tint , &
600	& ' climatology: ' , sdf(jf)%ln_clim , &
601	& ' weights : ' , TRIM( sdf(jf)%wgtname ), &
602	& ' pairing : ' , TRIM( sdf(jf)%vcomp ), &
603	& ' data type: ' , sdf(jf)%cltype
604	END DO
605	ENDIF
606
607	END SUBROUTINE fld_fill
608
609
610	SUBROUTINE wgt_list( sd, kwgt )
611	!!---------------------------------------------------------------------
612	!! * ROUTINE wgt_list *
613	!!
614	!! ** Purpose : search array of WGTs and find a weights file
615	!! entry, or return a new one adding it to the end
616	!! if it is a new entry, the weights data is read in and
617	!! restructured (fld_weight)
618	!!
619	!! ** Method :
620	!!----------------------------------------------------------------------
621	TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file
622	INTEGER, INTENT(inout) :: kwgt ! index of weights
623	!!
624	INTEGER :: kw
625	INTEGER :: nestid
626	LOGICAL :: found
627	!!----------------------------------------------------------------------
628	!
629	!! search down linked list
630	!! weights filename is either present or we hit the end of the list
631	found = .FALSE.
632
633	!! because agrif nest part of filenames are now added in iom_open
634	!! to distinguish between weights files on the different grids, need to track
635	!! nest number explicitly
636	nestid = 0
637	#if defined key_agrif
638	nestid = Agrif_Fixed()
639	#endif
640	DO kw = 1, nxt_wgt-1
641	IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
642	ref_wgts(kw)%nestid == nestid) THEN
643	kwgt = kw
644	found = .TRUE.
645	EXIT
646	ENDIF
647	END DO
648	IF( .NOT.found ) THEN
649	kwgt = nxt_wgt
650	CALL fld_weight( sd )
651	ENDIF
652
653	END SUBROUTINE wgt_list
654
655	SUBROUTINE wgt_print( )
656	!!---------------------------------------------------------------------
657	!! * ROUTINE wgt_print *
658	!!
659	!! ** Purpose : print the list of known weights
660	!!
661	!! ** Method :
662	!!----------------------------------------------------------------------
663	!!
664	INTEGER :: kw
665	!!----------------------------------------------------------------------
666	!
667
668	DO kw = 1, nxt_wgt-1
669	WRITE(numout,*) 'weight file: ',TRIM(ref_wgts(kw)%wgtname)
670	WRITE(numout,*) ' ddims: ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
671	WRITE(numout,*) ' numwgt: ',ref_wgts(kw)%numwgt
672	WRITE(numout,*) ' jpiwgt: ',ref_wgts(kw)%jpiwgt
673	WRITE(numout,*) ' jpjwgt: ',ref_wgts(kw)%jpjwgt
674	WRITE(numout,*) ' botleft: ',ref_wgts(kw)%botleft
675	WRITE(numout,*) ' topright: ',ref_wgts(kw)%topright
676	IF( ref_wgts(kw)%cyclic ) THEN
677	WRITE(numout,*) ' cyclical'
678	IF( ref_wgts(kw)%offset > 0 ) WRITE(numout,*) ' with offset'
679	ELSE
680	WRITE(numout,*) ' not cyclical'
681	ENDIF
682	IF( ASSOCIATED(ref_wgts(kw)%data_wgt) ) WRITE(numout,*) ' allocated'
683	END DO
684
685	END SUBROUTINE wgt_print
686
687	SUBROUTINE fld_weight( sd )
688	!!---------------------------------------------------------------------
689	!! * ROUTINE fld_weight *
690	!!
691	!! ** Purpose : create a new WGT structure and fill in data from
692	!! file, restructuring as required
693	!!
694	!! ** Method :
695	!!----------------------------------------------------------------------
696	TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file
697	!!
698	INTEGER :: jn ! dummy loop indices
699	INTEGER :: inum ! temporary logical unit
700	INTEGER :: id ! temporary variable id
701	CHARACTER (len=5) :: aname
702	INTEGER , DIMENSION(3) :: ddims
703	INTEGER , DIMENSION(jpi, jpj) :: data_src
704	REAL(wp), DIMENSION(jpi, jpj) :: data_tmp
705	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: line2, lines ! temporary array to read 2 lineumns
706	CHARACTER (len=34) :: lonvar
707	LOGICAL :: cyclical
708	REAL(wp) :: resid, dlon ! temporary array to read 2 lineumns
709	INTEGER :: offset ! temporary integer
710	!!----------------------------------------------------------------------
711	!
712	IF( nxt_wgt > tot_wgts ) THEN
713	CALL ctl_stop("fld_weights: weights array size exceeded, increase tot_wgts")
714	ENDIF
715	!
716	!! new weights file entry, add in extra information
717	!! a weights file represents a 2D grid of a certain shape, so we assume that the current
718	!! input data file is representative of all other files to be opened and processed with the
719	!! current weights file
720
721	!! open input data file (non-model grid)
722	CALL iom_open( sd%clname, inum, ldiof = LEN(TRIM(sd%wgtname)) > 0 )
723
724	!! get dimensions
725	id = iom_varid( inum, sd%clvar, ddims )
726
727	!! check for an east-west cyclic grid
728	!! try to guess name of longitude variable
729
730	lonvar = 'nav_lon'
731	id = iom_varid(inum, TRIM(lonvar), ldstop=.FALSE.)
732	IF( id <= 0 ) THEN
733	lonvar = 'lon'
734	id = iom_varid(inum, TRIM(lonvar), ldstop=.FALSE.)
735	ENDIF
736
737	offset = -1
738	cyclical = .FALSE.
739	IF( id > 0 ) THEN
740	!! found a longitude variable
741	!! now going to assume that grid is regular so we can read a single row
742
743	!! because input array is 2d, have to present iom with 2d array even though we only need 1d slice
744	!! worse, we cant pass line2(:,1) to iom_get since this is treated as a 1d array which doesnt match input file
745	ALLOCATE( lines(ddims(1),2) )
746	CALL iom_get(inum, jpdom_unknown, lonvar, lines(:,:), 1, kstart=(/1,1/), kcount=(/ddims(1),2/) )
747
748	!! find largest grid spacing
749	lines(1:ddims(1)-1,2) = lines(2:ddims(1),1) - lines(1:ddims(1)-1,1)
750	dlon = MAXVAL( lines(1:ddims(1)-1,2) )
751
752	resid = ABS(ABS(lines(ddims(1),1)-lines(1,1))-360.0)
753	IF( resid < rsmall ) THEN
754	!! end rows overlap in longitude
755	offset = 0
756	cyclical = .TRUE.
757	ELSEIF( resid < 2.0*dlon ) THEN
758	!! also call it cyclic if difference between end points is less than twice dlon from 360
759	offset = 1
760	cyclical = .TRUE.
761	ENDIF
762
763	DEALLOCATE( lines )
764
765	ELSE
766	!! guessing failed
767	!! read in first and last columns of data variable
768	!! since we dont know the name of the longitude variable (or even if there is one)
769	!! we assume that if these two columns are equal, file is cyclic east-west
770
771	!! because input array is 2d, have to present iom with 2d array even though we only need 1d slice
772	!! worse, we cant pass line2(1,:) to iom_get since this is treated as a 1d array which doesnt match input file
773	ALLOCATE( lines(2,ddims(2)), line2(2,ddims(2)) )
774	CALL iom_get(inum, jpdom_unknown, sd%clvar, line2(:,:), 1, kstart=(/1,1/), kcount=(/2,ddims(2)/) )
775	lines(2,:) = line2(1,:)
776
777	CALL iom_get(inum, jpdom_unknown, sd%clvar, line2(:,:), 1, kstart=(/ddims(1)-1,1/), kcount=(/2,ddims(2)/) )
778	lines(1,:) = line2(2,:)
779
780	resid = SUM( ABS(lines(1,:) - lines(2,:)) )
781	IF( resid < ddims(2)*rsmall ) THEN
782	offset = 0
783	cyclical = .TRUE.
784	ENDIF
785
786	DEALLOCATE( lines, line2 )
787	ENDIF
788
789	!! close it
790	CALL iom_close( inum )
791
792	!! now open the weights file
793
794	CALL iom_open ( sd%wgtname, inum ) ! interpolation weights
795	IF ( inum > 0 ) THEN
796
797	ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
798	ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
799	ref_wgts(nxt_wgt)%wgtname = sd%wgtname
800	ref_wgts(nxt_wgt)%offset = -1
801	ref_wgts(nxt_wgt)%cyclic = .FALSE.
802	IF( cyclical ) THEN
803	ref_wgts(nxt_wgt)%offset = offset
804	ref_wgts(nxt_wgt)%cyclic = .TRUE.
805	ENDIF
806	ref_wgts(nxt_wgt)%nestid = 0
807	#if defined key_agrif
808	ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
809	#endif
810	!! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
811	!! for each weight wgtNN there is an integer array srcNN which gives the point in
812	!! the input data grid which is to be multiplied by the weight
813	!! they are both arrays on the model grid so the result of the multiplication is
814	!! added into an output array on the model grid as a running sum
815
816	!! two possible cases: bilinear (4 weights) or bicubic (16 weights)
817	id = iom_varid(inum, 'src05', ldstop=.FALSE.)
818	IF( id <= 0) THEN
819	ref_wgts(nxt_wgt)%numwgt = 4
820	ELSE
821	ref_wgts(nxt_wgt)%numwgt = 16
822	ENDIF
823
824	ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
825	ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
826	ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
827
828	DO jn = 1,4
829	aname = ' '
830	WRITE(aname,'(a3,i2.2)') 'src',jn
831	data_tmp(:,:) = 0
832	CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
833	data_src(:,:) = INT(data_tmp(:,:))
834	ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
835	ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
836	END DO
837
838	DO jn = 1, ref_wgts(nxt_wgt)%numwgt
839	aname = ' '
840	WRITE(aname,'(a3,i2.2)') 'wgt',jn
841	ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
842	CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
843	END DO
844	CALL iom_close (inum)
845
846	! find min and max indices in grid
847	ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
848	ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
849	ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
850	ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
851
852	! and therefore dimensions of the input box
853	ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
854	ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
855
856	! shift indexing of source grid
857	ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
858	ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
859
860	! create input grid, give it a halo to allow gradient calculations
861	! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
862	! a more robust solution will be given in next release
863	ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3) )
864	IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col2(2,ref_wgts(nxt_wgt)%jpjwgt+3) )
865
866	nxt_wgt = nxt_wgt + 1
867
868	ELSE
869	CALL ctl_stop( ' fld_weight : unable to read the file ' )
870	ENDIF
871
872	END SUBROUTINE fld_weight
873
874	SUBROUTINE fld_interp(num, clvar, kw, dta, nrec)
875	!!---------------------------------------------------------------------
876	!! * ROUTINE fld_interp *
877	!!
878	!! ** Purpose : apply weights to input gridded data to create data
879	!! on model grid
880	!!
881	!! ** Method :
882	!!----------------------------------------------------------------------
883	INTEGER, INTENT(in) :: num ! stream number
884	CHARACTER(LEN=*), INTENT(in) :: clvar ! variable name
885	INTEGER, INTENT(in) :: kw ! weights number
886	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: dta ! output field on model grid
887	INTEGER, INTENT(in) :: nrec ! record number to read (ie time slice)
888	!!
889	INTEGER, DIMENSION(2) :: rec1,recn ! temporary arrays for start and length
890	INTEGER :: jk, jn, jm ! loop counters
891	INTEGER :: ni, nj ! lengths
892	INTEGER :: jpimin,jpiwid ! temporary indices
893	INTEGER :: jpjmin,jpjwid ! temporary indices
894	INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices
895	!!----------------------------------------------------------------------
896	!
897
898	!! for weighted interpolation we have weights at four corners of a box surrounding
899	!! a model grid point, each weight is multiplied by a grid value (bilinear case)
900	!! or by a grid value and gradients at the corner point (bicubic case)
901	!! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
902
903	!! sub grid where we already have weights
904	jpimin = ref_wgts(kw)%botleft(1)
905	jpjmin = ref_wgts(kw)%botleft(2)
906	jpiwid = ref_wgts(kw)%jpiwgt
907	jpjwid = ref_wgts(kw)%jpjwgt
908
909	!! what we need to read into sub grid in order to calculate gradients
910	rec1(1) = MAX( jpimin-1, 1 )
911	rec1(2) = MAX( jpjmin-1, 1 )
912	recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
913	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
914
915	!! where we need to read it to
916	jpi1 = 2 + rec1(1) - jpimin
917	jpj1 = 2 + rec1(2) - jpjmin
918	jpi2 = jpi1 + recn(1) - 1
919	jpj2 = jpj1 + recn(2) - 1
920
921	ref_wgts(kw)%fly_dta(:,:) = 0.0
922	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2), nrec, rec1, recn)
923
924	!! first four weights common to both bilinear and bicubic
925	!! note that we have to offset by 1 into fly_dta array because of halo
926	dta(:,:) = 0.0
927	DO jk = 1,4
928	DO jn = 1, jpj
929	DO jm = 1,jpi
930	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
931	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
932	dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1)
933	END DO
934	END DO
935	END DO
936
937	IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
938
939	!! fix up halo points that we couldnt read from file
940	IF( jpi1 == 2 ) THEN
941	ref_wgts(kw)%fly_dta(jpi1-1,:) = ref_wgts(kw)%fly_dta(jpi1,:)
942	ENDIF
943	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
944	ref_wgts(kw)%fly_dta(jpi2+1,:) = ref_wgts(kw)%fly_dta(jpi2,:)
945	ENDIF
946	IF( jpj1 == 2 ) THEN
947	ref_wgts(kw)%fly_dta(:,jpj1-1) = ref_wgts(kw)%fly_dta(:,jpj1)
948	ENDIF
949	IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
950	ref_wgts(kw)%fly_dta(:,jpj2+1) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2) - ref_wgts(kw)%fly_dta(:,jpj2-1)
951	ENDIF
952
953	!! if data grid is cyclic we can do better on east-west edges
954	!! but have to allow for whether first and last columns are coincident
955	IF( ref_wgts(kw)%cyclic ) THEN
956	rec1(2) = MAX( jpjmin-1, 1 )
957	recn(1) = 2
958	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
959	jpj1 = 2 + rec1(2) - jpjmin
960	jpj2 = jpj1 + recn(2) - 1
961	IF( jpi1 == 2 ) THEN
962	rec1(1) = ref_wgts(kw)%ddims(1) - 1
963	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2), nrec, rec1, recn)
964	ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2) = ref_wgts(kw)%col2(ref_wgts(kw)%offset+1,jpj1:jpj2)
965	ENDIF
966	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
967	rec1(1) = 1
968	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2), nrec, rec1, recn)
969	ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2) = ref_wgts(kw)%col2(2-ref_wgts(kw)%offset,jpj1:jpj2)
970	ENDIF
971	ENDIF
972
973	! gradient in the i direction
974	DO jk = 1,4
975	DO jn = 1, jpj
976	DO jm = 1,jpi
977	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
978	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
979	dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 * &
980	(ref_wgts(kw)%fly_dta(ni+2,nj+1) - ref_wgts(kw)%fly_dta(ni,nj+1))
981	END DO
982	END DO
983	END DO
984
985	! gradient in the j direction
986	DO jk = 1,4
987	DO jn = 1, jpj
988	DO jm = 1,jpi
989	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
990	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
991	dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 * &
992	(ref_wgts(kw)%fly_dta(ni+1,nj+2) - ref_wgts(kw)%fly_dta(ni+1,nj))
993	END DO
994	END DO
995	END DO
996
997	! gradient in the ij direction
998	DO jk = 1,4
999	DO jn = 1, jpj
1000	DO jm = 1,jpi
1001	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
1002	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
1003	dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
1004	(ref_wgts(kw)%fly_dta(ni+2,nj+2) - ref_wgts(kw)%fly_dta(ni ,nj+2)) - &
1005	(ref_wgts(kw)%fly_dta(ni+2,nj ) - ref_wgts(kw)%fly_dta(ni ,nj )))
1006	END DO
1007	END DO
1008	END DO
1009
1010	END IF
1011
1012	END SUBROUTINE fld_interp
1013
1014	END MODULE fldread

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

Download in other formats: