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_package/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90 @ 8400

Last change on this file since 8400 was 8400, checked in by timgraham, 7 years ago

GMED ticket 335:

  • Merge dev_r5518_GO6_package_inc_asm into package branch to make everything easier for data assimilation
  • No effect on configs without data assimilation
File size: 40.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
139   LOGICAL,  PRIVATE :: lalloc
140   REAL(wp), PRIVATE, DIMENSION(:,:,:), ALLOCATABLE :: &
141      & tbias_asm, &       !: Temperature bias field
142      & sbias_asm, &       !: Salinity bias field
143      & tbias_rlx, &       !: Temperature bias field
144      & sbias_rlx, &       !: Salinity bias field
145      & tbias_asm_out, &   !: Output temperature bias field
146      & sbias_asm_out, &   !: Output salinity bias field
147      & tbias_rlx_out, &   !: Output temperature bias field
148      & sbias_rlx_out, &   !: Output salinity bias field
149      & tbias_p_out,   &   !: Output temperature bias field for P correction
150      & sbias_p_out        !: Output salinity bias field for P correction
151
152   INTEGER, PRIVATE :: nn_lat_ramp     ! choice of latitude dependent ramp
153                                       ! for the pressure correction.
154                   ! 1:inertial ramp, 2:linear ramp, else:no ramp
155   LOGICAL, PRIVATE :: ln_bsyncro      ! syncronous or assincrous bias correction   
156   LOGICAL, PRIVATE :: ln_itdecay      ! evolve bias correction at every time step.   
157
158   REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef 
159
160   INTEGER, PRIVATE ::  &
161      & numbias_asm, &    ! netcdf id of bias file from assim
162      & numbias_tot, &    ! netcdf id of bias file with total bias
163      & nn_bias_itwrt     ! time step for outputting bias pressure corr
164         
165   CHARACTER(LEN=128), PRIVATE :: &
166      & cn_bias_asm,   &  ! name of bias file from assim
167      & cn_bias_tot       ! name of bias with total/rlx bias
168
169   ! Structure of input T and S bias offline (file informations, fields read)
170   TYPE(FLD), PRIVATE, ALLOCATABLE, DIMENSION(:) ::   sf_tbias_ofl 
171   TYPE(FLD), PRIVATE, ALLOCATABLE, DIMENSION(:) ::   sf_sbias_ofl
172
173   TYPE(FLD_N), PRIVATE ::&   ! information about the fields to be read
174      &  sn_tbias_ofl, sn_sbias_ofl 
175   
176CONTAINS
177
178   SUBROUTINE bias_init
179      !!----------------------------------------------------------------------
180      !!                    ***  ROUTINE bias_init  ***
181      !!
182      !! ** Purpose : Read in the bias namelist and read in the bias fields.
183      !!
184      !! ** Method  : Read in the bias namelist and read in the bias fields.
185      !!
186      !! ** Action  :
187      !!
188      !! History :
189      !!        !  08-05  (D. Lea)    Initial version
190      !!        !  08-10  (M. Martin) Tidy
191      !!        !  09-03  (M. Balmaseda). Generalize to estimate the bias
192      !!                                  from relax and offline bias term.
193      !!                                  Introduce parameters to control the
194      !!                                  model for the bias
195      !!                                  (variables and time evolution)
196      !!----------------------------------------------------------------------
197
198      IMPLICIT NONE
199     
200      !! * Local declarations
201
202      CHARACTER(len=100) ::  cn_dir       ! dir for location ofline bias
203      INTEGER            ::  ierror
204      INTEGER            ::  ios          ! Local integer output status for namelist read
205      REAL(wp)           ::  eft_rlx,  &  ! efolding time (bias memory) in days
206         &                   eft_asm,  &  !     "      "
207         &                   log2,     &
208         &                   lenscl_bias     !lengthscale of the pressure bias decay between minlat and maxlat.
209     
210      NAMELIST/nambias/ ln_bias, ln_bias_asm, ln_bias_rlx, ln_bias_ofl,   &
211         & ln_bias_ts_app, ln_bias_pc_app,                                &
212         & fb_t_asm, fb_t_rlx, fb_t_ofl, fb_p_asm, fb_p_rlx, fb_p_ofl,    &
213         & eft_rlx, t_rlx_upd, eft_asm, t_asm_upd, nn_lat_ramp,           &
214         & bias_time_unit_asm, bias_time_unit_rlx, bias_time_unit_ofl,    &
215         & cn_bias_tot, cn_bias_asm, cn_dir, sn_tbias_ofl, sn_sbias_ofl,  &
216         & ln_bsyncro, fctamp, rn_maxlat_bias, rn_minlat_bias,            &
217         & nn_bias_itwrt, ln_itdecay
218 
219
220      !-----------------------------------------------------------------------
221      ! Read Namelist : bias interface
222      !-----------------------------------------------------------------------
223
224
225      REWIND( numnam_ref )              ! Namelist nambias in reference namelist : Bias pressure correction
226      READ  ( numnam_ref, nambias, IOSTAT = ios, ERR = 901)
227901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambias in reference namelist', lwp )
228
229
230      ! Set additional default values (note that most values are set in the reference namelist)
231     
232      IF ( ln_asmiau ) nn_bias_itwrt = nitiaufin
233     
234      ! ... default values (NB: frequency positive => hours, negative => months)
235      !            !   file    ! frequency !  variable  ! time intep !  clim   ! 'yearly' or !
236      !            !   name    !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  !
237      sn_tbias_ofl = FLD_N( 'tbias_ofl'    ,    -1.    ,  'tbias'     ,  .TRUE.   , .FALSE. ,   'yearly', '', '', ''  )
238      sn_sbias_ofl = FLD_N( 'sbias_ofl'    ,    -1.    ,  'sbias'     ,  .TRUE.   , .FALSE. ,   'yearly', '', '', ''  )
239
240
241      REWIND( numnam_cfg )              ! Namelist nambias in configuration namelist : Bias pressure correction
242      READ  ( numnam_cfg, nambias, IOSTAT = ios, ERR = 902 )
243902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambias in configuration namelist', lwp )
244      IF(lwm) WRITE ( numond, nambias )
245
246
247      IF ( ( .NOT. ln_bias_asm ) .AND. ( .NOT. ln_bias_ofl ) .AND. ( .NOT. ln_bias_rlx ) ) THEN
248         ln_bias_ts_app = .FALSE.
249         ln_bias_pc_app = .FALSE.
250         ln_bias        = .FALSE.
251      ENDIF
252     
253      ! set up decay scales
254      log2           = LOG( 2.0_wp )
255      t_rlx_mem      = EXP( - log2 * rdt / ( eft_rlx * rday ) )
256      t_asm_mem      = EXP( - log2 * bias_time_unit_asm/ ( eft_asm * rday ) )
257     
258      ! Control print
259      IF(lwp) THEN
260         WRITE(numout,*)
261         WRITE(numout,*) 'bias_init : '
262         WRITE(numout,*) '~~~~~~~~~~~ '
263         WRITE(numout,*) '  Namelist nambias : '
264
265         WRITE(numout,*) '  Bias switches/options/variables '
266         WRITE(numout,*) '     bias main switch               ln_bias        = ',ln_bias
267         WRITE(numout,*) '     bias from assim                ln_bias_asm    = ',ln_bias_asm
268         WRITE(numout,*) '     bias from relax                ln_bias_rlx    = ',ln_bias_rlx
269         WRITE(numout,*) '     bias from offln                ln_bias_ofl    = ',ln_bias_ofl
270         WRITE(numout,*) '     bias T and S apply             ln_bias_ts_app = ',ln_bias_ts_app
271         WRITE(numout,*) '     bias pressure correctn apply   ln_bias_pc_app = ',ln_bias_pc_app
272         WRITE(numout,*) '     bias pressure correctn apply   ln_bias_pc_app = ',ln_bias_pc_app
273         WRITE(numout,*) '     lat ramp for bias correction   nn_lat_ramp    = ',nn_lat_ramp
274         WRITE(numout,*) '     time step for writing bias fld nn_bias_itwrt  = ',nn_bias_itwrt
275    WRITE(numout,*) '     evolve pcbias at each timestep ln_itdecay     = ',ln_itdecay
276         WRITE(numout,*) '  Parameters for parition of bias term '
277         WRITE(numout,*) '                                    fb_t_rlx       = ',fb_t_rlx
278         WRITE(numout,*) '                                    fb_t_asm       = ',fb_t_asm
279         WRITE(numout,*) '                                    fb_t_ofl       = ',fb_t_ofl
280         WRITE(numout,*) '                                    fb_p_rlx       = ',fb_p_rlx
281         WRITE(numout,*) '                                    fb_p_asm       = ',fb_p_asm
282         WRITE(numout,*) '                                    fb_p_ofl       = ',fb_p_ofl
283         WRITE(numout,*) '  Parameters for time evolution of bias '
284         WRITE(numout,*) '  Rlx   efolding time (mem)         eft_rlx,t_rlx_mem = ', eft_rlx, t_rlx_mem, 1. - log2 * rdt / (eft_rlx * rday)
285         WRITE(numout,*) '        uptdate factor              t_rlx_upd         = ',t_rlx_upd
286         WRITE(numout,*) '  Asm   efolding time (mem)         eft_asm,t_asm_mem = ', eft_asm, t_asm_mem, 1. - log2 * rdt / (eft_asm * rday)
287         WRITE(numout,*) '        uptdate factor              t_asm_upd         = ',t_asm_upd
288         WRITE(numout,*) '  Filenames and input structures'
289         WRITE(numout,*) '     bias_tot filename              cn_bias_to     = ',cn_bias_tot
290         WRITE(numout,*) '     bias_asm filename              cn_bias_asm    = ',cn_bias_asm
291         WRITE(numout,*) '     bias_asm time unit (secs)  bias_time_unit_asm = ',bias_time_unit_asm
292         WRITE(numout,*) '     structure Tem bias ofl         sn_tbias_ofl   = ',sn_tbias_ofl
293         WRITE(numout,*) '     structure Sal bias ofl         sn_sbias_ofl   = ',sn_sbias_ofl
294         
295         IF ( ( (.NOT. ln_tsd_tradmp) .OR. (.NOT. ln_tradmp) ) .AND. ln_bias_rlx ) &
296            &   CALL ctl_stop (' lk_dtatem, lk_dtasal and lk_tradmp need to be true with ln_bias_rlx' )
297
298      ENDIF
299      IF( .NOT. ln_bias ) RETURN
300
301      IF( .NOT. lalloc ) THEN
302
303         ALLOCATE( tbias(jpi,jpj,jpk)  , &
304            &      sbias(jpi,jpj,jpk)  , &
305            &      tbias_p(jpi,jpj,jpk), &
306            &      sbias_p(jpi,jpj,jpk), &
307            &      rhd_pc(jpi,jpj,jpk) , &
308            &      gru_pc(jpi,jpj)     , &
309            &      grv_pc(jpi,jpj)       )
310
311         ALLOCATE( fbcoef(jpi,jpj)       )
312         
313         IF( ln_bias_asm ) ALLOCATE(  tbias_asm(jpi,jpj,jpk),     &
314            &                         sbias_asm(jpi,jpj,jpk),     &
315                                      tbias_asm_out(jpi,jpj,jpk), & 
316                                      sbias_asm_out(jpi,jpj,jpk), & 
317                                      tbias_p_out(jpi,jpj,jpk),   &             
318                                      sbias_p_out(jpi,jpj,jpk)    )
319
320         IF( ln_bias_rlx ) ALLOCATE(  tbias_rlx(jpi,jpj,jpk),     &
321            &                         sbias_rlx(jpi,jpj,jpk),     &
322                                      tbias_rlx_out(jpi,jpj,jpk), & 
323                                      sbias_rlx_out(jpi,jpj,jpk)  )
324         lalloc = .TRUE.
325
326      ENDIF
327
328      IF( ln_bias_ofl ) THEN      ! set sf_tbias_ofl and sf_sbias_ofl strctrs
329         !
330         ! tbias
331         !
332         ALLOCATE( sf_tbias_ofl(1), STAT=ierror )
333         IF( ierror > 0 ) THEN
334            CALL ctl_stop( 'bias_init: unable to allocate sf_tbias_ofl structure' )   ;    RETURN
335         ENDIF
336         ALLOCATE( sf_tbias_ofl(1)%fnow(jpi,jpj,jpk) )
337         ALLOCATE( sf_tbias_ofl(1)%fdta(jpi,jpj,jpk,2) )
338         
339         ! fill structure with values and control print
340         CALL fld_fill( sf_tbias_ofl, (/ sn_tbias_ofl /), cn_dir, 'bias_init', 'Offline T bias term ', 'nam_tbias_ofl' )
341         !
342         ! salinity bias
343         !
344         ALLOCATE( sf_sbias_ofl(1), STAT=ierror )
345         IF( ierror > 0 ) THEN
346            CALL ctl_stop( 'bias_init: unable to allocate sf_sbias_ofl structure' )   ;   RETURN
347         ENDIF
348         ALLOCATE( sf_sbias_ofl(1)%fnow(jpi,jpj,jpk) )
349         ALLOCATE( sf_sbias_ofl(1)%fdta(jpi,jpj,jpk,2) )
350         
351         ! fill structure with values and control print
352         CALL fld_fill( sf_sbias_ofl, (/ sn_sbias_ofl /), cn_dir, 'bias_init', 'Offline S bias term ', 'nam_sbias_ofl' )
353         
354      ENDIF
355
356      ! Read total bias
357      IF ( ln_bias ) THEN
358         tbias(:,:,:)   = 0.0_wp
359         sbias(:,:,:)   = 0.0_wp
360         tbias_p(:,:,:) = 0.0_wp
361         sbias_p(:,:,:) = 0.0_wp
362         gru_pc(:,:)    = 0.0_wp
363         grv_pc(:,:)    = 0.0_wp
364         IF ( ln_bias_rlx ) THEN
365            tbias_rlx(:,:,:) = 0.0_wp
366            sbias_rlx(:,:,:) = 0.0_wp
367         ENDIF
368         IF ( ln_bias_asm ) THEN   !now rlx and asm bias in same file
369           tbias_asm(:,:,:) = 0.0_wp
370           sbias_asm(:,:,:) = 0.0_wp
371         ENDIF
372         numbias_tot    = 0
373         ! Get bias from file and prevent fail if the file does not exist
374         IF(lwp) WRITE(numout,*) 'Opening ',TRIM( cn_bias_tot ) 
375         CALL iom_open( cn_bias_tot, numbias_tot, ldstop=.FALSE. )
376         
377         IF ( numbias_tot > 0 ) THEN     
378            ! Could check validity time of bias fields here...
379            ! Get the T and S bias data
380            IF(lwp) WRITE(numout,*) 'Reading bias fields from tot...'
381
382            !Search for bias from relaxation term if needed. Use same file
383            IF ( ln_bias_rlx ) THEN
384               IF(lwp) WRITE(numout,*) 'Reading bias fields for bias rlx from file ',cn_bias_tot
385               IF( iom_varid( numbias_tot, 'tbias_rlx' ) > 0 ) THEN
386                  ! Get the T and S bias data
387                  CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_rlx', tbias_rlx )
388                  CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_rlx', sbias_rlx )
389               ELSE
390                  CALL ctl_stop( 'Bias relaxation variables not found in ',cn_bias_tot )
391               ENDIF
392            ENDIF
393
394
395            !Search for bias from assim term if needed. Use same file
396            IF ( ln_bias_asm .and. .not.ln_bsyncro ) THEN
397               IF(lwp) WRITE(numout,*) 'Reading a-syncro bias fields for bias asm from file ',cn_bias_tot
398               IF( iom_varid( numbias_tot, 'tbias_asm' ) > 0 ) THEN
399                  ! Get the T and S bias data
400                  CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_asm', tbias_asm )
401                  CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_asm', sbias_asm )
402               ELSE
403                  CALL ctl_stop( 'Bias assim variables not found in ',cn_bias_tot )
404               ENDIF
405            ENDIF
406
407            ! Close the file
408            CALL iom_close(numbias_tot)
409
410         ELSE
411            IF(lwp) WRITE(numout,*) 'No bias file found so T and S bias fields are set to zero'
412         ENDIF
413
414      ENDIF
415
416     ! for the time being, the bias_asm is read in the same file as
417     ! bias_rlx
418     ! Implications: bias_asm is estimated/evolved in time in the second outer
419     !               loop only, when the assimilation increments are ready.
420     !               bias_asm is kept cte during the first outer loop.
421     !              => Assyncronous bias correction.
422     ! Alternative: Syncronous bias correction:
423     !              bias_asm estimated/evolved in the first outer loop
424     !              with the asm incrments of the previous cycle.
425     !              bias_asm kept cte during the second outer loop.
426     !              Implication: bias_asm should be estimated really in the
427     !              inner loop.
428      IF ( ln_bsyncro ) THEN
429      ! Read bias from assimilation from a separate file
430      IF ( ln_bias_asm ) THEN
431         tbias_asm(:,:,:) = 0.0_wp
432         sbias_asm(:,:,:) = 0.0_wp
433         numbias_asm      = 0
434         ! Get bias from file and prevent fail if the file does not exist
435         IF(lwp) WRITE(numout,*) 'Opening file for syncro assim bias ',TRIM( cn_bias_asm ) 
436         CALL iom_open( cn_bias_asm, numbias_asm, ldstop=.FALSE. )
437         
438         IF ( numbias_asm > 0 ) THEN     
439            ! Could check validity time of bias fields here...
440           
441            ! Get the T and S bias data
442            IF(lwp) WRITE(numout,*) 'Reading syncro bias fields from asm from file ',cn_bias_asm
443            CALL iom_get( numbias_asm, jpdom_autoglo, 'tbias_asm', tbias_asm )
444            CALL iom_get( numbias_asm, jpdom_autoglo, 'sbias_asm', sbias_asm )
445           
446!  this is only applicable if tbias_asm were to be calculated in the inner loop
447            tbias_asm(:,:,:) = tbias_asm(:,:,:) * rdt / bias_time_unit_asm
448            sbias_asm(:,:,:) = sbias_asm(:,:,:) * rdt / bias_time_unit_asm
449           
450            ! Close the file
451            CALL iom_close(numbias_asm)
452           
453         ELSE
454            IF(lwp) WRITE(numout,*) 'No bias file found from asm so T and S bias fields are set to zero'
455         ENDIF
456         
457      ENDIF
458
459      ENDIF
460
461      !latitudinal dependence of partition coeficients. Adhoc
462      IF ( nn_lat_ramp == 1 ) THEN
463         ! Use the inertial ramp.
464         lenscl_bias = ( rn_maxlat_bias - rn_minlat_bias )*2._wp
465         WHERE ( abs( gphit(:,:) ) <= rn_minlat_bias )         
466            fbcoef(:,:) = 0._wp         
467         ELSEWHERE ( abs( gphit(:,:) ) >= rn_maxlat_bias )         
468            fbcoef(:,:) = 1._wp                   
469         ELSEWHERE       
470            fbcoef(:,:) = 1._wp - exp( -( abs( gphit(:,:) ) - rn_minlat_bias ) &
471                           * ( abs( gphit(:,:) ) - rn_minlat_bias ) / lenscl_bias )                         
472         ENDWHERE 
473      ELSEIF ( nn_lat_ramp == 2 ) THEN   
474         ! Use a linear ramp consist with the geostrophic velocity balance ramp in NEMOVAR
475     
476         WHERE ( abs( gphit(:,:) ) <= rn_minlat_bias )
477            fbcoef(:,:) = 0._wp
478         ELSEWHERE ( abs( gphit(:,:) ) >= rn_maxlat_bias ) 
479            fbcoef(:,:) = 1._wp
480         ELSEWHERE
481            fbcoef(:,:) = 1._wp - ((rn_maxlat_bias - abs( gphit(:,:)))/(rn_maxlat_bias - rn_minlat_bias))
482         ENDWHERE
483      ELSE
484         fbcoef(:,:) = 0.0_wp
485         fctamp      = 0.0_wp
486      ENDIF
487
488
489   END SUBROUTINE bias_init
490
491   SUBROUTINE tra_bias ( kt )
492      !!----------------------------------------------------------------------
493      !!                    ***  ROUTINE tra_bias  ***
494      !!
495      !! ** Purpose : Update bias field every time step
496      !!
497      !! ** Method  : add contributions to bias from 3 terms
498      !!
499      !! ** Action  : Bias from assimilation (read in bias_init)
500      !!              Bias from relaxation term is estimated according to
501      !!              the prescribed time evolution of the bias
502      !!              Bias from ofl term is read from external file
503      !!              The difference contributions are added and the partition
504      !!              into direct bias in T/S and pressure perfomed.
505      !!
506      !! History :  09-03  (M. Balmaseda)
507      !!----------------------------------------------------------------------
508      !! called every timestep after dta_sst if ln_bias true.
509
510      IMPLICIT NONE
511
512      !! * Arguments
513      INTEGER, INTENT(IN) ::   kt             ! ocean time-step index
514      !! * Local variables
515      INTEGER             ::   ji,jj,jk, it   ! local loop index
516      REAL(wp)            ::   tsclf          ! time saling factor
517      REAL(wp)            ::   fb_t_asm_max, fb_t_rlx_max, fb_t_ofl_max
518      REAL(wp)            ::   ztfrac, ztsday
519      REAL(wp), DIMENSION(jpi,jpj) :: zcof1, zcof2
520
521      IF ( .NOT. ln_bias ) RETURN
522      fb_t_rlx_max   = MIN(fb_t_rlx*fctamp,1.0_wp)
523      fb_t_asm_max   = MIN(fb_t_asm*fctamp,1.0_wp)
524      fb_t_ofl_max   = MIN(fb_t_ofl*fctamp,1.0_wp)
525
526      tbias(:,:,:)   = 0.0_wp
527      sbias(:,:,:)   = 0.0_wp
528      tbias_p(:,:,:) = 0.0_wp
529      sbias_p(:,:,:) = 0.0_wp
530      IF ( ln_bias_asm ) THEN
531     
532         tsclf = 1
533         IF ( .NOT.ln_bsyncro ) tsclf = rdt / bias_time_unit_asm 
534         zcof1(:,:) = tsclf * ( ( 1.0_wp - fbcoef(:,:) ) * fb_t_asm + &
535            &                              fbcoef(:,:)   * fb_t_asm_max )
536         zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_asm
537     
538         IF ( ln_itdecay ) THEN   !decay the pressure correction at each time step
539   
540       ztsday  = rday / real(rdt)
541
542            IF( ln_asmiau .and. ln_trainc ) THEN  ! nudge in increments and decay historical contribution
543
544               
545               IF ( kt <= nitiaufin ) THEN  ! During IAU calculate the fraction of increments to apply at each time step
546
547                  ztfrac = real(kt) / real(nitiaufin)  ! nudging factor during the IAU
548
549         
550                  IF (lwp) THEN
551                     WRITE(numout,*) 'tra_bias : bias weights'
552                     WRITE(numout,*) '~~~~~~~~~~~~'
553                     WRITE(numout,* ) "proportion of  increment applied in pcbias ",ztfrac
554                     WRITE(numout,* ) "proportion of  historical pcbias applied ",t_asm_mem**(real(kt)/ztsday)
555                  ENDIF
556           
557                  DO jk = 1, jpkm1
558                     tbias(:,:,jk) = tbias(:,:,jk) +                            &
559                     &                ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk)  +                    &
560                     &                ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:)
561                     sbias(:,:,jk) = sbias(:,:,jk) +                            &
562                     &               ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk)  +                     &
563                     &               ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:)
564                     tbias_p(:,:,jk) = tbias_p(:,:,jk) +                        &
565                     &               ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk)  +                     &
566                     &               ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:)
567                     sbias_p(:,:,jk) = sbias_p(:,:,jk) +                        & 
568                     &               ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk)  +                     &
569                     &               ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:)
570                  ENDDO
571
572                  IF ( .not.ln_bsyncro ) THEN
573                     IF ( kt == nn_bias_itwrt ) THEN
574                        DO jk = 1, jpkm1
575                           tbias_asm_out(:,:,jk) =  t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk)  +       &
576                           &                     ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
577                           sbias_asm_out(:,:,jk) =  t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk)  +       &
578                           &                     ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
579                         END DO
580                     ENDIF
581                  ENDIF
582
583                  ! update the historical component with the increments at the end of the IAU
584                  IF ( kt == nitiaufin ) THEN
585                     DO jk = 1, jpkm1
586                        tbias_asm(:,:,jk) =  t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk)  +       &
587                        &                     ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
588                        sbias_asm(:,:,jk) =  t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk)  +       &
589                        &                     ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
590                     END DO
591                  ENDIF
592               
593               ELSE ! decay pressure correction from combined historical component and increments after IAU
594
595                  ztfrac=t_asm_mem**(real(kt-nitiaufin)/real(nitiaufin))  ! decay from end of IAU
596
597                  DO jk = 1, jpkm1
598                     tbias(:,:,jk) = tbias(:,:,jk) +                            &
599                     &                ( ztfrac * tbias_asm(:,:,jk) ) * zcof1(:,:)
600                     sbias(:,:,jk) = sbias(:,:,jk) +                            &
601                     &               (  ztfrac * sbias_asm(:,:,jk) ) * zcof1(:,:)
602                     tbias_p(:,:,jk) = tbias_p(:,:,jk) +                        &
603                     &               (  ztfrac * tbias_asm(:,:,jk) ) * zcof2(:,:)
604                     sbias_p(:,:,jk) = sbias_p(:,:,jk) +                        & 
605                     &               ( ztfrac * sbias_asm(:,:,jk) ) * zcof2(:,:)
606                  ENDDO
607
608                  IF (lwp) THEN
609                     WRITE(numout,*) 'tra_bias : bias weights'
610                     WRITE(numout,*) '~~~~~~~~~~~~'
611                     WRITE(numout,* ) "proportion of increments + historical pcbias applied ",ztfrac
612                  ENDIF
613
614                  IF ( ln_trainc .and. .not.ln_bsyncro ) THEN
615                     IF ( kt == nn_bias_itwrt ) THEN
616                        DO jk = 1, jpkm1
617                           tbias_asm_out(:,:,jk) =  ztfrac * tbias_asm(:,:,jk) 
618                           sbias_asm_out(:,:,jk) =  ztfrac * sbias_asm(:,:,jk) 
619                         END DO
620                     ENDIF
621                  ENDIF
622
623               ENDIF
624
625
626            ELSE ! no assimilation increments, simply decay pressure correction (e.g for forecasts)
627
628               DO jk = 1, jpkm1
629                  tbias(:,:,jk) = tbias(:,:,jk) +                                                         &
630                  &               ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) ) * zcof1(:,:)
631                  sbias(:,:,jk) = sbias(:,:,jk) +                                                         &
632                  &               ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) ) * zcof1(:,:)
633                  tbias_p(:,:,jk) = tbias_p(:,:,jk) +                                                     &
634                  &               ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) ) * zcof2(:,:)
635                  sbias_p(:,:,jk) = sbias_p(:,:,jk) +                                                     & 
636                  &               ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) ) * zcof2(:,:)
637               ENDDO
638
639               IF (lwp) THEN
640                  WRITE(numout,*) 'tra_bias : bias weights'
641                  WRITE(numout,*) '~~~~~~~~~~~~'
642                  WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(real(kt)/ztsday)
643               ENDIF
644
645            ENDIF
646 
647         ELSE ! maintain a constant pressure correction 
648
649            DO jk = 1, jpkm1
650               tbias(:,:,jk) = tbias(:,:,jk) + tbias_asm(:,:,jk) * zcof1(:,:)
651               sbias(:,:,jk) = sbias(:,:,jk) + sbias_asm(:,:,jk) * zcof1(:,:)
652               tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_asm(:,:,jk) * zcof2(:,:)
653               sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_asm(:,:,jk) * zcof2(:,:)
654            END DO     
655
656            IF( ln_asmiau .and. ln_trainc .and. .not.ln_bsyncro ) THEN   
657            ! if last outer loop (ln_asmiau=true and ln_trainc=true). t/sbias_asm
658            ! is updated, only once (end of run) taking into account units.
659               IF ( kt == nn_bias_itwrt ) THEN
660                 IF(lwp) WRITE(numout,*)' estimating asm bias at timestep: ',kt
661                 DO jk = 1, jpkm1
662                   tbias_asm_out(:,:,jk) = t_asm_mem * tbias_asm(:,:,jk)  +             &
663                   &                      t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
664                   sbias_asm_out(:,:,jk) = t_asm_mem * sbias_asm(:,:,jk) +             &
665                   &                      t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)               
666                 END DO
667               ENDIF
668             ENDIF
669 
670         ENDIF
671
672      ENDIF
673
674
675#if   defined key_tradmp 
676      ! Time evolution of bias from relaxation
677      IF ( ln_bias_rlx ) THEN
678         tbias_rlx(:,:,:) = t_rlx_mem * tbias_rlx(:,:,:) + &
679            &               t_rlx_upd * ttrdmp(:,:,:) * rdt
680         sbias_rlx(:,:,:) = t_rlx_mem * sbias_rlx(:,:,:) + &
681            &               t_rlx_upd * strdmp(:,:,:) * rdt
682         zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_rlx +      &
683            &                    fbcoef(:,:)   * fb_t_rlx_max
684         zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_rlx
685         DO jk = 1, jpkm1
686            tbias(:,:,jk)   = tbias(:,:,jk) + tbias_rlx(:,:,jk) * zcof1(:,:)
687            sbias(:,:,jk)   = sbias(:,:,jk) + sbias_rlx(:,:,jk) * zcof1(:,:)
688            tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_rlx(:,:,jk) * zcof2(:,:)
689            sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_rlx(:,:,jk) * zcof2(:,:)
690         ENDDO
691         IF ( kt == nn_bias_itwrt ) THEN
692            tbias_rlx_out(:,:,:) = tbias_rlx(:,:,:)
693            sbias_rlx_out(:,:,:) = sbias_rlx(:,:,:)
694         ENDIF
695      ENDIF
696#endif
697      ! ofline bias
698      IF ( kt == nit000 ) THEN
699         IF(lwp) WRITE(numout,*) ' tra_bias: ln_bias_ofl = ',ln_bias_ofl
700         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~'
701      ENDIF
702      IF ( ln_bias_ofl ) THEN
703         IF(lwp) WRITE(numout,*) 'reading offline bias'
704         CALL fld_read( kt, 1, sf_tbias_ofl ) 
705         CALL fld_read( kt, 1, sf_sbias_ofl ) 
706
707         zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_ofl +           &
708            &                    fbcoef(:,:)   * fb_t_ofl_max
709         zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_ofl
710         DO jk = 1, jpkm1
711            tbias(:,:,jk)   = tbias(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:)
712            sbias(:,:,jk)   = sbias(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:)
713            tbias_p(:,:,jk) = tbias_p(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:)
714            sbias_p(:,:,jk) = sbias_p(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:)
715         ENDDO
716      ENDIF
717
718
719      !apply bias on tracers if needed     
720      IF ( ln_bias_ts_app ) THEN
721         
722         ! Add the bias directely to T/s
723         tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) + tmask(:,:,:) * tbias(:,:,:) / rdt
724         tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) + tmask(:,:,:) * sbias(:,:,:) / rdt
725
726         ! lateral boundary conditions (is this needed?)
727         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1.0_wp )
728         CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1.0_wp )   
729
730      ENDIF
731
732      IF ( kt == nn_bias_itwrt ) THEN
733         tbias_p_out(:,:,:)   = tbias_p(:,:,:)
734         sbias_p_out(:,:,:)   = sbias_p(:,:,:)
735      ENDIF
736
737   END SUBROUTINE tra_bias
738
739   SUBROUTINE dyn_bias( kt )
740      !!----------------------------------------------------------------------
741      !!                   ***  ROUTINE dyn_bias  ***
742      !!
743      !! ** Purpose :   Computes rhd_pc, gru/v_pc bias corrected
744      !!                for hydrostatic pressure gradient
745      !!                depending on time step (semi-implicit/centered)
746      !!                If partial steps computes bottom pressure gradient.
747      !!                These correction terms will affect only the dynamical
748      !!                component (pressure gradient calculation), but not
749      !!                the isopycnal calculation for the lateral mixing.
750      !!
751      !! ** Method  :   At this stage of the computation, ta and sa are the
752      !!                after temperature and salinity. If semi-implicit, these
753      !!                are used to compute rho and bottom pressure gradient.
754      !!                If centered, tb,sb are used instead.
755      !!                If bias key is activated, the temperature,salinity are
756      !!                bias corrected in the calculation of the density fields
757      !!                used in the pressure gradient calculation.
758      !!
759      !!
760      !! ** Action  : - rhd_pc ready. rhop will be overwriten later
761      !!              - if ln_zps, bottom density gradients gru/v_pc ready.
762      !!----------------------------------------------------------------------
763      !!
764      !! * Arguments
765      INTEGER, INTENT(IN) ::   kt    ! ocean time-step index
766      !! * Local variables
767      REAL(wp) :: tsw(jpi,jpj,jpk,jpts)
768      !!
769      !!----------------------------------------------------------------------
770      !
771      ! gtu,gsu,gtv,gsv rhop will be overwritten later in step.
772      !
773      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg
774         tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:)
775         tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:)
776      ELSE
777         tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:)
778         tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:)
779      ENDIF
780
781      CALL eos( tsw, rhd_pc, rhop , fsdept_n(:,:,:))
782     
783      ! is this needed?
784      CALL lbc_lnk( rhd_pc, 'T', 1.0_wp )
785
786      ! Partial steps: now horizontal gradient of t,s,rd
787      ! at the bottom ocean level
788      IF( ln_zps    )  THEN
789         CALL zps_hde( kt, jpts, tsw, gtsu, gtsv,  &
790            &          rhd_pc, gru_pc , grv_pc  )
791      ENDIF
792
793   END SUBROUTINE dyn_bias
794   
795   SUBROUTINE bias_opn( kt )
796      !!---------------------------------------------------------------------
797      !!                   ***  ROUTINE bias_opn  ***
798      !!                     
799      !! ** Purpose :  open bias restart file following the same logic as the
800      !!               standard restarts.
801      !!----------------------------------------------------------------------
802      !! * Arguments
803      INTEGER, INTENT(IN) ::   kt     ! ocean time-step
804      !! * Local variables
805      CHARACTER(LEN=20)   ::   clbkt    ! ocean time-step deine as a character
806      CHARACTER(LEN=50)   ::   clbias_tot   ! total bias restart file name
807      !!----------------------------------------------------------------------
808      !
809      IF( lrst_oce .AND. .NOT.lrst_bias ) THEN       ! create bias file
810         IF( nitend > 999999999 ) THEN   ;   WRITE(clbkt, *       ) nitend
811         ELSE                            ;   WRITE(clbkt, '(i8.8)') nitend
812         ENDIF
813         clbias_tot = TRIM(cexper)//"_"//TRIM(ADJUSTL(clbkt))//"_"//TRIM(cn_bias_tot)
814         IF(lwp) THEN
815            WRITE(numout,*)
816            SELECT CASE ( jprstlib )
817            CASE ( jprstdimg )   ;   WRITE(numout,*) '             open tot bias restart binary file: '//clbias_tot
818            CASE DEFAULT         ;   WRITE(numout,*) '             open tot bias restart NetCDF file: '//clbias_tot
819            END SELECT
820         ENDIF
821         CALL iom_open( clbias_tot, numbias_tot , ldwrt = .TRUE., kiolib = jprstlib )
822         lrst_bias=.TRUE.
823
824      ENDIF
825      !
826   END SUBROUTINE bias_opn
827
828   SUBROUTINE bias_wrt( kt )
829      !!---------------------------------------------------------------------
830      !!                   ***  ROUTINE bias_wrt  ***
831      !!                     
832      !! ** Purpose :   Write bias restart fields in the format corresponding to jprstlib
833      !!
834      !! ** Method  :   Write in numbias_tot when kt == nitend in output
835      !!                file, save fields which are necessary for restart.
836      !!
837      !! Changed the timestep for writing to nitend rather than nitrst as this causes a
838      !! problem if the nstock list is used to determine the restart output times and
839      !! means that the bias is not output at all. M. Martin. 08/2011.
840      !! Need to check with Magdalena that this is ok for ECMWF.
841      !!----------------------------------------------------------------------
842      !! * Arguments
843      INTEGER, INTENT(IN) ::   kt   ! ocean time-step
844      !!----------------------------------------------------------------------
845      !                                                                     ! the begining of the run [s]
846
847      IF ( ln_bias_rlx ) THEN
848         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_rlx' , tbias_rlx_out )   
849         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_rlx' , sbias_rlx_out )   
850      ENDIF
851     
852      IF ( ln_bias_asm ) THEN
853         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_asm' , tbias_asm_out )   
854         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_asm' , sbias_asm_out )   
855         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_p'   , tbias_p_out )
856         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_p'   , sbias_p_out )         
857      ENDIF
858     
859      IF( kt == nitend ) THEN
860         CALL iom_close( numbias_tot )     ! close the restart file (only at last time step)
861         lrst_bias = .FALSE.
862      ENDIF
863      !
864   END SUBROUTINE bias_wrt
865
866END MODULE bias
Note: See TracBrowser for help on using the repository browser.