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

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

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

Last change on this file since 8444 was 8444, checked in by anaguiar, 7 years ago

Changes to avoid merge conflicts with package branch

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