source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 20 months ago

GMED 450 add flush after prints

File size: 50.4 KB
Line 
1MODULE 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_mem*bias_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   !! * Substitutions needed to have the variable fsdept_n
116#  include "domzgr_substitute.h90"
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
135      & fctamp,               &  !: amplification factor for T if inertial
136      & rn_maxlat_bias,       &  !: Max lat for latitudinal ramp
137      & rn_minlat_bias,       &  !: Min lat for latitudinal ramp
138      & zwgt,                 &  !: weight for IPC
139      & ztscale                  !: decay rate for IPC
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
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
152      & sbias_p_out,   &   !: Output salinity bias field for P correction
153      & tbias_i_out,   &   !: Output temperature bias field for incremental P correction
154      & sbias_i_out,   &   !: Output salinity bias field for incremental P correction
155      & tbias_asm_stscale, &   !: Short time scale temperature bias field
156      & sbias_asm_stscale, &   !: Short time scale salinity bias field
157      & tbias_asm_stscale_out, &   !: Short time scale temperature bias output field
158      & sbias_asm_stscale_out   !: Short time scale salinity bias output field
159
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   
164   LOGICAL, PRIVATE :: ln_itdecay      ! evolve bias correction at every time step. 
165   LOGICAL, PRIVATE :: ln_incpc        ! incremental pressure correction plus pressure correction
166   LOGICAL, PRIVATE :: ln_incpc_only        ! incremental pressure correction only
167
168   REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef
169   REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef_stscale
170
171   INTEGER, PRIVATE ::  &
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
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   
187CONTAINS
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,     &
219         &                   lenscl_bias, &  ! lengthscale pressure bias decay between minlat and maxlat.
220         &                   minlat_bias, &  ! used in ipc
221         &                   maxlat_bias     ! used in ipc
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,    &
226         & eft_rlx, t_rlx_upd, eft_asm, t_asm_upd, nn_lat_ramp,           &
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,  &
229         & ln_bsyncro, fctamp, rn_maxlat_bias, rn_minlat_bias,            &
230         & nn_bias_itwrt, ln_itdecay, ln_incpc,ln_incpc_only, ztscale, zwgt
231 
232
233      !-----------------------------------------------------------------------
234      ! Read Namelist : bias interface
235      !-----------------------------------------------------------------------
236
237
238      REWIND( numnam_ref )              ! Namelist nambias in reference namelist : Bias pressure correction
239      READ  ( numnam_ref, nambias, IOSTAT = ios, ERR = 901)
240901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambias in reference namelist', lwp )
241
242
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     
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
253
254      REWIND( numnam_cfg )              ! Namelist nambias in configuration namelist : Bias pressure correction
255      READ  ( numnam_cfg, nambias, IOSTAT = ios, ERR = 902 )
256902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambias in configuration namelist', lwp )
257      IF(lwm .AND. nprint > 2) 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
265     
266      ! set up decay scales
267      log2           = LOG( 2.0_wp )
268      t_rlx_mem      = EXP( - log2 * rdt / ( eft_rlx * rday ) )
269      t_asm_mem      = EXP( - log2 * bias_time_unit_asm/ ( eft_asm * rday ) )
270     
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
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
289    WRITE(numout,*) '     incremental press. correction  ln_incpc       = ',ln_incpc
290    WRITE(numout,*) '     incremental press. correction only  ln_incpc_only       = ',ln_incpc_only
291    WRITE(numout,*) '     incremental press. correction  zwgt           = ',zwgt
292    WRITE(numout,*) '     incremental press. correction  ztscale        = ',ztscale
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
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
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' )
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 
322         IF(lflush) call flush(numout)
323      ENDIF
324      IF( .NOT. ln_bias ) RETURN
325
326      IF( .NOT. lalloc ) THEN
327
328         ALLOCATE( tbias(jpi,jpj,jpk)  , &
329            &      sbias(jpi,jpj,jpk)  , &
330            &      tbias_p(jpi,jpj,jpk), &
331            &      sbias_p(jpi,jpj,jpk), &
332            &      tbias_i(jpi,jpj,jpk), &
333            &      sbias_i(jpi,jpj,jpk), &
334            &      rhd_pc(jpi,jpj,jpk) , &
335            &      gru_pc(jpi,jpj)     , &
336            &      grv_pc(jpi,jpj)       )
337
338         ALLOCATE( fbcoef(jpi,jpj), fbcoef_stscale(jpi,jpj) )
339         
340         IF( ln_bias_asm ) ALLOCATE(  tbias_asm(jpi,jpj,jpk),     &
341            &                         sbias_asm(jpi,jpj,jpk),     &
342                                      tbias_asm_out(jpi,jpj,jpk), & 
343                                      sbias_asm_out(jpi,jpj,jpk), & 
344                                      tbias_p_out(jpi,jpj,jpk),   &             
345                                      sbias_p_out(jpi,jpj,jpk)    )
346
347         IF( ln_bias_rlx ) ALLOCATE(  tbias_rlx(jpi,jpj,jpk),     &
348            &                         sbias_rlx(jpi,jpj,jpk),     &
349                                      tbias_rlx_out(jpi,jpj,jpk), & 
350                                      sbias_rlx_out(jpi,jpj,jpk)  )
351
352         IF( ln_incpc )    ALLOCATE(  tbias_asm_stscale(jpi,jpj,jpk), &
353            &                         sbias_asm_stscale(jpi,jpj,jpk), &
354                                      tbias_asm_stscale_out(jpi,jpj,jpk), &
355                                      sbias_asm_stscale_out(jpi,jpj,jpk), &
356                                      tbias_i_out(jpi,jpj,jpk),   & 
357                                      sbias_i_out(jpi,jpj,jpk)   )
358
359         lalloc = .TRUE.
360
361      ENDIF
362
363      IF( ln_bias_ofl ) THEN      ! set sf_tbias_ofl and sf_sbias_ofl strctrs
364         !
365         ! tbias
366         !
367         ALLOCATE( sf_tbias_ofl(1), STAT=ierror )
368         IF( ierror > 0 ) THEN
369            CALL ctl_stop( 'bias_init: unable to allocate sf_tbias_ofl structure' )   ;    RETURN
370         ENDIF
371         ALLOCATE( sf_tbias_ofl(1)%fnow(jpi,jpj,jpk) )
372         ALLOCATE( sf_tbias_ofl(1)%fdta(jpi,jpj,jpk,2) )
373         
374         ! fill structure with values and control print
375         CALL fld_fill( sf_tbias_ofl, (/ sn_tbias_ofl /), cn_dir, 'bias_init', 'Offline T bias term ', 'nam_tbias_ofl' )
376         !
377         ! salinity bias
378         !
379         ALLOCATE( sf_sbias_ofl(1), STAT=ierror )
380         IF( ierror > 0 ) THEN
381            CALL ctl_stop( 'bias_init: unable to allocate sf_sbias_ofl structure' )   ;   RETURN
382         ENDIF
383         ALLOCATE( sf_sbias_ofl(1)%fnow(jpi,jpj,jpk) )
384         ALLOCATE( sf_sbias_ofl(1)%fdta(jpi,jpj,jpk,2) )
385         
386         ! fill structure with values and control print
387         CALL fld_fill( sf_sbias_ofl, (/ sn_sbias_ofl /), cn_dir, 'bias_init', 'Offline S bias term ', 'nam_sbias_ofl' )
388         
389      ENDIF
390
391      ! Read total bias
392      IF ( ln_bias ) THEN
393         tbias(:,:,:)   = 0.0_wp
394         sbias(:,:,:)   = 0.0_wp
395         tbias_p(:,:,:) = 0.0_wp
396         sbias_p(:,:,:) = 0.0_wp
397         tbias_i(:,:,:) = 0.0_wp
398         sbias_i(:,:,:) = 0.0_wp
399         gru_pc(:,:)    = 0.0_wp
400         grv_pc(:,:)    = 0.0_wp
401         
402         IF ( ln_bias_rlx ) THEN
403            tbias_rlx(:,:,:) = 0.0_wp
404            sbias_rlx(:,:,:) = 0.0_wp
405         ENDIF
406
407         IF ( ln_bias_asm ) THEN   !now rlx and asm bias in same file
408            tbias_asm(:,:,:) = 0.0_wp
409            sbias_asm(:,:,:) = 0.0_wp
410            tbias_asm_out(:,:,:) = 0.0_wp
411            sbias_asm_out(:,:,:) = 0.0_wp
412         ENDIF
413
414         IF ( ln_incpc ) THEN   !incr pressure correction
415            tbias_asm_stscale(:,:,:) = 0.0_wp
416            sbias_asm_stscale(:,:,:) = 0.0_wp
417            tbias_asm_stscale_out(:,:,:) = 0.0_wp
418            sbias_asm_stscale_out(:,:,:) = 0.0_wp
419         ENDIF
420
421
422         numbias_tot    = 0
423         ! Get bias from file and prevent fail if the file does not exist
424         IF(lwp) WRITE(numout,*) 'Opening ',TRIM( cn_bias_tot ) 
425         CALL iom_open( cn_bias_tot, numbias_tot, ldstop=.FALSE. )
426         
427         IF ( numbias_tot > 0 ) THEN     
428            ! Could check validity time of bias fields here...
429            ! Get the T and S bias data
430            IF(lwp) WRITE(numout,*) 'Reading bias fields from tot...'
431
432            !Search for bias from relaxation term if needed. Use same file
433            IF ( ln_bias_rlx ) THEN
434               IF(lwp) WRITE(numout,*) 'Reading bias fields for bias rlx from file ',cn_bias_tot
435               IF( iom_varid( numbias_tot, 'tbias_rlx' ) > 0 ) THEN
436                  ! Get the T and S bias data
437                  CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_rlx', tbias_rlx )
438                  CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_rlx', sbias_rlx )
439               ELSE
440                  CALL ctl_stop( 'Bias relaxation variables not found in ',cn_bias_tot )
441               ENDIF
442            ENDIF
443
444
445            !Search for bias from assim term if needed. Use same file
446            IF ( ln_bias_asm .and. .not.ln_bsyncro ) THEN
447               IF(lwp) WRITE(numout,*) 'Reading a-syncro bias fields for bias asm from file ',cn_bias_tot
448               IF( iom_varid( numbias_tot, 'tbias_asm' ) > 0 ) THEN
449                  ! Get the T and S bias data
450                  CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_asm', tbias_asm )
451                  CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_asm', sbias_asm )
452               ELSE
453                  CALL ctl_stop( 'Bias assim variables not found in ',cn_bias_tot )
454               ENDIF
455            ENDIF
456
457            IF ( ln_incpc .and. .not.ln_bsyncro ) THEN
458               IF(lwp) WRITE(numout,*) 'Reading short time scale bias correction fields for bias asm from file ',cn_bias_tot
459               IF( iom_varid( numbias_tot, 'tbias_asm_stscale' ) > 0 ) THEN
460                  ! Get the T and S bias data
461                  CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_asm_stscale', tbias_asm_stscale )
462                  CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_asm_stscale', sbias_asm_stscale )
463               ELSE
464                  CALL ctl_stop( 'Short time scale bias assim variables not found in ',cn_bias_tot )
465               ENDIF
466            ENDIF 
467
468
469
470            ! Close the file
471            CALL iom_close(numbias_tot)
472
473         ELSE
474            IF(lwp) WRITE(numout,*) 'No bias file found so T and S bias fields are set to zero'
475         ENDIF
476
477      ENDIF
478
479     ! for the time being, the bias_asm is read in the same file as
480     ! bias_rlx
481     ! Implications: bias_asm is estimated/evolved in time in the second outer
482     !               loop only, when the assimilation increments are ready.
483     !               bias_asm is kept cte during the first outer loop.
484     !              => Assyncronous bias correction.
485     ! Alternative: Syncronous bias correction:
486     !              bias_asm estimated/evolved in the first outer loop
487     !              with the asm incrments of the previous cycle.
488     !              bias_asm kept cte during the second outer loop.
489     !              Implication: bias_asm should be estimated really in the
490     !              inner loop.
491      IF ( ln_bsyncro ) THEN
492      ! Read bias from assimilation from a separate file
493      IF ( ln_bias_asm ) THEN
494         tbias_asm(:,:,:) = 0.0_wp
495         sbias_asm(:,:,:) = 0.0_wp
496         numbias_asm      = 0
497         ! Get bias from file and prevent fail if the file does not exist
498         IF(lwp) WRITE(numout,*) 'Opening file for syncro assim bias ',TRIM( cn_bias_asm ) 
499         CALL iom_open( cn_bias_asm, numbias_asm, ldstop=.FALSE. )
500         
501         IF ( numbias_asm > 0 ) THEN     
502            ! Could check validity time of bias fields here...
503           
504            ! Get the T and S bias data
505            IF(lwp) WRITE(numout,*) 'Reading syncro bias fields from asm from file ',cn_bias_asm
506            CALL iom_get( numbias_asm, jpdom_autoglo, 'tbias_asm', tbias_asm )
507            CALL iom_get( numbias_asm, jpdom_autoglo, 'sbias_asm', sbias_asm )
508           
509!  this is only applicable if tbias_asm were to be calculated in the inner loop
510            tbias_asm(:,:,:) = tbias_asm(:,:,:) * rdt / bias_time_unit_asm
511            sbias_asm(:,:,:) = sbias_asm(:,:,:) * rdt / bias_time_unit_asm
512           
513            ! Close the file
514            CALL iom_close(numbias_asm)
515           
516         ELSE
517            IF(lwp) WRITE(numout,*) 'No bias file found from asm so T and S bias fields are set to zero'
518         ENDIF
519         
520      ENDIF
521
522      ENDIF
523
524      !latitudinal dependence of partition coeficients. Adhoc
525      IF ( nn_lat_ramp == 1 ) THEN
526         ! Use the inertial ramp.
527         lenscl_bias = ( rn_maxlat_bias - rn_minlat_bias )*2._wp
528         WHERE ( abs( gphit(:,:) ) <= rn_minlat_bias )         
529            fbcoef(:,:) = 0._wp         
530         ELSEWHERE ( abs( gphit(:,:) ) >= rn_maxlat_bias )         
531            fbcoef(:,:) = 1._wp                   
532         ELSEWHERE       
533            fbcoef(:,:) = 1._wp - exp( -( abs( gphit(:,:) ) - rn_minlat_bias ) &
534                           * ( abs( gphit(:,:) ) - rn_minlat_bias ) / lenscl_bias )                         
535         ENDWHERE 
536      ELSEIF ( nn_lat_ramp == 2 ) THEN   
537         ! Use a linear ramp consist with the geostrophic velocity balance ramp in NEMOVAR
538     
539         WHERE ( abs( gphit(:,:) ) <= rn_minlat_bias )
540            fbcoef(:,:) = 0._wp
541         ELSEWHERE ( abs( gphit(:,:) ) >= rn_maxlat_bias ) 
542            fbcoef(:,:) = 1._wp
543         ELSEWHERE
544            fbcoef(:,:) = 1._wp - ((rn_maxlat_bias - abs( gphit(:,:)))/(rn_maxlat_bias - rn_minlat_bias))
545         ENDWHERE
546      ELSE
547         fbcoef(:,:) = 0.0_wp
548         fctamp      = 0.0_wp
549         fbcoef_stscale(:,:) = 0.0_wp
550      ENDIF
551
552
553      IF ( ln_incpc) THEN
554         minlat_bias = 3.0_wp
555         maxlat_bias = 8.0_wp 
556         WHERE ( abs( gphit(:,:) ) <= minlat_bias )
557            fbcoef_stscale(:,:)=0._wp
558         ELSEWHERE ( abs( gphit(:,:) ) >= maxlat_bias ) 
559            fbcoef_stscale(:,:)=1._wp
560         ELSEWHERE
561            fbcoef_stscale(:,:)=1._wp - ((maxlat_bias - abs( gphit(:,:)))/(maxlat_bias-minlat_bias))
562         ENDWHERE   
563      ENDIF
564     
565      IF(lwp .AND. lflush) call flush(numout)
566
567   END SUBROUTINE bias_init
568
569   SUBROUTINE tra_bias ( kt )
570      !!----------------------------------------------------------------------
571      !!                    ***  ROUTINE tra_bias  ***
572      !!
573      !! ** Purpose : Update bias field every time step
574      !!
575      !! ** Method  : add contributions to bias from 3 terms
576      !!
577      !! ** Action  : Bias from assimilation (read in bias_init)
578      !!              Bias from relaxation term is estimated according to
579      !!              the prescribed time evolution of the bias
580      !!              Bias from ofl term is read from external file
581      !!              The difference contributions are added and the partition
582      !!              into direct bias in T/S and pressure perfomed.
583      !!
584      !! History :  09-03  (M. Balmaseda)
585      !!----------------------------------------------------------------------
586      !! called every timestep after dta_sst if ln_bias true.
587
588      IMPLICIT NONE
589
590      !! * Arguments
591      INTEGER, INTENT(IN) ::   kt             ! ocean time-step index
592      !! * Local variables
593      INTEGER             ::   ji,jj,jk, it   ! local loop index
594      REAL(wp)            ::   tsclf          ! time saling factor
595      REAL(wp)            ::   fb_t_asm_max, fb_t_rlx_max, fb_t_ofl_max
596      REAL(wp)            ::   ztfrac, ztsday
597      REAL(wp)            ::   zfrac, zfrac1 ! temporal weights for inst pcbias (names could be changed)
598      REAL(wp)            ::   zdecay        ! used in inst pcorr
599      REAL(wp), DIMENSION(jpi,jpj) :: zcof1, zcof2
600
601      IF ( .NOT. ln_bias ) RETURN
602      fb_t_rlx_max   = MIN(fb_t_rlx*fctamp,1.0_wp)
603      fb_t_asm_max   = MIN(fb_t_asm*fctamp,1.0_wp)
604      fb_t_ofl_max   = MIN(fb_t_ofl*fctamp,1.0_wp)
605
606      tbias(:,:,:)   = 0.0_wp
607      sbias(:,:,:)   = 0.0_wp
608      tbias_p(:,:,:) = 0.0_wp
609      sbias_p(:,:,:) = 0.0_wp
610      tbias_i(:,:,:) = 0.0_wp
611      sbias_i(:,:,:) = 0.0_wp
612
613      IF ( ln_bias_asm ) THEN
614     
615         tsclf = 1
616         IF ( .NOT.ln_bsyncro ) tsclf = rdt / bias_time_unit_asm 
617         zcof1(:,:) = tsclf * ( ( 1.0_wp - fbcoef(:,:) ) * fb_t_asm + &
618            &                              fbcoef(:,:)   * fb_t_asm_max )
619         zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_asm
620
621         IF ( ln_itdecay ) THEN   !decay the pressure correction at each time step
622   
623       ztsday  = rday / REAL(rdt,wp)
624
625            zdecay = (1-ztscale)**(1/REAL(ztsday,wp)) ! used in ipc
626            zfrac1 = zdecay**REAL(kt,wp) ! used in ipc
627            IF ( zfrac1 <= 0.0 ) zfrac1 = 0.0_wp
628
629            IF( ln_asmiau .and. ln_trainc ) THEN  ! nudge in increments and decay historical contribution
630               
631               IF ( kt <= nitiaufin ) THEN  ! During IAU calculate the fraction of increments to apply at each time step
632
633                  ztfrac = REAL(kt,wp) / REAL(nitiaufin,wp)  ! nudging factor during the IAU
634         
635                  IF (lwp) THEN
636                     WRITE(numout,*) 'tra_bias : bias weights'
637                     WRITE(numout,*) '~~~~~~~~~~~~'
638                     WRITE(numout,* ) "proportion of  increment applied in pcbias ",ztfrac
639                     WRITE(numout,* ) "proportion of  historical pcbias applied ",t_asm_mem**(REAL(kt,wp)/ztsday)
640                  ENDIF
641
642                  DO jk = 1, jpkm1 
643                     tbias(:,:,jk) = tbias(:,:,jk) +                            &
644                     &                ( t_asm_mem**(REAL(kt,wp)/ztsday) * tbias_asm(:,:,jk)  +                    &
645                     &                ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:)
646                     sbias(:,:,jk) = sbias(:,:,jk) +                            &
647                     &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk)  +                     &
648                     &               ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:)
649
650                     tbias_p(:,:,jk) = tbias_p(:,:,jk) +                        &
651                     &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * tbias_asm(:,:,jk)  +                     &
652                     &               ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:)
653                     sbias_p(:,:,jk) = sbias_p(:,:,jk) +                        & 
654                     &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk)  +                     &
655                     &               ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:)
656                  ENDDO
657
658                  IF (ln_incpc) THEN
659
660                     IF (lwp) THEN
661                       WRITE(numout,*) 'tra_bias : bias weights'
662                       WRITE(numout,*) '~~~~~~~~~~~~'
663                       WRITE(numout,* ) "IPC - proportion of  increment applied in pcbias ",ztfrac
664                       WRITE(numout,* ) "IPC - proportion of  historical pcbias applied ",zfrac1
665                     ENDIF
666
667                     DO jk = 1, jpkm1
668
669                        tbias_i(:,:,jk) =  ( t_bkginc(:,:,jk) * tmask(:,:,jk)* zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) )         &
670                        &                + ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
671                        sbias_i(:,:,jk) =  ( s_bkginc(:,:,jk) * tmask(:,:,jk)* zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) )         &
672                        &                + ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
673
674                     ENDDO
675
676                  ENDIF
677
678                  IF ( .not.ln_bsyncro ) THEN
679
680                     IF ( kt == nn_bias_itwrt ) THEN
681
682                        DO jk = 1, jpkm1
683                           tbias_asm_out(:,:,jk) =  t_asm_mem**(REAL(kt,wp)/ztsday) * tbias_asm(:,:,jk)  +       &
684                           &                     ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
685                           sbias_asm_out(:,:,jk) =  t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk)  +       &
686                           &                     ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
687                         END DO
688
689                        IF ( ln_incpc) THEN
690                           DO jk = 1, jpkm1
691                              tbias_asm_stscale_out(:,:,jk) = ( t_bkginc(:,:,jk) * tmask(:,:,jk) *  zwgt * ztfrac ) + ( tbias_asm_stscale(:,:,jk) * zfrac1 )
692                              sbias_asm_stscale_out(:,:,jk) = ( s_bkginc(:,:,jk) * tmask(:,:,jk) *  zwgt * ztfrac ) + ( sbias_asm_stscale(:,:,jk) * zfrac1 )
693                           ENDDO
694                        ENDIF
695
696                     ENDIF
697
698                  ENDIF
699
700                  ! update the historical component with the increments at the end of the IAU
701                  IF ( kt == nitiaufin ) THEN
702                     DO jk = 1, jpkm1
703                        tbias_asm(:,:,jk) =  t_asm_mem**(REAL(kt,wp)/ztsday) * tbias_asm(:,:,jk)  +       &
704                        &                     ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
705                        sbias_asm(:,:,jk) =  t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk)  +       &
706                        &                     ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
707                     END DO
708
709                     IF (ln_incpc) THEN
710                        DO jk = 1, jpkm1
711                           tbias_asm_stscale(:,:,jk) = ( t_bkginc(:,:,jk) * tmask(:,:,jk) * zwgt * ztfrac ) + ( tbias_asm_stscale(:,:,jk) * zfrac1 )
712                           sbias_asm_stscale(:,:,jk) = ( s_bkginc(:,:,jk) * tmask(:,:,jk) * zwgt * ztfrac ) + ( sbias_asm_stscale(:,:,jk) * zfrac1 )
713                        ENDDO
714                     ENDIF
715     
716                  ENDIF
717               
718               ELSE ! decay pressure correction from combined historical component and increments after IAU
719
720                  ztfrac=t_asm_mem**(REAL(kt-nitiaufin,wp)/REAL(nitiaufin,wp))  ! decay from end of IAU
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                  ENDDO
732
733                 IF (ln_incpc) THEN
734
735                   zfrac  = zdecay**REAL(kt-nitiaufin,wp) 
736                   IF ( zfrac <= 0.0 ) zfrac = 0.0_wp
737   
738                   DO jk = 1, jpkm1
739                      tbias_i(:,:,jk) =  tbias_asm_stscale(:,:,jk) * zfrac * (1.0 - fbcoef_stscale(:,:))
740                      sbias_i(:,:,jk) =  sbias_asm_stscale(:,:,jk) * zfrac * (1.0 - fbcoef_stscale(:,:))
741                   ENDDO
742
743                   IF (lwp) THEN
744                      WRITE(numout,*) 'tra_bias : bias weights'
745                      WRITE(numout,*) '~~~~~~~~~~~~'
746                      WRITE(numout,* ) "IPC - proportion of increments + historical pcbias applied ",zfrac
747                   ENDIF
748
749                 ENDIF
750
751                 IF (lwp) THEN
752                    WRITE(numout,*) 'tra_bias : bias weights'
753                    WRITE(numout,*) '~~~~~~~~~~~~'
754                    WRITE(numout,* ) "proportion of increments + historical pcbias applied ",ztfrac
755                 ENDIF
756
757                 IF ( ln_trainc .and. .not.ln_bsyncro ) THEN
758                    IF ( kt == nn_bias_itwrt ) THEN
759                       DO jk = 1, jpkm1
760                          tbias_asm_out(:,:,jk) =  ztfrac * tbias_asm(:,:,jk) 
761                          sbias_asm_out(:,:,jk) =  ztfrac * sbias_asm(:,:,jk) 
762                       END DO
763
764                       IF ( ln_incpc ) THEN
765                          IF (lwp) WRITE(numout,*) 'after end of IAU - IPC - saving s/tbias_asm_stscale'
766                             DO jk = 1, jpkm1
767                                tbias_asm_stscale_out(:,:,jk) =  tbias_asm_stscale(:,:,jk) * zfrac
768                                sbias_asm_stscale_out(:,:,jk) =  sbias_asm_stscale(:,:,jk) * zfrac
769                             ENDDO
770                          ENDIF
771                       ENDIF
772                    ENDIF
773                 ENDIF
774
775            ELSE ! no assimilation increments, simply decay pressure correction (e.g for forecasts)
776                 ! this occurs at obsoper step as well ( ln_asmiau is false )
777
778               DO jk = 1, jpkm1
779                  tbias(:,:,jk) = tbias(:,:,jk) +                                                         &
780                  &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * tbias_asm(:,:,jk) ) * zcof1(:,:)
781                  sbias(:,:,jk) = sbias(:,:,jk) +                                                         &
782                  &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk) ) * zcof1(:,:)
783                  tbias_p(:,:,jk) = tbias_p(:,:,jk) +                                                     &
784                  &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * tbias_asm(:,:,jk) ) * zcof2(:,:)
785                  sbias_p(:,:,jk) = sbias_p(:,:,jk) +                                                     & 
786                  &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk) ) * zcof2(:,:)
787               ENDDO
788
789               IF (lwp) THEN
790                  WRITE(numout,*) 'tra_bias : bias weights'
791                  WRITE(numout,*) '~~~~~~~~~~~~'
792                  WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(REAL(kt,wp)/ztsday)
793               ENDIF
794
795               IF (ln_incpc) THEN
796                   IF (lwp) WRITE(numout,*) 'obsoper or forecast mode - IPC - computing tbias_i and sbias_i'
797                   DO jk = 1, jpkm1
798                      tbias_i(:,:,jk) = tbias_i(:,:,jk) + ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
799                      sbias_i(:,:,jk) = sbias_i(:,:,jk) + ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
800                   ENDDO
801                   tbias_i(:,:,:) =tbias_i(:,:,:)*tmask(:,:,:)
802                   sbias_i(:,:,:) =sbias_i(:,:,:)*tmask(:,:,:)
803               ENDIF
804            ENDIF
805 
806         ELSE ! maintain a constant pressure correction 
807
808            DO jk = 1, jpkm1
809               tbias(:,:,jk) = tbias(:,:,jk) + tbias_asm(:,:,jk) * zcof1(:,:)
810               sbias(:,:,jk) = sbias(:,:,jk) + sbias_asm(:,:,jk) * zcof1(:,:)
811               tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_asm(:,:,jk) * zcof2(:,:)
812               sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_asm(:,:,jk) * zcof2(:,:)
813            END DO     
814
815            IF( ln_asmiau .and. ln_trainc .and. .not.ln_bsyncro ) THEN   
816            ! if last outer loop (ln_asmiau=true and ln_trainc=true). t/sbias_asm
817            ! is updated, only once (end of run) taking into account units.
818               IF ( kt == nn_bias_itwrt ) THEN
819                 IF(lwp) WRITE(numout,*)' estimating asm bias at timestep: ',kt
820                 DO jk = 1, jpkm1
821                   tbias_asm_out(:,:,jk) = t_asm_mem * tbias_asm(:,:,jk)  +             &
822                   &                      t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
823                   sbias_asm_out(:,:,jk) = t_asm_mem * sbias_asm(:,:,jk) +             &
824                   &                      t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)               
825                 END DO
826               ENDIF
827             ENDIF
828 
829         ENDIF
830
831      ENDIF
832
833
834#if   defined key_tradmp 
835      ! Time evolution of bias from relaxation
836      IF ( ln_bias_rlx ) THEN
837         tbias_rlx(:,:,:) = t_rlx_mem * tbias_rlx(:,:,:) + &
838            &               t_rlx_upd * ttrdmp(:,:,:) * rdt
839         sbias_rlx(:,:,:) = t_rlx_mem * sbias_rlx(:,:,:) + &
840            &               t_rlx_upd * strdmp(:,:,:) * rdt
841         zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_rlx +      &
842            &                    fbcoef(:,:)   * fb_t_rlx_max
843         zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_rlx
844         DO jk = 1, jpkm1
845            tbias(:,:,jk)   = tbias(:,:,jk) + tbias_rlx(:,:,jk) * zcof1(:,:)
846            sbias(:,:,jk)   = sbias(:,:,jk) + sbias_rlx(:,:,jk) * zcof1(:,:)
847            tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_rlx(:,:,jk) * zcof2(:,:)
848            sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_rlx(:,:,jk) * zcof2(:,:)
849         ENDDO
850         IF ( kt == nn_bias_itwrt ) THEN
851            tbias_rlx_out(:,:,:) = tbias_rlx(:,:,:)
852            sbias_rlx_out(:,:,:) = sbias_rlx(:,:,:)
853         ENDIF
854      ENDIF
855#endif
856      ! ofline bias
857      IF ( kt == nit000 ) THEN
858         IF(lwp) WRITE(numout,*) ' tra_bias: ln_bias_ofl = ',ln_bias_ofl
859         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~'
860      ENDIF
861      IF ( ln_bias_ofl ) THEN
862         IF(lwp) WRITE(numout,*) 'reading offline bias'
863         CALL fld_read( kt, 1, sf_tbias_ofl ) 
864         CALL fld_read( kt, 1, sf_sbias_ofl ) 
865
866         zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_ofl +           &
867            &                    fbcoef(:,:)   * fb_t_ofl_max
868         zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_ofl
869         DO jk = 1, jpkm1
870            tbias(:,:,jk)   = tbias(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:)
871            sbias(:,:,jk)   = sbias(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:)
872            tbias_p(:,:,jk) = tbias_p(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:)
873            sbias_p(:,:,jk) = sbias_p(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:)
874         ENDDO
875      ENDIF
876
877
878      !apply bias on tracers if needed     
879      IF ( ln_bias_ts_app ) THEN
880         
881         ! Add the bias directely to T/s
882         tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) + tmask(:,:,:) * tbias(:,:,:) / rdt
883         tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) + tmask(:,:,:) * sbias(:,:,:) / rdt
884
885         ! lateral boundary conditions (is this needed?)
886         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1.0_wp )
887         CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1.0_wp )   
888
889      ENDIF
890
891      IF ( kt == nn_bias_itwrt ) THEN
892            tbias_p_out(:,:,:)   = tbias_p(:,:,:)
893            sbias_p_out(:,:,:)   = sbias_p(:,:,:)
894            IF (ln_incpc) THEN
895               tbias_i_out(:,:,:)   = tbias_i(:,:,:)
896               sbias_i_out(:,:,:)   = sbias_i(:,:,:)           
897            ENDIF
898      ENDIF
899      IF(lwp .AND. lflush) call flush(numout)
900   END SUBROUTINE tra_bias
901
902   SUBROUTINE dyn_bias( kt )
903      !!----------------------------------------------------------------------
904      !!                   ***  ROUTINE dyn_bias  ***
905      !!
906      !! ** Purpose :   Computes rhd_pc, gru/v_pc bias corrected
907      !!                for hydrostatic pressure gradient
908      !!                depending on time step (semi-implicit/centered)
909      !!                If partial steps computes bottom pressure gradient.
910      !!                These correction terms will affect only the dynamical
911      !!                component (pressure gradient calculation), but not
912      !!                the isopycnal calculation for the lateral mixing.
913      !!
914      !! ** Method  :   At this stage of the computation, ta and sa are the
915      !!                after temperature and salinity. If semi-implicit, these
916      !!                are used to compute rho and bottom pressure gradient.
917      !!                If centered, tb,sb are used instead.
918      !!                If bias key is activated, the temperature,salinity are
919      !!                bias corrected in the calculation of the density fields
920      !!                used in the pressure gradient calculation.
921      !!
922      !!
923      !! ** Action  : - rhd_pc ready. rhop will be overwriten later
924      !!              - if ln_zps, bottom density gradients gru/v_pc ready.
925      !!----------------------------------------------------------------------
926      !!
927      !! * Arguments
928      INTEGER, INTENT(IN) ::   kt    ! ocean time-step index
929      !! * Local variables
930      REAL(wp) :: tsw(jpi,jpj,jpk,jpts)
931      !!
932      !!----------------------------------------------------------------------
933      !
934      ! gtu,gsu,gtv,gsv rhop will be overwritten later in step.
935      !
936      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg
937         tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:)
938         tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:)
939         IF ( ln_incpc ) THEN
940            tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:) - tbias_i(:,:,:)
941            tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:) - sbias_i(:,:,:)
942            IF ( ln_incpc_only ) THEN
943               IF(lwp) WRITE(numout,*) 'ln_incpc_only =', ln_incpc_only, 'tsw updated with IPC only and ln_dynhpg_imp = ',ln_dynhpg_imp
944               tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_i(:,:,:)
945               tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_i(:,:,:)
946            ENDIF
947         ENDIF
948      ELSE
949         tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:)
950         tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:)
951         IF ( ln_incpc ) THEN
952            tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:) - tbias_i(:,:,:)
953            tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:) - sbias_i(:,:,:)
954            IF ( ln_incpc_only ) THEN
955               IF(lwp) WRITE(numout,*) 'ln_incpc_only =', ln_incpc_only, 'tsw updated with IPC only and ln_dynhpg_imp = ',ln_dynhpg_imp
956               tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_i(:,:,:)
957               tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_i(:,:,:)
958            ENDIF
959         ENDIF
960      ENDIF
961
962      IF(lwp) WRITE(numout,*) 'dyn_bias( kt ) calculating rhd_pc, kt =', kt
963      CALL eos( tsw, rhd_pc, rhop, fsdept_n(:,:,:) )
964     
965      CALL lbc_lnk( rhd_pc, 'T', 1.0_wp )
966
967
968      ! Partial steps: now horizontal gradient of t,s,rd
969      ! at the bottom ocean level
970
971      IF( ln_zps    )  THEN
972         CALL zps_hde( kt, jpts, tsw, gtsu, gtsv,  &
973            &          rhd_pc, gru_pc , grv_pc  )
974      ENDIF
975
976      IF(lwp .AND. lflush) call flush(numout)
977
978   END SUBROUTINE dyn_bias
979   
980   SUBROUTINE bias_opn( kt )
981      !!---------------------------------------------------------------------
982      !!                   ***  ROUTINE bias_opn  ***
983      !!                     
984      !! ** Purpose :  open bias restart file following the same logic as the
985      !!               standard restarts.
986      !!----------------------------------------------------------------------
987      !! * Arguments
988      INTEGER, INTENT(IN) ::   kt     ! ocean time-step
989      !! * Local variables
990      CHARACTER(LEN=20)   ::   clbkt    ! ocean time-step deine as a character
991      CHARACTER(LEN=50)   ::   clbias_tot   ! total bias restart file name
992      !!----------------------------------------------------------------------
993      !
994      IF( lrst_oce .AND. .NOT.lrst_bias ) THEN       ! create bias file
995         IF( nitend > 999999999 ) THEN   ;   WRITE(clbkt, *       ) nitend
996         ELSE                            ;   WRITE(clbkt, '(i8.8)') nitend
997         ENDIF
998         clbias_tot = TRIM(cexper)//"_"//TRIM(ADJUSTL(clbkt))//"_"//TRIM(cn_bias_tot)
999         IF(lwp) THEN
1000            WRITE(numout,*)
1001            SELECT CASE ( jprstlib )
1002            CASE ( jprstdimg )   ;   WRITE(numout,*) '             open tot bias restart binary file: '//clbias_tot
1003            CASE DEFAULT         ;   WRITE(numout,*) '             open tot bias restart NetCDF file: '//clbias_tot
1004            END SELECT
1005         ENDIF
1006         CALL iom_open( clbias_tot, numbias_tot , ldwrt = .TRUE., kiolib = jprstlib )
1007         lrst_bias=.TRUE.
1008
1009      ENDIF
1010      !
1011      IF(lwp .AND. lflush) call flush(numout)
1012
1013   END SUBROUTINE bias_opn
1014
1015   SUBROUTINE bias_wrt( kt )
1016      !!---------------------------------------------------------------------
1017      !!                   ***  ROUTINE bias_wrt  ***
1018      !!                     
1019      !! ** Purpose :   Write bias restart fields in the format corresponding to jprstlib
1020      !!
1021      !! ** Method  :   Write in numbias_tot when kt == nitend in output
1022      !!                file, save fields which are necessary for restart.
1023      !!
1024      !! Changed the timestep for writing to nitend rather than nitrst as this causes a
1025      !! problem if the nstock list is used to determine the restart output times and
1026      !! means that the bias is not output at all. M. Martin. 08/2011.
1027      !! Need to check with Magdalena that this is ok for ECMWF.
1028      !!----------------------------------------------------------------------
1029      !! * Arguments
1030      INTEGER, INTENT(IN) ::   kt   ! ocean time-step
1031      !!----------------------------------------------------------------------
1032      !                                                                     ! the begining of the run [s]
1033
1034      IF ( ln_bias_rlx ) THEN
1035         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_rlx' , tbias_rlx_out )   
1036         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_rlx' , sbias_rlx_out )   
1037      ENDIF
1038     
1039      IF ( ln_bias_asm ) THEN
1040         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_asm' , tbias_asm_out )   
1041         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_asm' , sbias_asm_out )   
1042         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_p'   , tbias_p_out )
1043         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_p'   , sbias_p_out )         
1044      ENDIF
1045
1046      IF ( ln_incpc ) THEN
1047         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_asm_stscale' , tbias_asm_stscale_out )   
1048         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_asm_stscale' , sbias_asm_stscale_out ) 
1049         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_i'   , tbias_i_out )
1050         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_i'   , sbias_i_out )   
1051      ENDIF
1052     
1053      IF( kt == nitend ) THEN
1054         CALL iom_close( numbias_tot )     ! close the restart file (only at last time step)
1055         lrst_bias = .FALSE.
1056      ENDIF
1057      !
1058   END SUBROUTINE bias_wrt
1059
1060END MODULE bias
Note: See TracBrowser for help on using the repository browser.