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 tags/nemo_v3_2/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: tags/nemo_v3_2/NEMO/OPA_SRC/SBC/fldread.F90 @ 1817

Last change on this file since 1817 was 1817, checked in by rblod, 14 years ago
Correct pathologic case in fld_init, see ticket #653
Property svn:keywords set to `Id`
File size: 51.8 KB

Rev	Line
[888]	1	MODULE fldread
	2	!!======================================================================
	3	!! * MODULE fldread *
	4	!! Ocean forcing: read input field for surface boundary condition
	5	!!=====================================================================
	6	!! History : 9.0 ! 06-06 (G. Madec) Original code
[1275]	7	!! ! 05-08 (S. Alderson) Modified for Interpolation in memory
	8	!! ! from input grid to model grid
[888]	9	!!----------------------------------------------------------------------
	10
	11	!!----------------------------------------------------------------------
	12	!! fld_read : read input fields used for the computation of the
	13	!! surface boundary condition
	14	!!----------------------------------------------------------------------
	15	USE oce ! ocean dynamics and tracers
	16	USE dom_oce ! ocean space and time domain
	17	USE phycst ! ???
	18	USE in_out_manager ! I/O manager
	19	USE iom ! I/O manager library
[1275]	20	USE geo2ocean ! for vector rotation on to model grid
[888]	21
	22	IMPLICIT NONE
	23	PRIVATE
	24
	25	TYPE, PUBLIC :: FLD_N !: Namelist field informations
[1730]	26	CHARACTER(len = 256) :: clname ! generic name of the NetCDF flux file
	27	INTEGER :: nfreqh ! frequency of each flux file
	28	CHARACTER(len = 34) :: clvar ! generic name of the variable in the NetCDF flux file
	29	LOGICAL :: ln_tint ! time interpolation or not (T/F)
	30	LOGICAL :: ln_clim ! climatology or not (T/F)
	31	CHARACTER(len = 7) :: cltype ! type of data file 'daily', 'monthly' or yearly'
	32	CHARACTER(len = 34) :: wname ! generic name of a NetCDF weights file to be used, blank if not
	33	CHARACTER(len = 34) :: vcomp ! symbolic component name if a vector that needs rotation
[1275]	34	! a string starting with "U" or "V" for each component
	35	! chars 2 onwards identify which components go together
[888]	36	END TYPE FLD_N
	37
	38	TYPE, PUBLIC :: FLD !: Input field related variables
	39	CHARACTER(len = 256) :: clrootname ! generic name of the NetCDF file
	40	CHARACTER(len = 256) :: clname ! current name of the NetCDF file
[1730]	41	INTEGER :: nfreqh ! frequency of each flux file
[888]	42	CHARACTER(len = 34) :: clvar ! generic name of the variable in the NetCDF flux file
	43	LOGICAL :: ln_tint ! time interpolation or not (T/F)
[1132]	44	LOGICAL :: ln_clim ! climatology or not (T/F)
[1628]	45	CHARACTER(len = 7) :: cltype ! type of data file 'daily', 'monthly' or yearly'
[1132]	46	INTEGER :: num ! iom id of the jpfld files to be read
[1730]	47	INTEGER :: nswap_sec ! swapping time in second since Jan. 1st 00h of nit000 year
	48	INTEGER , DIMENSION(2) :: nrec_b ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year)
	49	INTEGER , DIMENSION(2) :: nrec_a ! after record (1: index, 2: second since Jan. 1st 00h of nit000 year)
[1200]	50	REAL(wp) , ALLOCATABLE, DIMENSION(:,:) :: fnow ! input fields interpolated to now time step
	51	REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:) :: fdta ! 2 consecutive record of input fields
[1275]	52	CHARACTER(len = 256) :: wgtname ! current name of the NetCDF weight file acting as a key
	53	! into the WGTLIST structure
	54	CHARACTER(len = 34) :: vcomp ! symbolic name for a vector component that needs rotation
	55	LOGICAL , DIMENSION(2) :: rotn ! flag to indicate whether field has been rotated
[888]	56	END TYPE FLD
	57
[1275]	58	!$AGRIF_DO_NOT_TREAT
	59
	60	!! keep list of all weights variables so they're only read in once
	61	!! need to add AGRIF directives not to process this structure
	62	!! also need to force wgtname to include AGRIF nest number
	63	TYPE :: WGT !: Input weights related variables
	64	CHARACTER(len = 256) :: wgtname ! current name of the NetCDF weight file
	65	INTEGER , DIMENSION(2) :: ddims ! shape of input grid
	66	INTEGER , DIMENSION(2) :: botleft ! top left corner of box in input grid containing
	67	! current processor grid
	68	INTEGER , DIMENSION(2) :: topright ! top right corner of box
	69	INTEGER :: jpiwgt ! width of box on input grid
	70	INTEGER :: jpjwgt ! height of box on input grid
	71	INTEGER :: numwgt ! number of weights (4=bilinear, 16=bicubic)
	72	INTEGER :: nestid ! for agrif, keep track of nest we're in
	73	INTEGER :: offset ! =0 when cyclic grid has coincident first/last columns,
	74	! =1 when they assumed to be one grid spacing apart
	75	! =-1 otherwise
	76	LOGICAL :: cyclic ! east-west cyclic or not
	77	INTEGER, DIMENSION(:,:,:), POINTER :: data_jpi ! array of source integers
	78	INTEGER, DIMENSION(:,:,:), POINTER :: data_jpj ! array of source integers
	79	REAL(wp), DIMENSION(:,:,:), POINTER :: data_wgt ! array of weights on model grid
	80	REAL(wp), DIMENSION(:,:), POINTER :: fly_dta ! array of values on input grid
	81	REAL(wp), DIMENSION(:,:), POINTER :: col2 ! temporary array for reading in columns
	82	END TYPE WGT
	83
	84	INTEGER, PARAMETER :: tot_wgts = 10
	85	TYPE( WGT ), DIMENSION(tot_wgts) :: ref_wgts ! array of wgts
	86	INTEGER :: nxt_wgt = 1 ! point to next available space in ref_wgts array
	87
	88	!$AGRIF_END_DO_NOT_TREAT
	89
[1132]	90	PUBLIC fld_read, fld_fill ! called by sbc... modules
[888]	91
	92	!!----------------------------------------------------------------------
	93	!! OPA 9.0 , LOCEAN-IPSL (2006)
[1156]	94	!! $Id$
[888]	95	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
	96	!!----------------------------------------------------------------------
	97
	98	CONTAINS
	99
	100	SUBROUTINE fld_read( kt, kn_fsbc, sd )
	101	!!---------------------------------------------------------------------
	102	!! * ROUTINE fld_read *
	103	!!
	104	!! ** Purpose : provide at each time step the surface ocean fluxes
	105	!! (momentum, heat, freshwater and runoff)
	106	!!
	107	!! ** Method : READ each input fields in NetCDF files using IOM
	108	!! and intepolate it to the model time-step.
	109	!! Several assumptions are made on the input file:
	110	!! blahblahblah....
	111	!!----------------------------------------------------------------------
	112	INTEGER , INTENT(in ) :: kt ! ocean time step
[1132]	113	INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step)
[888]	114	TYPE(FLD), INTENT(inout), DIMENSION(:) :: sd ! input field related variables
	115	!!
[1275]	116	CHARACTER (LEN=34) :: acomp ! dummy weight name
	117	INTEGER :: kf, nf ! dummy indices
	118	INTEGER :: imf ! size of the structure sd
	119	REAL(wp), DIMENSION(jpi,jpj) :: utmp, vtmp! temporary arrays for vector rotation
	120
[1132]	121	INTEGER :: jf ! dummy indices
[1275]	122	INTEGER :: kw ! index into wgts array
[1730]	123	INTEGER :: ireclast ! last record to be read in the current year file
	124	INTEGER :: isecend ! number of second since Jan. 1st 00h of nit000 year at nitend
[1628]	125	LOGICAL :: llnxtyr ! open next year file?
	126	LOGICAL :: llnxtmth ! open next month file?
	127	LOGICAL :: llstop ! stop is the file does not exist
[1132]	128	REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation
	129	REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation
[1191]	130	CHARACTER(LEN=1000) :: clfmt ! write format
[888]	131	!!---------------------------------------------------------------------
[1275]	132	!
	133	imf = SIZE( sd )
[888]	134	! ! ===================== !
[1275]	135	DO jf = 1, imf ! LOOP OVER FIELD !
[888]	136	! ! ===================== !
	137	!
[1132]	138	IF( kt == nit000 ) CALL fld_init( sd(jf) )
	139	!
	140	! read/update the after data?
[1730]	141	IF( nsec_year + nsec1jan000 > sd(jf)%nswap_sec ) THEN
[888]	142
[1132]	143	IF( sd(jf)%ln_tint ) THEN ! time interpolation: swap before record field
[888]	144	!CDIR COLLAPSE
[1132]	145	sd(jf)%fdta(:,:,1) = sd(jf)%fdta(:,:,2)
[1275]	146	sd(jf)%rotn(1) = sd(jf)%rotn(2)
[888]	147	ENDIF
[1132]	148
	149	! update record informations
	150	CALL fld_rec( sd(jf) )
	151
[1628]	152	! do we have to change the year/month/day of the forcing field??
[888]	153	IF( sd(jf)%ln_tint ) THEN
[1628]	154	! if we do time interpolation we will need to open next year/month/day file before the end of the current one
[1730]	155	! if so, we are still before the end of the year/month/day when calling fld_rec so sd(jf)%nrec_a(1) will be
[1628]	156	! larger than the record number that should be read for current year/month/day (for ex. 13 for monthly mean file)
[1132]	157
	158	! last record to be read in the current file
[1730]	159	IF( sd(jf)%nfreqh == -1 ) THEN ; ireclast = 12
[1132]	160	ELSE
[1730]	161	IF( sd(jf)%cltype == 'monthly' ) THEN ; ireclast = 24 * nmonth_len(nmonth) / sd(jf)%nfreqh
	162	ELSEIF( sd(jf)%cltype == 'daily' ) THEN ; ireclast = 24 / sd(jf)%nfreqh
	163	ELSE ; ireclast = 24 * nyear_len( 1 ) / sd(jf)%nfreqh
[1132]	164	ENDIF
	165	ENDIF
	166
[1628]	167	! do we need next file data?
[1730]	168	IF( sd(jf)%nrec_a(1) > ireclast ) THEN
[1132]	169
[1730]	170	sd(jf)%nrec_a(1) = 1 ! force to read the first record of the next file
[1132]	171
	172	IF( .NOT. sd(jf)%ln_clim ) THEN ! close the current file and open a new one.
[1628]	173
	174	llnxtmth = sd(jf)%cltype == 'monthly' .OR. nday == nmonth_len(nmonth) ! open next month file?
	175	llnxtyr = sd(jf)%cltype == 'yearly' .OR. (nmonth == 12 .AND. llnxtmth) ! open next year file?
	176
	177	! if the run finishes at the end of the current year/month/day, we will allow next year/month/day file to be
	178	! not present. If the run continue further than the current year/month/day, next year/month/day file must exist
[1730]	179	isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdttra(1)) ! second at the end of the run
	180	llstop = isecend > sd(jf)%nswap_sec ! read more than 1 record of next year
[1132]	181
[1628]	182	CALL fld_clopn( sd(jf), nyear + COUNT((/llnxtyr /)) , &
	183	& nmonth + COUNT((/llnxtmth/)) - 12 * COUNT((/llnxtyr /)), &
	184	& nday + 1 - nmonth_len(nmonth) * COUNT((/llnxtmth/)), llstop )
[1132]	185
[1628]	186	IF( sd(jf)%num == 0 .AND. .NOT. llstop ) THEN ! next year file does not exist
	187	CALL ctl_warn('next year/month/day file: '//TRIM(sd(jf)%clname)// &
	188	& ' not present -> back to current year/month/day')
	189	CALL fld_clopn( sd(jf), nyear, nmonth, nday ) ! back to the current year/month/day
[1730]	190	sd(jf)%nrec_a(1) = ireclast ! force to read the last record to be read in the current year file
[1132]	191	ENDIF
	192
	193	ENDIF
	194	ENDIF
	195
[888]	196	ELSE
[1628]	197	! if we are not doing time interpolation, we must change the year/month/day of the file just after switching
	198	! to the NEW year/month/day. If it is the case, we are at the beginning of the year/month/day when calling
[1730]	199	! fld_rec so sd(jf)%nrec_a(1) = 1
	200	IF( sd(jf)%nrec_a(1) == 1 .AND. .NOT. sd(jf)%ln_clim ) CALL fld_clopn( sd(jf), nyear, nmonth, nday )
[888]	201	ENDIF
[1132]	202
	203	! read after data
[1275]	204	IF( LEN(TRIM(sd(jf)%wgtname)) > 0 ) THEN
	205	CALL wgt_list( sd(jf), kw )
[1730]	206	CALL fld_interp( sd(jf)%num, sd(jf)%clvar, kw, sd(jf)%fdta(:,:,2), sd(jf)%nrec_a(1) )
[1275]	207	ELSE
[1730]	208	CALL iom_get( sd(jf)%num, jpdom_data, sd(jf)%clvar, sd(jf)%fdta(:,:,2), sd(jf)%nrec_a(1) )
[1275]	209	ENDIF
	210	sd(jf)%rotn(2) = .FALSE.
[1132]	211
[888]	212	ENDIF
[1275]	213	! ! ===================== !
	214	END DO ! END LOOP OVER FIELD !
	215	! ! ===================== !
[888]	216
[1275]	217	IF( kt == nit000 .AND. lwp ) CALL wgt_print()
	218
	219	!! Vector fields may need to be rotated onto the local grid direction
	220	!! This has to happen before the time interpolations
	221	!! (sga: following code should be modified so that pairs arent searched for each time
	222
	223	DO jf = 1, imf
	224	!! find vector rotations required
	225	IF( LEN(TRIM(sd(jf)%vcomp)) > 0 ) THEN
	226	!! east-west component has symbolic name starting with 'U'
	227	IF( sd(jf)%vcomp(1:1) == 'U' ) THEN
	228	!! found an east-west component, look for the north-south component
	229	!! which has same symbolic name but with 'U' replaced with 'V'
	230	nf = LEN_TRIM( sd(jf)%vcomp )
	231	IF( nf == 1) THEN
	232	acomp = 'V'
	233	ELSE
	234	acomp = 'V' // sd(jf)%vcomp(2:nf)
	235	ENDIF
	236	kf = -1
	237	DO nf = 1, imf
	238	IF( TRIM(sd(nf)%vcomp) == TRIM(acomp) ) kf = nf
	239	END DO
	240	IF( kf > 0 ) THEN
	241	!! fields jf,kf are two components which need to be rotated together
	242	DO nf = 1,2
	243	!! check each time level of this pair
	244	IF( .NOT. sd(jf)%rotn(nf) .AND. .NOT. sd(kf)%rotn(nf) ) THEN
	245	utmp(:,:) = 0.0
	246	vtmp(:,:) = 0.0
	247	CALL rot_rep( sd(jf)%fdta(:,:,nf), sd(kf)%fdta(:,:,nf), 'T', 'en->i', utmp(:,:) )
	248	CALL rot_rep( sd(jf)%fdta(:,:,nf), sd(kf)%fdta(:,:,nf), 'T', 'en->j', vtmp(:,:) )
	249	sd(jf)%fdta(:,:,nf) = utmp(:,:)
	250	sd(kf)%fdta(:,:,nf) = vtmp(:,:)
	251	sd(jf)%rotn(nf) = .TRUE.
	252	sd(kf)%rotn(nf) = .TRUE.
	253	IF( lwp .AND. kt == nit000 ) &
	254	WRITE(numout,*) 'fld_read: vector pair (', &
	255	TRIM(sd(jf)%clvar),',',TRIM(sd(kf)%clvar), &
	256	') rotated on to model grid'
	257	ENDIF
	258	END DO
	259	ENDIF
	260	ENDIF
	261	ENDIF
	262	END DO
	263
	264	! ! ===================== !
	265	DO jf = 1, imf ! LOOP OVER FIELD !
	266	! ! ===================== !
	267	!
[1132]	268	! update field at each kn_fsbc time-step
	269	IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
[888]	270	!
[1132]	271	IF( sd(jf)%ln_tint ) THEN
[1191]	272	IF(lwp .AND. kt - nit000 <= 100 ) THEN
	273	clfmt = "('fld_read: var ', a, ' kt = ', i8,' Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," // &
	274	& "' records b/a: ', i4.4, '/', i4.4, ' (', f7.2,'/', f7.2, ' days)')"
	275	WRITE(numout, clfmt) TRIM( sd(jf)%clvar ), kt, nyear, nmonth, nday, &
[1730]	276	& sd(jf)%nrec_b(1), sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday
[1191]	277	ENDIF
[1132]	278	!
[1730]	279	ztinta = REAL( nsec_year + nsec1jan000 - sd(jf)%nrec_b(2), wp ) / REAL( sd(jf)%nrec_a(2) - sd(jf)%nrec_b(2), wp )
[1132]	280	ztintb = 1. - ztinta
[888]	281	!CDIR COLLAPSE
[1132]	282	sd(jf)%fnow(:,:) = ztintb * sd(jf)%fdta(:,:,1) + ztinta * sd(jf)%fdta(:,:,2)
[888]	283	ELSE
[1191]	284	IF(lwp .AND. kt - nit000 <= 100 ) THEN
	285	clfmt = "('fld_read: var ', a, ' kt = ', i8,' Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," // &
	286	& "' record: ', i4.4, ' at ', f7.2, ' day')"
[1730]	287	WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, nyear, nmonth, nday, sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_a(2),wp)/rday
[1191]	288	ENDIF
[888]	289	!CDIR COLLAPSE
[1132]	290	sd(jf)%fnow(:,:) = sd(jf)%fdta(:,:,2) ! piecewise constant field
	291
[888]	292	ENDIF
	293	!
	294	ENDIF
[1132]	295
	296	IF( kt == nitend ) CALL iom_close( sd(jf)%num ) ! Close the input files
	297
[888]	298	! ! ===================== !
	299	END DO ! END LOOP OVER FIELD !
	300	! ! ===================== !
	301	END SUBROUTINE fld_read
	302
	303
[1132]	304	SUBROUTINE fld_init( sdjf )
[888]	305	!!---------------------------------------------------------------------
[1132]	306	!! * ROUTINE fld_init *
	307	!!
	308	!! ** Purpose : - if time interpolation, read before data
	309	!! - open current year file
	310	!!
	311	!! ** Method :
	312	!!----------------------------------------------------------------------
	313	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
	314	!!
	315	LOGICAL :: llprevyr ! are we reading previous year file?
	316	LOGICAL :: llprevmth ! are we reading previous month file?
[1628]	317	LOGICAL :: llprevday ! are we reading previous day file?
	318	LOGICAL :: llprev ! llprevyr .OR. llprevmth .OR. llprevday
[1132]	319	INTEGER :: idvar ! variable id
	320	INTEGER :: inrec ! number of record existing for this variable
[1275]	321	INTEGER :: kwgt
[1191]	322	CHARACTER(LEN=1000) :: clfmt ! write format
[1132]	323	!!---------------------------------------------------------------------
	324
	325	! some default definitions...
	326	sdjf%num = 0 ! default definition for non-opened file
	327	IF( sdjf%ln_clim ) sdjf%clname = TRIM( sdjf%clrootname ) ! file name defaut definition, never change in this case
	328	llprevyr = .FALSE.
	329	llprevmth = .FALSE.
[1628]	330	llprevday = .FALSE.
[1132]	331
	332	! define record informations
	333	CALL fld_rec( sdjf )
	334
	335	IF( sdjf%ln_tint ) THEN ! we need to read the previous record and we will put it in the current record structure
	336
[1730]	337	IF( sdjf%nrec_b(1) == 0 ) THEN ! we redefine record sdjf%nrec_b(1) with the last record of previous year file
	338	IF( sdjf%nfreqh == -1 ) THEN ! monthly mean
[1628]	339	IF( sdjf%cltype == 'monthly' ) THEN ! monthly file
[1730]	340	sdjf%nrec_b(1) = 1 ! force to read the unique record
[1628]	341	llprevmth = .NOT. sdjf%ln_clim ! use previous month file?
	342	llprevyr = llprevmth .AND. nmonth == 1 ! use previous year file?
	343	ELSE ! yearly file
[1730]	344	sdjf%nrec_b(1) = 12 ! force to read december mean
[1628]	345	llprevyr = .NOT. sdjf%ln_clim ! use previous year file?
	346	ENDIF
[1132]	347	ELSE
	348	IF( sdjf%cltype == 'monthly' ) THEN ! monthly file
[1730]	349	sdjf%nrec_b(1) = 24 * nmonth_len(nmonth-1) / sdjf%nfreqh ! last record of previous month
[1628]	350	llprevmth = .NOT. sdjf%ln_clim ! use previous month file?
	351	llprevyr = llprevmth .AND. nmonth == 1 ! use previous year file?
	352	ELSEIF( sdjf%cltype == 'daily' ) THEN ! daily file
[1730]	353	sdjf%nrec_b(1) = 24 / sdjf%nfreqh ! last record of previous day
[1628]	354	llprevday = .NOT. sdjf%ln_clim ! use previous day file?
	355	llprevmth = llprevday .AND. nday == 1 ! use previous month file?
	356	llprevyr = llprevmth .AND. nmonth == 1 ! use previous year file?
[1132]	357	ELSE ! yearly file
[1730]	358	sdjf%nrec_b(1) = 24 * nyear_len(0) / sdjf%nfreqh ! last record of year month
[1628]	359	llprevyr = .NOT. sdjf%ln_clim ! use previous year file?
[1132]	360	ENDIF
	361	ENDIF
	362	ENDIF
[1628]	363	llprev = llprevyr .OR. llprevmth .OR. llprevday
[1132]	364
[1628]	365	CALL fld_clopn( sdjf, nyear - COUNT((/llprevyr /)) , &
	366	& nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /)), &
	367	& nday - COUNT((/llprevday/)) + nmonth_len(nmonth-1) * COUNT((/llprevmth/)), .NOT. llprev )
	368
	369	! if previous year/month/day file does not exist, we switch to the current year/month/day
[1817]	370	IF( llprev .AND. sdjf%num <= 0 ) THEN
[1628]	371	CALL ctl_warn( 'previous year/month/day file: '//TRIM(sdjf%clname)//' not present -> back to current year/month/day')
	372	! we force to read the first record of the current year/month/day instead of last record of previous year/month/day
[1132]	373	llprev = .false.
[1730]	374	sdjf%nrec_b(1) = 1
[1628]	375	CALL fld_clopn( sdjf, nyear, nmonth, nday )
[1132]	376	ENDIF
	377
[1730]	378	IF( llprev ) THEN ! check if the last record sdjf%nrec_n(1) exists in the file
[1132]	379	idvar = iom_varid( sdjf%num, sdjf%clvar ) ! id of the variable sdjf%clvar
	380	IF( idvar <= 0 ) RETURN
	381	inrec = iom_file( sdjf%num )%dimsz( iom_file( sdjf%num )%ndims(idvar), idvar ) ! size of the last dim of idvar
[1730]	382	sdjf%nrec_b(1) = MIN( sdjf%nrec_b(1), inrec ) ! make sure we select an existing record
[1132]	383	ENDIF
	384
	385	! read before data into sdjf%fdta(:,:,2) because we will swap data in the following part of fld_read
[1275]	386	IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
	387	CALL wgt_list( sdjf, kwgt )
[1730]	388	CALL fld_interp( sdjf%num, sdjf%clvar, kwgt, sdjf%fdta(:,:,2), sdjf%nrec_b(1) )
[1275]	389	ELSE
[1730]	390	CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,2), sdjf%nrec_b(1) )
[1275]	391	ENDIF
	392	sdjf%rotn(2) = .FALSE.
[1132]	393
[1191]	394	clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i4, ' at time = ', f7.2, ' days')"
[1730]	395	IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_b(1), REAL(sdjf%nrec_b(2),wp)/rday
[1132]	396
	397	IF( llprev ) CALL iom_close( sdjf%num ) ! close previous year file (-> redefine sdjf%num to 0)
	398
	399	ENDIF
	400
[1628]	401	IF( sdjf%num == 0 ) CALL fld_clopn( sdjf, nyear, nmonth, nday ) ! make sure current year/month/day file is opened
[1132]	402
[1730]	403	sdjf%nswap_sec = nsec_year + nsec1jan000 - 1 ! force read/update the after data in the following part of fld_read
[1132]	404
	405	END SUBROUTINE fld_init
	406
	407
	408	SUBROUTINE fld_rec( sdjf )
	409	!!---------------------------------------------------------------------
[888]	410	!! * ROUTINE fld_rec *
	411	!!
[1730]	412	!! ** Purpose : compute nrec_a, nrec_b and nswap_sec
[888]	413	!!
	414	!! ** Method :
	415	!!----------------------------------------------------------------------
[1132]	416	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
[888]	417	!!
[1132]	418	INTEGER :: irec ! record number
[1730]	419	INTEGER :: isecd ! rday
[1132]	420	REAL(wp) :: ztmp ! temporary variable
[1730]	421	INTEGER :: ifreq_sec ! frequency mean (in seconds)
[888]	422	!!----------------------------------------------------------------------
	423	!
[1730]	424	IF( sdjf%nfreqh == -1 ) THEN ! monthly mean
[888]	425	!
[1132]	426	IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record
	427	!
	428	! INT( ztmp )
	429	! /\|\
	430	! 1 \| *----
	431	! 0 \|----(
	432	! \|----+----\|--> time
	433	! 0 /\|\ 1 (nday/nmonth_len(nmonth))
	434	! \|
	435	! \|
	436	! forcing record : nmonth
	437	!
	438	ztmp = REAL( nday, wp ) / REAL( nmonth_len(nmonth), wp ) + 0.5
[888]	439	ELSE
[1132]	440	ztmp = 0.e0
[888]	441	ENDIF
[1132]	442	irec = nmonth + INT( ztmp )
	443
[1730]	444	IF( sdjf%ln_tint ) THEN ; sdjf%nswap_sec = nmonth_half(irec) + nsec1jan000 ! swap at the middle of the month
	445	ELSE ; sdjf%nswap_sec = nmonth_end (irec) + nsec1jan000 ! swap at the end of the month
[1132]	446	ENDIF
	447
[1730]	448	sdjf%nrec_a(:) = (/ irec, nmonth_half(irec) + nsec1jan000 /) ! define after record number and time
	449	irec = irec - 1 ! move back to previous record
	450	sdjf%nrec_b(:) = (/ irec, nmonth_half(irec) + nsec1jan000 /) ! define before record number and time
[888]	451	!
[1730]	452	ELSE ! higher frequency mean (in hours)
[888]	453	!
[1730]	454	ifreq_sec = sdjf%nfreqh * 3600 ! frequency mean (in seconds)
[1132]	455	! number of second since the beginning of the file
[1730]	456	IF( sdjf%cltype == 'monthly' ) THEN ; ztmp = REAL(nsec_month,wp) ! since 00h on the 1st day of the current month
	457	ELSEIF( sdjf%cltype == 'daily' ) THEN ; ztmp = REAL(nsec_day ,wp) ! since 00h of the current day
	458	ELSE ; ztmp = REAL(nsec_year ,wp) ! since 00h on Jan 1 of the current year
[1132]	459	ENDIF
	460	IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record
	461	!
	462	! INT( ztmp )
	463	! /\|\
	464	! 2 \| *-----(
	465	! 1 \| *-----(
	466	! 0 \|--(
	467	! \|--+--\|--+--\|--+--\|--> time
[1730]	468	! 0 /\|\ 1 /\|\ 2 /\|\ 3 (nsec_year/ifreq_sec) or (nsec_month/ifreq_sec)
[1132]	469	! \| \| \|
	470	! \| \| \|
	471	! forcing record : 1 2 3
	472	!
[1730]	473	ztmp= ztmp / ifreq_sec + 0.5
[888]	474	ELSE
[1132]	475	!
	476	! INT( ztmp )
	477	! /\|\
	478	! 2 \| *-----(
	479	! 1 \| *-----(
	480	! 0 \|-----(
	481	! \|--+--\|--+--\|--+--\|--> time
[1730]	482	! 0 /\|\ 1 /\|\ 2 /\|\ 3 (nsec_year/ifreq_sec) or (nsec_month/ifreq_sec)
[1132]	483	! \| \| \|
	484	! \| \| \|
	485	! forcing record : 1 2 3
	486	!
[1730]	487	ztmp= ztmp / ifreq_sec
[1132]	488	ENDIF
[1730]	489	irec = 1 + INT( ztmp )
[1132]	490
[1730]	491	isecd = NINT(rday)
[1132]	492	! after record index and second since Jan. 1st 00h of nit000 year
[1730]	493	sdjf%nrec_a(:) = (/ irec, ifreq_sec * irec - ifreq_sec / 2 + nsec1jan000 /)
[1628]	494	IF( sdjf%cltype == 'monthly' ) & ! add the number of seconds between 00h Jan 1 and the end of previous month
[1730]	495	sdjf%nrec_a(2) = sdjf%nrec_a(2) + isecd * SUM(nmonth_len(1:nmonth -1)) ! ok if nmonth=1
[1628]	496	IF( sdjf%cltype == 'daily' ) & ! add the number of seconds between 00h Jan 1 and the end of previous day
[1730]	497	sdjf%nrec_a(2) = sdjf%nrec_a(2) + isecd * ( nday_year - 1 )
[1132]	498
	499	! before record index and second since Jan. 1st 00h of nit000 year
[1730]	500	irec = irec - 1. ! move back to previous record
	501	sdjf%nrec_b(:) = (/ irec, ifreq_sec * irec - ifreq_sec / 2 + nsec1jan000 /)
[1628]	502	IF( sdjf%cltype == 'monthly' ) & ! add the number of seconds between 00h Jan 1 and the end of previous month
[1730]	503	sdjf%nrec_b(2) = sdjf%nrec_b(2) + isecd * SUM(nmonth_len(1:nmonth -1)) ! ok if nmonth=1
[1628]	504	IF( sdjf%cltype == 'daily' ) & ! add the number of seconds between 00h Jan 1 and the end of previous day
[1730]	505	sdjf%nrec_b(2) = sdjf%nrec_b(2) + isecd * ( nday_year - 1 )
[1132]	506
	507	! swapping time in second since Jan. 1st 00h of nit000 year
[1730]	508	IF( sdjf%ln_tint ) THEN ; sdjf%nswap_sec = sdjf%nrec_a(2) ! swap at the middle of the record
	509	ELSE ; sdjf%nswap_sec = sdjf%nrec_a(2) + ifreq_sec / 2 ! swap at the end of the record
[1132]	510	ENDIF
[888]	511	!
	512	ENDIF
	513	!
[1132]	514	END SUBROUTINE fld_rec
	515
	516
[1628]	517	SUBROUTINE fld_clopn( sdjf, kyear, kmonth, kday, ldstop )
[1132]	518	!!---------------------------------------------------------------------
	519	!! * ROUTINE fld_clopn *
	520	!!
	521	!! ** Purpose : update the file name and open the file
	522	!!
	523	!! ** Method :
	524	!!----------------------------------------------------------------------
	525	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
	526	INTEGER , INTENT(in ) :: kyear ! year value
	527	INTEGER , INTENT(in ) :: kmonth ! month value
[1628]	528	INTEGER , INTENT(in ) :: kday ! day value
[1132]	529	LOGICAL , INTENT(in ), OPTIONAL :: ldstop ! stop if open to read a non-existing file (default = .TRUE.)
	530
	531	IF( sdjf%num /= 0 ) CALL iom_close( sdjf%num ) ! close file if already open
	532	! build the new filename if not climatological data
[1212]	533	IF( .NOT. sdjf%ln_clim ) THEN ; WRITE(sdjf%clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), kyear ! add year
[1628]	534	IF( sdjf%cltype /= 'yearly' ) WRITE(sdjf%clname, '(a,"m" ,i2.2)' ) TRIM( sdjf%clname ), kmonth ! add month
	535	IF( sdjf%cltype == 'daily' ) WRITE(sdjf%clname, '(a,"d" ,i2.2)' ) TRIM( sdjf%clname ), kday ! add day
[888]	536	ENDIF
[1319]	537	CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 )
[888]	538	!
[1132]	539	END SUBROUTINE fld_clopn
	540
	541
	542	SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )
	543	!!---------------------------------------------------------------------
	544	!! * ROUTINE fld_fill *
	545	!!
	546	!! ** Purpose : fill sdf with sdf_n and control print
	547	!!
	548	!! ** Method :
	549	!!----------------------------------------------------------------------
	550	TYPE(FLD) , DIMENSION(:), INTENT(inout) :: sdf ! structure of input fields (file informations, fields read)
	551	TYPE(FLD_N), DIMENSION(:), INTENT(in ) :: sdf_n ! array of namelist information structures
	552	CHARACTER(len=*) , INTENT(in ) :: cdir ! Root directory for location of flx files
	553	CHARACTER(len=*) , INTENT(in ) :: cdcaller !
	554	CHARACTER(len=*) , INTENT(in ) :: cdtitle !
	555	CHARACTER(len=*) , INTENT(in ) :: cdnam !
[888]	556	!
[1132]	557	INTEGER :: jf ! dummy indices
	558	!!---------------------------------------------------------------------
[888]	559
[1132]	560	DO jf = 1, SIZE(sdf)
	561	sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
[1730]	562	sdf(jf)%nfreqh = sdf_n(jf)%nfreqh
[1132]	563	sdf(jf)%clvar = sdf_n(jf)%clvar
	564	sdf(jf)%ln_tint = sdf_n(jf)%ln_tint
	565	sdf(jf)%ln_clim = sdf_n(jf)%ln_clim
[1730]	566	IF( sdf(jf)%nfreqh == -1. ) THEN ; sdf(jf)%cltype = 'yearly'
	567	ELSE ; sdf(jf)%cltype = sdf_n(jf)%cltype
[1132]	568	ENDIF
[1275]	569	sdf(jf)%wgtname = " "
[1730]	570	IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
[1275]	571	sdf(jf)%vcomp = sdf_n(jf)%vcomp
[1132]	572	END DO
	573
	574	IF(lwp) THEN ! control print
	575	WRITE(numout,*)
	576	WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
	577	WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
	578	WRITE(numout,*) ' '//TRIM( cdnam )//' Namelist'
	579	WRITE(numout,*) ' list of files and frequency (>0: in hours ; <0 in months)'
	580	DO jf = 1, SIZE(sdf)
	581	WRITE(numout,*) ' root filename: ' , TRIM( sdf(jf)%clrootname ), &
	582	& ' variable name: ' , TRIM( sdf(jf)%clvar )
[1730]	583	WRITE(numout,*) ' frequency: ' , sdf(jf)%nfreqh , &
[1132]	584	& ' time interp: ' , sdf(jf)%ln_tint , &
	585	& ' climatology: ' , sdf(jf)%ln_clim , &
[1275]	586	& ' weights : ' , TRIM( sdf(jf)%wgtname ), &
	587	& ' pairing : ' , TRIM( sdf(jf)%vcomp ), &
[1132]	588	& ' data type: ' , sdf(jf)%cltype
	589	END DO
	590	ENDIF
	591
	592	END SUBROUTINE fld_fill
	593
	594
[1275]	595	SUBROUTINE wgt_list( sd, kwgt )
	596	!!---------------------------------------------------------------------
	597	!! * ROUTINE wgt_list *
	598	!!
	599	!! ** Purpose : search array of WGTs and find a weights file
	600	!! entry, or return a new one adding it to the end
	601	!! if it is a new entry, the weights data is read in and
	602	!! restructured (fld_weight)
	603	!!
	604	!! ** Method :
	605	!!----------------------------------------------------------------------
	606	TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file
	607	INTEGER, INTENT(inout) :: kwgt ! index of weights
	608	!!
	609	INTEGER :: kw
	610	INTEGER :: nestid
	611	LOGICAL :: found
	612	!!----------------------------------------------------------------------
	613	!
	614	!! search down linked list
	615	!! weights filename is either present or we hit the end of the list
	616	found = .FALSE.
	617
	618	!! because agrif nest part of filenames are now added in iom_open
	619	!! to distinguish between weights files on the different grids, need to track
	620	!! nest number explicitly
	621	nestid = 0
	622	#if defined key_agrif
	623	nestid = Agrif_Fixed()
	624	#endif
	625	DO kw = 1, nxt_wgt-1
	626	IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
	627	ref_wgts(kw)%nestid == nestid) THEN
	628	kwgt = kw
	629	found = .TRUE.
	630	EXIT
	631	ENDIF
	632	END DO
	633	IF( .NOT.found ) THEN
	634	kwgt = nxt_wgt
	635	CALL fld_weight( sd )
	636	ENDIF
	637
	638	END SUBROUTINE wgt_list
	639
	640	SUBROUTINE wgt_print( )
	641	!!---------------------------------------------------------------------
	642	!! * ROUTINE wgt_print *
	643	!!
	644	!! ** Purpose : print the list of known weights
	645	!!
	646	!! ** Method :
	647	!!----------------------------------------------------------------------
	648	!!
	649	INTEGER :: kw
	650	!!----------------------------------------------------------------------
	651	!
	652
	653	DO kw = 1, nxt_wgt-1
	654	WRITE(numout,*) 'weight file: ',TRIM(ref_wgts(kw)%wgtname)
	655	WRITE(numout,*) ' ddims: ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
	656	WRITE(numout,*) ' numwgt: ',ref_wgts(kw)%numwgt
	657	WRITE(numout,*) ' jpiwgt: ',ref_wgts(kw)%jpiwgt
	658	WRITE(numout,*) ' jpjwgt: ',ref_wgts(kw)%jpjwgt
	659	WRITE(numout,*) ' botleft: ',ref_wgts(kw)%botleft
	660	WRITE(numout,*) ' topright: ',ref_wgts(kw)%topright
	661	IF( ref_wgts(kw)%cyclic ) THEN
	662	WRITE(numout,*) ' cyclical'
	663	IF( ref_wgts(kw)%offset > 0 ) WRITE(numout,*) ' with offset'
	664	ELSE
	665	WRITE(numout,*) ' not cyclical'
	666	ENDIF
	667	IF( ASSOCIATED(ref_wgts(kw)%data_wgt) ) WRITE(numout,*) ' allocated'
	668	END DO
	669
	670	END SUBROUTINE wgt_print
	671
	672	SUBROUTINE fld_weight( sd )
	673	!!---------------------------------------------------------------------
	674	!! * ROUTINE fld_weight *
	675	!!
	676	!! ** Purpose : create a new WGT structure and fill in data from
	677	!! file, restructuring as required
	678	!!
	679	!! ** Method :
	680	!!----------------------------------------------------------------------
	681	TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file
	682	!!
	683	INTEGER :: jn ! dummy loop indices
	684	INTEGER :: inum ! temporary logical unit
	685	INTEGER :: id ! temporary variable id
	686	CHARACTER (len=5) :: aname
	687	INTEGER , DIMENSION(3) :: ddims
	688	INTEGER , DIMENSION(jpi, jpj) :: data_src
	689	REAL(wp), DIMENSION(jpi, jpj) :: data_tmp
	690	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: line2, lines ! temporary array to read 2 lineumns
	691	CHARACTER (len=34) :: lonvar
	692	LOGICAL :: cyclical
	693	REAL(wp) :: resid, dlon ! temporary array to read 2 lineumns
	694	INTEGER :: offset ! temporary integer
	695	!!----------------------------------------------------------------------
	696	!
	697	IF( nxt_wgt > tot_wgts ) THEN
	698	CALL ctl_stop("fld_weights: weights array size exceeded, increase tot_wgts")
	699	ENDIF
	700	!
	701	!! new weights file entry, add in extra information
	702	!! a weights file represents a 2D grid of a certain shape, so we assume that the current
	703	!! input data file is representative of all other files to be opened and processed with the
	704	!! current weights file
	705
	706	!! open input data file (non-model grid)
[1319]	707	CALL iom_open( sd%clname, inum, ldiof = LEN(TRIM(sd%wgtname)) > 0 )
[1275]	708
	709	!! get dimensions
	710	id = iom_varid( inum, sd%clvar, ddims )
	711
	712	!! check for an east-west cyclic grid
	713	!! try to guess name of longitude variable
	714
	715	lonvar = 'nav_lon'
	716	id = iom_varid(inum, TRIM(lonvar), ldstop=.FALSE.)
	717	IF( id <= 0 ) THEN
	718	lonvar = 'lon'
	719	id = iom_varid(inum, TRIM(lonvar), ldstop=.FALSE.)
	720	ENDIF
	721
	722	offset = -1
	723	cyclical = .FALSE.
	724	IF( id > 0 ) THEN
	725	!! found a longitude variable
	726	!! now going to assume that grid is regular so we can read a single row
	727
	728	!! because input array is 2d, have to present iom with 2d array even though we only need 1d slice
	729	!! worse, we cant pass line2(:,1) to iom_get since this is treated as a 1d array which doesnt match input file
	730	ALLOCATE( lines(ddims(1),2) )
	731	CALL iom_get(inum, jpdom_unknown, lonvar, lines(:,:), 1, kstart=(/1,1/), kcount=(/ddims(1),2/) )
	732
	733	!! find largest grid spacing
	734	lines(1:ddims(1)-1,2) = lines(2:ddims(1),1) - lines(1:ddims(1)-1,1)
	735	dlon = MAXVAL( lines(1:ddims(1)-1,2) )
	736
	737	resid = ABS(ABS(lines(ddims(1),1)-lines(1,1))-360.0)
	738	IF( resid < rsmall ) THEN
	739	!! end rows overlap in longitude
	740	offset = 0
	741	cyclical = .TRUE.
	742	ELSEIF( resid < 2.0*dlon ) THEN
	743	!! also call it cyclic if difference between end points is less than twice dlon from 360
	744	offset = 1
	745	cyclical = .TRUE.
	746	ENDIF
	747
	748	DEALLOCATE( lines )
	749
	750	ELSE
	751	!! guessing failed
	752	!! read in first and last columns of data variable
	753	!! since we dont know the name of the longitude variable (or even if there is one)
	754	!! we assume that if these two columns are equal, file is cyclic east-west
	755
	756	!! because input array is 2d, have to present iom with 2d array even though we only need 1d slice
	757	!! worse, we cant pass line2(1,:) to iom_get since this is treated as a 1d array which doesnt match input file
	758	ALLOCATE( lines(2,ddims(2)), line2(2,ddims(2)) )
	759	CALL iom_get(inum, jpdom_unknown, sd%clvar, line2(:,:), 1, kstart=(/1,1/), kcount=(/2,ddims(2)/) )
	760	lines(2,:) = line2(1,:)
	761
	762	CALL iom_get(inum, jpdom_unknown, sd%clvar, line2(:,:), 1, kstart=(/ddims(1)-1,1/), kcount=(/2,ddims(2)/) )
	763	lines(1,:) = line2(2,:)
	764
	765	resid = SUM( ABS(lines(1,:) - lines(2,:)) )
	766	IF( resid < ddims(2)*rsmall ) THEN
	767	offset = 0
	768	cyclical = .TRUE.
	769	ENDIF
	770
	771	DEALLOCATE( lines, line2 )
	772	ENDIF
	773
	774	!! close it
	775	CALL iom_close( inum )
	776
	777	!! now open the weights file
	778
	779	CALL iom_open ( sd%wgtname, inum ) ! interpolation weights
	780	IF ( inum > 0 ) THEN
	781
	782	ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
	783	ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
	784	ref_wgts(nxt_wgt)%wgtname = sd%wgtname
	785	ref_wgts(nxt_wgt)%offset = -1
	786	ref_wgts(nxt_wgt)%cyclic = .FALSE.
	787	IF( cyclical ) THEN
	788	ref_wgts(nxt_wgt)%offset = offset
	789	ref_wgts(nxt_wgt)%cyclic = .TRUE.
	790	ENDIF
	791	ref_wgts(nxt_wgt)%nestid = 0
	792	#if defined key_agrif
	793	ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
	794	#endif
	795	!! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
	796	!! for each weight wgtNN there is an integer array srcNN which gives the point in
	797	!! the input data grid which is to be multiplied by the weight
	798	!! they are both arrays on the model grid so the result of the multiplication is
	799	!! added into an output array on the model grid as a running sum
	800
	801	!! two possible cases: bilinear (4 weights) or bicubic (16 weights)
	802	id = iom_varid(inum, 'src05', ldstop=.FALSE.)
	803	IF( id <= 0) THEN
	804	ref_wgts(nxt_wgt)%numwgt = 4
	805	ELSE
	806	ref_wgts(nxt_wgt)%numwgt = 16
	807	ENDIF
	808
	809	ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
	810	ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
	811	ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
	812
	813	DO jn = 1,4
	814	aname = ' '
	815	WRITE(aname,'(a3,i2.2)') 'src',jn
	816	data_tmp(:,:) = 0
	817	CALL iom_get ( inum, jpdom_unknown, aname, data_tmp(1:nlci,1:nlcj), &
	818	kstart=(/nimpp,njmpp/), kcount=(/nlci,nlcj/) )
	819	data_src(:,:) = INT(data_tmp(:,:))
	820	ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
	821	ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
	822	END DO
	823
	824	DO jn = 1, ref_wgts(nxt_wgt)%numwgt
	825	aname = ' '
	826	WRITE(aname,'(a3,i2.2)') 'wgt',jn
	827	ref_wgts(nxt_wgt)%data_wgt(1:nlci,1:nlcj,jn) = 0.0
	828	CALL iom_get ( inum, jpdom_unknown, aname, ref_wgts(nxt_wgt)%data_wgt(1:nlci,1:nlcj,jn), &
	829	kstart=(/nimpp,njmpp/), kcount=(/nlci,nlcj/) )
	830	END DO
	831	CALL iom_close (inum)
	832
	833	! find min and max indices in grid
	834	ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(1:nlci,1:nlcj,:))
	835	ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(1:nlci,1:nlcj,:))
	836	ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(1:nlci,1:nlcj,:))
	837	ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(1:nlci,1:nlcj,:))
	838
	839	! and therefore dimensions of the input box
	840	ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
	841	ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
	842
	843	! shift indexing of source grid
	844	ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
	845	ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
	846
	847	! create input grid, give it a halo to allow gradient calculations
[1702]	848	! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
	849	! a more robust solution will be given in next release
	850	ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3) )
	851	IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col2(2,ref_wgts(nxt_wgt)%jpjwgt+3) )
[1275]	852
	853	nxt_wgt = nxt_wgt + 1
	854
	855	ELSE
	856	CALL ctl_stop( ' fld_weight : unable to read the file ' )
	857	ENDIF
	858
	859	END SUBROUTINE fld_weight
	860
	861	SUBROUTINE fld_interp(num, clvar, kw, dta, nrec)
	862	!!---------------------------------------------------------------------
	863	!! * ROUTINE fld_interp *
	864	!!
	865	!! ** Purpose : apply weights to input gridded data to create data
	866	!! on model grid
	867	!!
	868	!! ** Method :
	869	!!----------------------------------------------------------------------
	870	INTEGER, INTENT(in) :: num ! stream number
	871	CHARACTER(LEN=*), INTENT(in) :: clvar ! variable name
	872	INTEGER, INTENT(in) :: kw ! weights number
	873	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: dta ! output field on model grid
	874	INTEGER, INTENT(in) :: nrec ! record number to read (ie time slice)
	875	!!
	876	INTEGER, DIMENSION(2) :: rec1,recn ! temporary arrays for start and length
	877	INTEGER :: jk, jn, jm ! loop counters
	878	INTEGER :: ni, nj ! lengths
	879	INTEGER :: jpimin,jpiwid ! temporary indices
	880	INTEGER :: jpjmin,jpjwid ! temporary indices
	881	INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices
	882	!!----------------------------------------------------------------------
	883	!
	884
	885	!! for weighted interpolation we have weights at four corners of a box surrounding
	886	!! a model grid point, each weight is multiplied by a grid value (bilinear case)
	887	!! or by a grid value and gradients at the corner point (bicubic case)
	888	!! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
	889
	890	!! sub grid where we already have weights
	891	jpimin = ref_wgts(kw)%botleft(1)
	892	jpjmin = ref_wgts(kw)%botleft(2)
	893	jpiwid = ref_wgts(kw)%jpiwgt
	894	jpjwid = ref_wgts(kw)%jpjwgt
	895
	896	!! what we need to read into sub grid in order to calculate gradients
	897	rec1(1) = MAX( jpimin-1, 1 )
	898	rec1(2) = MAX( jpjmin-1, 1 )
	899	recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
	900	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
	901
	902	!! where we need to read it to
	903	jpi1 = 2 + rec1(1) - jpimin
	904	jpj1 = 2 + rec1(2) - jpjmin
	905	jpi2 = jpi1 + recn(1) - 1
	906	jpj2 = jpj1 + recn(2) - 1
	907
	908	ref_wgts(kw)%fly_dta(:,:) = 0.0
	909	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2), nrec, rec1, recn)
	910
	911	!! first four weights common to both bilinear and bicubic
	912	!! note that we have to offset by 1 into fly_dta array because of halo
	913	dta(:,:) = 0.0
	914	DO jk = 1,4
	915	DO jn = 1, jpj
	916	DO jm = 1,jpi
	917	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
	918	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
	919	dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1)
	920	END DO
	921	END DO
	922	END DO
	923
	924	IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
	925
	926	!! fix up halo points that we couldnt read from file
	927	IF( jpi1 == 2 ) THEN
	928	ref_wgts(kw)%fly_dta(jpi1-1,:) = ref_wgts(kw)%fly_dta(jpi1,:)
	929	ENDIF
	930	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
	931	ref_wgts(kw)%fly_dta(jpi2+1,:) = ref_wgts(kw)%fly_dta(jpi2,:)
	932	ENDIF
	933	IF( jpj1 == 2 ) THEN
	934	ref_wgts(kw)%fly_dta(:,jpj1-1) = ref_wgts(kw)%fly_dta(:,jpj1)
	935	ENDIF
	936	IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
	937	ref_wgts(kw)%fly_dta(:,jpj2+1) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2) - ref_wgts(kw)%fly_dta(:,jpj2-1)
	938	ENDIF
	939
	940	!! if data grid is cyclic we can do better on east-west edges
	941	!! but have to allow for whether first and last columns are coincident
	942	IF( ref_wgts(kw)%cyclic ) THEN
	943	rec1(2) = MAX( jpjmin-1, 1 )
	944	recn(1) = 2
	945	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
	946	jpj1 = 2 + rec1(2) - jpjmin
	947	jpj2 = jpj1 + recn(2) - 1
	948	IF( jpi1 == 2 ) THEN
	949	rec1(1) = ref_wgts(kw)%ddims(1) - 1
	950	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2), nrec, rec1, recn)
	951	ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2) = ref_wgts(kw)%col2(ref_wgts(kw)%offset+1,jpj1:jpj2)
	952	ENDIF
	953	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
	954	rec1(1) = 1
	955	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2), nrec, rec1, recn)
	956	ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2) = ref_wgts(kw)%col2(2-ref_wgts(kw)%offset,jpj1:jpj2)
	957	ENDIF
	958	ENDIF
	959
	960	! gradient in the i direction
	961	DO jk = 1,4
	962	DO jn = 1, jpj
	963	DO jm = 1,jpi
	964	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
	965	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
	966	dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 * &
	967	(ref_wgts(kw)%fly_dta(ni+2,nj+1) - ref_wgts(kw)%fly_dta(ni,nj+1))
	968	END DO
	969	END DO
	970	END DO
	971
	972	! gradient in the j direction
	973	DO jk = 1,4
	974	DO jn = 1, jpj
	975	DO jm = 1,jpi
	976	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
	977	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
	978	dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 * &
	979	(ref_wgts(kw)%fly_dta(ni+1,nj+2) - ref_wgts(kw)%fly_dta(ni+1,nj))
	980	END DO
	981	END DO
	982	END DO
	983
	984	! gradient in the ij direction
	985	DO jk = 1,4
	986	DO jn = 1, jpj
	987	DO jm = 1,jpi
	988	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
	989	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
	990	dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
	991	(ref_wgts(kw)%fly_dta(ni+2,nj+2) - ref_wgts(kw)%fly_dta(ni ,nj+2)) - &
	992	(ref_wgts(kw)%fly_dta(ni+2,nj ) - ref_wgts(kw)%fly_dta(ni ,nj )))
	993	END DO
	994	END DO
	995	END DO
	996
	997	END IF
	998
	999	END SUBROUTINE fld_interp
	1000
[888]	1001	END MODULE fldread

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

Download in other formats: