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

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

fldread.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

File size: 61.6 KB

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

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

Download in other formats: