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

Last change on this file since 7660 was 7660, checked in by isabella, 4 years ago

multiply tbias_i by the tmask

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