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 NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/fldread.F90 @ 12077

Last change on this file since 12077 was 11822, checked in by acc, 5 years ago
Branch 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. Sette tested updates to branch to align with trunk changes between 10721 and 11740. Sette tests are passing but results differ from branch before these changes (except for GYRE_PISCES and VORTEX) and branch results already differed from trunk because of algorithmic fixes. Will need more checks to confirm correctness.
Property svn:keywords set to `Id`
File size: 101.7 KB

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

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

Download in other formats: