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

Context Navigation

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 2601

Last change on this file since 2601 was 2590, checked in by trackstand2, 13 years ago
Merge branch 'dynamic_memory' into master-svn-dyn
Property svn:keywords set to `Id`
File size: 62.2 KB

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

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

Download in other formats: