New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
bias.F90 in branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90 @ 10120

Last change on this file since 10120 was 10120, checked in by charris, 6 years ago

Added in obsoper updates (and commented out some lines in bias.F90 which don't work as intended currently).

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