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

Context Navigation

source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 5232

Last change on this file since 5232 was 5232, checked in by davestorkey, 9 years ago
Svn keywords deactivated using "svn propdel" in branch 2015/dev_r5021_UKMO1_CICE_coupling.
File size: 84.1 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'
[4663]	42	CHARACTER(len = 256) :: wname ! generic name of a NetCDF weights file to be used, blank if not
[1730]	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)
[5232]	110	!! $Id: fldread.F90 4784 2014-09-24 08:44:53Z jamesharle $
[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	!
[4784]	475	ztmp = REAL( nsec_year, wp ) / ( REAL( nyear_len(1), wp ) * rday ) + 0.5 &
	476	& + REAL( it_offset, wp ) / ( REAL( nyear_len(1), wp ) * rday )
[2528]	477	sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
	478	! swap at the middle of the year
[4784]	479	IF( llbefore ) THEN ; sdjf%nrec_a(2) = nsec1jan000 - (1 - INT(ztmp)) * NINT(0.5 * rday) * nyear_len(0) + &
	480	& INT(ztmp) * NINT( 0.5 * rday) * nyear_len(1)
	481	ELSE ; sdjf%nrec_a(2) = nsec1jan000 + (1 - INT(ztmp)) * NINT(0.5 * rday) * nyear_len(1) + &
	482	& INT(ztmp) * INT(rday) * nyear_len(1) + INT(ztmp) * NINT( 0.5 * rday) * nyear_len(2)
[2528]	483	ENDIF
	484	ELSE ! no time interpolation
	485	sdjf%nrec_a(1) = 1
	486	sdjf%nrec_a(2) = NINT(rday) * nyear_len(1) + nsec1jan000 ! swap at the end of the year
	487	sdjf%nrec_b(2) = nsec1jan000 ! beginning of the year (only for print)
	488	ENDIF
	489	!
	490	! ! ============ !
	491	ELSEIF( sdjf%nfreqh == -1 ) THEN ! monthly mean !
	492	! ! ============ !
	493	!
	494	IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record
	495	!
	496	! INT( ztmp )
	497	! /\|\
	498	! 1 \| *----
	499	! 0 \|----(
	500	! \|----+----\|--> time
[1132]	501	! 0 /\|\ 1 (nday/nmonth_len(nmonth))
	502	! \|
	503	! \|
	504	! forcing record : nmonth
	505	!
[4784]	506	ztmp = REAL( nsec_month, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) + 0.5 &
	507	& + REAL( it_offset, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday )
[2528]	508	imth = nmonth + INT( ztmp ) - COUNT((/llbefore/))
	509	IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
	510	ELSE ; sdjf%nrec_a(1) = imth
	511	ENDIF
	512	sdjf%nrec_a(2) = nmonth_half( imth ) + nsec1jan000 ! swap at the middle of the month
	513	ELSE ! no time interpolation
	514	IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nrec_a(1) = 1
	515	ELSE ; sdjf%nrec_a(1) = nmonth
	516	ENDIF
	517	sdjf%nrec_a(2) = nmonth_end(nmonth ) + nsec1jan000 ! swap at the end of the month
	518	sdjf%nrec_b(2) = nmonth_end(nmonth-1) + nsec1jan000 ! beginning of the month (only for print)
[888]	519	ENDIF
	520	!
[2528]	521	! ! ================================ !
	522	ELSE ! higher frequency mean (in hours)
	523	! ! ================================ !
[888]	524	!
[4245]	525	ifreq_sec = NINT( sdjf%nfreqh * 3600 ) ! frequency mean (in seconds)
[2528]	526	IF( sdjf%cltype(1:4) == 'week' ) isec_week = ksec_week( sdjf%cltype(6:8) ) ! since the first day of the current week
[1132]	527	! number of second since the beginning of the file
[2528]	528	IF( sdjf%cltype == 'monthly' ) THEN ; ztmp = REAL(nsec_month,wp) ! since the first day of the current month
	529	ELSEIF( sdjf%cltype(1:4) == 'week' ) THEN ; ztmp = REAL(isec_week ,wp) ! since the first day of the current week
	530	ELSEIF( sdjf%cltype == 'daily' ) THEN ; ztmp = REAL(nsec_day ,wp) ! since 00h of the current day
	531	ELSE ; ztmp = REAL(nsec_year ,wp) ! since 00h on Jan 1 of the current year
[1132]	532	ENDIF
[3851]	533	ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdttra(1) + REAL( it_offset, wp ) ! centrered in the middle of sbc time step
	534	ztmp = ztmp + 0.01 * rdttra(1) ! avoid truncation error
[1132]	535	IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record
	536	!
[3851]	537	! INT( ztmp/ifreq_sec + 0.5 )
[1132]	538	! /\|\
	539	! 2 \| *-----(
	540	! 1 \| *-----(
	541	! 0 \|--(
	542	! \|--+--\|--+--\|--+--\|--> time
[3851]	543	! 0 /\|\ 1 /\|\ 2 /\|\ 3 (ztmp/ifreq_sec)
[1132]	544	! \| \| \|
	545	! \| \| \|
	546	! forcing record : 1 2 3
	547	!
[2528]	548	ztmp= ztmp / REAL(ifreq_sec, wp) + 0.5
	549	ELSE ! no time interpolation
[1132]	550	!
[3851]	551	! INT( ztmp/ifreq_sec )
[1132]	552	! /\|\
	553	! 2 \| *-----(
	554	! 1 \| *-----(
	555	! 0 \|-----(
	556	! \|--+--\|--+--\|--+--\|--> time
[3851]	557	! 0 /\|\ 1 /\|\ 2 /\|\ 3 (ztmp/ifreq_sec)
[1132]	558	! \| \| \|
	559	! \| \| \|
	560	! forcing record : 1 2 3
	561	!
[2528]	562	ztmp= ztmp / REAL(ifreq_sec, wp)
[1132]	563	ENDIF
[3851]	564	sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/)) ! record number to be read
[1132]	565
[2528]	566	iendrec = ifreq_sec * sdjf%nrec_a(1) + nsec1jan000 ! end of this record (in second)
	567	! add the number of seconds between 00h Jan 1 and the end of previous month/week/day (ok if nmonth=1)
	568	IF( sdjf%cltype == 'monthly' ) iendrec = iendrec + NINT(rday) * SUM(nmonth_len(1:nmonth -1))
	569	IF( sdjf%cltype(1:4) == 'week' ) iendrec = iendrec + ( nsec_year - isec_week )
	570	IF( sdjf%cltype == 'daily' ) iendrec = iendrec + NINT(rday) * ( nday_year - 1 )
	571	IF( sdjf%ln_tint ) THEN
	572	sdjf%nrec_a(2) = iendrec - ifreq_sec / 2 ! swap at the middle of the record
	573	ELSE
	574	sdjf%nrec_a(2) = iendrec ! swap at the end of the record
	575	sdjf%nrec_b(2) = iendrec - ifreq_sec ! beginning of the record (only for print)
	576	ENDIF
[888]	577	!
	578	ENDIF
	579	!
[1132]	580	END SUBROUTINE fld_rec
	581
	582
[3294]	583	SUBROUTINE fld_get( sdjf, map )
[2528]	584	!!---------------------------------------------------------------------
[3294]	585	!! * ROUTINE fld_get *
[2528]	586	!!
	587	!! ** Purpose : read the data
	588	!!----------------------------------------------------------------------
[2715]	589	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
[3851]	590	TYPE(MAP_POINTER),INTENT(in) :: map ! global-to-local mapping indices
[2528]	591	!!
[2715]	592	INTEGER :: ipk ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
	593	INTEGER :: iw ! index into wgts array
[3680]	594	INTEGER :: ipdom ! index of the domain
[4245]	595	INTEGER :: idvar ! variable ID
	596	INTEGER :: idmspc ! number of spatial dimensions
	597	LOGICAL :: lmoor ! C1D case: point data
[2528]	598	!!---------------------------------------------------------------------
[3851]	599	!
[2528]	600	ipk = SIZE( sdjf%fnow, 3 )
[3680]	601	!
[3851]	602	IF( ASSOCIATED(map%ptr) ) THEN
	603	IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map%ptr )
	604	ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1), map%ptr )
[3294]	605	ENDIF
	606	ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
[2528]	607	CALL wgt_list( sdjf, iw )
[4230]	608	IF( sdjf%ln_tint ) THEN ; CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fdta(:,:,:,2), &
	609	& sdjf%nrec_a(1), sdjf%lsmname )
	610	ELSE ; CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fnow(:,:,: ), &
	611	& sdjf%nrec_a(1), sdjf%lsmname )
[2528]	612	ENDIF
	613	ELSE
[3680]	614	IF( SIZE(sdjf%fnow, 1) == jpi ) THEN ; ipdom = jpdom_data
	615	ELSE ; ipdom = jpdom_unknown
	616	ENDIF
[4245]	617	! C1D case: If product of spatial dimensions == ipk, then x,y are of
	618	! size 1 (point/mooring data): this must be read onto the central grid point
	619	idvar = iom_varid( sdjf%num, sdjf%clvar )
	620	idmspc = iom_file( sdjf%num )%ndims( idvar )
	621	IF( iom_file( sdjf%num )%luld( idvar ) ) idmspc = idmspc - 1
	622	lmoor = (idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk)
	623	!
[2528]	624	SELECT CASE( ipk )
[3680]	625	CASE(1)
[4245]	626	IF( lk_c1d .AND. lmoor ) THEN
	627	IF( sdjf%ln_tint ) THEN
	628	CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fdta(2,2,1,2), sdjf%nrec_a(1) )
	629	CALL lbc_lnk( sdjf%fdta(:,:,1,2),'Z',1. )
	630	ELSE
	631	CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fnow(2,2,1 ), sdjf%nrec_a(1) )
	632	CALL lbc_lnk( sdjf%fnow(:,:,1 ),'Z',1. )
	633	ENDIF
	634	ELSE
	635	IF( sdjf%ln_tint ) THEN ; CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,1,2), sdjf%nrec_a(1) )
	636	ELSE ; CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,1 ), sdjf%nrec_a(1) )
	637	ENDIF
[2528]	638	ENDIF
	639	CASE DEFAULT
[4245]	640	IF (lk_c1d .AND. lmoor ) THEN
	641	IF( sdjf%ln_tint ) THEN
	642	CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fdta(2,2,:,2), sdjf%nrec_a(1) )
	643	CALL lbc_lnk( sdjf%fdta(:,:,:,2),'Z',1. )
	644	ELSE
	645	CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fnow(2,2,: ), sdjf%nrec_a(1) )
	646	CALL lbc_lnk( sdjf%fnow(:,:,: ),'Z',1. )
	647	ENDIF
	648	ELSE
	649	IF( sdjf%ln_tint ) THEN ; CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) )
	650	ELSE ; CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1) )
	651	ENDIF
[2528]	652	ENDIF
	653	END SELECT
	654	ENDIF
	655	!
[3851]	656	sdjf%rotn(2) = .false. ! vector not yet rotated
[2528]	657
	658	END SUBROUTINE fld_get
	659
[3294]	660	SUBROUTINE fld_map( num, clvar, dta, nrec, map )
	661	!!---------------------------------------------------------------------
[3851]	662	!! * ROUTINE fld_map *
[3294]	663	!!
	664	!! ** Purpose : read global data from file and map onto local data
	665	!! using a general mapping (for open boundaries)
	666	!!----------------------------------------------------------------------
	667	#if defined key_bdy
[3651]	668	USE bdy_oce, ONLY: dta_global, dta_global2 ! workspace to read in global data arrays
[3294]	669	#endif
	670	INTEGER , INTENT(in ) :: num ! stream number
	671	CHARACTER(LEN=*) , INTENT(in ) :: clvar ! variable name
	672	REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional)
	673	INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice)
	674	INTEGER, DIMENSION(:) , INTENT(in ) :: map ! global-to-local mapping indices
	675	!!
	676	INTEGER :: ipi ! length of boundary data on local process
	677	INTEGER :: ipj ! length of dummy dimension ( = 1 )
	678	INTEGER :: ipk ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk )
	679	INTEGER :: ilendta ! length of data in file
	680	INTEGER :: idvar ! variable ID
[3651]	681	INTEGER :: ib, ik, ji, jj ! loop counters
[3294]	682	INTEGER :: ierr
[3651]	683	REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read ! work space for global data
[3294]	684	!!---------------------------------------------------------------------
	685
	686	ipi = SIZE( dta, 1 )
	687	ipj = 1
	688	ipk = SIZE( dta, 3 )
	689
	690	idvar = iom_varid( num, clvar )
	691	ilendta = iom_file(num)%dimsz(1,idvar)
[3651]	692
	693	#if defined key_bdy
	694	ipj = iom_file(num)%dimsz(2,idvar)
	695	IF (ipj == 1) THEN ! we assume that this is a structured open boundary file
	696	dta_read => dta_global
	697	ELSE
	698	dta_read => dta_global2
	699	ENDIF
	700	#endif
	701
[3294]	702	IF(lwp) WRITE(numout,*) 'Dim size for ',TRIM(clvar),' is ', ilendta
	703	IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk
	704
	705	SELECT CASE( ipk )
[3851]	706	CASE(1) ; CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1 ), nrec )
	707	CASE DEFAULT ; CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec )
[3294]	708	END SELECT
	709	!
[3651]	710	IF (ipj==1) THEN
	711	DO ib = 1, ipi
	712	DO ik = 1, ipk
	713	dta(ib,1,ik) = dta_read(map(ib),1,ik)
	714	END DO
[3294]	715	END DO
[3651]	716	ELSE ! we assume that this is a structured open boundary file
	717	DO ib = 1, ipi
	718	jj=1+floor(REAL(map(ib)-1)/REAL(ilendta))
	719	ji=map(ib)-(jj-1)*ilendta
	720	DO ik = 1, ipk
	721	dta(ib,1,ik) = dta_read(ji,jj,ik)
	722	END DO
	723	END DO
	724	ENDIF
[3294]	725
	726	END SUBROUTINE fld_map
	727
	728
[2528]	729	SUBROUTINE fld_rot( kt, sd )
	730	!!---------------------------------------------------------------------
[3294]	731	!! * ROUTINE fld_rot *
[2528]	732	!!
	733	!! ** Purpose : Vector fields may need to be rotated onto the local grid direction
[2715]	734	!!----------------------------------------------------------------------
[2528]	735	INTEGER , INTENT(in ) :: kt ! ocean time step
	736	TYPE(FLD), INTENT(inout), DIMENSION(:) :: sd ! input field related variables
	737	!!
[3851]	738	INTEGER :: ju,jv,jk,jn ! loop indices
[3294]	739	INTEGER :: imf ! size of the structure sd
	740	INTEGER :: ill ! character length
	741	INTEGER :: iv ! indice of V component
	742	REAL(wp), POINTER, DIMENSION(:,:) :: utmp, vtmp ! temporary arrays for vector rotation
	743	CHARACTER (LEN=100) :: clcomp ! dummy weight name
[2528]	744	!!---------------------------------------------------------------------
[2715]	745
[3294]	746	CALL wrk_alloc( jpi,jpj, utmp, vtmp )
[2715]	747
[2528]	748	!! (sga: following code should be modified so that pairs arent searched for each time
	749	!
	750	imf = SIZE( sd )
	751	DO ju = 1, imf
	752	ill = LEN_TRIM( sd(ju)%vcomp )
[3851]	753	DO jn = 2-COUNT((/sd(ju)%ln_tint/)), 2
	754	IF( ill > 0 .AND. .NOT. sd(ju)%rotn(jn) ) THEN ! find vector rotations required
	755	IF( sd(ju)%vcomp(1:1) == 'U' ) THEN ! east-west component has symbolic name starting with 'U'
	756	! look for the north-south component which has same symbolic name but with 'U' replaced with 'V'
	757	clcomp = 'V' // sd(ju)%vcomp(2:ill) ! works even if ill == 1
	758	iv = -1
	759	DO jv = 1, imf
	760	IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) ) iv = jv
	761	END DO
	762	IF( iv > 0 ) THEN ! fields ju and iv are two components which need to be rotated together
	763	DO jk = 1, SIZE( sd(ju)%fnow, 3 )
	764	IF( sd(ju)%ln_tint )THEN
	765	CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->i', utmp(:,:) )
	766	CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->j', vtmp(:,:) )
	767	sd(ju)%fdta(:,:,jk,jn) = utmp(:,:) ; sd(iv)%fdta(:,:,jk,jn) = vtmp(:,:)
	768	ELSE
	769	CALL rot_rep( sd(ju)%fnow(:,:,jk ), sd(iv)%fnow(:,:,jk ), 'T', 'en->i', utmp(:,:) )
	770	CALL rot_rep( sd(ju)%fnow(:,:,jk ), sd(iv)%fnow(:,:,jk ), 'T', 'en->j', vtmp(:,:) )
	771	sd(ju)%fnow(:,:,jk ) = utmp(:,:) ; sd(iv)%fnow(:,:,jk ) = vtmp(:,:)
	772	ENDIF
	773	END DO
	774	sd(ju)%rotn(jn) = .TRUE. ! vector was rotated
	775	IF( lwp .AND. kt == nit000 ) WRITE(numout,*) &
	776	& 'fld_read: vector pair ('//TRIM(sd(ju)%clvar)//', '//TRIM(sd(iv)%clvar)//') rotated on to model grid'
	777	ENDIF
	778	ENDIF
	779	ENDIF
	780	END DO
[2528]	781	END DO
[2715]	782	!
[3294]	783	CALL wrk_dealloc( jpi,jpj, utmp, vtmp )
[2715]	784	!
[2528]	785	END SUBROUTINE fld_rot
	786
	787
[1628]	788	SUBROUTINE fld_clopn( sdjf, kyear, kmonth, kday, ldstop )
[1132]	789	!!---------------------------------------------------------------------
	790	!! * ROUTINE fld_clopn *
	791	!!
	792	!! ** Purpose : update the file name and open the file
	793	!!----------------------------------------------------------------------
[2715]	794	TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables
[3851]	795	INTEGER, OPTIONAL, INTENT(in ) :: kyear ! year value
	796	INTEGER, OPTIONAL, INTENT(in ) :: kmonth ! month value
	797	INTEGER, OPTIONAL, INTENT(in ) :: kday ! day value
[2715]	798	LOGICAL, OPTIONAL, INTENT(in ) :: ldstop ! stop if open to read a non-existing file (default = .TRUE.)
[3851]	799	!!
	800	LOGICAL :: llprevyr ! are we reading previous year file?
	801	LOGICAL :: llprevmth ! are we reading previous month file?
	802	INTEGER :: iyear, imonth, iday ! first day of the current file in yyyy mm dd
	803	INTEGER :: isec_week ! number of seconds since start of the weekly file
	804	INTEGER :: indexyr ! year undex (O/1/2: previous/current/next)
	805	INTEGER :: iyear_len, imonth_len ! length (days) of iyear and imonth !
	806	CHARACTER(len = 256):: clname ! temporary file name
[2715]	807	!!----------------------------------------------------------------------
[3851]	808	IF( PRESENT(kyear) ) THEN ! use given values
	809	iyear = kyear
	810	imonth = kmonth
	811	iday = kday
	812	ELSE ! use current day values
	813	IF ( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of the week
	814	isec_week = ksec_week( sdjf%cltype(6:8) ) ! second since the beginning of the week
	815	llprevmth = isec_week > nsec_month ! longer time since beginning of the week than the month
	816	llprevyr = llprevmth .AND. nmonth == 1
	817	ELSE
	818	isec_week = 0
	819	llprevmth = .FALSE.
	820	llprevyr = .FALSE.
	821	ENDIF
	822	iyear = nyear - COUNT((/llprevyr /))
	823	imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
	824	iday = nday + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
	825	ENDIF
[1132]	826
	827	! build the new filename if not climatological data
[3851]	828	clname=TRIM(sdjf%clrootname)
[2528]	829	!
[3851]	830	! note that sdjf%ln_clim is is only acting on the presence of the year in the file name
[2528]	831	IF( .NOT. sdjf%ln_clim ) THEN
[3851]	832	WRITE(clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), iyear ! add year
	833	IF( sdjf%cltype /= 'yearly' ) WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname ), imonth ! add month
[2528]	834	ELSE
	835	! build the new filename if climatological data
[3851]	836	IF( sdjf%cltype /= 'yearly' ) WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), imonth ! add month
[888]	837	ENDIF
[2528]	838	IF( sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) &
[3851]	839	& WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname ), iday ! add day
[2528]	840	!
[3851]	841	IF( TRIM(clname) /= TRIM(sdjf%clname) .OR. sdjf%num == 0 ) THEN ! new file to be open
	842
	843	sdjf%clname = TRIM(clname)
	844	IF( sdjf%num /= 0 ) CALL iom_close( sdjf%num ) ! close file if already open
	845	CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 )
	846
	847	! find the last record to be read -> update sdjf%nreclast
	848	indexyr = iyear - nyear + 1
	849	iyear_len = nyear_len( indexyr )
	850	SELECT CASE ( indexyr )
	851	CASE ( 0 ) ; imonth_len = 31 ! previous year -> imonth = 12
	852	CASE ( 1 ) ; imonth_len = nmonth_len(imonth)
	853	CASE ( 2 ) ; imonth_len = 31 ! next year -> imonth = 1
	854	END SELECT
	855
	856	! last record to be read in the current file
	857	IF ( sdjf%nfreqh == -12 ) THEN ; sdjf%nreclast = 1 ! yearly mean
	858	ELSEIF( sdjf%nfreqh == -1 ) THEN ! monthly mean
	859	IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nreclast = 1
	860	ELSE ; sdjf%nreclast = 12
	861	ENDIF
	862	ELSE ! higher frequency mean (in hours)
[4245]	863	IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nreclast = NINT( 24 * imonth_len / sdjf%nfreqh )
	864	ELSEIF( sdjf%cltype(1:4) == 'week' ) THEN ; sdjf%nreclast = NINT( 24 * 7 / sdjf%nfreqh )
	865	ELSEIF( sdjf%cltype == 'daily' ) THEN ; sdjf%nreclast = NINT( 24 / sdjf%nfreqh )
	866	ELSE ; sdjf%nreclast = NINT( 24 * iyear_len / sdjf%nfreqh )
[3851]	867	ENDIF
	868	ENDIF
	869
	870	ENDIF
	871	!
[1132]	872	END SUBROUTINE fld_clopn
	873
	874
	875	SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )
	876	!!---------------------------------------------------------------------
	877	!! * ROUTINE fld_fill *
	878	!!
	879	!! ** Purpose : fill sdf with sdf_n and control print
	880	!!----------------------------------------------------------------------
	881	TYPE(FLD) , DIMENSION(:), INTENT(inout) :: sdf ! structure of input fields (file informations, fields read)
	882	TYPE(FLD_N), DIMENSION(:), INTENT(in ) :: sdf_n ! array of namelist information structures
	883	CHARACTER(len=*) , INTENT(in ) :: cdir ! Root directory for location of flx files
	884	CHARACTER(len=*) , INTENT(in ) :: cdcaller !
	885	CHARACTER(len=*) , INTENT(in ) :: cdtitle !
	886	CHARACTER(len=*) , INTENT(in ) :: cdnam !
[888]	887	!
[1132]	888	INTEGER :: jf ! dummy indices
	889	!!---------------------------------------------------------------------
[888]	890
[1132]	891	DO jf = 1, SIZE(sdf)
	892	sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
[3851]	893	sdf(jf)%clname = "not yet defined"
[1730]	894	sdf(jf)%nfreqh = sdf_n(jf)%nfreqh
[1132]	895	sdf(jf)%clvar = sdf_n(jf)%clvar
	896	sdf(jf)%ln_tint = sdf_n(jf)%ln_tint
	897	sdf(jf)%ln_clim = sdf_n(jf)%ln_clim
[2528]	898	sdf(jf)%cltype = sdf_n(jf)%cltype
[3851]	899	sdf(jf)%num = -1
	900	sdf(jf)%wgtname = " "
[1730]	901	IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
[4230]	902	sdf(jf)%lsmname = " "
	903	IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 ) sdf(jf)%lsmname = TRIM( cdir )//TRIM( sdf_n(jf)%lname )
[3851]	904	sdf(jf)%vcomp = sdf_n(jf)%vcomp
	905	sdf(jf)%rotn(:) = .TRUE. ! pretend to be rotated -> won't try to rotate data before the first call to fld_get
	906	IF( sdf(jf)%cltype(1:4) == 'week' .AND. nn_leapy == 0 ) &
	907	& CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs nn_leapy = 1')
	908	IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim ) &
	909	& CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.')
[4371]	910	sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
[1132]	911	END DO
	912
	913	IF(lwp) THEN ! control print
	914	WRITE(numout,*)
	915	WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
	916	WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
	917	WRITE(numout,*) ' '//TRIM( cdnam )//' Namelist'
	918	WRITE(numout,*) ' list of files and frequency (>0: in hours ; <0 in months)'
	919	DO jf = 1, SIZE(sdf)
	920	WRITE(numout,*) ' root filename: ' , TRIM( sdf(jf)%clrootname ), &
	921	& ' variable name: ' , TRIM( sdf(jf)%clvar )
[1730]	922	WRITE(numout,*) ' frequency: ' , sdf(jf)%nfreqh , &
[1132]	923	& ' time interp: ' , sdf(jf)%ln_tint , &
	924	& ' climatology: ' , sdf(jf)%ln_clim , &
[1275]	925	& ' weights : ' , TRIM( sdf(jf)%wgtname ), &
	926	& ' pairing : ' , TRIM( sdf(jf)%vcomp ), &
[4230]	927	& ' data type: ' , sdf(jf)%cltype , &
	928	& ' land/sea mask:' , TRIM( sdf(jf)%lsmname )
[2528]	929	call flush(numout)
[1132]	930	END DO
	931	ENDIF
	932
	933	END SUBROUTINE fld_fill
	934
	935
[1275]	936	SUBROUTINE wgt_list( sd, kwgt )
	937	!!---------------------------------------------------------------------
	938	!! * ROUTINE wgt_list *
	939	!!
	940	!! ** Purpose : search array of WGTs and find a weights file
	941	!! entry, or return a new one adding it to the end
	942	!! if it is a new entry, the weights data is read in and
	943	!! restructured (fld_weight)
	944	!!----------------------------------------------------------------------
[2715]	945	TYPE( FLD ), INTENT(in ) :: sd ! field with name of weights file
	946	INTEGER , INTENT(inout) :: kwgt ! index of weights
[1275]	947	!!
[2715]	948	INTEGER :: kw, nestid ! local integer
	949	LOGICAL :: found ! local logical
[1275]	950	!!----------------------------------------------------------------------
	951	!
	952	!! search down linked list
	953	!! weights filename is either present or we hit the end of the list
	954	found = .FALSE.
	955
	956	!! because agrif nest part of filenames are now added in iom_open
	957	!! to distinguish between weights files on the different grids, need to track
	958	!! nest number explicitly
	959	nestid = 0
	960	#if defined key_agrif
	961	nestid = Agrif_Fixed()
	962	#endif
	963	DO kw = 1, nxt_wgt-1
	964	IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
	965	ref_wgts(kw)%nestid == nestid) THEN
	966	kwgt = kw
	967	found = .TRUE.
	968	EXIT
	969	ENDIF
	970	END DO
	971	IF( .NOT.found ) THEN
	972	kwgt = nxt_wgt
	973	CALL fld_weight( sd )
	974	ENDIF
[2715]	975	!
[1275]	976	END SUBROUTINE wgt_list
	977
[2715]	978
[1275]	979	SUBROUTINE wgt_print( )
	980	!!---------------------------------------------------------------------
	981	!! * ROUTINE wgt_print *
	982	!!
	983	!! ** Purpose : print the list of known weights
	984	!!----------------------------------------------------------------------
[2715]	985	INTEGER :: kw !
[1275]	986	!!----------------------------------------------------------------------
	987	!
	988	DO kw = 1, nxt_wgt-1
	989	WRITE(numout,*) 'weight file: ',TRIM(ref_wgts(kw)%wgtname)
	990	WRITE(numout,*) ' ddims: ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
	991	WRITE(numout,*) ' numwgt: ',ref_wgts(kw)%numwgt
	992	WRITE(numout,*) ' jpiwgt: ',ref_wgts(kw)%jpiwgt
	993	WRITE(numout,*) ' jpjwgt: ',ref_wgts(kw)%jpjwgt
	994	WRITE(numout,*) ' botleft: ',ref_wgts(kw)%botleft
	995	WRITE(numout,*) ' topright: ',ref_wgts(kw)%topright
	996	IF( ref_wgts(kw)%cyclic ) THEN
	997	WRITE(numout,*) ' cyclical'
[2528]	998	IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) ' with overlap of ', ref_wgts(kw)%overlap
[1275]	999	ELSE
	1000	WRITE(numout,*) ' not cyclical'
	1001	ENDIF
	1002	IF( ASSOCIATED(ref_wgts(kw)%data_wgt) ) WRITE(numout,*) ' allocated'
	1003	END DO
[2715]	1004	!
[1275]	1005	END SUBROUTINE wgt_print
	1006
[2715]	1007
[1275]	1008	SUBROUTINE fld_weight( sd )
	1009	!!---------------------------------------------------------------------
	1010	!! * ROUTINE fld_weight *
	1011	!!
	1012	!! ** Purpose : create a new WGT structure and fill in data from
	1013	!! file, restructuring as required
	1014	!!----------------------------------------------------------------------
[2715]	1015	TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file
	1016	!!
[3294]	1017	INTEGER :: jn ! dummy loop indices
	1018	INTEGER :: inum ! temporary logical unit
	1019	INTEGER :: id ! temporary variable id
	1020	INTEGER :: ipk ! temporary vertical dimension
	1021	CHARACTER (len=5) :: aname
	1022	INTEGER , DIMENSION(3) :: ddims
	1023	INTEGER , POINTER, DIMENSION(:,:) :: data_src
	1024	REAL(wp), POINTER, DIMENSION(:,:) :: data_tmp
	1025	LOGICAL :: cyclical
	1026	INTEGER :: zwrap ! local integer
[1275]	1027	!!----------------------------------------------------------------------
	1028	!
[3294]	1029	CALL wrk_alloc( jpi,jpj, data_src ) ! integer
	1030	CALL wrk_alloc( jpi,jpj, data_tmp )
[2715]	1031	!
[1275]	1032	IF( nxt_wgt > tot_wgts ) THEN
[2777]	1033	CALL ctl_stop("fld_weight: weights array size exceeded, increase tot_wgts")
[1275]	1034	ENDIF
	1035	!
	1036	!! new weights file entry, add in extra information
	1037	!! a weights file represents a 2D grid of a certain shape, so we assume that the current
	1038	!! input data file is representative of all other files to be opened and processed with the
	1039	!! current weights file
	1040
	1041	!! open input data file (non-model grid)
[1319]	1042	CALL iom_open( sd%clname, inum, ldiof = LEN(TRIM(sd%wgtname)) > 0 )
[1275]	1043
	1044	!! get dimensions
	1045	id = iom_varid( inum, sd%clvar, ddims )
	1046
[2528]	1047	!! close it
	1048	CALL iom_close( inum )
[1275]	1049
[2528]	1050	!! now open the weights file
[1275]	1051
[2528]	1052	CALL iom_open ( sd%wgtname, inum ) ! interpolation weights
	1053	IF ( inum > 0 ) THEN
[1275]	1054
[2528]	1055	!! determine whether we have an east-west cyclic grid
	1056	!! from global attribute called "ew_wrap" in the weights file
	1057	!! note that if not found, iom_getatt returns -999 and cyclic with no overlap is assumed
	1058	!! since this is the most common forcing configuration
[1275]	1059
[2528]	1060	CALL iom_getatt(inum, 'ew_wrap', zwrap)
	1061	IF( zwrap >= 0 ) THEN
[1275]	1062	cyclical = .TRUE.
[2528]	1063	ELSE IF( zwrap == -999 ) THEN
[1275]	1064	cyclical = .TRUE.
[2528]	1065	zwrap = 0
	1066	ELSE
	1067	cyclical = .FALSE.
[1275]	1068	ENDIF
	1069
	1070	ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
	1071	ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
	1072	ref_wgts(nxt_wgt)%wgtname = sd%wgtname
[2528]	1073	ref_wgts(nxt_wgt)%overlap = zwrap
	1074	ref_wgts(nxt_wgt)%cyclic = cyclical
[1275]	1075	ref_wgts(nxt_wgt)%nestid = 0
	1076	#if defined key_agrif
	1077	ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
	1078	#endif
	1079	!! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
	1080	!! for each weight wgtNN there is an integer array srcNN which gives the point in
	1081	!! the input data grid which is to be multiplied by the weight
	1082	!! they are both arrays on the model grid so the result of the multiplication is
	1083	!! added into an output array on the model grid as a running sum
	1084
	1085	!! two possible cases: bilinear (4 weights) or bicubic (16 weights)
	1086	id = iom_varid(inum, 'src05', ldstop=.FALSE.)
	1087	IF( id <= 0) THEN
	1088	ref_wgts(nxt_wgt)%numwgt = 4
	1089	ELSE
	1090	ref_wgts(nxt_wgt)%numwgt = 16
	1091	ENDIF
	1092
	1093	ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
	1094	ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
	1095	ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
	1096
	1097	DO jn = 1,4
	1098	aname = ' '
	1099	WRITE(aname,'(a3,i2.2)') 'src',jn
	1100	data_tmp(:,:) = 0
[1955]	1101	CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
[1275]	1102	data_src(:,:) = INT(data_tmp(:,:))
	1103	ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
	1104	ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
	1105	END DO
	1106
	1107	DO jn = 1, ref_wgts(nxt_wgt)%numwgt
	1108	aname = ' '
	1109	WRITE(aname,'(a3,i2.2)') 'wgt',jn
[1955]	1110	ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
	1111	CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
[1275]	1112	END DO
	1113	CALL iom_close (inum)
	1114
	1115	! find min and max indices in grid
[1955]	1116	ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
	1117	ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
	1118	ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
	1119	ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
[1275]	1120
	1121	! and therefore dimensions of the input box
	1122	ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
	1123	ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
	1124
	1125	! shift indexing of source grid
	1126	ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
	1127	ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
	1128
	1129	! create input grid, give it a halo to allow gradient calculations
[1702]	1130	! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
	1131	! a more robust solution will be given in next release
[2528]	1132	ipk = SIZE(sd%fnow, 3)
	1133	ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
	1134	IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
[1275]	1135
	1136	nxt_wgt = nxt_wgt + 1
	1137
	1138	ELSE
	1139	CALL ctl_stop( ' fld_weight : unable to read the file ' )
	1140	ENDIF
	1141
[3294]	1142	CALL wrk_dealloc( jpi,jpj, data_src ) ! integer
	1143	CALL wrk_dealloc( jpi,jpj, data_tmp )
[2715]	1144	!
[1275]	1145	END SUBROUTINE fld_weight
	1146
[2715]	1147
[4230]	1148	SUBROUTINE apply_seaoverland(clmaskfile,zfieldo,jpi1_lsm,jpi2_lsm,jpj1_lsm, &
	1149	& jpj2_lsm,itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
[1275]	1150	!!---------------------------------------------------------------------
[4230]	1151	!! * ROUTINE apply_seaoverland *
	1152	!!
	1153	!! ** Purpose : avoid spurious fluxes in coastal or near-coastal areas
	1154	!! due to the wrong usage of "land" values from the coarse
	1155	!! atmospheric model when spatial interpolation is required
	1156	!! D. Delrosso INGV
	1157	!!----------------------------------------------------------------------
	1158	INTEGER :: inum,jni,jnj,jnz,jc ! temporary indices
	1159	INTEGER, INTENT(in) :: itmpi,itmpj,itmpz ! lengths
	1160	INTEGER, INTENT(in) :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices
	1161	INTEGER, DIMENSION(3), INTENT(in) :: rec1_lsm,recn_lsm ! temporary arrays for start and length
	1162	REAL(wp),DIMENSION (:,:,:),INTENT(inout) :: zfieldo ! input/output array for seaoverland application
	1163	REAL(wp),DIMENSION (:,:,:),ALLOCATABLE :: zslmec1 ! temporary array for land point detection
	1164	REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfieldn ! array of forcing field with undeff for land points
	1165	REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfield ! array of forcing field
	1166	CHARACTER (len=100), INTENT(in) :: clmaskfile ! land/sea mask file name
	1167	!!---------------------------------------------------------------------
	1168	ALLOCATE ( zslmec1(itmpi,itmpj,itmpz) )
	1169	ALLOCATE ( zfieldn(itmpi,itmpj) )
	1170	ALLOCATE ( zfield(itmpi,itmpj) )
	1171
	1172	! Retrieve the land sea mask data
	1173	CALL iom_open( clmaskfile, inum )
	1174	SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
	1175	CASE(1)
	1176	CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm)
	1177	CASE DEFAULT
	1178	CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm)
	1179	END SELECT
	1180	CALL iom_close( inum )
	1181
	1182	DO jnz=1,rec1_lsm(3) !! Loop over k dimension
	1183
	1184	DO jni=1,itmpi !! copy the original field into a tmp array
	1185	DO jnj=1,itmpj !! substituting undeff over land points
	1186	zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz)
	1187	IF ( zslmec1(jni,jnj,jnz) == 1. ) THEN
	1188	zfieldn(jni,jnj) = undeff_lsm
	1189	ENDIF
	1190	END DO
	1191	END DO
	1192
	1193	CALL seaoverland(zfieldn,itmpi,itmpj,zfield)
	1194	DO jc=1,nn_lsm
	1195	CALL seaoverland(zfield,itmpi,itmpj,zfield)
	1196	END DO
	1197
	1198	! Check for Undeff and substitute original values
	1199	IF(ANY(zfield==undeff_lsm)) THEN
	1200	DO jni=1,itmpi
	1201	DO jnj=1,itmpj
	1202	IF (zfield(jni,jnj)==undeff_lsm) THEN
	1203	zfield(jni,jnj) = zfieldo(jni,jnj,jnz)
	1204	ENDIF
	1205	ENDDO
	1206	ENDDO
	1207	ENDIF
	1208
	1209	zfieldo(:,:,jnz)=zfield(:,:)
	1210
	1211	END DO !! End Loop over k dimension
	1212
	1213	DEALLOCATE ( zslmec1 )
	1214	DEALLOCATE ( zfieldn )
	1215	DEALLOCATE ( zfield )
	1216
	1217	END SUBROUTINE apply_seaoverland
	1218
	1219
	1220	SUBROUTINE seaoverland(zfieldn,ileni,ilenj,zfield)
	1221	!!---------------------------------------------------------------------
	1222	!! * ROUTINE seaoverland *
	1223	!!
	1224	!! ** Purpose : create shifted matrices for seaoverland application
	1225	!! D. Delrosso INGV
	1226	!!----------------------------------------------------------------------
	1227	INTEGER,INTENT(in) :: ileni,ilenj ! lengths
	1228	REAL,DIMENSION (ileni,ilenj),INTENT(in) :: zfieldn ! array of forcing field with undeff for land points
	1229	REAL,DIMENSION (ileni,ilenj),INTENT(out) :: zfield ! array of forcing field
	1230	REAL,DIMENSION (ileni,ilenj) :: zmat1,zmat2,zmat3,zmat4 ! temporary arrays for seaoverland application
	1231	REAL,DIMENSION (ileni,ilenj) :: zmat5,zmat6,zmat7,zmat8 ! temporary arrays for seaoverland application
	1232	REAL,DIMENSION (ileni,ilenj) :: zlsm2d ! temporary arrays for seaoverland application
	1233	REAL,DIMENSION (ileni,ilenj,8) :: zlsm3d ! temporary arrays for seaoverland application
	1234	LOGICAL,DIMENSION (ileni,ilenj,8) :: ll_msknan3d ! logical mask for undeff detection
	1235	LOGICAL,DIMENSION (ileni,ilenj) :: ll_msknan2d ! logical mask for undeff detection
	1236	!!----------------------------------------------------------------------
	1237	zmat8 = eoshift(zfieldn , SHIFT=-1, BOUNDARY = (/zfieldn(:,1)/) ,DIM=2)
	1238	zmat1 = eoshift(zmat8 , SHIFT=-1, BOUNDARY = (/zmat8(1,:)/) ,DIM=1)
	1239	zmat2 = eoshift(zfieldn , SHIFT=-1, BOUNDARY = (/zfieldn(1,:)/) ,DIM=1)
	1240	zmat4 = eoshift(zfieldn , SHIFT= 1, BOUNDARY = (/zfieldn(:,ilenj)/),DIM=2)
	1241	zmat3 = eoshift(zmat4 , SHIFT=-1, BOUNDARY = (/zmat4(1,:)/) ,DIM=1)
	1242	zmat5 = eoshift(zmat4 , SHIFT= 1, BOUNDARY = (/zmat4(ileni,:)/) ,DIM=1)
	1243	zmat6 = eoshift(zfieldn , SHIFT= 1, BOUNDARY = (/zfieldn(ileni,:)/),DIM=1)
	1244	zmat7 = eoshift(zmat8 , SHIFT= 1, BOUNDARY = (/zmat8(ileni,:)/) ,DIM=1)
	1245
	1246	zlsm3d = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /))
	1247	ll_msknan3d = .not.(zlsm3d==undeff_lsm)
	1248	ll_msknan2d = .not.(zfieldn==undeff_lsm) ! FALSE where is Undeff (land)
	1249	zlsm2d = (SUM ( zlsm3d, 3 , ll_msknan3d ) )/(MAX(1,(COUNT( ll_msknan3d , 3 )) ))
	1250	WHERE ((COUNT( ll_msknan3d , 3 )) == 0.0_wp) zlsm2d = undeff_lsm
	1251	zfield = MERGE (zfieldn,zlsm2d,ll_msknan2d)
	1252	END SUBROUTINE seaoverland
	1253
	1254
	1255	SUBROUTINE fld_interp( num, clvar, kw, kk, dta, &
	1256	& nrec, lsmfile)
	1257	!!---------------------------------------------------------------------
[1275]	1258	!! * ROUTINE fld_interp *
	1259	!!
	1260	!! ** Purpose : apply weights to input gridded data to create data
	1261	!! on model grid
	1262	!!----------------------------------------------------------------------
[2715]	1263	INTEGER , INTENT(in ) :: num ! stream number
	1264	CHARACTER(LEN=*) , INTENT(in ) :: clvar ! variable name
	1265	INTEGER , INTENT(in ) :: kw ! weights number
	1266	INTEGER , INTENT(in ) :: kk ! vertical dimension of kk
	1267	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: dta ! output field on model grid
	1268	INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice)
[4230]	1269	CHARACTER(LEN=*) , INTENT(in ) :: lsmfile ! land sea mask file name
[1275]	1270	!!
[4230]	1271	REAL(wp),DIMENSION(:,:,:),ALLOCATABLE :: ztmp_fly_dta,zfieldo ! temporary array of values on input grid
	1272	INTEGER, DIMENSION(3) :: rec1,recn ! temporary arrays for start and length
	1273	INTEGER, DIMENSION(3) :: rec1_lsm,recn_lsm ! temporary arrays for start and length in case of seaoverland
	1274	INTEGER :: ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2 ! temporary indices
	1275	INTEGER :: jk, jn, jm, jir, jjr ! loop counters
	1276	INTEGER :: ni, nj ! lengths
	1277	INTEGER :: jpimin,jpiwid ! temporary indices
	1278	INTEGER :: jpimin_lsm,jpiwid_lsm ! temporary indices
	1279	INTEGER :: jpjmin,jpjwid ! temporary indices
	1280	INTEGER :: jpjmin_lsm,jpjwid_lsm ! temporary indices
	1281	INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices
	1282	INTEGER :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices
	1283	INTEGER :: itmpi,itmpj,itmpz ! lengths
	1284
[1275]	1285	!!----------------------------------------------------------------------
	1286	!
	1287	!! for weighted interpolation we have weights at four corners of a box surrounding
	1288	!! a model grid point, each weight is multiplied by a grid value (bilinear case)
	1289	!! or by a grid value and gradients at the corner point (bicubic case)
	1290	!! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
	1291
[2528]	1292	!! sub grid from non-model input grid which encloses all grid points in this nemo process
[1275]	1293	jpimin = ref_wgts(kw)%botleft(1)
	1294	jpjmin = ref_wgts(kw)%botleft(2)
	1295	jpiwid = ref_wgts(kw)%jpiwgt
	1296	jpjwid = ref_wgts(kw)%jpjwgt
	1297
[2528]	1298	!! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients
[1275]	1299	rec1(1) = MAX( jpimin-1, 1 )
	1300	rec1(2) = MAX( jpjmin-1, 1 )
[2528]	1301	rec1(3) = 1
[1275]	1302	recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
	1303	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
[2528]	1304	recn(3) = kk
[1275]	1305
[2528]	1306	!! where we need to put it in the non-nemo grid fly_dta
	1307	!! note that jpi1 and jpj1 only differ from 1 when jpimin and jpjmin are 1
	1308	!! (ie at the extreme west or south of the whole input grid) and similarly for jpi2 and jpj2
[1275]	1309	jpi1 = 2 + rec1(1) - jpimin
	1310	jpj1 = 2 + rec1(2) - jpjmin
	1311	jpi2 = jpi1 + recn(1) - 1
	1312	jpj2 = jpj1 + recn(2) - 1
	1313
	1314
[4230]	1315	IF( LEN( TRIM(lsmfile) ) > 0 ) THEN
	1316	!! indeces for ztmp_fly_dta
	1317	! --------------------------
	1318	rec1_lsm(1)=MAX(rec1(1)-nn_lsm,1) ! starting index for enlarged external data, x direction
	1319	rec1_lsm(2)=MAX(rec1(2)-nn_lsm,1) ! starting index for enlarged external data, y direction
	1320	rec1_lsm(3) = 1 ! vertical dimension
	1321	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
	1322	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
	1323	recn_lsm(3) = kk ! number of vertical levels in the input file
	1324
	1325	! Avoid out of bound
	1326	jpimin_lsm = MAX( rec1_lsm(1)+1, 1 )
	1327	jpjmin_lsm = MAX( rec1_lsm(2)+1, 1 )
	1328	jpiwid_lsm = MIN( recn_lsm(1)-2,ref_wgts(kw)%ddims(1)-rec1(1)+1)
	1329	jpjwid_lsm = MIN( recn_lsm(2)-2,ref_wgts(kw)%ddims(2)-rec1(2)+1)
	1330
	1331	jpi1_lsm = 2+rec1_lsm(1)-jpimin_lsm
	1332	jpj1_lsm = 2+rec1_lsm(2)-jpjmin_lsm
	1333	jpi2_lsm = jpi1_lsm + recn_lsm(1) - 1
	1334	jpj2_lsm = jpj1_lsm + recn_lsm(2) - 1
	1335
	1336
	1337	itmpi=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),1)
	1338	itmpj=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),2)
	1339	itmpz=kk
	1340	ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz))
	1341	ztmp_fly_dta(:,:,:) = 0.0
	1342	SELECT CASE( SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
	1343	CASE(1)
	1344	CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), &
	1345	& nrec, rec1_lsm, recn_lsm)
	1346	CASE DEFAULT
	1347	CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), &
	1348	& nrec, rec1_lsm, recn_lsm)
	1349	END SELECT
	1350	CALL apply_seaoverland(lsmfile,ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), &
	1351	& jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm, &
	1352	& itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
	1353
	1354
	1355	! Relative indeces for remapping
	1356	ii_lsm1 = (rec1(1)-rec1_lsm(1))+1
	1357	ii_lsm2 = (ii_lsm1+recn(1))-1
	1358	ij_lsm1 = (rec1(2)-rec1_lsm(2))+1
	1359	ij_lsm2 = (ij_lsm1+recn(2))-1
	1360
	1361	ref_wgts(kw)%fly_dta(:,:,:) = 0.0
	1362	ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:) = ztmp_fly_dta(ii_lsm1:ii_lsm2,ij_lsm1:ij_lsm2,:)
	1363	DEALLOCATE(ztmp_fly_dta)
	1364
	1365	ELSE
	1366
	1367	ref_wgts(kw)%fly_dta(:,:,:) = 0.0
	1368	SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) )
	1369	CASE(1)
	1370	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn)
	1371	CASE DEFAULT
	1372	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn)
	1373	END SELECT
	1374	ENDIF
	1375
	1376
[1275]	1377	!! first four weights common to both bilinear and bicubic
[2528]	1378	!! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft
[1275]	1379	!! note that we have to offset by 1 into fly_dta array because of halo
[2528]	1380	dta(:,:,:) = 0.0
[1275]	1381	DO jk = 1,4
	1382	DO jn = 1, jpj
	1383	DO jm = 1,jpi
	1384	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
	1385	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
[2528]	1386	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,:)
[1275]	1387	END DO
	1388	END DO
	1389	END DO
	1390
	1391	IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
	1392
	1393	!! fix up halo points that we couldnt read from file
	1394	IF( jpi1 == 2 ) THEN
[2528]	1395	ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:)
[1275]	1396	ENDIF
	1397	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
[2528]	1398	ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:)
[1275]	1399	ENDIF
	1400	IF( jpj1 == 2 ) THEN
[2528]	1401	ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:)
[1275]	1402	ENDIF
	1403	IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
[2528]	1404	ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:)
[1275]	1405	ENDIF
	1406
	1407	!! if data grid is cyclic we can do better on east-west edges
	1408	!! but have to allow for whether first and last columns are coincident
	1409	IF( ref_wgts(kw)%cyclic ) THEN
	1410	rec1(2) = MAX( jpjmin-1, 1 )
[2528]	1411	recn(1) = 1
[1275]	1412	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
	1413	jpj1 = 2 + rec1(2) - jpjmin
	1414	jpj2 = jpj1 + recn(2) - 1
	1415	IF( jpi1 == 2 ) THEN
[2528]	1416	rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap
	1417	SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
	1418	CASE(1)
	1419	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
	1420	CASE DEFAULT
	1421	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
	1422	END SELECT
	1423	ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
[1275]	1424	ENDIF
	1425	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
[2528]	1426	rec1(1) = 1 + ref_wgts(kw)%overlap
	1427	SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
	1428	CASE(1)
	1429	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
	1430	CASE DEFAULT
	1431	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
	1432	END SELECT
	1433	ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
[1275]	1434	ENDIF
	1435	ENDIF
	1436
	1437	! gradient in the i direction
	1438	DO jk = 1,4
	1439	DO jn = 1, jpj
	1440	DO jm = 1,jpi
	1441	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
	1442	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
[2528]	1443	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 * &
	1444	(ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:))
[1275]	1445	END DO
	1446	END DO
	1447	END DO
	1448
	1449	! gradient in the j direction
	1450	DO jk = 1,4
	1451	DO jn = 1, jpj
	1452	DO jm = 1,jpi
	1453	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
	1454	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
[2528]	1455	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 * &
	1456	(ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:))
[1275]	1457	END DO
	1458	END DO
	1459	END DO
	1460
[2715]	1461	! gradient in the ij direction
	1462	DO jk = 1,4
	1463	DO jn = 1, jpj
	1464	DO jm = 1,jpi
	1465	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
	1466	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
	1467	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
[2528]	1468	(ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni ,nj+2,:)) - &
	1469	(ref_wgts(kw)%fly_dta(ni+2,nj ,:) - ref_wgts(kw)%fly_dta(ni ,nj ,:)))
[2715]	1470	END DO
[1275]	1471	END DO
[2715]	1472	END DO
	1473	!
[1275]	1474	END IF
[2715]	1475	!
[1275]	1476	END SUBROUTINE fld_interp
[2528]	1477
	1478
	1479	FUNCTION ksec_week( cdday )
	1480	!!---------------------------------------------------------------------
	1481	!! * FUNCTION kshift_week *
	1482	!!
	1483	!! ** Purpose :
	1484	!!---------------------------------------------------------------------
	1485	CHARACTER(len=*), INTENT(in) :: cdday !3 first letters of the first day of the weekly file
	1486	!!
	1487	INTEGER :: ksec_week ! output variable
	1488	INTEGER :: ijul !temp variable
	1489	INTEGER :: ishift !temp variable
	1490	CHARACTER(len=3),DIMENSION(7) :: cl_week
	1491	!!----------------------------------------------------------------------
	1492	cl_week = (/"sun","sat","fri","thu","wed","tue","mon"/)
	1493	DO ijul = 1, 7
	1494	IF( cl_week(ijul) == TRIM(cdday) ) EXIT
[2715]	1495	END DO
[2528]	1496	IF( ijul .GT. 7 ) CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): '//TRIM(cdday) )
	1497	!
	1498	ishift = ijul * NINT(rday)
	1499	!
	1500	ksec_week = nsec_week + ishift
	1501	ksec_week = MOD( ksec_week, 7*NINT(rday) )
	1502	!
	1503	END FUNCTION ksec_week
	1504
[2715]	1505	!!======================================================================
[888]	1506	END MODULE fldread

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

Download in other formats: