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 @ 7666

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

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

Download in other formats: