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

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 4371

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

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

Download in other formats: