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

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

fldread.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 4161

Last change on this file since 4161 was 3851, checked in by smasson, 11 years ago
trunk: bugfix in fldread when nn_fsbc * rn_rdt > frcing file period, see #958
Property svn:keywords set to `Id`
File size: 70.9 KB

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

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

Download in other formats: