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

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

fldread.F90 in branches/2013/dev_CMCC_INGV_2013/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2013/dev_CMCC_INGV_2013/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 4217

Last change on this file since 4217 was 4217, checked in by poddo, 11 years ago
Solved problems with namelist parameter land/sea mask
Property svn:keywords set to `Id`
File size: 81.8 KB

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

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

Download in other formats: