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 @ 12248

Last change on this file since 12248 was 12248, checked in by smasson, 4 years ago
trunk : clean useless warnings in fldread
Property svn:keywords set to `Id`
File size: 100.0 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	SELECT CASE( kgrd )
	868	CASE(1)
	869	! depth of T points:
	870	zdepth(:) = gdept_n(ji,jj,:)
	871	CASE(2)
	872	! depth of U points: we must not use gdept_n as we don't want to do a communication
	873	! --> copy what is done for gdept_n in domvvl...
	874	zdhalf(1) = 0.0_wp
	875	zdepth(1) = 0.5_wp * e3uw_n(ji,jj,1)
	876	DO jk = 2, jpk ! vertical sum
	877	! zcoef = umask - wumask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
	878	! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
	879	! ! 0.5 where jk = mikt
	880	!!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ??
	881	zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) )
	882	zdhalf(jk) = zdhalf(jk-1) + e3u_n(ji,jj,jk-1)
	883	zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5 * e3uw_n(ji,jj,jk)) &
	884	& + (1-zcoef) * ( zdepth(jk-1) + e3uw_n(ji,jj,jk))
	885	END DO
	886	CASE(3)
	887	! depth of V points: we must not use gdept_n as we don't want to do a communication
	888	! --> copy what is done for gdept_n in domvvl...
	889	zdhalf(1) = 0.0_wp
	890	zdepth(1) = 0.5_wp * e3vw_n(ji,jj,1)
	891	DO jk = 2, jpk ! vertical sum
	892	! zcoef = vmask - wvmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
	893	! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
	894	! ! 0.5 where jk = mikt
	895	!!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ??
	896	zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) )
	897	zdhalf(jk) = zdhalf(jk-1) + e3v_n(ji,jj,jk-1)
	898	zdepth(jk) = zcoef * ( zdhalf(jk ) + 0.5 * e3vw_n(ji,jj,jk)) &
	899	& + (1-zcoef) * ( zdepth(jk-1) + e3vw_n(ji,jj,jk))
	900	END DO
	901	END SELECT
	902	!
	903	DO jk = 1, jpk
	904	IF( zdepth(jk) < pdta_read_z(jb,1, 1) ) THEN ! above the first level of external data
	905	pdta(jb,1,jk) = pdta_read(jb,1,1)
	906	ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkb) ) THEN ! below the last level of external data
	907	pdta(jb,1,jk) = pdta_read(jb,1,MAXLOC(pdta_read_z(jb,1,:),1))
	908	ELSE ! inbetween: vertical interpolation between jkb & jkb+1
	909	DO jkb = 1, ipkb-1 ! when gdept_n(jkb) < zdepth(jk) < gdept_n(jkb+1)
	910	IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) <= 0._wp ) &
	911	& .AND. ( pdta_read_z(jb,1,jkb+1) /= zfv_alt) ) THEN ! linear interpolation between 2 levels
	912	zi = ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) / ( pdta_read_z(jb,1,jkb+1) - pdta_read_z(jb,1,jkb) )
	913	pdta(jb,1,jk) = pdta_read(jb,1,jkb) + ( pdta_read (jb,1,jkb+1) - pdta_read (jb,1,jkb) ) * zi
[7646]	914	ENDIF
[11536]	915	END DO
	916	ENDIF
	917	END DO ! jpk
	918	!
	919	END DO ! ipi
	920
	921	IF(kgrd == 2) THEN ! do we need to adjust the transport term?
	922	DO jb = 1, ipi
	923	ji = idx_bdy(kbdy)%nbi(jb,kgrd)
	924	jj = idx_bdy(kbdy)%nbj(jb,kgrd)
	925	zh = SUM(pdta_read_dz(jb,1,:) )
	926	ztrans = 0._wp
	927	ztrans_new = 0._wp
	928	DO jkb = 1, ipkb ! calculate transport on input grid
	929	ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)
[7646]	930	ENDDO
[11536]	931	DO jk = 1, jpk ! calculate transport on model grid
	932	ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3u_n(ji,jj,jk ) * umask(ji,jj,jk)
[7646]	933	ENDDO
[11536]	934	DO jk = 1, jpk ! make transport correction
	935	IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data
	936	pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu_n(ji,jj) ) * umask(ji,jj,jk)
	937	ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
	938	pdta(jb,1,jk) = pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hu_n(ji,jj) * umask(ji,jj,jk)
[7646]	939	ENDIF
	940	ENDDO
[11536]	941	ENDDO
	942	ENDIF
	943
	944	IF(kgrd == 3) THEN ! do we need to adjust the transport term?
	945	DO jb = 1, ipi
	946	ji = idx_bdy(kbdy)%nbi(jb,kgrd)
	947	jj = idx_bdy(kbdy)%nbj(jb,kgrd)
	948	zh = SUM(pdta_read_dz(jb,1,:) )
	949	ztrans = 0._wp
	950	ztrans_new = 0._wp
	951	DO jkb = 1, ipkb ! calculate transport on input grid
	952	ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)
	953	ENDDO
	954	DO jk = 1, jpk ! calculate transport on model grid
	955	ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3v_n(ji,jj,jk ) * vmask(ji,jj,jk)
	956	ENDDO
	957	DO jk = 1, jpk ! make transport correction
	958	IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data
	959	pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv_n(ji,jj) ) * vmask(ji,jj,jk)
	960	ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
	961	pdta(jb,1,jk) = pdta(jb,1,jk) + ( 0._wp - ztrans_new ) * r1_hv_n(ji,jj) * vmask(ji,jj,jk)
[7646]	962	ENDIF
	963	ENDDO
[11536]	964	ENDDO
	965	ENDIF
	966
[7646]	967	END SUBROUTINE fld_bdy_interp
	968
	969
[2528]	970	SUBROUTINE fld_rot( kt, sd )
	971	!!---------------------------------------------------------------------
[3294]	972	!! * ROUTINE fld_rot *
[2528]	973	!!
	974	!! ** Purpose : Vector fields may need to be rotated onto the local grid direction
[2715]	975	!!----------------------------------------------------------------------
[6140]	976	INTEGER , INTENT(in ) :: kt ! ocean time step
	977	TYPE(FLD), DIMENSION(:), INTENT(inout) :: sd ! input field related variables
	978	!
	979	INTEGER :: ju, jv, jk, jn ! loop indices
	980	INTEGER :: imf ! size of the structure sd
	981	INTEGER :: ill ! character length
	982	INTEGER :: iv ! indice of V component
[9125]	983	CHARACTER (LEN=100) :: clcomp ! dummy weight name
	984	REAL(wp), DIMENSION(jpi,jpj) :: utmp, vtmp ! temporary arrays for vector rotation
[2528]	985	!!---------------------------------------------------------------------
[6140]	986	!
[2528]	987	!! (sga: following code should be modified so that pairs arent searched for each time
	988	!
	989	imf = SIZE( sd )
	990	DO ju = 1, imf
[11536]	991	IF( TRIM(sd(ju)%clrootname) == 'NOT USED' ) CYCLE
[2528]	992	ill = LEN_TRIM( sd(ju)%vcomp )
[3851]	993	DO jn = 2-COUNT((/sd(ju)%ln_tint/)), 2
	994	IF( ill > 0 .AND. .NOT. sd(ju)%rotn(jn) ) THEN ! find vector rotations required
	995	IF( sd(ju)%vcomp(1:1) == 'U' ) THEN ! east-west component has symbolic name starting with 'U'
	996	! look for the north-south component which has same symbolic name but with 'U' replaced with 'V'
	997	clcomp = 'V' // sd(ju)%vcomp(2:ill) ! works even if ill == 1
	998	iv = -1
	999	DO jv = 1, imf
[11536]	1000	IF( TRIM(sd(jv)%clrootname) == 'NOT USED' ) CYCLE
[3851]	1001	IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) ) iv = jv
	1002	END DO
	1003	IF( iv > 0 ) THEN ! fields ju and iv are two components which need to be rotated together
	1004	DO jk = 1, SIZE( sd(ju)%fnow, 3 )
	1005	IF( sd(ju)%ln_tint )THEN
	1006	CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->i', utmp(:,:) )
	1007	CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->j', vtmp(:,:) )
	1008	sd(ju)%fdta(:,:,jk,jn) = utmp(:,:) ; sd(iv)%fdta(:,:,jk,jn) = vtmp(:,:)
	1009	ELSE
	1010	CALL rot_rep( sd(ju)%fnow(:,:,jk ), sd(iv)%fnow(:,:,jk ), 'T', 'en->i', utmp(:,:) )
	1011	CALL rot_rep( sd(ju)%fnow(:,:,jk ), sd(iv)%fnow(:,:,jk ), 'T', 'en->j', vtmp(:,:) )
	1012	sd(ju)%fnow(:,:,jk ) = utmp(:,:) ; sd(iv)%fnow(:,:,jk ) = vtmp(:,:)
	1013	ENDIF
	1014	END DO
	1015	sd(ju)%rotn(jn) = .TRUE. ! vector was rotated
	1016	IF( lwp .AND. kt == nit000 ) WRITE(numout,*) &
	1017	& 'fld_read: vector pair ('//TRIM(sd(ju)%clvar)//', '//TRIM(sd(iv)%clvar)//') rotated on to model grid'
	1018	ENDIF
	1019	ENDIF
	1020	ENDIF
	1021	END DO
[2528]	1022	END DO
[2715]	1023	!
[2528]	1024	END SUBROUTINE fld_rot
	1025
	1026
[1628]	1027	SUBROUTINE fld_clopn( sdjf, kyear, kmonth, kday, ldstop )
[1132]	1028	!!---------------------------------------------------------------------
	1029	!! * ROUTINE fld_clopn *
	1030	!!
[7646]	1031	!! ** Purpose : update the file name and close/open the files
[1132]	1032	!!----------------------------------------------------------------------
[2715]	1033	TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables
[3851]	1034	INTEGER, OPTIONAL, INTENT(in ) :: kyear ! year value
	1035	INTEGER, OPTIONAL, INTENT(in ) :: kmonth ! month value
	1036	INTEGER, OPTIONAL, INTENT(in ) :: kday ! day value
[2715]	1037	LOGICAL, OPTIONAL, INTENT(in ) :: ldstop ! stop if open to read a non-existing file (default = .TRUE.)
[6140]	1038	!
[11536]	1039	LOGICAL :: llprevyr ! are we reading previous year file?
	1040	LOGICAL :: llprevmth ! are we reading previous month file?
	1041	INTEGER :: iyear, imonth, iday ! first day of the current file in yyyy mm dd
	1042	INTEGER :: isec_week ! number of seconds since start of the weekly file
	1043	INTEGER :: indexyr ! year undex (O/1/2: previous/current/next)
	1044	REAL(wp) :: zyear_len, zmonth_len ! length (days) of iyear and imonth !
	1045	CHARACTER(len = 256) :: clname ! temporary file name
[2715]	1046	!!----------------------------------------------------------------------
[3851]	1047	IF( PRESENT(kyear) ) THEN ! use given values
	1048	iyear = kyear
	1049	imonth = kmonth
	1050	iday = kday
[5751]	1051	IF ( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of the week
[5749]	1052	isec_week = ksec_week( sdjf%cltype(6:8) )- (86400 * 8 )
[5751]	1053	llprevmth = isec_week > nsec_month ! longer time since beginning of the week than the month
[5749]	1054	llprevyr = llprevmth .AND. nmonth == 1
	1055	iyear = nyear - COUNT((/llprevyr /))
	1056	imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
[5751]	1057	iday = nday + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
[5749]	1058	ENDIF
[3851]	1059	ELSE ! use current day values
	1060	IF ( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of the week
	1061	isec_week = ksec_week( sdjf%cltype(6:8) ) ! second since the beginning of the week
	1062	llprevmth = isec_week > nsec_month ! longer time since beginning of the week than the month
	1063	llprevyr = llprevmth .AND. nmonth == 1
	1064	ELSE
	1065	isec_week = 0
	1066	llprevmth = .FALSE.
	1067	llprevyr = .FALSE.
	1068	ENDIF
	1069	iyear = nyear - COUNT((/llprevyr /))
	1070	imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
	1071	iday = nday + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
	1072	ENDIF
[1132]	1073
	1074	! build the new filename if not climatological data
[3851]	1075	clname=TRIM(sdjf%clrootname)
[2528]	1076	!
[3851]	1077	! note that sdjf%ln_clim is is only acting on the presence of the year in the file name
[2528]	1078	IF( .NOT. sdjf%ln_clim ) THEN
[3851]	1079	WRITE(clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), iyear ! add year
	1080	IF( sdjf%cltype /= 'yearly' ) WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname ), imonth ! add month
[2528]	1081	ELSE
	1082	! build the new filename if climatological data
[3851]	1083	IF( sdjf%cltype /= 'yearly' ) WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), imonth ! add month
[888]	1084	ENDIF
[2528]	1085	IF( sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) &
[3851]	1086	& WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname ), iday ! add day
[2528]	1087	!
[3851]	1088	IF( TRIM(clname) /= TRIM(sdjf%clname) .OR. sdjf%num == 0 ) THEN ! new file to be open
[6140]	1089	!
[3851]	1090	sdjf%clname = TRIM(clname)
	1091	IF( sdjf%num /= 0 ) CALL iom_close( sdjf%num ) ! close file if already open
	1092	CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 )
[6140]	1093	!
[3851]	1094	! find the last record to be read -> update sdjf%nreclast
	1095	indexyr = iyear - nyear + 1
[11536]	1096	zyear_len = REAL(nyear_len( indexyr ), wp)
[3851]	1097	SELECT CASE ( indexyr )
[11536]	1098	CASE ( 0 ) ; zmonth_len = 31. ! previous year -> imonth = 12
	1099	CASE ( 1 ) ; zmonth_len = REAL(nmonth_len(imonth), wp)
	1100	CASE ( 2 ) ; zmonth_len = 31. ! next year -> imonth = 1
[3851]	1101	END SELECT
[6140]	1102	!
[3851]	1103	! last record to be read in the current file
[11536]	1104	IF ( sdjf%freqh == -12. ) THEN ; sdjf%nreclast = 1 ! yearly mean
	1105	ELSEIF( sdjf%freqh == -1. ) THEN ! monthly mean
[3851]	1106	IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nreclast = 1
	1107	ELSE ; sdjf%nreclast = 12
	1108	ENDIF
	1109	ELSE ! higher frequency mean (in hours)
[11536]	1110	IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nreclast = NINT( 24. * zmonth_len / sdjf%freqh )
	1111	ELSEIF( sdjf%cltype(1:4) == 'week' ) THEN ; sdjf%nreclast = NINT( 24. * 7. / sdjf%freqh )
	1112	ELSEIF( sdjf%cltype == 'daily' ) THEN ; sdjf%nreclast = NINT( 24. / sdjf%freqh )
	1113	ELSE ; sdjf%nreclast = NINT( 24. * zyear_len / sdjf%freqh )
[3851]	1114	ENDIF
	1115	ENDIF
[6140]	1116	!
[3851]	1117	ENDIF
	1118	!
[1132]	1119	END SUBROUTINE fld_clopn
	1120
	1121
[7646]	1122	SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam, knoprint )
[1132]	1123	!!---------------------------------------------------------------------
	1124	!! * ROUTINE fld_fill *
	1125	!!
[7646]	1126	!! ** Purpose : fill the data structure (sdf) with the associated information
	1127	!! read in namelist (sdf_n) and control print
[1132]	1128	!!----------------------------------------------------------------------
[7646]	1129	TYPE(FLD) , DIMENSION(:) , INTENT(inout) :: sdf ! structure of input fields (file informations, fields read)
	1130	TYPE(FLD_N), DIMENSION(:) , INTENT(in ) :: sdf_n ! array of namelist information structures
	1131	CHARACTER(len=*) , INTENT(in ) :: cdir ! Root directory for location of flx files
	1132	CHARACTER(len=*) , INTENT(in ) :: cdcaller ! name of the calling routine
	1133	CHARACTER(len=*) , INTENT(in ) :: cdtitle ! description of the calling routine
	1134	CHARACTER(len=*) , INTENT(in ) :: cdnam ! name of the namelist from which sdf_n comes
	1135	INTEGER , OPTIONAL, INTENT(in ) :: knoprint ! no calling routine information printed
[888]	1136	!
[7646]	1137	INTEGER :: jf ! dummy indices
[1132]	1138	!!---------------------------------------------------------------------
[6140]	1139	!
[1132]	1140	DO jf = 1, SIZE(sdf)
[11536]	1141	sdf(jf)%clrootname = sdf_n(jf)%clname
	1142	IF( TRIM(sdf_n(jf)%clname) /= 'NOT USED' ) sdf(jf)%clrootname = TRIM( cdir )//sdf(jf)%clrootname
[3851]	1143	sdf(jf)%clname = "not yet defined"
[11536]	1144	sdf(jf)%freqh = sdf_n(jf)%freqh
[1132]	1145	sdf(jf)%clvar = sdf_n(jf)%clvar
	1146	sdf(jf)%ln_tint = sdf_n(jf)%ln_tint
	1147	sdf(jf)%ln_clim = sdf_n(jf)%ln_clim
[2528]	1148	sdf(jf)%cltype = sdf_n(jf)%cltype
[3851]	1149	sdf(jf)%num = -1
	1150	sdf(jf)%wgtname = " "
[11536]	1151	IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) sdf(jf)%wgtname = TRIM( cdir )//sdf_n(jf)%wname
[4230]	1152	sdf(jf)%lsmname = " "
[11536]	1153	IF( LEN( TRIM(sdf_n(jf)%lname) ) > 0 ) sdf(jf)%lsmname = TRIM( cdir )//sdf_n(jf)%lname
[3851]	1154	sdf(jf)%vcomp = sdf_n(jf)%vcomp
	1155	sdf(jf)%rotn(:) = .TRUE. ! pretend to be rotated -> won't try to rotate data before the first call to fld_get
	1156	IF( sdf(jf)%cltype(1:4) == 'week' .AND. nn_leapy == 0 ) &
	1157	& CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs nn_leapy = 1')
	1158	IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim ) &
	1159	& CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.')
[11536]	1160	sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
	1161	sdf(jf)%igrd = 0
	1162	sdf(jf)%ibdy = 0
	1163	sdf(jf)%imap => NULL()
	1164	sdf(jf)%ltotvel = .FALSE.
	1165	sdf(jf)%lzint = .FALSE.
[1132]	1166	END DO
[6140]	1167	!
[1132]	1168	IF(lwp) THEN ! control print
	1169	WRITE(numout,*)
[7646]	1170	IF( .NOT.PRESENT( knoprint) ) THEN
	1171	WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
	1172	WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
	1173	ENDIF
	1174	WRITE(numout,*) ' fld_fill : fill data structure with information from namelist '//TRIM( cdnam )
	1175	WRITE(numout,*) ' ~~~~~~~~'
	1176	WRITE(numout,*) ' list of files and frequency (>0: in hours ; <0 in months)'
[1132]	1177	DO jf = 1, SIZE(sdf)
[7646]	1178	WRITE(numout,*) ' root filename: ' , TRIM( sdf(jf)%clrootname ), ' variable name: ', TRIM( sdf(jf)%clvar )
[11536]	1179	WRITE(numout,*) ' frequency: ' , sdf(jf)%freqh , &
[7646]	1180	& ' time interp: ' , sdf(jf)%ln_tint , &
	1181	& ' climatology: ' , sdf(jf)%ln_clim
	1182	WRITE(numout,*) ' weights: ' , TRIM( sdf(jf)%wgtname ), &
	1183	& ' pairing: ' , TRIM( sdf(jf)%vcomp ), &
	1184	& ' data type: ' , sdf(jf)%cltype , &
	1185	& ' land/sea mask:' , TRIM( sdf(jf)%lsmname )
[2528]	1186	call flush(numout)
[1132]	1187	END DO
	1188	ENDIF
[6140]	1189	!
[1132]	1190	END SUBROUTINE fld_fill
	1191
	1192
[1275]	1193	SUBROUTINE wgt_list( sd, kwgt )
	1194	!!---------------------------------------------------------------------
	1195	!! * ROUTINE wgt_list *
	1196	!!
[7646]	1197	!! ** Purpose : search array of WGTs and find a weights file entry,
	1198	!! or return a new one adding it to the end if new entry.
	1199	!! the weights data is read in and restructured (fld_weight)
[1275]	1200	!!----------------------------------------------------------------------
[2715]	1201	TYPE( FLD ), INTENT(in ) :: sd ! field with name of weights file
	1202	INTEGER , INTENT(inout) :: kwgt ! index of weights
[6140]	1203	!
[2715]	1204	INTEGER :: kw, nestid ! local integer
	1205	LOGICAL :: found ! local logical
[1275]	1206	!!----------------------------------------------------------------------
	1207	!
	1208	!! search down linked list
	1209	!! weights filename is either present or we hit the end of the list
	1210	found = .FALSE.
[6140]	1211	!
[1275]	1212	!! because agrif nest part of filenames are now added in iom_open
	1213	!! to distinguish between weights files on the different grids, need to track
	1214	!! nest number explicitly
	1215	nestid = 0
	1216	#if defined key_agrif
	1217	nestid = Agrif_Fixed()
	1218	#endif
	1219	DO kw = 1, nxt_wgt-1
	1220	IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
	1221	ref_wgts(kw)%nestid == nestid) THEN
	1222	kwgt = kw
	1223	found = .TRUE.
	1224	EXIT
	1225	ENDIF
	1226	END DO
	1227	IF( .NOT.found ) THEN
	1228	kwgt = nxt_wgt
	1229	CALL fld_weight( sd )
	1230	ENDIF
[2715]	1231	!
[1275]	1232	END SUBROUTINE wgt_list
	1233
[2715]	1234
[1275]	1235	SUBROUTINE wgt_print( )
	1236	!!---------------------------------------------------------------------
	1237	!! * ROUTINE wgt_print *
	1238	!!
	1239	!! ** Purpose : print the list of known weights
	1240	!!----------------------------------------------------------------------
[2715]	1241	INTEGER :: kw !
[1275]	1242	!!----------------------------------------------------------------------
	1243	!
	1244	DO kw = 1, nxt_wgt-1
	1245	WRITE(numout,*) 'weight file: ',TRIM(ref_wgts(kw)%wgtname)
	1246	WRITE(numout,*) ' ddims: ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
	1247	WRITE(numout,*) ' numwgt: ',ref_wgts(kw)%numwgt
	1248	WRITE(numout,*) ' jpiwgt: ',ref_wgts(kw)%jpiwgt
	1249	WRITE(numout,*) ' jpjwgt: ',ref_wgts(kw)%jpjwgt
	1250	WRITE(numout,*) ' botleft: ',ref_wgts(kw)%botleft
	1251	WRITE(numout,*) ' topright: ',ref_wgts(kw)%topright
	1252	IF( ref_wgts(kw)%cyclic ) THEN
	1253	WRITE(numout,*) ' cyclical'
[2528]	1254	IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) ' with overlap of ', ref_wgts(kw)%overlap
[1275]	1255	ELSE
	1256	WRITE(numout,*) ' not cyclical'
	1257	ENDIF
	1258	IF( ASSOCIATED(ref_wgts(kw)%data_wgt) ) WRITE(numout,*) ' allocated'
	1259	END DO
[2715]	1260	!
[1275]	1261	END SUBROUTINE wgt_print
	1262
[2715]	1263
[1275]	1264	SUBROUTINE fld_weight( sd )
	1265	!!---------------------------------------------------------------------
	1266	!! * ROUTINE fld_weight *
	1267	!!
[7646]	1268	!! ** Purpose : create a new WGT structure and fill in data from file,
	1269	!! restructuring as required
[1275]	1270	!!----------------------------------------------------------------------
[2715]	1271	TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file
	1272	!!
[6140]	1273	INTEGER :: jn ! dummy loop indices
	1274	INTEGER :: inum ! local logical unit
	1275	INTEGER :: id ! local variable id
	1276	INTEGER :: ipk ! local vertical dimension
	1277	INTEGER :: zwrap ! local integer
	1278	LOGICAL :: cyclical !
	1279	CHARACTER (len=5) :: aname !
[5399]	1280	INTEGER , DIMENSION(:), ALLOCATABLE :: ddims
[9125]	1281	INTEGER, DIMENSION(jpi,jpj) :: data_src
	1282	REAL(wp), DIMENSION(jpi,jpj) :: data_tmp
[1275]	1283	!!----------------------------------------------------------------------
	1284	!
	1285	IF( nxt_wgt > tot_wgts ) THEN
[2777]	1286	CALL ctl_stop("fld_weight: weights array size exceeded, increase tot_wgts")
[1275]	1287	ENDIF
	1288	!
	1289	!! new weights file entry, add in extra information
	1290	!! a weights file represents a 2D grid of a certain shape, so we assume that the current
	1291	!! input data file is representative of all other files to be opened and processed with the
	1292	!! current weights file
	1293
	1294	!! open input data file (non-model grid)
[1319]	1295	CALL iom_open( sd%clname, inum, ldiof = LEN(TRIM(sd%wgtname)) > 0 )
[1275]	1296
	1297	!! get dimensions
[5399]	1298	IF ( SIZE(sd%fnow, 3) > 1 ) THEN
	1299	ALLOCATE( ddims(4) )
	1300	ELSE
	1301	ALLOCATE( ddims(3) )
	1302	ENDIF
[1275]	1303	id = iom_varid( inum, sd%clvar, ddims )
	1304
[2528]	1305	!! close it
	1306	CALL iom_close( inum )
[1275]	1307
[2528]	1308	!! now open the weights file
[1275]	1309
[2528]	1310	CALL iom_open ( sd%wgtname, inum ) ! interpolation weights
	1311	IF ( inum > 0 ) THEN
[1275]	1312
[2528]	1313	!! determine whether we have an east-west cyclic grid
	1314	!! from global attribute called "ew_wrap" in the weights file
	1315	!! note that if not found, iom_getatt returns -999 and cyclic with no overlap is assumed
	1316	!! since this is the most common forcing configuration
[1275]	1317
[2528]	1318	CALL iom_getatt(inum, 'ew_wrap', zwrap)
	1319	IF( zwrap >= 0 ) THEN
[1275]	1320	cyclical = .TRUE.
[2528]	1321	ELSE IF( zwrap == -999 ) THEN
[1275]	1322	cyclical = .TRUE.
[2528]	1323	zwrap = 0
	1324	ELSE
	1325	cyclical = .FALSE.
[1275]	1326	ENDIF
	1327
	1328	ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
	1329	ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
	1330	ref_wgts(nxt_wgt)%wgtname = sd%wgtname
[2528]	1331	ref_wgts(nxt_wgt)%overlap = zwrap
	1332	ref_wgts(nxt_wgt)%cyclic = cyclical
[1275]	1333	ref_wgts(nxt_wgt)%nestid = 0
	1334	#if defined key_agrif
	1335	ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
	1336	#endif
	1337	!! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
	1338	!! for each weight wgtNN there is an integer array srcNN which gives the point in
	1339	!! the input data grid which is to be multiplied by the weight
	1340	!! they are both arrays on the model grid so the result of the multiplication is
	1341	!! added into an output array on the model grid as a running sum
	1342
	1343	!! two possible cases: bilinear (4 weights) or bicubic (16 weights)
	1344	id = iom_varid(inum, 'src05', ldstop=.FALSE.)
	1345	IF( id <= 0) THEN
	1346	ref_wgts(nxt_wgt)%numwgt = 4
	1347	ELSE
	1348	ref_wgts(nxt_wgt)%numwgt = 16
	1349	ENDIF
	1350
	1351	ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
	1352	ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
	1353	ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
	1354
	1355	DO jn = 1,4
	1356	aname = ' '
	1357	WRITE(aname,'(a3,i2.2)') 'src',jn
	1358	data_tmp(:,:) = 0
[1955]	1359	CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
[1275]	1360	data_src(:,:) = INT(data_tmp(:,:))
	1361	ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
	1362	ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
	1363	END DO
	1364
	1365	DO jn = 1, ref_wgts(nxt_wgt)%numwgt
	1366	aname = ' '
	1367	WRITE(aname,'(a3,i2.2)') 'wgt',jn
[1955]	1368	ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
	1369	CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
[1275]	1370	END DO
	1371	CALL iom_close (inum)
	1372
	1373	! find min and max indices in grid
[1955]	1374	ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
	1375	ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
	1376	ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
	1377	ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
[1275]	1378
	1379	! and therefore dimensions of the input box
	1380	ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
	1381	ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
	1382
	1383	! shift indexing of source grid
	1384	ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
	1385	ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
	1386
	1387	! create input grid, give it a halo to allow gradient calculations
[1702]	1388	! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
	1389	! a more robust solution will be given in next release
[2528]	1390	ipk = SIZE(sd%fnow, 3)
	1391	ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
	1392	IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
[6140]	1393	!
[1275]	1394	nxt_wgt = nxt_wgt + 1
[6140]	1395	!
[1275]	1396	ELSE
	1397	CALL ctl_stop( ' fld_weight : unable to read the file ' )
	1398	ENDIF
	1399
[5399]	1400	DEALLOCATE (ddims )
[2715]	1401	!
[1275]	1402	END SUBROUTINE fld_weight
	1403
[2715]	1404
[6140]	1405	SUBROUTINE apply_seaoverland( clmaskfile, zfieldo, jpi1_lsm, jpi2_lsm, jpj1_lsm, &
[7646]	1406	& jpj2_lsm, itmpi, itmpj, itmpz, rec1_lsm, recn_lsm )
[1275]	1407	!!---------------------------------------------------------------------
[4230]	1408	!! * ROUTINE apply_seaoverland *
	1409	!!
	1410	!! ** Purpose : avoid spurious fluxes in coastal or near-coastal areas
	1411	!! due to the wrong usage of "land" values from the coarse
	1412	!! atmospheric model when spatial interpolation is required
	1413	!! D. Delrosso INGV
	1414	!!----------------------------------------------------------------------
[6140]	1415	INTEGER, INTENT(in ) :: itmpi,itmpj,itmpz ! lengths
	1416	INTEGER, INTENT(in ) :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices
	1417	INTEGER, DIMENSION(3), INTENT(in ) :: rec1_lsm,recn_lsm ! temporary arrays for start and length
	1418	REAL(wp),DIMENSION (:,:,:),INTENT(inout) :: zfieldo ! input/output array for seaoverland application
	1419	CHARACTER (len=100), INTENT(in ) :: clmaskfile ! land/sea mask file name
	1420	!
	1421	INTEGER :: inum,jni,jnj,jnz,jc ! local indices
	1422	REAL(wp),DIMENSION (:,:,:),ALLOCATABLE :: zslmec1 ! local array for land point detection
	1423	REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfieldn ! array of forcing field with undeff for land points
	1424	REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfield ! array of forcing field
[4230]	1425	!!---------------------------------------------------------------------
[6140]	1426	!
	1427	ALLOCATE ( zslmec1(itmpi,itmpj,itmpz), zfieldn(itmpi,itmpj), zfield(itmpi,itmpj) )
	1428	!
[4230]	1429	! Retrieve the land sea mask data
	1430	CALL iom_open( clmaskfile, inum )
	1431	SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
	1432	CASE(1)
[6140]	1433	CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm)
[4230]	1434	CASE DEFAULT
[6140]	1435	CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm)
[4230]	1436	END SELECT
	1437	CALL iom_close( inum )
[6140]	1438	!
	1439	DO jnz=1,rec1_lsm(3) !! Loop over k dimension
	1440	!
	1441	DO jni = 1, itmpi !! copy the original field into a tmp array
	1442	DO jnj = 1, itmpj !! substituting undeff over land points
	1443	zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz)
	1444	IF( zslmec1(jni,jnj,jnz) == 1. ) zfieldn(jni,jnj) = undeff_lsm
[4230]	1445	END DO
	1446	END DO
[6140]	1447	!
	1448	CALL seaoverland( zfieldn, itmpi, itmpj, zfield )
	1449	DO jc = 1, nn_lsm
	1450	CALL seaoverland( zfield, itmpi, itmpj, zfield )
	1451	END DO
	1452	!
	1453	! Check for Undeff and substitute original values
	1454	IF( ANY(zfield==undeff_lsm) ) THEN
	1455	DO jni = 1, itmpi
	1456	DO jnj = 1, itmpj
	1457	IF( zfield(jni,jnj)==undeff_lsm ) zfield(jni,jnj) = zfieldo(jni,jnj,jnz)
	1458	END DO
	1459	END DO
	1460	ENDIF
	1461	!
	1462	zfieldo(:,:,jnz) = zfield(:,:)
	1463	!
	1464	END DO !! End Loop over k dimension
	1465	!
	1466	DEALLOCATE ( zslmec1, zfieldn, zfield )
	1467	!
[4230]	1468	END SUBROUTINE apply_seaoverland
	1469
	1470
[6140]	1471	SUBROUTINE seaoverland( zfieldn, ileni, ilenj, zfield )
[4230]	1472	!!---------------------------------------------------------------------
	1473	!! * ROUTINE seaoverland *
	1474	!!
	1475	!! ** Purpose : create shifted matrices for seaoverland application
	1476	!! D. Delrosso INGV
	1477	!!----------------------------------------------------------------------
[6140]	1478	INTEGER , INTENT(in ) :: ileni,ilenj ! lengths
	1479	REAL, DIMENSION (ileni,ilenj), INTENT(in ) :: zfieldn ! array of forcing field with undeff for land points
	1480	REAL, DIMENSION (ileni,ilenj), INTENT( out) :: zfield ! array of forcing field
	1481	!
	1482	REAL , DIMENSION (ileni,ilenj) :: zmat1, zmat2, zmat3, zmat4 ! local arrays
	1483	REAL , DIMENSION (ileni,ilenj) :: zmat5, zmat6, zmat7, zmat8 ! - -
	1484	REAL , DIMENSION (ileni,ilenj) :: zlsm2d ! - -
	1485	REAL , DIMENSION (ileni,ilenj,8) :: zlsm3d ! - -
	1486	LOGICAL, DIMENSION (ileni,ilenj,8) :: ll_msknan3d ! logical mask for undeff detection
	1487	LOGICAL, DIMENSION (ileni,ilenj) :: ll_msknan2d ! logical mask for undeff detection
[4230]	1488	!!----------------------------------------------------------------------
[6140]	1489	zmat8 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(:,1)/) , DIM=2 )
	1490	zmat1 = eoshift( zmat8 , SHIFT=-1 , BOUNDARY = (/zmat8(1,:)/) , DIM=1 )
	1491	zmat2 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(1,:)/) , DIM=1 )
	1492	zmat4 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(:,ilenj)/) , DIM=2 )
	1493	zmat3 = eoshift( zmat4 , SHIFT=-1 , BOUNDARY = (/zmat4(1,:)/) , DIM=1 )
	1494	zmat5 = eoshift( zmat4 , SHIFT= 1 , BOUNDARY = (/zmat4(ileni,:)/) , DIM=1 )
	1495	zmat6 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(ileni,:)/) , DIM=1 )
	1496	zmat7 = eoshift( zmat8 , SHIFT= 1 , BOUNDARY = (/zmat8(ileni,:)/) , DIM=1 )
	1497	!
[4230]	1498	zlsm3d = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /))
[6140]	1499	ll_msknan3d = .NOT.( zlsm3d == undeff_lsm )
	1500	ll_msknan2d = .NOT.( zfieldn == undeff_lsm ) ! FALSE where is Undeff (land)
	1501	zlsm2d = SUM( zlsm3d, 3 , ll_msknan3d ) / MAX( 1 , COUNT( ll_msknan3d , 3 ) )
	1502	WHERE( COUNT( ll_msknan3d , 3 ) == 0._wp ) zlsm2d = undeff_lsm
	1503	zfield = MERGE( zfieldn, zlsm2d, ll_msknan2d )
	1504	!
[4230]	1505	END SUBROUTINE seaoverland
	1506
	1507
	1508	SUBROUTINE fld_interp( num, clvar, kw, kk, dta, &
	1509	& nrec, lsmfile)
	1510	!!---------------------------------------------------------------------
[1275]	1511	!! * ROUTINE fld_interp *
	1512	!!
	1513	!! ** Purpose : apply weights to input gridded data to create data
	1514	!! on model grid
	1515	!!----------------------------------------------------------------------
[2715]	1516	INTEGER , INTENT(in ) :: num ! stream number
	1517	CHARACTER(LEN=*) , INTENT(in ) :: clvar ! variable name
	1518	INTEGER , INTENT(in ) :: kw ! weights number
	1519	INTEGER , INTENT(in ) :: kk ! vertical dimension of kk
	1520	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: dta ! output field on model grid
	1521	INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice)
[4230]	1522	CHARACTER(LEN=*) , INTENT(in ) :: lsmfile ! land sea mask file name
[6140]	1523	!
	1524	INTEGER, DIMENSION(3) :: rec1, recn ! temporary arrays for start and length
	1525	INTEGER, DIMENSION(3) :: rec1_lsm, recn_lsm ! temporary arrays for start and length in case of seaoverland
	1526	INTEGER :: ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2 ! temporary indices
	1527	INTEGER :: jk, jn, jm, jir, jjr ! loop counters
	1528	INTEGER :: ni, nj ! lengths
	1529	INTEGER :: jpimin,jpiwid ! temporary indices
	1530	INTEGER :: jpimin_lsm,jpiwid_lsm ! temporary indices
	1531	INTEGER :: jpjmin,jpjwid ! temporary indices
	1532	INTEGER :: jpjmin_lsm,jpjwid_lsm ! temporary indices
	1533	INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices
	1534	INTEGER :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices
	1535	INTEGER :: itmpi,itmpj,itmpz ! lengths
	1536	REAL(wp),DIMENSION(:,:,:), ALLOCATABLE :: ztmp_fly_dta ! local array of values on input grid
[1275]	1537	!!----------------------------------------------------------------------
	1538	!
	1539	!! for weighted interpolation we have weights at four corners of a box surrounding
	1540	!! a model grid point, each weight is multiplied by a grid value (bilinear case)
	1541	!! or by a grid value and gradients at the corner point (bicubic case)
	1542	!! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
	1543
[2528]	1544	!! sub grid from non-model input grid which encloses all grid points in this nemo process
[1275]	1545	jpimin = ref_wgts(kw)%botleft(1)
	1546	jpjmin = ref_wgts(kw)%botleft(2)
	1547	jpiwid = ref_wgts(kw)%jpiwgt
	1548	jpjwid = ref_wgts(kw)%jpjwgt
	1549
[2528]	1550	!! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients
[1275]	1551	rec1(1) = MAX( jpimin-1, 1 )
	1552	rec1(2) = MAX( jpjmin-1, 1 )
[2528]	1553	rec1(3) = 1
[1275]	1554	recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
	1555	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
[2528]	1556	recn(3) = kk
[1275]	1557
[2528]	1558	!! where we need to put it in the non-nemo grid fly_dta
	1559	!! note that jpi1 and jpj1 only differ from 1 when jpimin and jpjmin are 1
	1560	!! (ie at the extreme west or south of the whole input grid) and similarly for jpi2 and jpj2
[1275]	1561	jpi1 = 2 + rec1(1) - jpimin
	1562	jpj1 = 2 + rec1(2) - jpjmin
	1563	jpi2 = jpi1 + recn(1) - 1
	1564	jpj2 = jpj1 + recn(2) - 1
	1565
	1566
[4230]	1567	IF( LEN( TRIM(lsmfile) ) > 0 ) THEN
	1568	!! indeces for ztmp_fly_dta
	1569	! --------------------------
	1570	rec1_lsm(1)=MAX(rec1(1)-nn_lsm,1) ! starting index for enlarged external data, x direction
	1571	rec1_lsm(2)=MAX(rec1(2)-nn_lsm,1) ! starting index for enlarged external data, y direction
	1572	rec1_lsm(3) = 1 ! vertical dimension
	1573	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
	1574	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
	1575	recn_lsm(3) = kk ! number of vertical levels in the input file
	1576
	1577	! Avoid out of bound
	1578	jpimin_lsm = MAX( rec1_lsm(1)+1, 1 )
	1579	jpjmin_lsm = MAX( rec1_lsm(2)+1, 1 )
	1580	jpiwid_lsm = MIN( recn_lsm(1)-2,ref_wgts(kw)%ddims(1)-rec1(1)+1)
	1581	jpjwid_lsm = MIN( recn_lsm(2)-2,ref_wgts(kw)%ddims(2)-rec1(2)+1)
	1582
	1583	jpi1_lsm = 2+rec1_lsm(1)-jpimin_lsm
	1584	jpj1_lsm = 2+rec1_lsm(2)-jpjmin_lsm
	1585	jpi2_lsm = jpi1_lsm + recn_lsm(1) - 1
	1586	jpj2_lsm = jpj1_lsm + recn_lsm(2) - 1
	1587
	1588
[6140]	1589	itmpi=jpi2_lsm-jpi1_lsm+1
	1590	itmpj=jpj2_lsm-jpj1_lsm+1
[4230]	1591	itmpz=kk
	1592	ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz))
	1593	ztmp_fly_dta(:,:,:) = 0.0
	1594	SELECT CASE( SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
	1595	CASE(1)
	1596	CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), &
	1597	& nrec, rec1_lsm, recn_lsm)
	1598	CASE DEFAULT
	1599	CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), &
	1600	& nrec, rec1_lsm, recn_lsm)
	1601	END SELECT
	1602	CALL apply_seaoverland(lsmfile,ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), &
	1603	& jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm, &
	1604	& itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
	1605
	1606
	1607	! Relative indeces for remapping
	1608	ii_lsm1 = (rec1(1)-rec1_lsm(1))+1
	1609	ii_lsm2 = (ii_lsm1+recn(1))-1
	1610	ij_lsm1 = (rec1(2)-rec1_lsm(2))+1
	1611	ij_lsm2 = (ij_lsm1+recn(2))-1
	1612
	1613	ref_wgts(kw)%fly_dta(:,:,:) = 0.0
	1614	ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:) = ztmp_fly_dta(ii_lsm1:ii_lsm2,ij_lsm1:ij_lsm2,:)
	1615	DEALLOCATE(ztmp_fly_dta)
	1616
	1617	ELSE
	1618
	1619	ref_wgts(kw)%fly_dta(:,:,:) = 0.0
	1620	SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) )
	1621	CASE(1)
	1622	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn)
	1623	CASE DEFAULT
	1624	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn)
	1625	END SELECT
	1626	ENDIF
	1627
	1628
[1275]	1629	!! first four weights common to both bilinear and bicubic
[2528]	1630	!! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft
[1275]	1631	!! note that we have to offset by 1 into fly_dta array because of halo
[2528]	1632	dta(:,:,:) = 0.0
[1275]	1633	DO jk = 1,4
	1634	DO jn = 1, jpj
	1635	DO jm = 1,jpi
	1636	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
	1637	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
[2528]	1638	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,:)
[1275]	1639	END DO
	1640	END DO
	1641	END DO
	1642
	1643	IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
	1644
	1645	!! fix up halo points that we couldnt read from file
	1646	IF( jpi1 == 2 ) THEN
[2528]	1647	ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:)
[1275]	1648	ENDIF
	1649	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
[2528]	1650	ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:)
[1275]	1651	ENDIF
	1652	IF( jpj1 == 2 ) THEN
[2528]	1653	ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:)
[1275]	1654	ENDIF
	1655	IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
[2528]	1656	ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:)
[1275]	1657	ENDIF
	1658
	1659	!! if data grid is cyclic we can do better on east-west edges
	1660	!! but have to allow for whether first and last columns are coincident
	1661	IF( ref_wgts(kw)%cyclic ) THEN
	1662	rec1(2) = MAX( jpjmin-1, 1 )
[2528]	1663	recn(1) = 1
[1275]	1664	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
	1665	jpj1 = 2 + rec1(2) - jpjmin
	1666	jpj2 = jpj1 + recn(2) - 1
	1667	IF( jpi1 == 2 ) THEN
[2528]	1668	rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap
	1669	SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
	1670	CASE(1)
	1671	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
	1672	CASE DEFAULT
	1673	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
	1674	END SELECT
	1675	ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
[1275]	1676	ENDIF
	1677	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
[2528]	1678	rec1(1) = 1 + ref_wgts(kw)%overlap
	1679	SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
	1680	CASE(1)
	1681	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
	1682	CASE DEFAULT
	1683	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
	1684	END SELECT
	1685	ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
[1275]	1686	ENDIF
	1687	ENDIF
	1688
	1689	! gradient in the i direction
	1690	DO jk = 1,4
	1691	DO jn = 1, jpj
	1692	DO jm = 1,jpi
	1693	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
	1694	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
[2528]	1695	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 * &
	1696	(ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:))
[1275]	1697	END DO
	1698	END DO
	1699	END DO
	1700
	1701	! gradient in the j direction
	1702	DO jk = 1,4
	1703	DO jn = 1, jpj
	1704	DO jm = 1,jpi
	1705	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
	1706	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
[2528]	1707	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 * &
	1708	(ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:))
[1275]	1709	END DO
	1710	END DO
	1711	END DO
	1712
[2715]	1713	! gradient in the ij direction
	1714	DO jk = 1,4
	1715	DO jn = 1, jpj
	1716	DO jm = 1,jpi
	1717	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
	1718	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
	1719	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
[2528]	1720	(ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni ,nj+2,:)) - &
	1721	(ref_wgts(kw)%fly_dta(ni+2,nj ,:) - ref_wgts(kw)%fly_dta(ni ,nj ,:)))
[2715]	1722	END DO
[1275]	1723	END DO
[2715]	1724	END DO
	1725	!
[1275]	1726	END IF
[2715]	1727	!
[1275]	1728	END SUBROUTINE fld_interp
[2528]	1729
	1730
	1731	FUNCTION ksec_week( cdday )
	1732	!!---------------------------------------------------------------------
	1733	!! * FUNCTION kshift_week *
	1734	!!
[7646]	1735	!! ** Purpose : return the first 3 letters of the first day of the weekly file
[2528]	1736	!!---------------------------------------------------------------------
[7646]	1737	CHARACTER(len=*), INTENT(in) :: cdday ! first 3 letters of the first day of the weekly file
[2528]	1738	!!
[7646]	1739	INTEGER :: ksec_week ! output variable
	1740	INTEGER :: ijul, ishift ! local integer
[2528]	1741	CHARACTER(len=3),DIMENSION(7) :: cl_week
	1742	!!----------------------------------------------------------------------
	1743	cl_week = (/"sun","sat","fri","thu","wed","tue","mon"/)
	1744	DO ijul = 1, 7
	1745	IF( cl_week(ijul) == TRIM(cdday) ) EXIT
[2715]	1746	END DO
[2528]	1747	IF( ijul .GT. 7 ) CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): '//TRIM(cdday) )
	1748	!
	1749	ishift = ijul * NINT(rday)
	1750	!
	1751	ksec_week = nsec_week + ishift
	1752	ksec_week = MOD( ksec_week, 7*NINT(rday) )
	1753	!
	1754	END FUNCTION ksec_week
	1755
[2715]	1756	!!======================================================================
[888]	1757	END MODULE fldread

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

Download in other formats: