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.

bias.F90 in branches/UKMO/dev_r5518_pcbias_ipc/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_pcbias_ipc/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90 @ 7172

Last change on this file since 7172 was 7172, checked in by isabella, 7 years ago
adding updates to variables tbias_asm_stscale, sbias_asm_stscale
File size: 46.0 KB

Rev	Line
[6279]	1	MODULE bias
	2	!! Is used by OPA and STEP
	3	!!======================================================================
	4	!! * Module bias *
	5	!! Code to estimate and apply bias correction.
	6	!! The bias is in T/S and Pressure. It follows the generalized
	7	!! bias algorithm presented in Balmaseda et al 2007.
	8	!!
	9	!! It can be read from a file offline, estimated online from relaxation
	10	!! terms or from assimilation increments (this latter estimtd in inner loop)
	11	!!
	12	!! The parameters for the time evolution of the bias can be different
	13	!! from the relaxation terms and for the assim increments. Only the
	14	!! parameter for the relaxtion terms are considered here.
	15	!!
	16	!! The offline bias term can contain the seasonal cycle.
	17	!!
	18	!! The time evolution of the bias for relaxtion is estimated as followed
	19	!! bias_rlx(t+1)=t_rlx_membias_rlx(t)+t_rlx_upd(t/s)trdmp.
	20	!!
	21	!! The total bias in T/S is partion between the correction to T/S only
	22	!! (fb_t) and the correction applied to the pressure gradient (fb_p).
	23	!! We impose that (fb_p=1.-fb_t). These factors can be different for the
	24	!! different terms(fb_t_rxl,fb_t_asm,fb_t_ofl)
	25	!!
	26	!! (t/s)bias = fb_t_ofl * (t/s)bias_ofl +
	27	!! fb_t_rlx * (t/s)bias_rlx +
	28	!! fb_t_asm * (t/s)bias_asm
	29	!! (t/s)bias_p =fb_p_ofl * (t/s)bias_ofl+
	30	!! fb_p_rlx * (t/s)bias_rlx_p +
	31	!! fb_p_asm * (t/s)bias_asm_p
	32	!! (t/s)bias is applied directely to correct T and S
	33	!! (t/s)bias_p is used to compute rhd_pc and gru/v_pc
	34	!!
	35	!! Note: the above is an adhoc /simple way of making the partition
	36	!! between bias in pressure and bias in T/S. It would be better
	37	!! if the partition is done at the time of estimating the time
	38	!! evolution of the bias. That would mean doubling the number of
	39	!! 3D arrays.
	40	!!
	41	!! New addtion: (not in Balmaseda et al 2007):
	42	!! A more physical alternative to the partition of the bias can be
	43	!! done in terms of the inertial frequency: when the time scales of
	44	!! adjustment are shorter that >1/f (Equator), the bias correction should
	45	!! be in the the pressure term. Otherwise, it can act directly in T/S.
	46	!! NOTE that the bias correction in the pressure term here (following
	47	!! (Bell et al 2007) is just the "opposite" as the semi-prognostic method
	48	!! in Greatbatch et al 2004.
	49	!! The use of this partition is controlled by ln_inertial=.true.
	50	!!
	51	!!
	52	!! 2009-03 (M.A. Balmaseda ECMWF)
	53	!!======================================================================
	54
	55	!!----------------------------------------------------------------------
	56	!! bias_init : Read in the bias namelist and the bias arrays
	57	!! tra_bias : Apply the bias fields on T/S directly
	58	!! dyn_bias : Compute density correction for dynamic hpg
	59	!! bias_opn : open bias files for restart capabilities
	60	!! bias_wrt : write bias fies " " "
	61	!!----------------------------------------------------------------------
	62	!! * Modules used
	63	USE par_kind, ONLY: &
	64	& wp
	65	USE par_oce, ONLY: &
	66	& jpi, &
	67	& jpj, &
	68	& jpk
	69	USE dom_oce, ONLY: &
	70	& rdt, &
	71	& ln_zps, &
	72	& gphit
	73	USE phycst, ONLY: &
	74	& rday, &
	75	& rad
	76	USE oce, ONLY: &
	77	& tsb, tsn, tsa, &
	78	& rhop, &
	79	& gtsu, gtsv
	80	USE dynhpg, ONLY: &
	81	& ln_dynhpg_imp
	82	USE tradmp
	83	USE dtatsd, ONLY: &
	84	& ln_tsd_tradmp
	85	USE in_out_manager, ONLY: &
	86	& lwp, &
	87	& numnam_ref, &
	88	& numnam_cfg, &
	89	& numond, &
	90	& numout, &
	91	& lrst_oce, &
	92	& nit000
	93	USE iom
	94	USE eosbn2
	95	USE zpshde ! partial step: hor. derivative (zps_hde routine)
	96	USE biaspar
	97	USE fldread ! read input fields
	98	USE lbclnk ! lateral boundary conditions (or mpp link)
	99	USE asmpar
	100	USE asminc
	101	USE lib_mpp, ONLY: &
	102	& ctl_stop, &
	103	& ctl_nam
	104
	105	IMPLICIT NONE
	106
	107	!! * Routine accessibility
	108	PRIVATE
	109	PUBLIC bias_init, & !: Read in the bias namelist and the bias arrays
	110	& tra_bias, & !: Estimate/apply bias on traces
	111	& dyn_bias, & !: " density correction for pressure gradient.
	112	& bias_opn, &
	113	& bias_wrt
	114
	115
	116	!! * Shared variables
	117	!! * Private module variables
	118	REAL(wp), PRIVATE :: &
	119	& bias_time_unit_asm, & !: bias_asm units in s ( per day = 86400 s)
	120	& bias_time_unit_rlx, & !: bias_rlx units in s ( 1 second)
	121	& bias_time_unit_ofl, & !: bias_ofl units in s ( 1 second)
	122	& t_rlx_mem, & !: time param for mem in bias_rlx model
	123	& t_rlx_upd, & !: time param for update in bias_rlx model
	124	!: (pct of incr for computation of bias)
	125	& t_asm_mem, & !: time param for mem in bias_asm model
	126	& t_asm_upd, & !: time param for update in bias_asm model
	127	!: (pct of incr for computation of bias)
	128	& fb_t_rlx, & !: parition of bias in T for rlx bias term
	129	& fb_t_asm, & !: parition of bias in T for asm bias term
	130	& fb_t_ofl, & !: parition of bias in T for ofl bias term
	131	& fb_p_rlx, & !: parition of bias in P for rlx bias term
	132	& fb_p_asm, & !: parition of bias in P for asm bias term
	133	& fb_p_ofl, & !: parition of bias in P for ofl bias term
[6300]	134	& fctamp, & !: amplification factor for T if inertial
[7081]	135	& rn_maxlat_bias, & !: Max lat for latitudinal ramp
	136	& rn_minlat_bias !: Min lat for latitudinal ramp
[6279]	137
	138	LOGICAL, PRIVATE :: lalloc
	139	REAL(wp), PRIVATE, DIMENSION(:,:,:), ALLOCATABLE :: &
	140	& tbias_asm, & !: Temperature bias field
	141	& sbias_asm, & !: Salinity bias field
	142	& tbias_rlx, & !: Temperature bias field
[6300]	143	& sbias_rlx, & !: Salinity bias field
	144	& tbias_asm_out, & !: Output temperature bias field
	145	& sbias_asm_out, & !: Output salinity bias field
	146	& tbias_rlx_out, & !: Output temperature bias field
	147	& sbias_rlx_out, & !: Output salinity bias field
	148	& tbias_p_out, & !: Output temperature bias field for P correction
[7092]	149	& sbias_p_out, & !: Output salinity bias field for P correction
[7154]	150	& tbias_i_out, & !: Output temperature bias field for incremental P correction
	151	& sbias_i_out, & !: Output salinity bias field for incremental P correction
[7092]	152	& tbias_asm_stscale, & !: Short time scale temperature bias field
	153	& sbias_asm_stscale !: Short time scale salinity bias field
[6279]	154
[6300]	155	INTEGER, PRIVATE :: nn_lat_ramp ! choice of latitude dependent ramp
	156	! for the pressure correction.
	157	! 1:inertial ramp, 2:linear ramp, else:no ramp
	158	LOGICAL, PRIVATE :: ln_bsyncro ! syncronous or assincrous bias correction
[7081]	159	LOGICAL, PRIVATE :: ln_itdecay ! evolve bias correction at every time step.
	160	LOGICAL, PRIVATE :: ln_incpc ! incremental pressure correction
[6279]	161
[7081]	162	REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef
	163	REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef_stscale
[6279]	164
	165	INTEGER, PRIVATE :: &
[6300]	166	& numbias_asm, & ! netcdf id of bias file from assim
	167	& numbias_tot, & ! netcdf id of bias file with total bias
	168	& nn_bias_itwrt ! time step for outputting bias pressure corr
[6279]	169
	170	CHARACTER(LEN=128), PRIVATE :: &
	171	& cn_bias_asm, & ! name of bias file from assim
	172	& cn_bias_tot ! name of bias with total/rlx bias
	173
	174	! Structure of input T and S bias offline (file informations, fields read)
	175	TYPE(FLD), PRIVATE, ALLOCATABLE, DIMENSION(:) :: sf_tbias_ofl
	176	TYPE(FLD), PRIVATE, ALLOCATABLE, DIMENSION(:) :: sf_sbias_ofl
	177
	178	TYPE(FLD_N), PRIVATE ::& ! information about the fields to be read
	179	& sn_tbias_ofl, sn_sbias_ofl
	180
	181	CONTAINS
	182
	183	SUBROUTINE bias_init
	184	!!----------------------------------------------------------------------
	185	!! * ROUTINE bias_init *
	186	!!
	187	!! ** Purpose : Read in the bias namelist and read in the bias fields.
	188	!!
	189	!! ** Method : Read in the bias namelist and read in the bias fields.
	190	!!
	191	!! ** Action :
	192	!!
	193	!! History :
	194	!! ! 08-05 (D. Lea) Initial version
	195	!! ! 08-10 (M. Martin) Tidy
	196	!! ! 09-03 (M. Balmaseda). Generalize to estimate the bias
	197	!! from relax and offline bias term.
	198	!! Introduce parameters to control the
	199	!! model for the bias
	200	!! (variables and time evolution)
	201	!!----------------------------------------------------------------------
	202
	203	IMPLICIT NONE
	204
	205	!! * Local declarations
	206
	207	CHARACTER(len=100) :: cn_dir ! dir for location ofline bias
	208	INTEGER :: ierror
	209	INTEGER :: ios ! Local integer output status for namelist read
	210	REAL(wp) :: eft_rlx, & ! efolding time (bias memory) in days
	211	& eft_asm, & ! " "
	212	& log2, &
[7125]	213	& lenscl_bias, & !lengthscale of the pressure bias decay between minlat and maxlat.
	214	& minlat_bias, & !used in ipc
	215	& maxlat_bias !used in ipc
[6279]	216
	217	NAMELIST/nambias/ ln_bias, ln_bias_asm, ln_bias_rlx, ln_bias_ofl, &
	218	& ln_bias_ts_app, ln_bias_pc_app, &
	219	& fb_t_asm, fb_t_rlx, fb_t_ofl, fb_p_asm, fb_p_rlx, fb_p_ofl, &
[6300]	220	& eft_rlx, t_rlx_upd, eft_asm, t_asm_upd, nn_lat_ramp, &
[6279]	221	& bias_time_unit_asm, bias_time_unit_rlx, bias_time_unit_ofl, &
	222	& cn_bias_tot, cn_bias_asm, cn_dir, sn_tbias_ofl, sn_sbias_ofl, &
[6328]	223	& ln_bsyncro, fctamp, rn_maxlat_bias, rn_minlat_bias, &
[7089]	224	& nn_bias_itwrt, ln_itdecay, ln_incpc
[6279]	225
	226
	227	!-----------------------------------------------------------------------
[6328]	228	! Read Namelist : bias interface
[6279]	229	!-----------------------------------------------------------------------
	230
	231
[6328]	232	REWIND( numnam_ref ) ! Namelist nambias in reference namelist : Bias pressure correction
	233	READ ( numnam_ref, nambias, IOSTAT = ios, ERR = 901)
	234	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambias in reference namelist', lwp )
[6279]	235
	236
[6328]	237	! Set additional default values (note that most values are set in the reference namelist)
	238
	239	IF ( ln_asmiau ) nn_bias_itwrt = nitiaufin
	240
[6279]	241	! ... default values (NB: frequency positive => hours, negative => months)
	242	! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or !
	243	! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' !
	244	sn_tbias_ofl = FLD_N( 'tbias_ofl' , -1. , 'tbias' , .TRUE. , .FALSE. , 'yearly', '', '', '' )
	245	sn_sbias_ofl = FLD_N( 'sbias_ofl' , -1. , 'sbias' , .TRUE. , .FALSE. , 'yearly', '', '', '' )
	246
[6328]	247
[6279]	248	REWIND( numnam_cfg ) ! Namelist nambias in configuration namelist : Bias pressure correction
	249	READ ( numnam_cfg, nambias, IOSTAT = ios, ERR = 902 )
	250	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambias in configuration namelist', lwp )
	251	IF(lwm) WRITE ( numond, nambias )
	252
	253
	254	IF ( ( .NOT. ln_bias_asm ) .AND. ( .NOT. ln_bias_ofl ) .AND. ( .NOT. ln_bias_rlx ) ) THEN
	255	ln_bias_ts_app = .FALSE.
	256	ln_bias_pc_app = .FALSE.
	257	ln_bias = .FALSE.
	258	ENDIF
[6328]	259
	260	! set up decay scales
	261	log2 = LOG( 2.0_wp )
[6279]	262	t_rlx_mem = EXP( - log2 * rdt / ( eft_rlx * rday ) )
	263	t_asm_mem = EXP( - log2 * bias_time_unit_asm/ ( eft_asm * rday ) )
[6328]	264
[6279]	265	! Control print
	266	IF(lwp) THEN
	267	WRITE(numout,*)
	268	WRITE(numout,*) 'bias_init : '
	269	WRITE(numout,*) '~~~~~~~~~~~ '
	270	WRITE(numout,*) ' Namelist nambias : '
	271
	272	WRITE(numout,*) ' Bias switches/options/variables '
	273	WRITE(numout,*) ' bias main switch ln_bias = ',ln_bias
	274	WRITE(numout,*) ' bias from assim ln_bias_asm = ',ln_bias_asm
	275	WRITE(numout,*) ' bias from relax ln_bias_rlx = ',ln_bias_rlx
	276	WRITE(numout,*) ' bias from offln ln_bias_ofl = ',ln_bias_ofl
	277	WRITE(numout,*) ' bias T and S apply ln_bias_ts_app = ',ln_bias_ts_app
	278	WRITE(numout,*) ' bias pressure correctn apply ln_bias_pc_app = ',ln_bias_pc_app
	279	WRITE(numout,*) ' bias pressure correctn apply ln_bias_pc_app = ',ln_bias_pc_app
[6300]	280	WRITE(numout,*) ' lat ramp for bias correction nn_lat_ramp = ',nn_lat_ramp
	281	WRITE(numout,*) ' time step for writing bias fld nn_bias_itwrt = ',nn_bias_itwrt
	282	WRITE(numout,*) ' evolve pcbias at each timestep ln_itdecay = ',ln_itdecay
[7081]	283	WRITE(numout,*) ' incremental press. correction ln_incpc = ',ln_incpc
[6279]	284	WRITE(numout,*) ' Parameters for parition of bias term '
	285	WRITE(numout,*) ' fb_t_rlx = ',fb_t_rlx
	286	WRITE(numout,*) ' fb_t_asm = ',fb_t_asm
	287	WRITE(numout,*) ' fb_t_ofl = ',fb_t_ofl
	288	WRITE(numout,*) ' fb_p_rlx = ',fb_p_rlx
	289	WRITE(numout,*) ' fb_p_asm = ',fb_p_asm
	290	WRITE(numout,*) ' fb_p_ofl = ',fb_p_ofl
	291	WRITE(numout,*) ' Parameters for time evolution of bias '
	292	WRITE(numout,) ' Rlx efolding time (mem) eft_rlx,t_rlx_mem = ', eft_rlx, t_rlx_mem, 1. - log2 rdt / (eft_rlx * rday)
	293	WRITE(numout,*) ' uptdate factor t_rlx_upd = ',t_rlx_upd
	294	WRITE(numout,) ' Asm efolding time (mem) eft_asm,t_asm_mem = ', eft_asm, t_asm_mem, 1. - log2 rdt / (eft_asm * rday)
	295	WRITE(numout,*) ' uptdate factor t_asm_upd = ',t_asm_upd
	296	WRITE(numout,*) ' Filenames and input structures'
	297	WRITE(numout,*) ' bias_tot filename cn_bias_to = ',cn_bias_tot
	298	WRITE(numout,*) ' bias_asm filename cn_bias_asm = ',cn_bias_asm
	299	WRITE(numout,*) ' bias_asm time unit (secs) bias_time_unit_asm = ',bias_time_unit_asm
	300	WRITE(numout,*) ' structure Tem bias ofl sn_tbias_ofl = ',sn_tbias_ofl
	301	WRITE(numout,*) ' structure Sal bias ofl sn_sbias_ofl = ',sn_sbias_ofl
	302
	303	IF ( ( (.NOT. ln_tsd_tradmp) .OR. (.NOT. ln_tradmp) ) .AND. ln_bias_rlx ) &
	304	& CALL ctl_stop (' lk_dtatem, lk_dtasal and lk_tradmp need to be true with ln_bias_rlx' )
	305
	306	ENDIF
	307	IF( .NOT. ln_bias ) RETURN
	308
	309	IF( .NOT. lalloc ) THEN
	310
	311	ALLOCATE( tbias(jpi,jpj,jpk) , &
	312	& sbias(jpi,jpj,jpk) , &
	313	& tbias_p(jpi,jpj,jpk), &
	314	& sbias_p(jpi,jpj,jpk), &
[7154]	315	& tbias_i(jpi,jpj,jpk), &
	316	& sbias_i(jpi,jpj,jpk), &
[6279]	317	& rhd_pc(jpi,jpj,jpk) , &
	318	& gru_pc(jpi,jpj) , &
	319	& grv_pc(jpi,jpj) )
	320
[7081]	321	ALLOCATE( fbcoef(jpi,jpj), fbcoef_stscale(jpi,jpj) )
[6279]	322
[6300]	323	IF( ln_bias_asm ) ALLOCATE( tbias_asm(jpi,jpj,jpk), &
	324	& sbias_asm(jpi,jpj,jpk), &
	325	tbias_asm_out(jpi,jpj,jpk), &
	326	sbias_asm_out(jpi,jpj,jpk), &
	327	tbias_p_out(jpi,jpj,jpk), &
[7154]	328	sbias_p_out(jpi,jpj,jpk), &
	329	tbias_i_out(jpi,jpj,jpk), &
	330	sbias_i_out(jpi,jpj,jpk) )
[6279]	331
[6300]	332	IF( ln_bias_rlx ) ALLOCATE( tbias_rlx(jpi,jpj,jpk), &
	333	& sbias_rlx(jpi,jpj,jpk), &
	334	tbias_rlx_out(jpi,jpj,jpk), &
	335	sbias_rlx_out(jpi,jpj,jpk) )
[7081]	336
	337	IF( ln_incpc ) ALLOCATE( tbias_asm_stscale(jpi,jpj,jpk), &
	338	& sbias_asm_stscale(jpi,jpj,jpk))
	339
[6279]	340	lalloc = .TRUE.
	341
	342	ENDIF
	343
	344	IF( ln_bias_ofl ) THEN ! set sf_tbias_ofl and sf_sbias_ofl strctrs
	345	!
	346	! tbias
	347	!
	348	ALLOCATE( sf_tbias_ofl(1), STAT=ierror )
	349	IF( ierror > 0 ) THEN
	350	CALL ctl_stop( 'bias_init: unable to allocate sf_tbias_ofl structure' ) ; RETURN
	351	ENDIF
	352	ALLOCATE( sf_tbias_ofl(1)%fnow(jpi,jpj,jpk) )
	353	ALLOCATE( sf_tbias_ofl(1)%fdta(jpi,jpj,jpk,2) )
	354
	355	! fill structure with values and control print
	356	CALL fld_fill( sf_tbias_ofl, (/ sn_tbias_ofl /), cn_dir, 'bias_init', 'Offline T bias term ', 'nam_tbias_ofl' )
	357	!
	358	! salinity bias
	359	!
	360	ALLOCATE( sf_sbias_ofl(1), STAT=ierror )
	361	IF( ierror > 0 ) THEN
	362	CALL ctl_stop( 'bias_init: unable to allocate sf_sbias_ofl structure' ) ; RETURN
	363	ENDIF
	364	ALLOCATE( sf_sbias_ofl(1)%fnow(jpi,jpj,jpk) )
	365	ALLOCATE( sf_sbias_ofl(1)%fdta(jpi,jpj,jpk,2) )
	366
	367	! fill structure with values and control print
	368	CALL fld_fill( sf_sbias_ofl, (/ sn_sbias_ofl /), cn_dir, 'bias_init', 'Offline S bias term ', 'nam_sbias_ofl' )
	369
	370	ENDIF
	371
	372	! Read total bias
	373	IF ( ln_bias ) THEN
	374	tbias(:,:,:) = 0.0_wp
	375	sbias(:,:,:) = 0.0_wp
	376	tbias_p(:,:,:) = 0.0_wp
	377	sbias_p(:,:,:) = 0.0_wp
[7154]	378	tbias_i(:,:,:) = 0.0_wp
	379	sbias_i(:,:,:) = 0.0_wp
[6279]	380	gru_pc(:,:) = 0.0_wp
	381	grv_pc(:,:) = 0.0_wp
	382	IF ( ln_bias_rlx ) THEN
	383	tbias_rlx(:,:,:) = 0.0_wp
	384	sbias_rlx(:,:,:) = 0.0_wp
	385	ENDIF
	386	IF ( ln_bias_asm ) THEN !now rlx and asm bias in same file
	387	tbias_asm(:,:,:) = 0.0_wp
	388	sbias_asm(:,:,:) = 0.0_wp
	389	ENDIF
[7081]	390
	391	IF ( ln_incpc ) THEN !incr pressure correction
	392	tbias_asm_stscale(:,:,:) = 0.0_wp
	393	sbias_asm_stscale(:,:,:) = 0.0_wp
	394	ENDIF
	395
	396
[6279]	397	numbias_tot = 0
	398	! Get bias from file and prevent fail if the file does not exist
	399	IF(lwp) WRITE(numout,*) 'Opening ',TRIM( cn_bias_tot )
	400	CALL iom_open( cn_bias_tot, numbias_tot, ldstop=.FALSE. )
	401
	402	IF ( numbias_tot > 0 ) THEN
	403	! Could check validity time of bias fields here...
	404	! Get the T and S bias data
	405	IF(lwp) WRITE(numout,*) 'Reading bias fields from tot...'
	406
	407	!Search for bias from relaxation term if needed. Use same file
	408	IF ( ln_bias_rlx ) THEN
	409	IF(lwp) WRITE(numout,*) 'Reading bias fields for bias rlx from file ',cn_bias_tot
	410	IF( iom_varid( numbias_tot, 'tbias_rlx' ) > 0 ) THEN
	411	! Get the T and S bias data
	412	CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_rlx', tbias_rlx )
	413	CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_rlx', sbias_rlx )
	414	ELSE
	415	CALL ctl_stop( 'Bias relaxation variables not found in ',cn_bias_tot )
	416	ENDIF
	417	ENDIF
	418
	419
	420	!Search for bias from assim term if needed. Use same file
	421	IF ( ln_bias_asm .and. .not.ln_bsyncro ) THEN
	422	IF(lwp) WRITE(numout,*) 'Reading a-syncro bias fields for bias asm from file ',cn_bias_tot
	423	IF( iom_varid( numbias_tot, 'tbias_asm' ) > 0 ) THEN
	424	! Get the T and S bias data
	425	CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_asm', tbias_asm )
	426	CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_asm', sbias_asm )
	427	ELSE
	428	CALL ctl_stop( 'Bias assim variables not found in ',cn_bias_tot )
	429	ENDIF
[7081]	430
	431
	432	IF ( ln_incpc ) THEN
	433	IF(lwp) WRITE(numout,*) 'Reading short time scale bias correction fields for bias asm from file ',cn_bias_tot
	434	IF( iom_varid( numbias_tot, 'tbias_asm_stscale' ) > 0 ) THEN
	435	! Get the T and S bias data
	436	CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_asm_stscale', tbias_asm_stscale )
	437	CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_asm_stscale', sbias_asm_stscale )
	438	ELSE
	439	CALL ctl_stop( 'Short time scale bias assim variables not found in ',cn_bias_tot )
	440	ENDIF
	441	ENDIF
	442
[6279]	443	ENDIF
	444
[7081]	445
	446
[6279]	447	! Close the file
	448	CALL iom_close(numbias_tot)
	449
	450	ELSE
	451	IF(lwp) WRITE(numout,*) 'No bias file found so T and S bias fields are set to zero'
	452	ENDIF
	453
	454	ENDIF
	455
	456	! for the time being, the bias_asm is read in the same file as
	457	! bias_rlx
	458	! Implications: bias_asm is estimated/evolved in time in the second outer
	459	! loop only, when the assimilation increments are ready.
	460	! bias_asm is kept cte during the first outer loop.
	461	! => Assyncronous bias correction.
	462	! Alternative: Syncronous bias correction:
	463	! bias_asm estimated/evolved in the first outer loop
	464	! with the asm incrments of the previous cycle.
	465	! bias_asm kept cte during the second outer loop.
	466	! Implication: bias_asm should be estimated really in the
	467	! inner loop.
	468	IF ( ln_bsyncro ) THEN
	469	! Read bias from assimilation from a separate file
	470	IF ( ln_bias_asm ) THEN
	471	tbias_asm(:,:,:) = 0.0_wp
	472	sbias_asm(:,:,:) = 0.0_wp
	473	numbias_asm = 0
	474	! Get bias from file and prevent fail if the file does not exist
	475	IF(lwp) WRITE(numout,*) 'Opening file for syncro assim bias ',TRIM( cn_bias_asm )
	476	CALL iom_open( cn_bias_asm, numbias_asm, ldstop=.FALSE. )
	477
	478	IF ( numbias_asm > 0 ) THEN
	479	! Could check validity time of bias fields here...
	480
	481	! Get the T and S bias data
	482	IF(lwp) WRITE(numout,*) 'Reading syncro bias fields from asm from file ',cn_bias_asm
	483	CALL iom_get( numbias_asm, jpdom_autoglo, 'tbias_asm', tbias_asm )
	484	CALL iom_get( numbias_asm, jpdom_autoglo, 'sbias_asm', sbias_asm )
	485
	486	! this is only applicable if tbias_asm were to be calculated in the inner loop
	487	tbias_asm(:,:,:) = tbias_asm(:,:,:) * rdt / bias_time_unit_asm
	488	sbias_asm(:,:,:) = sbias_asm(:,:,:) * rdt / bias_time_unit_asm
	489
	490	! Close the file
	491	CALL iom_close(numbias_asm)
	492
	493	ELSE
	494	IF(lwp) WRITE(numout,*) 'No bias file found from asm so T and S bias fields are set to zero'
	495	ENDIF
	496
	497	ENDIF
	498
	499	ENDIF
	500
	501	!latitudinal dependence of partition coeficients. Adhoc
[6300]	502	IF ( nn_lat_ramp == 1 ) THEN
[6328]	503	! Use the inertial ramp.
	504	lenscl_bias = ( rn_maxlat_bias - rn_minlat_bias )*2._wp
	505	WHERE ( abs( gphit(:,:) ) <= rn_minlat_bias )
[6279]	506	fbcoef(:,:) = 0._wp
[6328]	507	ELSEWHERE ( abs( gphit(:,:) ) >= rn_maxlat_bias )
[6279]	508	fbcoef(:,:) = 1._wp
	509	ELSEWHERE
[6328]	510	fbcoef(:,:) = 1._wp - exp( -( abs( gphit(:,:) ) - rn_minlat_bias ) &
	511	* ( abs( gphit(:,:) ) - rn_minlat_bias ) / lenscl_bias )
[6300]	512	ENDWHERE
	513	ELSEIF ( nn_lat_ramp == 2 ) THEN
	514	! Use a linear ramp consist with the geostrophic velocity balance ramp in NEMOVAR
	515
[6328]	516	WHERE ( abs( gphit(:,:) ) <= rn_minlat_bias )
[6300]	517	fbcoef(:,:) = 0._wp
[6328]	518	ELSEWHERE ( abs( gphit(:,:) ) >= rn_maxlat_bias )
[6300]	519	fbcoef(:,:) = 1._wp
	520	ELSEWHERE
[6328]	521	fbcoef(:,:) = 1._wp - ((rn_maxlat_bias - abs( gphit(:,:)))/(rn_maxlat_bias - rn_minlat_bias))
[6300]	522	ENDWHERE
[6279]	523	ELSE
	524	fbcoef(:,:) = 0.0_wp
	525	fctamp = 0.0_wp
[7081]	526	fbcoef_stscale(:,:) = 0.0_wp
[6279]	527	ENDIF
	528
[6300]	529
[7089]	530	IF ( ln_incpc) THEN
[7081]	531	! not sure if this should be a special case of nn_lat_ramp == 2
	532	minlat_bias = 3.0_wp
	533	maxlat_bias = 8.0_wp
	534	WHERE ( abs( gphit(:,:) ) <= minlat_bias )
	535	fbcoef_stscale(:,:)=0._wp
	536	ELSEWHERE ( abs( gphit(:,:) ) >= maxlat_bias )
	537	fbcoef_stscale(:,:)=1._wp
	538	ELSEWHERE
	539	fbcoef_stscale(:,:)=1._wp - ((maxlat_bias - abs( gphit(:,:)))/(maxlat_bias-minlat_bias))
	540	ENDWHERE
	541	ENDIF
	542
	543
[6279]	544	END SUBROUTINE bias_init
	545
	546	SUBROUTINE tra_bias ( kt )
	547	!!----------------------------------------------------------------------
	548	!! * ROUTINE tra_bias *
	549	!!
	550	!! ** Purpose : Update bias field every time step
	551	!!
	552	!! ** Method : add contributions to bias from 3 terms
	553	!!
	554	!! ** Action : Bias from assimilation (read in bias_init)
	555	!! Bias from relaxation term is estimated according to
	556	!! the prescribed time evolution of the bias
	557	!! Bias from ofl term is read from external file
	558	!! The difference contributions are added and the partition
	559	!! into direct bias in T/S and pressure perfomed.
	560	!!
	561	!! History : 09-03 (M. Balmaseda)
	562	!!----------------------------------------------------------------------
	563	!! called every timestep after dta_sst if ln_bias true.
	564
	565	IMPLICIT NONE
	566
	567	!! * Arguments
	568	INTEGER, INTENT(IN) :: kt ! ocean time-step index
	569	!! * Local variables
	570	INTEGER :: ji,jj,jk, it ! local loop index
	571	REAL(wp) :: tsclf ! time saling factor
	572	REAL(wp) :: fb_t_asm_max, fb_t_rlx_max, fb_t_ofl_max
[6300]	573	REAL(wp) :: ztfrac, ztsday
[7081]	574	REAL(wp) :: zfrac, zfrac1 ! temporal weights for inst pcbias (names could be changed)
	575	REAL(wp) :: ztscale ! reduce the inst pcbias by this amount per 24 hours
	576	REAL(wp) :: zwgt ! Weight for the inst pcorr term
	577	REAL(wp) :: zdecay ! used in inst pcorr
[6279]	578	REAL(wp), DIMENSION(jpi,jpj) :: zcof1, zcof2
	579
	580	IF ( .NOT. ln_bias ) RETURN
	581	fb_t_rlx_max = MIN(fb_t_rlx*fctamp,1.0_wp)
	582	fb_t_asm_max = MIN(fb_t_asm*fctamp,1.0_wp)
	583	fb_t_ofl_max = MIN(fb_t_ofl*fctamp,1.0_wp)
	584
	585	tbias(:,:,:) = 0.0_wp
	586	sbias(:,:,:) = 0.0_wp
	587	tbias_p(:,:,:) = 0.0_wp
	588	sbias_p(:,:,:) = 0.0_wp
[7154]	589	tbias_i(:,:,:) = 0.0_wp
	590	sbias_i(:,:,:) = 0.0_wp
[7081]	591
	592	ztscale = 0.1_wp
	593	zwgt = 1.0_wp
	594
[6279]	595	IF ( ln_bias_asm ) THEN
[6300]	596
[6279]	597	tsclf = 1
	598	IF ( .NOT.ln_bsyncro ) tsclf = rdt / bias_time_unit_asm
	599	zcof1(:,:) = tsclf * ( ( 1.0_wp - fbcoef(:,:) ) * fb_t_asm + &
	600	& fbcoef(:,:) * fb_t_asm_max )
	601	zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_asm
[6300]	602
	603	IF ( ln_itdecay ) THEN !decay the pressure correction at each time step
	604
	605	ztsday = rday / real(rdt)
	606
[7081]	607	zdecay = (1-ztscale)**(1/real(ztsday)) ! used in ipc
	608	zfrac1 = max(0.0_wp, zdecay**real(kt)) ! used in ipc
	609
	610
[6337]	611	IF( ln_asmiau .and. ln_trainc ) THEN ! nudge in increments and decay historical contribution
[6300]	612
	613
	614	IF ( kt <= nitiaufin ) THEN ! During IAU calculate the fraction of increments to apply at each time step
	615
	616	ztfrac = real(kt) / real(nitiaufin) ! nudging factor during the IAU
[7081]	617
[6300]	618
	619	IF (lwp) THEN
	620	WRITE(numout,*) 'tra_bias : bias weights'
	621	WRITE(numout,*) '~~~~~~~~~~~~'
	622	WRITE(numout,* ) "proportion of increment applied in pcbias ",ztfrac
	623	WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(real(kt)/ztsday)
	624	ENDIF
	625
	626	DO jk = 1, jpkm1
	627	tbias(:,:,jk) = tbias(:,:,jk) + &
	628	& ( t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) + &
	629	& ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:)
	630	sbias(:,:,jk) = sbias(:,:,jk) + &
	631	& ( t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) + &
	632	& ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:)
	633	tbias_p(:,:,jk) = tbias_p(:,:,jk) + &
	634	& ( t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) + &
	635	& ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:)
	636	sbias_p(:,:,jk) = sbias_p(:,:,jk) + &
	637	& ( t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) + &
	638	& ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:)
	639	ENDDO
	640
[7081]	641
	642	IF (ln_incpc) THEN
	643
	644	DO jk = 1, jpkm1
[7154]	645	tbias_i(:,:,jk) = ( t_bkginc(:,:,jk) * zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) ) &
[7081]	646	& - ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
[7154]	647	sbias_i(:,:,jk) = ( s_bkginc(:,:,jk) * zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) ) &
[7081]	648	& - ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
	649	ENDDO
	650
	651	ENDIF
	652
[6337]	653	IF ( .not.ln_bsyncro ) THEN
[6300]	654	IF ( kt == nn_bias_itwrt ) THEN
	655	DO jk = 1, jpkm1
	656	tbias_asm_out(:,:,jk) = t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) + &
	657	& ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
	658	sbias_asm_out(:,:,jk) = t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) + &
	659	& ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
	660	END DO
	661	ENDIF
	662	ENDIF
	663
	664	! update the historical component with the increments at the end of the IAU
	665	IF ( kt == nitiaufin ) THEN
	666	DO jk = 1, jpkm1
	667	tbias_asm(:,:,jk) = t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) + &
	668	& ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
	669	sbias_asm(:,:,jk) = t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) + &
	670	& ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
	671	END DO
[7172]	672
	673	IF (ln_incpc) THEN
	674	DO jk = 1, jpk
	675	tbias_asm_stscale(:,:,jk) = ( t_bkginc(:,:,jk) * zwgt * ztfrac ) + ( tbias_asm_stscale(:,:,jk) * zfrac1 )
	676	sbias_asm_stscale(:,:,jk) = ( s_bkginc(:,:,jk) * zwgt * ztfrac ) + ( sbias_asm_stscale(:,:,jk) * zfrac1 )
	677	ENDDO
	678	ENDIF
[6300]	679	ENDIF
	680
	681	ELSE ! decay pressure correction from combined historical component and increments after IAU
	682
	683	ztfrac=t_asm_mem**(real(kt-nitiaufin)/real(nitiaufin)) ! decay from end of IAU
[7081]	684
	685
[6300]	686	DO jk = 1, jpkm1
	687	tbias(:,:,jk) = tbias(:,:,jk) + &
	688	& ( ztfrac * tbias_asm(:,:,jk) ) * zcof1(:,:)
	689	sbias(:,:,jk) = sbias(:,:,jk) + &
	690	& ( ztfrac * sbias_asm(:,:,jk) ) * zcof1(:,:)
	691	tbias_p(:,:,jk) = tbias_p(:,:,jk) + &
	692	& ( ztfrac * tbias_asm(:,:,jk) ) * zcof2(:,:)
	693	sbias_p(:,:,jk) = sbias_p(:,:,jk) + &
	694	& ( ztfrac * sbias_asm(:,:,jk) ) * zcof2(:,:)
[7081]	695
[6300]	696	ENDDO
	697
[7081]	698	IF (ln_incpc) THEN
	699
	700	zfrac = max(0.0_wp, zdecay**real(kt-nitiaufin) )
	701
	702	DO jk = 1, jpkm1
[7154]	703	tbias_i(:,:,jk) = ( t_bkginc(:,:,jk) * zwgt * zfrac * (1.0 - fbcoef_stscale(:,:)) ) &
[7081]	704	& - ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
[7154]	705	sbias_i(:,:,jk) = ( s_bkginc(:,:,jk) * zwgt * zfrac * (1.0 - fbcoef_stscale(:,:)) ) &
[7081]	706	& - ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
[7089]	707	ENDDO
[7081]	708
	709	ENDIF
	710
[6300]	711	IF (lwp) THEN
	712	WRITE(numout,*) 'tra_bias : bias weights'
	713	WRITE(numout,*) '~~~~~~~~~~~~'
	714	WRITE(numout,* ) "proportion of increments + historical pcbias applied ",ztfrac
	715	ENDIF
	716
	717	IF ( ln_trainc .and. .not.ln_bsyncro ) THEN
	718	IF ( kt == nn_bias_itwrt ) THEN
	719	DO jk = 1, jpkm1
	720	tbias_asm_out(:,:,jk) = ztfrac * tbias_asm(:,:,jk)
	721	sbias_asm_out(:,:,jk) = ztfrac * sbias_asm(:,:,jk)
	722	END DO
	723	ENDIF
	724	ENDIF
	725
	726	ENDIF
	727
	728
	729	ELSE ! no assimilation increments, simply decay pressure correction (e.g for forecasts)
	730
	731	DO jk = 1, jpkm1
	732	tbias(:,:,jk) = tbias(:,:,jk) + &
	733	& ( t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) ) * zcof1(:,:)
	734	sbias(:,:,jk) = sbias(:,:,jk) + &
	735	& ( t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) ) * zcof1(:,:)
	736	tbias_p(:,:,jk) = tbias_p(:,:,jk) + &
	737	& ( t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) ) * zcof2(:,:)
	738	sbias_p(:,:,jk) = sbias_p(:,:,jk) + &
	739	& ( t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) ) * zcof2(:,:)
	740	ENDDO
	741
	742	IF (lwp) THEN
	743	WRITE(numout,*) 'tra_bias : bias weights'
	744	WRITE(numout,*) '~~~~~~~~~~~~'
	745	WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(real(kt)/ztsday)
	746	ENDIF
	747
[6279]	748	ENDIF
[6300]	749
	750	ELSE ! maintain a constant pressure correction
	751
	752	DO jk = 1, jpkm1
	753	tbias(:,:,jk) = tbias(:,:,jk) + tbias_asm(:,:,jk) * zcof1(:,:)
	754	sbias(:,:,jk) = sbias(:,:,jk) + sbias_asm(:,:,jk) * zcof1(:,:)
	755	tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_asm(:,:,jk) * zcof2(:,:)
	756	sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_asm(:,:,jk) * zcof2(:,:)
	757	END DO
	758
[6337]	759	IF( ln_asmiau .and. ln_trainc .and. .not.ln_bsyncro ) THEN
	760	! if last outer loop (ln_asmiau=true and ln_trainc=true). t/sbias_asm
[6300]	761	! is updated, only once (end of run) taking into account units.
[6328]	762	IF ( kt == nn_bias_itwrt ) THEN
	763	IF(lwp) WRITE(numout,*)' estimating asm bias at timestep: ',kt
[6300]	764	DO jk = 1, jpkm1
	765	tbias_asm_out(:,:,jk) = t_asm_mem * tbias_asm(:,:,jk) + &
	766	& t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
	767	sbias_asm_out(:,:,jk) = t_asm_mem * sbias_asm(:,:,jk) + &
	768	& t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
	769	END DO
	770	ENDIF
	771	ENDIF
	772
	773	ENDIF
	774
[6279]	775	ENDIF
	776
	777
	778	#if defined key_tradmp
	779	! Time evolution of bias from relaxation
	780	IF ( ln_bias_rlx ) THEN
	781	tbias_rlx(:,:,:) = t_rlx_mem * tbias_rlx(:,:,:) + &
	782	& t_rlx_upd * ttrdmp(:,:,:) * rdt
	783	sbias_rlx(:,:,:) = t_rlx_mem * sbias_rlx(:,:,:) + &
	784	& t_rlx_upd * strdmp(:,:,:) * rdt
	785	zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_rlx + &
	786	& fbcoef(:,:) * fb_t_rlx_max
	787	zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_rlx
	788	DO jk = 1, jpkm1
	789	tbias(:,:,jk) = tbias(:,:,jk) + tbias_rlx(:,:,jk) * zcof1(:,:)
	790	sbias(:,:,jk) = sbias(:,:,jk) + sbias_rlx(:,:,jk) * zcof1(:,:)
	791	tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_rlx(:,:,jk) * zcof2(:,:)
	792	sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_rlx(:,:,jk) * zcof2(:,:)
	793	ENDDO
[6300]	794	IF ( kt == nn_bias_itwrt ) THEN
	795	tbias_rlx_out(:,:,:) = tbias_rlx(:,:,:)
	796	sbias_rlx_out(:,:,:) = sbias_rlx(:,:,:)
	797	ENDIF
[6279]	798	ENDIF
	799	#endif
	800	! ofline bias
	801	IF ( kt == nit000 ) THEN
	802	IF(lwp) WRITE(numout,*) ' tra_bias: ln_bias_ofl = ',ln_bias_ofl
	803	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~'
	804	ENDIF
	805	IF ( ln_bias_ofl ) THEN
	806	IF(lwp) WRITE(numout,*) 'reading offline bias'
	807	CALL fld_read( kt, 1, sf_tbias_ofl )
	808	CALL fld_read( kt, 1, sf_sbias_ofl )
	809
	810	zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_ofl + &
	811	& fbcoef(:,:) * fb_t_ofl_max
	812	zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_ofl
	813	DO jk = 1, jpkm1
	814	tbias(:,:,jk) = tbias(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:)
	815	sbias(:,:,jk) = sbias(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:)
	816	tbias_p(:,:,jk) = tbias_p(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:)
	817	sbias_p(:,:,jk) = sbias_p(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:)
	818	ENDDO
	819	ENDIF
	820
	821
	822	!apply bias on tracers if needed
	823	IF ( ln_bias_ts_app ) THEN
	824
	825	! Add the bias directely to T/s
	826	tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) + tmask(:,:,:) * tbias(:,:,:) / rdt
	827	tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) + tmask(:,:,:) * sbias(:,:,:) / rdt
	828
	829	! lateral boundary conditions (is this needed?)
	830	CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1.0_wp )
	831	CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1.0_wp )
	832
	833	ENDIF
	834
[6300]	835	IF ( kt == nn_bias_itwrt ) THEN
	836	tbias_p_out(:,:,:) = tbias_p(:,:,:)
	837	sbias_p_out(:,:,:) = sbias_p(:,:,:)
[7154]	838	tbias_i_out(:,:,:) = tbias_i(:,:,:)
	839	sbias_i_out(:,:,:) = sbias_i(:,:,:)
[6300]	840	ENDIF
	841
[6279]	842	END SUBROUTINE tra_bias
	843
	844	SUBROUTINE dyn_bias( kt )
	845	!!----------------------------------------------------------------------
	846	!! * ROUTINE dyn_bias *
	847	!!
	848	!! ** Purpose : Computes rhd_pc, gru/v_pc bias corrected
	849	!! for hydrostatic pressure gradient
	850	!! depending on time step (semi-implicit/centered)
	851	!! If partial steps computes bottom pressure gradient.
	852	!! These correction terms will affect only the dynamical
	853	!! component (pressure gradient calculation), but not
	854	!! the isopycnal calculation for the lateral mixing.
	855	!!
	856	!! ** Method : At this stage of the computation, ta and sa are the
	857	!! after temperature and salinity. If semi-implicit, these
	858	!! are used to compute rho and bottom pressure gradient.
	859	!! If centered, tb,sb are used instead.
	860	!! If bias key is activated, the temperature,salinity are
	861	!! bias corrected in the calculation of the density fields
	862	!! used in the pressure gradient calculation.
	863	!!
	864	!!
	865	!! ** Action : - rhd_pc ready. rhop will be overwriten later
	866	!! - if ln_zps, bottom density gradients gru/v_pc ready.
	867	!!----------------------------------------------------------------------
	868	!!
	869	!! * Arguments
	870	INTEGER, INTENT(IN) :: kt ! ocean time-step index
	871	!! * Local variables
	872	REAL(wp) :: tsw(jpi,jpj,jpk,jpts)
	873	!!
	874	!!----------------------------------------------------------------------
	875	!
	876	! gtu,gsu,gtv,gsv rhop will be overwritten later in step.
	877	!
	878	IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg
	879	tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:)
	880	tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:)
[7154]	881	IF ( ln_incpc ) THEN
[7160]	882	tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:) - tbias_i(:,:,:)
[7154]	883	tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:) - sbias_i(:,:,:)
	884	ENDIF
[6279]	885	ELSE
	886	tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:)
	887	tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:)
[7160]	888	IF ( ln_incpc ) THEN
	889	tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:) - tbias_i(:,:,:)
	890	tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:) - sbias_i(:,:,:)
	891	ENDIF
[6279]	892	ENDIF
	893
	894	CALL eos( tsw, rhd_pc, rhop )
	895
	896	! is this needed?
	897	CALL lbc_lnk( rhd_pc, 'T', 1.0_wp )
	898
	899	! Partial steps: now horizontal gradient of t,s,rd
	900	! at the bottom ocean level
	901	IF( ln_zps ) THEN
	902	CALL zps_hde( kt, jpts, tsw, gtsu, gtsv, &
	903	& rhd_pc, gru_pc , grv_pc )
	904	ENDIF
	905
	906	END SUBROUTINE dyn_bias
	907
	908	SUBROUTINE bias_opn( kt )
	909	!!---------------------------------------------------------------------
	910	!! * ROUTINE bias_opn *
	911	!!
	912	!! ** Purpose : open bias restart file following the same logic as the
	913	!! standard restarts.
	914	!!----------------------------------------------------------------------
	915	!! * Arguments
	916	INTEGER, INTENT(IN) :: kt ! ocean time-step
	917	!! * Local variables
	918	CHARACTER(LEN=20) :: clbkt ! ocean time-step deine as a character
	919	CHARACTER(LEN=50) :: clbias_tot ! total bias restart file name
	920	!!----------------------------------------------------------------------
	921	!
	922	IF( lrst_oce .AND. .NOT.lrst_bias ) THEN ! create bias file
	923	IF( nitend > 999999999 ) THEN ; WRITE(clbkt, * ) nitend
	924	ELSE ; WRITE(clbkt, '(i8.8)') nitend
	925	ENDIF
	926	clbias_tot = TRIM(cexper)//"_"//TRIM(ADJUSTL(clbkt))//"_"//TRIM(cn_bias_tot)
	927	IF(lwp) THEN
	928	WRITE(numout,*)
	929	SELECT CASE ( jprstlib )
	930	CASE ( jprstdimg ) ; WRITE(numout,*) ' open tot bias restart binary file: '//clbias_tot
	931	CASE DEFAULT ; WRITE(numout,*) ' open tot bias restart NetCDF file: '//clbias_tot
	932	END SELECT
	933	ENDIF
	934	CALL iom_open( clbias_tot, numbias_tot , ldwrt = .TRUE., kiolib = jprstlib )
	935	lrst_bias=.TRUE.
	936
	937	ENDIF
	938	!
	939	END SUBROUTINE bias_opn
	940
	941	SUBROUTINE bias_wrt( kt )
	942	!!---------------------------------------------------------------------
	943	!! * ROUTINE bias_wrt *
	944	!!
	945	!! ** Purpose : Write bias restart fields in the format corresponding to jprstlib
	946	!!
	947	!! ** Method : Write in numbias_tot when kt == nitend in output
	948	!! file, save fields which are necessary for restart.
	949	!!
	950	!! Changed the timestep for writing to nitend rather than nitrst as this causes a
	951	!! problem if the nstock list is used to determine the restart output times and
	952	!! means that the bias is not output at all. M. Martin. 08/2011.
	953	!! Need to check with Magdalena that this is ok for ECMWF.
	954	!!----------------------------------------------------------------------
	955	!! * Arguments
	956	INTEGER, INTENT(IN) :: kt ! ocean time-step
	957	!!----------------------------------------------------------------------
	958	! ! the begining of the run [s]
	959
	960	IF ( ln_bias_rlx ) THEN
[6300]	961	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_rlx' , tbias_rlx_out )
	962	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_rlx' , sbias_rlx_out )
[6279]	963	ENDIF
	964
	965	IF ( ln_bias_asm ) THEN
[6300]	966	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_asm' , tbias_asm_out )
	967	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_asm' , sbias_asm_out )
	968	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_p' , tbias_p_out )
[7154]	969	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_p' , sbias_p_out )
	970	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_i' , tbias_i_out )
	971	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_i' , sbias_i_out )
[6279]	972	ENDIF
[7093]	973
	974	IF ( ln_incpc ) THEN
	975	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_asm_stscale' , tbias_asm_stscale )
[7094]	976	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_asm_stscale' , sbias_asm_stscale )
[7093]	977	ENDIF
[6279]	978
	979	IF( kt == nitend ) THEN
	980	CALL iom_close( numbias_tot ) ! close the restart file (only at last time step)
	981	lrst_bias = .FALSE.
	982	ENDIF
	983	!
	984	END SUBROUTINE bias_wrt
	985
	986	END MODULE bias

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

Download in other formats: