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/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 4694

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

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

Download in other formats: