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

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

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

Context Navigation

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

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

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

Download in other formats: