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

Context Navigation

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 2800

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

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

Download in other formats: