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

Last change on this file since 7553 was 7553, checked in by isabella, 7 years ago
added switch data type
File size: 49.4 KB

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

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

Download in other formats: