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/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/fldread.F90 @ 12779

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

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

Download in other formats: