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

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

fldread.F90 in NEMO/trunk/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/SBC/fldread.F90 @ 11536

Last change on this file since 11536 was 11536, checked in by smasson, 5 years ago
trunk: merge dev_r10984_HPC-13 into the trunk
Property svn:keywords set to `Id`
File size: 101.3 KB

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

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

Download in other formats: