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

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

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

Context Navigation

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

Last change on this file since 14205 was 13546, checked in by smasson, 4 years ago
trunk: supress the use of undefined values, see #2535
Property svn:keywords set to `Id`
File size: 88.3 KB

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

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

Download in other formats: