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

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

fldread.F90 in branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 @ 2601

Last change on this file since 2601 was 2590, checked in by trackstand2, 13 years ago
Merge branch 'dynamic_memory' into master-svn-dyn
Property svn:keywords set to `Id`
File size: 62.2 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
[2528]	17	USE ioipsl, ONLY : ymds2ju, ju2ymds ! for calendar
[888]	18	USE phycst ! ???
	19	USE in_out_manager ! I/O manager
	20	USE iom ! I/O manager library
[1275]	21	USE geo2ocean ! for vector rotation on to model grid
[888]	22
	23	IMPLICIT NONE
	24	PRIVATE
	25
	26	TYPE, PUBLIC :: FLD_N !: Namelist field informations
[1730]	27	CHARACTER(len = 256) :: clname ! generic name of the NetCDF flux file
	28	INTEGER :: nfreqh ! frequency of each flux file
	29	CHARACTER(len = 34) :: clvar ! generic name of the variable in the NetCDF flux file
	30	LOGICAL :: ln_tint ! time interpolation or not (T/F)
	31	LOGICAL :: ln_clim ! climatology or not (T/F)
[2528]	32	CHARACTER(len = 8) :: cltype ! type of data file 'daily', 'monthly' or yearly'
[1730]	33	CHARACTER(len = 34) :: wname ! generic name of a NetCDF weights file to be used, blank if not
	34	CHARACTER(len = 34) :: vcomp ! symbolic component name if a vector that needs rotation
[2528]	35	! a string starting with "U" or "V" for each component
	36	! chars 2 onwards identify which components go together
[888]	37	END TYPE FLD_N
	38
	39	TYPE, PUBLIC :: FLD !: Input field related variables
	40	CHARACTER(len = 256) :: clrootname ! generic name of the NetCDF file
	41	CHARACTER(len = 256) :: clname ! current name of the NetCDF file
[1730]	42	INTEGER :: nfreqh ! frequency of each flux file
[888]	43	CHARACTER(len = 34) :: clvar ! generic name of the variable in the NetCDF flux file
	44	LOGICAL :: ln_tint ! time interpolation or not (T/F)
[1132]	45	LOGICAL :: ln_clim ! climatology or not (T/F)
[2528]	46	CHARACTER(len = 8) :: cltype ! type of data file 'daily', 'monthly' or yearly'
[1132]	47	INTEGER :: num ! iom id of the jpfld files to be read
[1730]	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)
[2528]	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
[2528]	55	LOGICAL :: 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
[2528]	73	INTEGER :: overlap ! =0 when cyclic grid has no overlapping EW columns
	74	! =>1 when they have one or more overlapping columns
	75	! =-1 not cyclic
[1275]	76	LOGICAL :: cyclic ! east-west cyclic or not
[2528]	77	INTEGER, DIMENSION(:,:,:), POINTER :: data_jpi ! array of source integers
	78	INTEGER, DIMENSION(:,:,:), POINTER :: data_jpj ! array of source integers
[1275]	79	REAL(wp), DIMENSION(:,:,:), POINTER :: data_wgt ! array of weights on model grid
[2528]	80	REAL(wp), DIMENSION(:,:,:), POINTER :: fly_dta ! array of values on input grid
	81	REAL(wp), DIMENSION(:,:,:), POINTER :: col ! temporary array for reading in columns
[1275]	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	!!----------------------------------------------------------------------
[2528]	93	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1156]	94	!! $Id$
[2528]	95	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[888]	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	!!
[2528]	116	INTEGER :: imf ! size of the structure sd
[1132]	117	INTEGER :: jf ! dummy indices
[1730]	118	INTEGER :: ireclast ! last record to be read in the current year file
	119	INTEGER :: isecend ! number of second since Jan. 1st 00h of nit000 year at nitend
[2323]	120	INTEGER :: isecsbc ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step
[1628]	121	LOGICAL :: llnxtyr ! open next year file?
	122	LOGICAL :: llnxtmth ! open next month file?
	123	LOGICAL :: llstop ! stop is the file does not exist
[1132]	124	REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation
	125	REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation
[1191]	126	CHARACTER(LEN=1000) :: clfmt ! write format
[888]	127	!!---------------------------------------------------------------------
[2528]	128	! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
	129	isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdttra(1)) ! middle of sbc time step
[1275]	130	imf = SIZE( sd )
[2323]	131	!
[2528]	132	IF( kt == nit000 ) THEN ! initialization
	133	DO jf = 1, imf
	134	CALL fld_init( kn_fsbc, sd(jf) ) ! read each before field (put them in after as they will be swapped)
	135	END DO
	136	IF( lwp ) CALL wgt_print() ! control print
	137	CALL fld_rot( kt, sd ) ! rotate vector fiels if needed
	138	ENDIF
	139	! ! ====================================== !
	140	IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN ! update field at each kn_fsbc time-step !
	141	! ! ====================================== !
[888]	142	!
[2528]	143	DO jf = 1, imf ! --- loop over field --- !
	144
	145	IF( isecsbc > sd(jf)%nrec_a(2) .OR. kt == nit000 ) THEN ! read/update the after data?
[888]	146
[2528]	147	IF( sd(jf)%ln_tint ) THEN ! swap before record field and informations
	148	sd(jf)%nrec_b(:) = sd(jf)%nrec_a(:)
[888]	149	!CDIR COLLAPSE
[2528]	150	sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2)
	151	ENDIF
[1132]	152
[2528]	153	CALL fld_rec( kn_fsbc, sd(jf) ) ! update record informations
[1132]	154
[2528]	155	! do we have to change the year/month/week/day of the forcing field??
	156	IF( sd(jf)%ln_tint ) THEN
	157	! if we do time interpolation we will need to open next year/month/week/day file before the end of the current
	158	! one. If so, we are still before the end of the year/month/week/day when calling fld_rec so sd(jf)%nrec_a(1)
	159	! will be larger than the record number that should be read for current year/month/week/day
[1132]	160
[2528]	161	! last record to be read in the current file
	162	IF ( sd(jf)%nfreqh == -12 ) THEN ; ireclast = 1 ! yearly mean
	163	ELSEIF( sd(jf)%nfreqh == -1 ) THEN ! monthly mean
	164	IF( sd(jf)%cltype == 'monthly' ) THEN ; ireclast = 1
	165	ELSE ; ireclast = 12
	166	ENDIF
	167	ELSE ! higher frequency mean (in hours)
	168	IF( sd(jf)%cltype == 'monthly' ) THEN ; ireclast = 24 * nmonth_len(nmonth) / sd(jf)%nfreqh
	169	ELSEIF( sd(jf)%cltype(1:4) == 'week' ) THEN ; ireclast = 24 * 7 / sd(jf)%nfreqh
	170	ELSEIF( sd(jf)%cltype == 'daily' ) THEN ; ireclast = 24 / sd(jf)%nfreqh
	171	ELSE ; ireclast = 24 * nyear_len( 1 ) / sd(jf)%nfreqh
	172	ENDIF
[1132]	173	ENDIF
	174
[2528]	175	! do we need next file data?
	176	IF( sd(jf)%nrec_a(1) > ireclast ) THEN
[1132]	177
[2528]	178	sd(jf)%nrec_a(1) = 1 ! force to read the first record of the next file
[1628]	179
[2528]	180	IF( .NOT. sd(jf)%ln_clim ) THEN ! close the current file and open a new one.
[1628]	181
[2528]	182	llnxtmth = sd(jf)%cltype == 'monthly' .OR. nday == nmonth_len(nmonth) ! open next month file?
	183	llnxtyr = sd(jf)%cltype == 'yearly' .OR. (nmonth == 12 .AND. llnxtmth) ! open next year file?
[1132]	184
[2528]	185	! if the run finishes at the end of the current year/month/week/day, we will allow next
	186	! year/month/week/day file to be not present. If the run continue further than the current
	187	! year/month/week/day, next year/month/week/day file must exist
	188	isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdttra(1)) ! second at the end of the run
	189	llstop = isecend > sd(jf)%nrec_a(2) ! read more than 1 record of next year
[1132]	190
[2528]	191	CALL fld_clopn( sd(jf), nyear + COUNT((/llnxtyr /)) , &
	192	& nmonth + COUNT((/llnxtmth/)) - 12 * COUNT((/llnxtyr /)), &
	193	& nday + 1 - nmonth_len(nmonth) * COUNT((/llnxtmth/)), llstop )
	194
	195	IF( sd(jf)%num <= 0 .AND. .NOT. llstop ) THEN ! next year file does not exist
	196	CALL ctl_warn('next year/month/week/day file: '//TRIM(sd(jf)%clname)// &
	197	& ' not present -> back to current year/month/day')
	198	CALL fld_clopn( sd(jf), nyear, nmonth, nday ) ! back to the current year/month/day
	199	sd(jf)%nrec_a(1) = ireclast ! force to read the last record to be read in the current year file
	200	ENDIF
	201
[1132]	202	ENDIF
[2528]	203	ENDIF
[1132]	204
[2528]	205	ELSE
	206	! if we are not doing time interpolation, we must change the year/month/week/day of the file just after
	207	! switching to the NEW year/month/week/day. If it is the case, we are at the beginning of the
	208	! year/month/week/day when calling fld_rec so sd(jf)%nrec_a(1) = 1
	209	IF( sd(jf)%nrec_a(1) == 1 .AND. .NOT. ( sd(jf)%ln_clim .AND. sd(jf)%cltype == 'yearly' ) ) &
	210	& CALL fld_clopn( sd(jf), nyear, nmonth, nday )
[1132]	211	ENDIF
	212
[2528]	213	! read after data
	214	CALL fld_get( sd(jf) )
	215
[1275]	216	ENDIF
[2528]	217	END DO ! --- end loop over field --- !
[1132]	218
[2528]	219	CALL fld_rot( kt, sd ) ! rotate vector fiels if needed
[888]	220
[2528]	221	DO jf = 1, imf ! --- loop over field --- !
[888]	222	!
[2528]	223	IF( sd(jf)%ln_tint ) THEN ! temporal interpolation
[1191]	224	IF(lwp .AND. kt - nit000 <= 100 ) THEN
[2528]	225	clfmt = "('fld_read: var ', a, ' kt = ', i8, ' (', f7.2,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," // &
	226	& "', records b/a: ', i4.4, '/', i4.4, ' (days ', f7.2,'/', f7.2, ')')"
	227	WRITE(numout, clfmt) TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday, &
[1730]	228	& 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]	229	ENDIF
[2528]	230	! temporal interpolation weights
[2323]	231	ztinta = REAL( isecsbc - sd(jf)%nrec_b(2), wp ) / REAL( sd(jf)%nrec_a(2) - sd(jf)%nrec_b(2), wp )
[1132]	232	ztintb = 1. - ztinta
[888]	233	!CDIR COLLAPSE
[2528]	234	sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:,1) + ztinta * sd(jf)%fdta(:,:,:,2)
	235	ELSE ! nothing to do...
[1191]	236	IF(lwp .AND. kt - nit000 <= 100 ) THEN
[2528]	237	clfmt = "('fld_read: var ', a, ' kt = ', i8,' (', f7.2,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," // &
	238	& "', record: ', i4.4, ' (days ', f7.2, ' <-> ', f7.2, ')')"
	239	WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday, &
	240	& sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday
[1191]	241	ENDIF
[888]	242	ENDIF
	243	!
[2528]	244	IF( kt == nitend - kn_fsbc + 1 ) CALL iom_close( sd(jf)%num ) ! Close the input files
[1132]	245
[2528]	246	END DO ! --- end loop over field --- !
	247	!
	248	! ! ====================================== !
	249	ENDIF ! update field at each kn_fsbc time-step !
	250	! ! ====================================== !
	251	!
[888]	252	END SUBROUTINE fld_read
	253
	254
[2323]	255	SUBROUTINE fld_init( kn_fsbc, sdjf )
[888]	256	!!---------------------------------------------------------------------
[1132]	257	!! * ROUTINE fld_init *
	258	!!
	259	!! ** Purpose : - if time interpolation, read before data
	260	!! - open current year file
	261	!!
	262	!! ** Method :
	263	!!----------------------------------------------------------------------
[2528]	264	INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step)
	265	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
[1132]	266	!!
[2528]	267	LOGICAL :: llprevyr ! are we reading previous year file?
	268	LOGICAL :: llprevmth ! are we reading previous month file?
	269	LOGICAL :: llprevweek ! are we reading previous week file?
	270	LOGICAL :: llprevday ! are we reading previous day file?
	271	LOGICAL :: llprev ! llprevyr .OR. llprevmth .OR. llprevweek .OR. llprevday
	272	INTEGER :: idvar ! variable id
	273	INTEGER :: inrec ! number of record existing for this variable
	274	INTEGER :: iyear, imonth, iday ! first day of the current file in yyyy mm dd
	275	INTEGER :: isec_week ! number of seconds since start of the weekly file
[1191]	276	CHARACTER(LEN=1000) :: clfmt ! write format
[1132]	277	!!---------------------------------------------------------------------
[2528]	278
[1132]	279	! some default definitions...
	280	sdjf%num = 0 ! default definition for non-opened file
	281	IF( sdjf%ln_clim ) sdjf%clname = TRIM( sdjf%clrootname ) ! file name defaut definition, never change in this case
[2528]	282	llprevyr = .FALSE.
	283	llprevmth = .FALSE.
	284	llprevweek = .FALSE.
	285	llprevday = .FALSE.
	286	isec_week = 0
[1132]	287
[2528]	288	IF( sdjf%cltype(1:4) == 'week' .AND. nn_leapy == 0 ) &
	289	& CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdjf%clrootname)//') needs nn_leapy = 1')
	290	IF( sdjf%cltype(1:4) == 'week' .AND. sdjf%ln_clim ) &
	291	& CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdjf%clrootname)//') needs ln_clim = .FALSE.')
	292
[1132]	293	! define record informations
[2528]	294	CALL fld_rec( kn_fsbc, sdjf, ldbefore = .TRUE. ) ! return before values in sdjf%nrec_a (as we will swap it later)
[1132]	295
[2528]	296	! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
[2323]	297
[1132]	298	IF( sdjf%ln_tint ) THEN ! we need to read the previous record and we will put it in the current record structure
[2528]	299
	300	IF( sdjf%nrec_a(1) == 0 ) THEN ! we redefine record sdjf%nrec_a(1) with the last record of previous year file
	301	IF ( sdjf%nfreqh == -12 ) THEN ! yearly mean
	302	IF( sdjf%cltype == 'yearly' ) THEN ! yearly file
	303	sdjf%nrec_a(1) = 1 ! force to read the unique record
	304	llprevyr = .NOT. sdjf%ln_clim ! use previous year file?
	305	ELSE
	306	CALL ctl_stop( "fld_init: yearly mean file must be in a yearly type of file: "//TRIM(sdjf%clname) )
	307	ENDIF
	308	ELSEIF( sdjf%nfreqh == -1 ) THEN ! monthly mean
	309	IF( sdjf%cltype == 'monthly' ) THEN ! monthly file
	310	sdjf%nrec_a(1) = 1 ! force to read the unique record
	311	llprevmth = .TRUE. ! use previous month file?
[1628]	312	llprevyr = llprevmth .AND. nmonth == 1 ! use previous year file?
[2528]	313	ELSE ! yearly file
	314	sdjf%nrec_a(1) = 12 ! force to read december mean
[1628]	315	llprevyr = .NOT. sdjf%ln_clim ! use previous year file?
	316	ENDIF
[2528]	317	ELSE ! higher frequency mean (in hours)
	318	IF ( sdjf%cltype == 'monthly' ) THEN ! monthly file
	319	sdjf%nrec_a(1) = 24 * nmonth_len(nmonth-1) / sdjf%nfreqh ! last record of previous month
	320	llprevmth = .TRUE. ! use previous month file?
[1628]	321	llprevyr = llprevmth .AND. nmonth == 1 ! use previous year file?
[2528]	322	ELSEIF( sdjf%cltype(1:4) == 'week' ) THEN ! weekly file
	323	llprevweek = .TRUE. ! use previous week file?
	324	sdjf%nrec_a(1) = 24 * 7 / sdjf%nfreqh ! last record of previous week
	325	isec_week = NINT(rday) * 7 ! add a shift toward previous week
	326	ELSEIF( sdjf%cltype == 'daily' ) THEN ! daily file
	327	sdjf%nrec_a(1) = 24 / sdjf%nfreqh ! last record of previous day
	328	llprevday = .TRUE. ! use previous day file?
[1628]	329	llprevmth = llprevday .AND. nday == 1 ! use previous month file?
	330	llprevyr = llprevmth .AND. nmonth == 1 ! use previous year file?
[2528]	331	ELSE ! yearly file
	332	sdjf%nrec_a(1) = 24 * nyear_len(0) / sdjf%nfreqh ! last record of previous year
[1628]	333	llprevyr = .NOT. sdjf%ln_clim ! use previous year file?
[1132]	334	ENDIF
	335	ENDIF
	336	ENDIF
[2528]	337	IF ( sdjf%cltype(1:4) == 'week' ) THEN
	338	isec_week = isec_week + ksec_week( sdjf%cltype(6:8) ) ! second since the beginning of the week
	339	llprevmth = isec_week > nsec_month ! longer time since the beginning of the week than the month
	340	llprevyr = llprevmth .AND. nmonth == 1
	341	ENDIF
	342	llprev = llprevyr .OR. llprevmth .OR. llprevweek .OR. llprevday
	343	!
	344	iyear = nyear - COUNT((/llprevyr /))
	345	imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
	346	iday = nday - COUNT((/llprevday/)) + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
	347	!
	348	CALL fld_clopn( sdjf, iyear, imonth, iday, .NOT. llprev )
[1132]	349
[1628]	350	! if previous year/month/day file does not exist, we switch to the current year/month/day
[1818]	351	IF( llprev .AND. sdjf%num <= 0 ) THEN
[2528]	352	CALL ctl_warn( 'previous year/month/week/day file: '//TRIM(sdjf%clname)// &
	353	& ' not present -> back to current year/month/week/day' )
[1628]	354	! we force to read the first record of the current year/month/day instead of last record of previous year/month/day
[2528]	355	llprev = .FALSE.
	356	sdjf%nrec_a(1) = 1
[1628]	357	CALL fld_clopn( sdjf, nyear, nmonth, nday )
[1132]	358	ENDIF
	359
[1730]	360	IF( llprev ) THEN ! check if the last record sdjf%nrec_n(1) exists in the file
[1132]	361	idvar = iom_varid( sdjf%num, sdjf%clvar ) ! id of the variable sdjf%clvar
	362	IF( idvar <= 0 ) RETURN
	363	inrec = iom_file( sdjf%num )%dimsz( iom_file( sdjf%num )%ndims(idvar), idvar ) ! size of the last dim of idvar
[2528]	364	sdjf%nrec_a(1) = MIN( sdjf%nrec_a(1), inrec ) ! make sure we select an existing record
[1132]	365	ENDIF
	366
[2528]	367	! read before data
	368	CALL fld_get( sdjf ) ! read before values in after arrays(as we will swap it later)
[1132]	369
[1191]	370	clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i4, ' at time = ', f7.2, ' days')"
[2528]	371	IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_a(1), REAL(sdjf%nrec_a(2),wp)/rday
[1132]	372
[2528]	373	IF( llprev ) CALL iom_close( sdjf%num ) ! force to close previous year file (-> redefine sdjf%num to 0)
[1132]	374
	375	ENDIF
	376
[2528]	377	! make sure current year/month/day file is opened
	378	IF( sdjf%num <= 0 ) THEN
	379	!
	380	IF ( sdjf%cltype(1:4) == 'week' ) THEN
	381	isec_week = ksec_week( sdjf%cltype(6:8) ) ! second since the beginning of the week
	382	llprevmth = isec_week > nsec_month ! longer time since beginning of the week than the month
	383	llprevyr = llprevmth .AND. nmonth == 1
	384	ELSE
	385	isec_week = 0
	386	llprevmth = .FALSE.
	387	llprevyr = .FALSE.
	388	ENDIF
	389	!
	390	iyear = nyear - COUNT((/llprevyr /))
	391	imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
	392	iday = nday + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
	393	!
	394	CALL fld_clopn( sdjf, iyear, imonth, iday )
	395	ENDIF
[1132]	396
	397	END SUBROUTINE fld_init
	398
	399
[2528]	400	SUBROUTINE fld_rec( kn_fsbc, sdjf, ldbefore )
[1132]	401	!!---------------------------------------------------------------------
[888]	402	!! * ROUTINE fld_rec *
	403	!!
[2528]	404	!! ** Purpose : Compute
	405	!! if sdjf%ln_tint = .TRUE.
	406	!! nrec_a: record number and its time (nrec_b is obtained from nrec_a when swapping)
	407	!! if sdjf%ln_tint = .FALSE.
	408	!! nrec_a(1): record number
	409	!! nrec_b(2) and nrec_a(2): time of the beginning and end of the record (for print only)
[888]	410	!!
	411	!! ** Method :
	412	!!----------------------------------------------------------------------
[2528]	413	INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step)
	414	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
	415	LOGICAL , INTENT(in ), OPTIONAL :: ldbefore ! sent back before record values (default = .FALSE.)
	416	! used only if sdjf%ln_tint = .TRUE.
[888]	417	!!
[2528]	418	LOGICAL :: llbefore ! local definition of ldbefore
	419	INTEGER :: iendrec ! end of this record (in seconds)
	420	INTEGER :: imth ! month number
	421	INTEGER :: ifreq_sec ! frequency mean (in seconds)
	422	INTEGER :: isec_week ! number of seconds since the start of the weekly file
[1132]	423	REAL(wp) :: ztmp ! temporary variable
[888]	424	!!----------------------------------------------------------------------
	425	!
[2528]	426	! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar
[2323]	427	!
[2528]	428	IF( PRESENT(ldbefore) ) THEN ; llbefore = ldbefore .AND. sdjf%ln_tint ! needed only if sdjf%ln_tint = .TRUE.
	429	ELSE ; llbefore = .FALSE.
	430	ENDIF
	431	!
	432	! ! =========== !
	433	IF ( sdjf%nfreqh == -12 ) THEN ! yearly mean
	434	! ! =========== !
[888]	435	!
[1132]	436	IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record
	437	!
	438	! INT( ztmp )
	439	! /\|\
	440	! 1 \| *----
	441	! 0 \|----(
	442	! \|----+----\|--> time
[2528]	443	! 0 /\|\ 1 (nday/nyear_len(1))
	444	! \|
	445	! \|
	446	! forcing record : 1
	447	!
	448	ztmp = REAL( nday, wp ) / REAL( nyear_len(1), wp ) + 0.5
	449	sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
	450	! swap at the middle of the year
	451	IF( llbefore ) THEN ; sdjf%nrec_a(2) = nsec1jan000 - NINT(0.5 * rday) * nyear_len(0)
	452	ELSE ; sdjf%nrec_a(2) = nsec1jan000 + NINT(0.5 * rday) * nyear_len(1)
	453	ENDIF
	454	ELSE ! no time interpolation
	455	sdjf%nrec_a(1) = 1
	456	sdjf%nrec_a(2) = NINT(rday) * nyear_len(1) + nsec1jan000 ! swap at the end of the year
	457	sdjf%nrec_b(2) = nsec1jan000 ! beginning of the year (only for print)
	458	ENDIF
	459	!
	460	! ! ============ !
	461	ELSEIF( sdjf%nfreqh == -1 ) THEN ! monthly mean !
	462	! ! ============ !
	463	!
	464	IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record
	465	!
	466	! INT( ztmp )
	467	! /\|\
	468	! 1 \| *----
	469	! 0 \|----(
	470	! \|----+----\|--> time
[1132]	471	! 0 /\|\ 1 (nday/nmonth_len(nmonth))
	472	! \|
	473	! \|
	474	! forcing record : nmonth
	475	!
[2528]	476	ztmp = REAL( nday, wp ) / REAL( nmonth_len(nmonth), wp ) + 0.5
	477	imth = nmonth + INT( ztmp ) - COUNT((/llbefore/))
	478	IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
	479	ELSE ; sdjf%nrec_a(1) = imth
	480	ENDIF
	481	sdjf%nrec_a(2) = nmonth_half( imth ) + nsec1jan000 ! swap at the middle of the month
	482	ELSE ! no time interpolation
	483	IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nrec_a(1) = 1
	484	ELSE ; sdjf%nrec_a(1) = nmonth
	485	ENDIF
	486	sdjf%nrec_a(2) = nmonth_end(nmonth ) + nsec1jan000 ! swap at the end of the month
	487	sdjf%nrec_b(2) = nmonth_end(nmonth-1) + nsec1jan000 ! beginning of the month (only for print)
[888]	488	ENDIF
	489	!
[2528]	490	! ! ================================ !
	491	ELSE ! higher frequency mean (in hours)
	492	! ! ================================ !
[888]	493	!
[2528]	494	ifreq_sec = sdjf%nfreqh * 3600 ! frequency mean (in seconds)
	495	IF( sdjf%cltype(1:4) == 'week' ) isec_week = ksec_week( sdjf%cltype(6:8) ) ! since the first day of the current week
[1132]	496	! number of second since the beginning of the file
[2528]	497	IF( sdjf%cltype == 'monthly' ) THEN ; ztmp = REAL(nsec_month,wp) ! since the first day of the current month
	498	ELSEIF( sdjf%cltype(1:4) == 'week' ) THEN ; ztmp = REAL(isec_week ,wp) ! since the first day of the current week
	499	ELSEIF( sdjf%cltype == 'daily' ) THEN ; ztmp = REAL(nsec_day ,wp) ! since 00h of the current day
	500	ELSE ; ztmp = REAL(nsec_year ,wp) ! since 00h on Jan 1 of the current year
[1132]	501	ENDIF
[2323]	502	ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdttra(1) ! shift time to be centrered in the middle of sbc time step
	503	ztmp = ztmp + 0.01 * rdttra(1) ! add 0.01 time step to avoid truncation error
[1132]	504	IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record
	505	!
	506	! INT( ztmp )
	507	! /\|\
	508	! 2 \| *-----(
	509	! 1 \| *-----(
	510	! 0 \|--(
	511	! \|--+--\|--+--\|--+--\|--> time
[1730]	512	! 0 /\|\ 1 /\|\ 2 /\|\ 3 (nsec_year/ifreq_sec) or (nsec_month/ifreq_sec)
[1132]	513	! \| \| \|
	514	! \| \| \|
	515	! forcing record : 1 2 3
	516	!
[2528]	517	ztmp= ztmp / REAL(ifreq_sec, wp) + 0.5
	518	ELSE ! no time interpolation
[1132]	519	!
	520	! INT( ztmp )
	521	! /\|\
	522	! 2 \| *-----(
	523	! 1 \| *-----(
	524	! 0 \|-----(
	525	! \|--+--\|--+--\|--+--\|--> time
[1730]	526	! 0 /\|\ 1 /\|\ 2 /\|\ 3 (nsec_year/ifreq_sec) or (nsec_month/ifreq_sec)
[1132]	527	! \| \| \|
	528	! \| \| \|
	529	! forcing record : 1 2 3
	530	!
[2528]	531	ztmp= ztmp / REAL(ifreq_sec, wp)
[1132]	532	ENDIF
[2528]	533	sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/)) ! record nomber to be read
[1132]	534
[2528]	535	iendrec = ifreq_sec * sdjf%nrec_a(1) + nsec1jan000 ! end of this record (in second)
	536	! add the number of seconds between 00h Jan 1 and the end of previous month/week/day (ok if nmonth=1)
	537	IF( sdjf%cltype == 'monthly' ) iendrec = iendrec + NINT(rday) * SUM(nmonth_len(1:nmonth -1))
	538	IF( sdjf%cltype(1:4) == 'week' ) iendrec = iendrec + ( nsec_year - isec_week )
	539	IF( sdjf%cltype == 'daily' ) iendrec = iendrec + NINT(rday) * ( nday_year - 1 )
	540	IF( sdjf%ln_tint ) THEN
	541	sdjf%nrec_a(2) = iendrec - ifreq_sec / 2 ! swap at the middle of the record
	542	ELSE
	543	sdjf%nrec_a(2) = iendrec ! swap at the end of the record
	544	sdjf%nrec_b(2) = iendrec - ifreq_sec ! beginning of the record (only for print)
	545	ENDIF
[888]	546	!
	547	ENDIF
	548	!
[1132]	549	END SUBROUTINE fld_rec
	550
	551
[2528]	552	SUBROUTINE fld_get( sdjf )
	553	!!---------------------------------------------------------------------
	554	!! * ROUTINE fld_clopn *
	555	!!
	556	!! ** Purpose : read the data
	557	!!
	558	!! ** Method :
	559	!!----------------------------------------------------------------------
	560	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
	561	!!
	562	INTEGER :: ipk ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
	563	INTEGER :: iw ! index into wgts array
	564	!!---------------------------------------------------------------------
	565
	566	ipk = SIZE( sdjf%fnow, 3 )
	567	IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
	568	CALL wgt_list( sdjf, iw )
	569	IF( sdjf%ln_tint ) THEN ; CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) )
	570	ELSE ; CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fnow(:,:,: ), sdjf%nrec_a(1) )
	571	ENDIF
	572	ELSE
	573	SELECT CASE( ipk )
	574	CASE(1)
	575	IF( sdjf%ln_tint ) THEN ; CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,1,2), sdjf%nrec_a(1) )
	576	ELSE ; CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fnow(:,:,1 ), sdjf%nrec_a(1) )
	577	ENDIF
	578	CASE DEFAULT
	579	IF( sdjf%ln_tint ) THEN ; CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) )
	580	ELSE ; CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1) )
	581	ENDIF
	582	END SELECT
	583	ENDIF
	584	!
	585	sdjf%rotn = .false. ! vector not yet rotated
	586
	587	END SUBROUTINE fld_get
	588
	589
	590	SUBROUTINE fld_rot( kt, sd )
	591	!!---------------------------------------------------------------------
	592	!! * ROUTINE fld_clopn *
	593	!!
	594	!! ** Purpose : Vector fields may need to be rotated onto the local grid direction
	595	!!
	596	!! ** Method :
	597	!!----------------------------------------------------------------------
[2590]	598	USE wrk_nemo, ONLY: wrk_use, wrk_release
	599	USE wrk_nemo, ONLY: utmp => wrk_2d_4, vtmp => wrk_2d_5
	600	!!
[2528]	601	INTEGER , INTENT(in ) :: kt ! ocean time step
	602	TYPE(FLD), INTENT(inout), DIMENSION(:) :: sd ! input field related variables
	603	!!
	604	INTEGER :: ju, jv, jk ! loop indices
	605	INTEGER :: imf ! size of the structure sd
	606	INTEGER :: ill ! character length
	607	INTEGER :: iv ! indice of V component
	608	CHARACTER (LEN=100) :: clcomp ! dummy weight name
	609	!!---------------------------------------------------------------------
[2590]	610
	611	IF(.not. wrk_use(2, 4,5))THEN
	612	CALL ctl_stop('fld_rot: ERROR: requested workspace arrays are unavailable.')
	613	RETURN
	614	END IF
	615
[2528]	616	!! (sga: following code should be modified so that pairs arent searched for each time
	617	!
	618	imf = SIZE( sd )
	619	DO ju = 1, imf
	620	ill = LEN_TRIM( sd(ju)%vcomp )
	621	IF( ill > 0 .AND. .NOT. sd(ju)%rotn ) THEN ! find vector rotations required
	622	IF( sd(ju)%vcomp(1:1) == 'U' ) THEN ! east-west component has symbolic name starting with 'U'
	623	! look for the north-south component which has same symbolic name but with 'U' replaced with 'V'
	624	clcomp = 'V' // sd(ju)%vcomp(2:ill) ! works even if ill == 1
	625	iv = -1
	626	DO jv = 1, imf
	627	IF( TRIM(sd(jv)%vcomp) == TRIM(clcomp) ) iv = jv
	628	END DO
	629	IF( iv > 0 ) THEN ! fields ju and iv are two components which need to be rotated together
	630	DO jk = 1, SIZE( sd(ju)%fnow, 3 )
	631	IF( sd(ju)%ln_tint )THEN
	632	CALL rot_rep( sd(ju)%fdta(:,:,jk,2), sd(iv)%fdta(:,:,jk,2), 'T', 'en->i', utmp(:,:) )
	633	CALL rot_rep( sd(ju)%fdta(:,:,jk,2), sd(iv)%fdta(:,:,jk,2), 'T', 'en->j', vtmp(:,:) )
	634	sd(ju)%fdta(:,:,jk,2) = utmp(:,:) ; sd(iv)%fdta(:,:,jk,2) = vtmp(:,:)
	635	ELSE
	636	CALL rot_rep( sd(ju)%fnow(:,:,jk ), sd(iv)%fnow(:,:,jk ), 'T', 'en->i', utmp(:,:) )
	637	CALL rot_rep( sd(ju)%fnow(:,:,jk ), sd(iv)%fnow(:,:,jk ), 'T', 'en->j', vtmp(:,:) )
	638	sd(ju)%fnow(:,:,jk ) = utmp(:,:) ; sd(iv)%fnow(:,:,jk ) = vtmp(:,:)
	639	ENDIF
	640	END DO
	641	sd(ju)%rotn = .TRUE. ! vector was rotated
	642	IF( lwp .AND. kt == nit000 ) WRITE(numout,*) &
	643	& 'fld_read: vector pair ('//TRIM(sd(ju)%clvar)//', '//TRIM(sd(iv)%clvar)//') rotated on to model grid'
	644	ENDIF
	645	ENDIF
	646	ENDIF
	647	END DO
[2590]	648
	649	IF(.not. wrk_release(2, 4,5))THEN
	650	CALL ctl_stop('fld_rot: ERROR: failed to release workspace arrays.')
	651	END IF
	652
[2528]	653	END SUBROUTINE fld_rot
	654
	655
[1628]	656	SUBROUTINE fld_clopn( sdjf, kyear, kmonth, kday, ldstop )
[1132]	657	!!---------------------------------------------------------------------
	658	!! * ROUTINE fld_clopn *
	659	!!
	660	!! ** Purpose : update the file name and open the file
	661	!!
	662	!! ** Method :
	663	!!----------------------------------------------------------------------
[2528]	664	TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables
	665	INTEGER , INTENT(in ) :: kyear ! year value
	666	INTEGER , INTENT(in ) :: kmonth ! month value
	667	INTEGER , INTENT(in ) :: kday ! day value
	668	LOGICAL , INTENT(in ), OPTIONAL :: ldstop ! stop if open to read a non-existing file (default = .TRUE.)
[1132]	669
	670	IF( sdjf%num /= 0 ) CALL iom_close( sdjf%num ) ! close file if already open
	671	! build the new filename if not climatological data
[2528]	672	sdjf%clname=TRIM(sdjf%clrootname)
	673	!
	674	! note that sdjf%ln_clim is is only acting on presence of the year in the file
	675	IF( .NOT. sdjf%ln_clim ) THEN
	676	WRITE(sdjf%clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), kyear ! add year
	677	IF( sdjf%cltype /= 'yearly' ) WRITE(sdjf%clname, '(a,"m" ,i2.2)' ) TRIM( sdjf%clname ), kmonth ! add month
	678	ELSE
	679	! build the new filename if climatological data
	680	IF( sdjf%cltype /= 'yearly' ) WRITE(sdjf%clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), kmonth ! add month
[888]	681	ENDIF
[2528]	682	IF( sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) &
	683	& WRITE(sdjf%clname, '(a,"d" ,i2.2)' ) TRIM( sdjf%clname ), kday ! add day
	684	!
[1319]	685	CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 )
[888]	686	!
[1132]	687	END SUBROUTINE fld_clopn
	688
	689
	690	SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )
	691	!!---------------------------------------------------------------------
	692	!! * ROUTINE fld_fill *
	693	!!
	694	!! ** Purpose : fill sdf with sdf_n and control print
	695	!!
	696	!! ** Method :
	697	!!----------------------------------------------------------------------
	698	TYPE(FLD) , DIMENSION(:), INTENT(inout) :: sdf ! structure of input fields (file informations, fields read)
	699	TYPE(FLD_N), DIMENSION(:), INTENT(in ) :: sdf_n ! array of namelist information structures
	700	CHARACTER(len=*) , INTENT(in ) :: cdir ! Root directory for location of flx files
	701	CHARACTER(len=*) , INTENT(in ) :: cdcaller !
	702	CHARACTER(len=*) , INTENT(in ) :: cdtitle !
	703	CHARACTER(len=*) , INTENT(in ) :: cdnam !
[888]	704	!
[1132]	705	INTEGER :: jf ! dummy indices
	706	!!---------------------------------------------------------------------
[888]	707
[1132]	708	DO jf = 1, SIZE(sdf)
	709	sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname )
[1730]	710	sdf(jf)%nfreqh = sdf_n(jf)%nfreqh
[1132]	711	sdf(jf)%clvar = sdf_n(jf)%clvar
	712	sdf(jf)%ln_tint = sdf_n(jf)%ln_tint
	713	sdf(jf)%ln_clim = sdf_n(jf)%ln_clim
[2528]	714	sdf(jf)%cltype = sdf_n(jf)%cltype
[1275]	715	sdf(jf)%wgtname = " "
[1730]	716	IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 ) sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
[1275]	717	sdf(jf)%vcomp = sdf_n(jf)%vcomp
[1132]	718	END DO
	719
	720	IF(lwp) THEN ! control print
	721	WRITE(numout,*)
	722	WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
	723	WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
	724	WRITE(numout,*) ' '//TRIM( cdnam )//' Namelist'
	725	WRITE(numout,*) ' list of files and frequency (>0: in hours ; <0 in months)'
	726	DO jf = 1, SIZE(sdf)
	727	WRITE(numout,*) ' root filename: ' , TRIM( sdf(jf)%clrootname ), &
	728	& ' variable name: ' , TRIM( sdf(jf)%clvar )
[1730]	729	WRITE(numout,*) ' frequency: ' , sdf(jf)%nfreqh , &
[1132]	730	& ' time interp: ' , sdf(jf)%ln_tint , &
	731	& ' climatology: ' , sdf(jf)%ln_clim , &
[1275]	732	& ' weights : ' , TRIM( sdf(jf)%wgtname ), &
	733	& ' pairing : ' , TRIM( sdf(jf)%vcomp ), &
[1132]	734	& ' data type: ' , sdf(jf)%cltype
[2528]	735	call flush(numout)
[1132]	736	END DO
	737	ENDIF
	738
	739	END SUBROUTINE fld_fill
	740
	741
[1275]	742	SUBROUTINE wgt_list( sd, kwgt )
	743	!!---------------------------------------------------------------------
	744	!! * ROUTINE wgt_list *
	745	!!
	746	!! ** Purpose : search array of WGTs and find a weights file
	747	!! entry, or return a new one adding it to the end
	748	!! if it is a new entry, the weights data is read in and
	749	!! restructured (fld_weight)
	750	!!
	751	!! ** Method :
	752	!!----------------------------------------------------------------------
	753	TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file
	754	INTEGER, INTENT(inout) :: kwgt ! index of weights
	755	!!
	756	INTEGER :: kw
	757	INTEGER :: nestid
	758	LOGICAL :: found
	759	!!----------------------------------------------------------------------
	760	!
	761	!! search down linked list
	762	!! weights filename is either present or we hit the end of the list
	763	found = .FALSE.
	764
	765	!! because agrif nest part of filenames are now added in iom_open
	766	!! to distinguish between weights files on the different grids, need to track
	767	!! nest number explicitly
	768	nestid = 0
	769	#if defined key_agrif
	770	nestid = Agrif_Fixed()
	771	#endif
	772	DO kw = 1, nxt_wgt-1
	773	IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. &
	774	ref_wgts(kw)%nestid == nestid) THEN
	775	kwgt = kw
	776	found = .TRUE.
	777	EXIT
	778	ENDIF
	779	END DO
	780	IF( .NOT.found ) THEN
	781	kwgt = nxt_wgt
	782	CALL fld_weight( sd )
	783	ENDIF
	784
	785	END SUBROUTINE wgt_list
	786
	787	SUBROUTINE wgt_print( )
	788	!!---------------------------------------------------------------------
	789	!! * ROUTINE wgt_print *
	790	!!
	791	!! ** Purpose : print the list of known weights
	792	!!
	793	!! ** Method :
	794	!!----------------------------------------------------------------------
	795	!!
	796	INTEGER :: kw
	797	!!----------------------------------------------------------------------
	798	!
	799
	800	DO kw = 1, nxt_wgt-1
	801	WRITE(numout,*) 'weight file: ',TRIM(ref_wgts(kw)%wgtname)
	802	WRITE(numout,*) ' ddims: ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
	803	WRITE(numout,*) ' numwgt: ',ref_wgts(kw)%numwgt
	804	WRITE(numout,*) ' jpiwgt: ',ref_wgts(kw)%jpiwgt
	805	WRITE(numout,*) ' jpjwgt: ',ref_wgts(kw)%jpjwgt
	806	WRITE(numout,*) ' botleft: ',ref_wgts(kw)%botleft
	807	WRITE(numout,*) ' topright: ',ref_wgts(kw)%topright
	808	IF( ref_wgts(kw)%cyclic ) THEN
	809	WRITE(numout,*) ' cyclical'
[2528]	810	IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) ' with overlap of ', ref_wgts(kw)%overlap
[1275]	811	ELSE
	812	WRITE(numout,*) ' not cyclical'
	813	ENDIF
	814	IF( ASSOCIATED(ref_wgts(kw)%data_wgt) ) WRITE(numout,*) ' allocated'
	815	END DO
	816
	817	END SUBROUTINE wgt_print
	818
	819	SUBROUTINE fld_weight( sd )
	820	!!---------------------------------------------------------------------
	821	!! * ROUTINE fld_weight *
	822	!!
	823	!! ** Purpose : create a new WGT structure and fill in data from
	824	!! file, restructuring as required
	825	!!
	826	!! ** Method :
	827	!!----------------------------------------------------------------------
[2590]	828	USE wrk_nemo, ONLY: wrk_use, wrk_release, iwrk_use, iwrk_release
	829	USE wrk_nemo, ONLY: data_tmp => wrk_2d_1
	830	USE wrk_nemo, ONLY: data_src => iwrk_2d_1
	831	!!
[1275]	832	TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file
	833	!!
	834	INTEGER :: jn ! dummy loop indices
	835	INTEGER :: inum ! temporary logical unit
	836	INTEGER :: id ! temporary variable id
[2528]	837	INTEGER :: ipk ! temporary vertical dimension
[1275]	838	CHARACTER (len=5) :: aname
	839	INTEGER , DIMENSION(3) :: ddims
	840	LOGICAL :: cyclical
[2528]	841	INTEGER :: zwrap ! temporary integer
[1275]	842	!!----------------------------------------------------------------------
	843	!
[2590]	844	IF( (.NOT. wrk_use(2, 1)) .OR. (.NOT. iwrk_use(2,1)) )THEN
	845	CALL ctl_stop('fld_weights: requested workspace arrays are unavailable.')
	846	RETURN
	847	END IF
	848	!
[1275]	849	IF( nxt_wgt > tot_wgts ) THEN
	850	CALL ctl_stop("fld_weights: weights array size exceeded, increase tot_wgts")
	851	ENDIF
	852	!
	853	!! new weights file entry, add in extra information
	854	!! a weights file represents a 2D grid of a certain shape, so we assume that the current
	855	!! input data file is representative of all other files to be opened and processed with the
	856	!! current weights file
	857
	858	!! open input data file (non-model grid)
[1319]	859	CALL iom_open( sd%clname, inum, ldiof = LEN(TRIM(sd%wgtname)) > 0 )
[1275]	860
	861	!! get dimensions
	862	id = iom_varid( inum, sd%clvar, ddims )
	863
[2528]	864	!! close it
	865	CALL iom_close( inum )
[1275]	866
[2528]	867	!! now open the weights file
[1275]	868
[2528]	869	CALL iom_open ( sd%wgtname, inum ) ! interpolation weights
	870	IF ( inum > 0 ) THEN
[1275]	871
[2528]	872	!! determine whether we have an east-west cyclic grid
	873	!! from global attribute called "ew_wrap" in the weights file
	874	!! note that if not found, iom_getatt returns -999 and cyclic with no overlap is assumed
	875	!! since this is the most common forcing configuration
[1275]	876
[2528]	877	CALL iom_getatt(inum, 'ew_wrap', zwrap)
	878	IF( zwrap >= 0 ) THEN
[1275]	879	cyclical = .TRUE.
[2528]	880	ELSE IF( zwrap == -999 ) THEN
[1275]	881	cyclical = .TRUE.
[2528]	882	zwrap = 0
	883	ELSE
	884	cyclical = .FALSE.
[1275]	885	ENDIF
	886
	887	ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
	888	ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
	889	ref_wgts(nxt_wgt)%wgtname = sd%wgtname
[2528]	890	ref_wgts(nxt_wgt)%overlap = zwrap
	891	ref_wgts(nxt_wgt)%cyclic = cyclical
[1275]	892	ref_wgts(nxt_wgt)%nestid = 0
	893	#if defined key_agrif
	894	ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
	895	#endif
	896	!! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
	897	!! for each weight wgtNN there is an integer array srcNN which gives the point in
	898	!! the input data grid which is to be multiplied by the weight
	899	!! they are both arrays on the model grid so the result of the multiplication is
	900	!! added into an output array on the model grid as a running sum
	901
	902	!! two possible cases: bilinear (4 weights) or bicubic (16 weights)
	903	id = iom_varid(inum, 'src05', ldstop=.FALSE.)
	904	IF( id <= 0) THEN
	905	ref_wgts(nxt_wgt)%numwgt = 4
	906	ELSE
	907	ref_wgts(nxt_wgt)%numwgt = 16
	908	ENDIF
	909
	910	ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) )
	911	ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) )
	912	ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) )
	913
	914	DO jn = 1,4
	915	aname = ' '
	916	WRITE(aname,'(a3,i2.2)') 'src',jn
	917	data_tmp(:,:) = 0
[1955]	918	CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) )
[1275]	919	data_src(:,:) = INT(data_tmp(:,:))
	920	ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1)
	921	ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1)
	922	END DO
	923
	924	DO jn = 1, ref_wgts(nxt_wgt)%numwgt
	925	aname = ' '
	926	WRITE(aname,'(a3,i2.2)') 'wgt',jn
[1955]	927	ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0
	928	CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) )
[1275]	929	END DO
	930	CALL iom_close (inum)
	931
	932	! find min and max indices in grid
[1955]	933	ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
	934	ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
	935	ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
	936	ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
[1275]	937
	938	! and therefore dimensions of the input box
	939	ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
	940	ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
	941
	942	! shift indexing of source grid
	943	ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
	944	ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
	945
	946	! create input grid, give it a halo to allow gradient calculations
[1702]	947	! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
	948	! a more robust solution will be given in next release
[2528]	949	ipk = SIZE(sd%fnow, 3)
	950	ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
	951	IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
[1275]	952
	953	nxt_wgt = nxt_wgt + 1
	954
	955	ELSE
	956	CALL ctl_stop( ' fld_weight : unable to read the file ' )
	957	ENDIF
	958
[2590]	959	IF( (.NOT. wrk_release(2, 1)) .OR. (.NOT. iwrk_release(2,1)) )THEN
	960	CALL ctl_stop('fld_weights: failed to release workspace arrays.')
	961	END IF
	962
[1275]	963	END SUBROUTINE fld_weight
	964
[2528]	965	SUBROUTINE fld_interp(num, clvar, kw, kk, dta, nrec)
[1275]	966	!!---------------------------------------------------------------------
	967	!! * ROUTINE fld_interp *
	968	!!
	969	!! ** Purpose : apply weights to input gridded data to create data
	970	!! on model grid
	971	!!
	972	!! ** Method :
	973	!!----------------------------------------------------------------------
	974	INTEGER, INTENT(in) :: num ! stream number
	975	CHARACTER(LEN=*), INTENT(in) :: clvar ! variable name
	976	INTEGER, INTENT(in) :: kw ! weights number
[2528]	977	INTEGER, INTENT(in) :: kk ! vertical dimension of kk
[2590]	978	REAL(wp), INTENT(inout), DIMENSION(:,:,:) :: dta ! output field on model grid
[1275]	979	INTEGER, INTENT(in) :: nrec ! record number to read (ie time slice)
	980	!!
[2528]	981	INTEGER, DIMENSION(3) :: rec1,recn ! temporary arrays for start and length
[1275]	982	INTEGER :: jk, jn, jm ! loop counters
	983	INTEGER :: ni, nj ! lengths
	984	INTEGER :: jpimin,jpiwid ! temporary indices
	985	INTEGER :: jpjmin,jpjwid ! temporary indices
	986	INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices
	987	!!----------------------------------------------------------------------
	988	!
	989
	990	!! for weighted interpolation we have weights at four corners of a box surrounding
	991	!! a model grid point, each weight is multiplied by a grid value (bilinear case)
	992	!! or by a grid value and gradients at the corner point (bicubic case)
	993	!! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
	994
[2528]	995	!! sub grid from non-model input grid which encloses all grid points in this nemo process
[1275]	996	jpimin = ref_wgts(kw)%botleft(1)
	997	jpjmin = ref_wgts(kw)%botleft(2)
	998	jpiwid = ref_wgts(kw)%jpiwgt
	999	jpjwid = ref_wgts(kw)%jpjwgt
	1000
[2528]	1001	!! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients
[1275]	1002	rec1(1) = MAX( jpimin-1, 1 )
	1003	rec1(2) = MAX( jpjmin-1, 1 )
[2528]	1004	rec1(3) = 1
[1275]	1005	recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
	1006	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
[2528]	1007	recn(3) = kk
[1275]	1008
[2528]	1009	!! where we need to put it in the non-nemo grid fly_dta
	1010	!! note that jpi1 and jpj1 only differ from 1 when jpimin and jpjmin are 1
	1011	!! (ie at the extreme west or south of the whole input grid) and similarly for jpi2 and jpj2
[1275]	1012	jpi1 = 2 + rec1(1) - jpimin
	1013	jpj1 = 2 + rec1(2) - jpjmin
	1014	jpi2 = jpi1 + recn(1) - 1
	1015	jpj2 = jpj1 + recn(2) - 1
	1016
[2528]	1017	ref_wgts(kw)%fly_dta(:,:,:) = 0.0
	1018	SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) )
	1019	CASE(1)
	1020	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn)
	1021	CASE DEFAULT
	1022	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn)
	1023	END SELECT
[1275]	1024
	1025	!! first four weights common to both bilinear and bicubic
[2528]	1026	!! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft
[1275]	1027	!! note that we have to offset by 1 into fly_dta array because of halo
[2528]	1028	dta(:,:,:) = 0.0
[1275]	1029	DO jk = 1,4
	1030	DO jn = 1, jpj
	1031	DO jm = 1,jpi
	1032	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
	1033	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
[2528]	1034	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,:)
[1275]	1035	END DO
	1036	END DO
	1037	END DO
	1038
	1039	IF (ref_wgts(kw)%numwgt .EQ. 16) THEN
	1040
	1041	!! fix up halo points that we couldnt read from file
	1042	IF( jpi1 == 2 ) THEN
[2528]	1043	ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:)
[1275]	1044	ENDIF
	1045	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
[2528]	1046	ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:)
[1275]	1047	ENDIF
	1048	IF( jpj1 == 2 ) THEN
[2528]	1049	ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:)
[1275]	1050	ENDIF
	1051	IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN
[2528]	1052	ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:)
[1275]	1053	ENDIF
	1054
	1055	!! if data grid is cyclic we can do better on east-west edges
	1056	!! but have to allow for whether first and last columns are coincident
	1057	IF( ref_wgts(kw)%cyclic ) THEN
	1058	rec1(2) = MAX( jpjmin-1, 1 )
[2528]	1059	recn(1) = 1
[1275]	1060	recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
	1061	jpj1 = 2 + rec1(2) - jpjmin
	1062	jpj2 = jpj1 + recn(2) - 1
	1063	IF( jpi1 == 2 ) THEN
[2528]	1064	rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap
	1065	SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
	1066	CASE(1)
	1067	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
	1068	CASE DEFAULT
	1069	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
	1070	END SELECT
	1071	ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
[1275]	1072	ENDIF
	1073	IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
[2528]	1074	rec1(1) = 1 + ref_wgts(kw)%overlap
	1075	SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) )
	1076	CASE(1)
	1077	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn)
	1078	CASE DEFAULT
	1079	CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn)
	1080	END SELECT
	1081	ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
[1275]	1082	ENDIF
	1083	ENDIF
	1084
	1085	! gradient in the i direction
	1086	DO jk = 1,4
	1087	DO jn = 1, jpj
	1088	DO jm = 1,jpi
	1089	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
	1090	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
[2528]	1091	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 * &
	1092	(ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:))
[1275]	1093	END DO
	1094	END DO
	1095	END DO
	1096
	1097	! gradient in the j direction
	1098	DO jk = 1,4
	1099	DO jn = 1, jpj
	1100	DO jm = 1,jpi
	1101	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
	1102	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
[2528]	1103	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 * &
	1104	(ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:))
[1275]	1105	END DO
	1106	END DO
	1107	END DO
	1108
	1109	! gradient in the ij direction
	1110	DO jk = 1,4
	1111	DO jn = 1, jpj
	1112	DO jm = 1,jpi
	1113	ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
	1114	nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
[2528]	1115	dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
	1116	(ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni ,nj+2,:)) - &
	1117	(ref_wgts(kw)%fly_dta(ni+2,nj ,:) - ref_wgts(kw)%fly_dta(ni ,nj ,:)))
[1275]	1118	END DO
	1119	END DO
	1120	END DO
	1121
	1122	END IF
	1123
	1124	END SUBROUTINE fld_interp
[2528]	1125
	1126
	1127	FUNCTION ksec_week( cdday )
	1128	!!---------------------------------------------------------------------
	1129	!! * FUNCTION kshift_week *
	1130	!!
	1131	!! ** Purpose :
	1132	!!
	1133	!! ** Method :
	1134	!!---------------------------------------------------------------------
	1135	CHARACTER(len=*), INTENT(in) :: cdday !3 first letters of the first day of the weekly file
	1136	!!
	1137	INTEGER :: ksec_week ! output variable
	1138	INTEGER :: ijul !temp variable
	1139	INTEGER :: ishift !temp variable
	1140	CHARACTER(len=3),DIMENSION(7) :: cl_week
	1141	!!----------------------------------------------------------------------
	1142	cl_week = (/"sun","sat","fri","thu","wed","tue","mon"/)
	1143	DO ijul = 1, 7
	1144	IF( cl_week(ijul) == TRIM(cdday) ) EXIT
	1145	ENDDO
	1146	IF( ijul .GT. 7 ) CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): '//TRIM(cdday) )
	1147	!
	1148	ishift = ijul * NINT(rday)
	1149	!
	1150	ksec_week = nsec_week + ishift
	1151	ksec_week = MOD( ksec_week, 7*NINT(rday) )
	1152	!
	1153	END FUNCTION ksec_week
	1154
	1155
[888]	1156	END MODULE fldread

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

Download in other formats: