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

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

including domzgr_substitute.h90 so that fsdept_n is known to bias

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