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 @ 1730

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

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

Download in other formats: