source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/ASM/asminc.F90 @ 11099

Last change on this file since 11099 was 11099, checked in by davestorkey, 18 months ago

dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap : Updates and bug fixes. This version still doesn't pass all SETTE tests.

  • Property svn:keywords set to Id
File size: 49.0 KB
Line 
1MODULE asminc
2   !!======================================================================
3   !!                       ***  MODULE asminc  ***
4   !! Assimilation increment : Apply an increment generated by data
5   !!                          assimilation
6   !!======================================================================
7   !! History :       ! 2007-03  (M. Martin)  Met Office version
8   !!                 ! 2007-04  (A. Weaver)  calc_date original code
9   !!                 ! 2007-04  (A. Weaver)  Merge with OPAVAR/NEMOVAR
10   !!   NEMO     3.3  ! 2010-05  (D. Lea)  Update to work with NEMO v3.2
11   !!             -   ! 2010-05  (D. Lea)  add calc_month_len routine based on day_init
12   !!            3.4  ! 2012-10  (A. Weaver and K. Mogensen) Fix for direct initialization
13   !!                 ! 2014-09  (D. Lea)  Local calc_date removed use routine from OBS
14   !!                 ! 2015-11  (D. Lea)  Handle non-zero initial time of day
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!   asm_inc_init   : Initialize the increment arrays and IAU weights
19   !!   tra_asm_inc    : Apply the tracer (T and S) increments
20   !!   dyn_asm_inc    : Apply the dynamic (u and v) increments
21   !!   ssh_asm_inc    : Apply the SSH increment
22   !!   ssh_asm_div    : Apply divergence associated with SSH increment
23   !!   seaice_asm_inc : Apply the seaice increment
24   !!----------------------------------------------------------------------
25   USE oce             ! Dynamics and active tracers defined in memory
26   USE par_oce         ! Ocean space and time domain variables
27   USE dom_oce         ! Ocean space and time domain
28   USE domvvl          ! domain: variable volume level
29   USE ldfdyn          ! lateral diffusion: eddy viscosity coefficients
30   USE eosbn2          ! Equation of state - in situ and potential density
31   USE zpshde          ! Partial step : Horizontal Derivative
32   USE asmpar          ! Parameters for the assmilation interface
33   USE asmbkg          !
34   USE c1d             ! 1D initialization
35   USE sbc_oce         ! Surface boundary condition variables.
36   USE diaobs   , ONLY : calc_date     ! Compute the calendar date on a given step
37#if defined key_si3
38   USE ice      , ONLY : hm_i, at_i, at_i_b
39#endif
40   !
41   USE in_out_manager  ! I/O manager
42   USE iom             ! Library to read input files
43   USE lib_mpp         ! MPP library
44
45   IMPLICIT NONE
46   PRIVATE
47   
48   PUBLIC   asm_inc_init   !: Initialize the increment arrays and IAU weights
49   PUBLIC   tra_asm_inc    !: Apply the tracer (T and S) increments
50   PUBLIC   dyn_asm_inc    !: Apply the dynamic (u and v) increments
51   PUBLIC   ssh_asm_inc    !: Apply the SSH increment
52   PUBLIC   ssh_asm_div    !: Apply the SSH divergence
53   PUBLIC   seaice_asm_inc !: Apply the seaice increment
54
55#if defined key_asminc
56    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .TRUE.   !: Logical switch for assimilation increment interface
57#else
58    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE.  !: No assimilation increments
59#endif
60   LOGICAL, PUBLIC :: ln_bkgwri     !: No output of the background state fields
61   LOGICAL, PUBLIC :: ln_asmiau     !: No applying forcing with an assimilation increment
62   LOGICAL, PUBLIC :: ln_asmdin     !: No direct initialization
63   LOGICAL, PUBLIC :: ln_trainc     !: No tracer (T and S) assimilation increments
64   LOGICAL, PUBLIC :: ln_dyninc     !: No dynamics (u and v) assimilation increments
65   LOGICAL, PUBLIC :: ln_sshinc     !: No sea surface height assimilation increment
66   LOGICAL, PUBLIC :: ln_seaiceinc  !: No sea ice concentration increment
67   LOGICAL, PUBLIC :: ln_salfix     !: Apply minimum salinity check
68   LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
69   INTEGER, PUBLIC :: nn_divdmp     !: Apply divergence damping filter nn_divdmp times
70
71   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkg   , s_bkg      !: Background temperature and salinity
72   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkg   , v_bkg      !: Background u- & v- velocity components
73   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkginc, s_bkginc   !: Increment to the background T & S
74   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components
75   REAL(wp), PUBLIC, DIMENSION(:)    , ALLOCATABLE ::   wgtiau               !: IAU weights for each time step
76#if defined key_asminc
77   REAL(wp), PUBLIC, DIMENSION(:,:)  , ALLOCATABLE ::   ssh_iau              !: IAU-weighted sea surface height increment
78#endif
79   !                                !!! time steps relative to the cycle interval [0,nitend-nit000-1]
80   INTEGER , PUBLIC ::   nitbkg      !: Time step of the background state used in the Jb term
81   INTEGER , PUBLIC ::   nitdin      !: Time step of the background state for direct initialization
82   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval
83   INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval
84   !
85   INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting
86   !                                 !: = 1   Linear hat-like, centred in middle of IAU interval
87   REAL(wp), PUBLIC ::   salfixmin   !: Ensure that the salinity is larger than this  value if (ln_salfix)
88
89   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment
90   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc
91#if defined key_cice && defined key_asminc
92   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ndaice_da             ! ice increment tendency into CICE
93#endif
94
95   !! * Substitutions
96#  include "vectopt_loop_substitute.h90"
97   !!----------------------------------------------------------------------
98   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
99   !! $Id$
100   !! Software governed by the CeCILL license (see ./LICENSE)
101   !!----------------------------------------------------------------------
102CONTAINS
103
104   SUBROUTINE asm_inc_init( Kbb, Kmm, Krhs )
105      !!----------------------------------------------------------------------
106      !!                    ***  ROUTINE asm_inc_init  ***
107      !!         
108      !! ** Purpose : Initialize the assimilation increment and IAU weights.
109      !!
110      !! ** Method  : Initialize the assimilation increment and IAU weights.
111      !!
112      !! ** Action  :
113      !!----------------------------------------------------------------------
114      INTEGER, INTENT(in) ::  Kbb, Kmm, Krhs  ! time level indices
115      !
116      INTEGER :: ji, jj, jk, jt  ! dummy loop indices
117      INTEGER :: imid, inum      ! local integers
118      INTEGER :: ios             ! Local integer output status for namelist read
119      INTEGER :: iiauper         ! Number of time steps in the IAU period
120      INTEGER :: icycper         ! Number of time steps in the cycle
121      REAL(KIND=dp) :: ditend_date     ! Date YYYYMMDD.HHMMSS of final time step
122      REAL(KIND=dp) :: ditbkg_date     ! Date YYYYMMDD.HHMMSS of background time step for Jb term
123      REAL(KIND=dp) :: ditdin_date     ! Date YYYYMMDD.HHMMSS of background time step for DI
124      REAL(KIND=dp) :: ditiaustr_date  ! Date YYYYMMDD.HHMMSS of IAU interval start time step
125      REAL(KIND=dp) :: ditiaufin_date  ! Date YYYYMMDD.HHMMSS of IAU interval final time step
126
127      REAL(wp) :: znorm        ! Normalization factor for IAU weights
128      REAL(wp) :: ztotwgt      ! Value of time-integrated IAU weights (should be equal to one)
129      REAL(wp) :: z_inc_dateb  ! Start date of interval on which increment is valid
130      REAL(wp) :: z_inc_datef  ! End date of interval on which increment is valid
131      REAL(wp) :: zdate_bkg    ! Date in background state file for DI
132      REAL(wp) :: zdate_inc    ! Time axis in increments file
133      !
134      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zhdiv   ! 2D workspace
135      !!
136      NAMELIST/nam_asminc/ ln_bkgwri,                                      &
137         &                 ln_trainc, ln_dyninc, ln_sshinc,                &
138         &                 ln_asmdin, ln_asmiau,                           &
139         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
140         &                 ln_salfix, salfixmin, nn_divdmp
141      !!----------------------------------------------------------------------
142
143      !-----------------------------------------------------------------------
144      ! Read Namelist nam_asminc : assimilation increment interface
145      !-----------------------------------------------------------------------
146      ln_seaiceinc   = .FALSE.
147      ln_temnofreeze = .FALSE.
148
149      REWIND( numnam_ref )              ! Namelist nam_asminc in reference namelist : Assimilation increment
150      READ  ( numnam_ref, nam_asminc, IOSTAT = ios, ERR = 901)
151901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_asminc in reference namelist', lwp )
152      REWIND( numnam_cfg )              ! Namelist nam_asminc in configuration namelist : Assimilation increment
153      READ  ( numnam_cfg, nam_asminc, IOSTAT = ios, ERR = 902 )
154902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_asminc in configuration namelist', lwp )
155      IF(lwm) WRITE ( numond, nam_asminc )
156
157      ! Control print
158      IF(lwp) THEN
159         WRITE(numout,*)
160         WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :'
161         WRITE(numout,*) '~~~~~~~~~~~~'
162         WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters'
163         WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri
164         WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc
165         WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc
166         WRITE(numout,*) '      Logical switch for applying SSH increments               ln_sshinc = ', ln_sshinc
167         WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin
168         WRITE(numout,*) '      Logical switch for applying sea ice increments        ln_seaiceinc = ', ln_seaiceinc
169         WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau
170         WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg
171         WRITE(numout,*) '      Timestep of background for DI in [0,nitend-nit000-1]     nitdin    = ', nitdin
172         WRITE(numout,*) '      Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr
173         WRITE(numout,*) '      Timestep of end of IAU interval in [0,nitend-nit000-1]   nitiaufin = ', nitiaufin
174         WRITE(numout,*) '      Type of IAU weighting function                           niaufn    = ', niaufn
175         WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix
176         WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin
177      ENDIF
178
179      nitbkg_r    = nitbkg    + nit000 - 1            ! Background time referenced to nit000
180      nitdin_r    = nitdin    + nit000 - 1            ! Background time for DI referenced to nit000
181      nitiaustr_r = nitiaustr + nit000 - 1            ! Start of IAU interval referenced to nit000
182      nitiaufin_r = nitiaufin + nit000 - 1            ! End of IAU interval referenced to nit000
183
184      iiauper     = nitiaufin_r - nitiaustr_r + 1     ! IAU interval length
185      icycper     = nitend      - nit000      + 1     ! Cycle interval length
186
187      CALL calc_date( nitend     , ditend_date    )   ! Date of final time step
188      CALL calc_date( nitbkg_r   , ditbkg_date    )   ! Background time for Jb referenced to ndate0
189      CALL calc_date( nitdin_r   , ditdin_date    )   ! Background time for DI referenced to ndate0
190      CALL calc_date( nitiaustr_r, ditiaustr_date )   ! IAU start time referenced to ndate0
191      CALL calc_date( nitiaufin_r, ditiaufin_date )   ! IAU end time referenced to ndate0
192
193      IF(lwp) THEN
194         WRITE(numout,*)
195         WRITE(numout,*) '   Time steps referenced to current cycle:'
196         WRITE(numout,*) '       iitrst      = ', nit000 - 1
197         WRITE(numout,*) '       nit000      = ', nit000
198         WRITE(numout,*) '       nitend      = ', nitend
199         WRITE(numout,*) '       nitbkg_r    = ', nitbkg_r
200         WRITE(numout,*) '       nitdin_r    = ', nitdin_r
201         WRITE(numout,*) '       nitiaustr_r = ', nitiaustr_r
202         WRITE(numout,*) '       nitiaufin_r = ', nitiaufin_r
203         WRITE(numout,*)
204         WRITE(numout,*) '   Dates referenced to current cycle:'
205         WRITE(numout,*) '       ndastp         = ', ndastp
206         WRITE(numout,*) '       ndate0         = ', ndate0
207         WRITE(numout,*) '       nn_time0       = ', nn_time0
208         WRITE(numout,*) '       ditend_date    = ', ditend_date
209         WRITE(numout,*) '       ditbkg_date    = ', ditbkg_date
210         WRITE(numout,*) '       ditdin_date    = ', ditdin_date
211         WRITE(numout,*) '       ditiaustr_date = ', ditiaustr_date
212         WRITE(numout,*) '       ditiaufin_date = ', ditiaufin_date
213      ENDIF
214
215
216      IF ( ( ln_asmdin ).AND.( ln_asmiau ) )   &
217         & CALL ctl_stop( ' ln_asmdin and ln_asmiau :', &
218         &                ' Choose Direct Initialization OR Incremental Analysis Updating')
219
220      IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
221           .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) &
222         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', &
223         &                ' but ln_asmdin and ln_asmiau are both set to .false. :', &
224         &                ' Inconsistent options')
225
226      IF ( ( niaufn /= 0 ).AND.( niaufn /= 1 ) ) &
227         & CALL ctl_stop( ' niaufn /= 0 or niaufn /=1 :',  &
228         &                ' Type IAU weighting function is invalid')
229
230      IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) &
231         &                     )  &
232         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', &
233         &                ' The assimilation increments are not applied')
234
235      IF ( ( ln_asmiau ).AND.( nitiaustr == nitiaufin ) ) &
236         & CALL ctl_stop( ' nitiaustr = nitiaufin :',  &
237         &                ' IAU interval is of zero length')
238
239      IF ( ( ln_asmiau ).AND.( ( nitiaustr_r < nit000 ).OR.( nitiaufin_r > nitend ) ) ) &
240         & CALL ctl_stop( ' nitiaustr or nitiaufin :',  &
241         &                ' IAU starting or final time step is outside the cycle interval', &
242         &                 ' Valid range nit000 to nitend')
243
244      IF ( ( nitbkg_r < nit000 - 1 ).OR.( nitbkg_r > nitend ) ) &
245         & CALL ctl_stop( ' nitbkg :',  &
246         &                ' Background time step is outside the cycle interval')
247
248      IF ( ( nitdin_r < nit000 - 1 ).OR.( nitdin_r > nitend ) ) &
249         & CALL ctl_stop( ' nitdin :',  &
250         &                ' Background time step for Direct Initialization is outside', &
251         &                ' the cycle interval')
252
253      IF ( nstop > 0 ) RETURN       ! if there are any errors then go no further
254
255      !--------------------------------------------------------------------
256      ! Initialize the Incremental Analysis Updating weighting function
257      !--------------------------------------------------------------------
258
259      IF( ln_asmiau ) THEN
260         !
261         ALLOCATE( wgtiau( icycper ) )
262         !
263         wgtiau(:) = 0._wp
264         !
265         !                                !---------------------------------------------------------
266         IF( niaufn == 0 ) THEN           ! Constant IAU forcing
267            !                             !---------------------------------------------------------
268            DO jt = 1, iiauper
269               wgtiau(jt+nitiaustr-1) = 1.0 / REAL( iiauper )
270            END DO
271            !                             !---------------------------------------------------------
272         ELSEIF ( niaufn == 1 ) THEN      ! Linear hat-like, centred in middle of IAU interval
273            !                             !---------------------------------------------------------
274            ! Compute the normalization factor
275            znorm = 0._wp
276            IF( MOD( iiauper, 2 ) == 0 ) THEN   ! Even number of time steps in IAU interval
277               imid = iiauper / 2 
278               DO jt = 1, imid
279                  znorm = znorm + REAL( jt )
280               END DO
281               znorm = 2.0 * znorm
282            ELSE                                ! Odd number of time steps in IAU interval
283               imid = ( iiauper + 1 ) / 2       
284               DO jt = 1, imid - 1
285                  znorm = znorm + REAL( jt )
286               END DO
287               znorm = 2.0 * znorm + REAL( imid )
288            ENDIF
289            znorm = 1.0 / znorm
290            !
291            DO jt = 1, imid - 1
292               wgtiau(jt+nitiaustr-1) = REAL( jt ) * znorm
293            END DO
294            DO jt = imid, iiauper
295               wgtiau(jt+nitiaustr-1) = REAL( iiauper - jt + 1 ) * znorm
296            END DO
297            !
298         ENDIF
299
300         ! Test that the integral of the weights over the weighting interval equals 1
301          IF(lwp) THEN
302             WRITE(numout,*)
303             WRITE(numout,*) 'asm_inc_init : IAU weights'
304             WRITE(numout,*) '~~~~~~~~~~~~'
305             WRITE(numout,*) '             time step         IAU  weight'
306             WRITE(numout,*) '             =========     ====================='
307             ztotwgt = 0.0
308             DO jt = 1, icycper
309                ztotwgt = ztotwgt + wgtiau(jt)
310                WRITE(numout,*) '         ', jt, '       ', wgtiau(jt) 
311             END DO   
312             WRITE(numout,*) '         ==================================='
313             WRITE(numout,*) '         Time-integrated weight = ', ztotwgt
314             WRITE(numout,*) '         ==================================='
315          ENDIF
316         
317      ENDIF
318
319      !--------------------------------------------------------------------
320      ! Allocate and initialize the increment arrays
321      !--------------------------------------------------------------------
322
323      ALLOCATE( t_bkginc     (jpi,jpj,jpk) )   ;   t_bkginc     (:,:,:) = 0._wp
324      ALLOCATE( s_bkginc     (jpi,jpj,jpk) )   ;   s_bkginc     (:,:,:) = 0._wp
325      ALLOCATE( u_bkginc     (jpi,jpj,jpk) )   ;   u_bkginc     (:,:,:) = 0._wp
326      ALLOCATE( v_bkginc     (jpi,jpj,jpk) )   ;   v_bkginc     (:,:,:) = 0._wp
327      ALLOCATE( ssh_bkginc   (jpi,jpj)     )   ;   ssh_bkginc   (:,:)   = 0._wp
328      ALLOCATE( seaice_bkginc(jpi,jpj)     )   ;   seaice_bkginc(:,:)   = 0._wp
329#if defined key_asminc
330      ALLOCATE( ssh_iau      (jpi,jpj)     )   ;   ssh_iau      (:,:)   = 0._wp
331#endif
332#if defined key_cice && defined key_asminc
333      ALLOCATE( ndaice_da    (jpi,jpj)     )   ;   ndaice_da    (:,:)   = 0._wp
334#endif
335      !
336      IF ( ln_trainc .OR. ln_dyninc .OR.   &       !--------------------------------------
337         & ln_sshinc .OR. ln_seaiceinc   ) THEN    ! Read the increments from file
338         !                                         !--------------------------------------
339         CALL iom_open( c_asminc, inum )
340         !
341         CALL iom_get( inum, 'time'       , zdate_inc   ) 
342         CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb )
343         CALL iom_get( inum, 'z_inc_datef', z_inc_datef )
344         z_inc_dateb = zdate_inc
345         z_inc_datef = zdate_inc
346         !
347         IF(lwp) THEN
348            WRITE(numout,*) 
349            WRITE(numout,*) 'asm_inc_init : Assimilation increments valid between dates ', z_inc_dateb,' and ', z_inc_datef
350            WRITE(numout,*) '~~~~~~~~~~~~'
351         ENDIF
352         !
353         IF ( ( z_inc_dateb < ndastp + nn_time0*0.0001_wp ) .OR.   &
354            & ( z_inc_datef > ditend_date ) ) &
355            &    CALL ctl_warn( ' Validity time of assimilation increments is ', &
356            &                   ' outside the assimilation interval' )
357
358         IF ( ( ln_asmdin ).AND.( zdate_inc /= ditdin_date ) ) &
359            & CALL ctl_warn( ' Validity time of assimilation increments does ', &
360            &                ' not agree with Direct Initialization time' )
361
362         IF ( ln_trainc ) THEN   
363            CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 )
364            CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 )
365            ! Apply the masks
366            t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:)
367            s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:)
368            ! Set missing increments to 0.0 rather than 1e+20
369            ! to allow for differences in masks
370            WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0
371            WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0
372         ENDIF
373
374         IF ( ln_dyninc ) THEN   
375            CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )             
376            CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )             
377            ! Apply the masks
378            u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:)
379            v_bkginc(:,:,:) = v_bkginc(:,:,:) * vmask(:,:,:)
380            ! Set missing increments to 0.0 rather than 1e+20
381            ! to allow for differences in masks
382            WHERE( ABS( u_bkginc(:,:,:) ) > 1.0e+10 ) u_bkginc(:,:,:) = 0.0
383            WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0
384         ENDIF
385       
386         IF ( ln_sshinc ) THEN
387            CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 )
388            ! Apply the masks
389            ssh_bkginc(:,:) = ssh_bkginc(:,:) * tmask(:,:,1)
390            ! Set missing increments to 0.0 rather than 1e+20
391            ! to allow for differences in masks
392            WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0
393         ENDIF
394
395         IF ( ln_seaiceinc ) THEN
396            CALL iom_get( inum, jpdom_autoglo, 'bckinseaice', seaice_bkginc, 1 )
397            ! Apply the masks
398            seaice_bkginc(:,:) = seaice_bkginc(:,:) * tmask(:,:,1)
399            ! Set missing increments to 0.0 rather than 1e+20
400            ! to allow for differences in masks
401            WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0
402         ENDIF
403         !
404         CALL iom_close( inum )
405         !
406      ENDIF
407      !
408      !                                            !--------------------------------------
409      IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN    ! Apply divergence damping filter
410         !                                         !--------------------------------------
411         ALLOCATE( zhdiv(jpi,jpj) ) 
412         !
413         DO jt = 1, nn_divdmp
414            !
415            DO jk = 1, jpkm1           ! zhdiv = e1e1 * div
416               zhdiv(:,:) = 0._wp
417               DO jj = 2, jpjm1
418                  DO ji = fs_2, fs_jpim1   ! vector opt.
419                     zhdiv(ji,jj) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * u_bkginc(ji  ,jj,jk)    &
420                        &            - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * u_bkginc(ji-1,jj,jk)    &
421                        &            + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * v_bkginc(ji,jj  ,jk)    &
422                        &            - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk)  ) / e3t(ji,jj,jk,Kmm)
423                  END DO
424               END DO
425               CALL lbc_lnk( 'asminc', zhdiv, 'T', 1. )   ! lateral boundary cond. (no sign change)
426               !
427               DO jj = 2, jpjm1
428                  DO ji = fs_2, fs_jpim1   ! vector opt.
429                     u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk)                         &
430                        &               + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji  ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
431                     v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk)                         &
432                        &               + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj  ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 
433                  END DO
434               END DO
435            END DO
436            !
437         END DO
438         !
439         DEALLOCATE( zhdiv ) 
440         !
441      ENDIF
442      !
443      !                             !-----------------------------------------------------
444      IF ( ln_asmdin ) THEN         ! Allocate and initialize the background state arrays
445         !                          !-----------------------------------------------------
446         !
447         ALLOCATE( t_bkg  (jpi,jpj,jpk) )   ;   t_bkg  (:,:,:) = 0._wp
448         ALLOCATE( s_bkg  (jpi,jpj,jpk) )   ;   s_bkg  (:,:,:) = 0._wp
449         ALLOCATE( u_bkg  (jpi,jpj,jpk) )   ;   u_bkg  (:,:,:) = 0._wp
450         ALLOCATE( v_bkg  (jpi,jpj,jpk) )   ;   v_bkg  (:,:,:) = 0._wp
451         ALLOCATE( ssh_bkg(jpi,jpj)     )   ;   ssh_bkg(:,:)   = 0._wp
452         !
453         !
454         !--------------------------------------------------------------------
455         ! Read from file the background state at analysis time
456         !--------------------------------------------------------------------
457         !
458         CALL iom_open( c_asmdin, inum )
459         !
460         CALL iom_get( inum, 'rdastp', zdate_bkg ) 
461         !
462         IF(lwp) THEN
463            WRITE(numout,*) 
464            WRITE(numout,*) '   ==>>>  Assimilation background state valid at : ', zdate_bkg
465            WRITE(numout,*)
466         ENDIF
467         !
468         IF ( zdate_bkg /= ditdin_date )   &
469            & CALL ctl_warn( ' Validity time of assimilation background state does', &
470            &                ' not agree with Direct Initialization time' )
471         !
472         IF ( ln_trainc ) THEN   
473            CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg )
474            CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg )
475            t_bkg(:,:,:) = t_bkg(:,:,:) * tmask(:,:,:)
476            s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:)
477         ENDIF
478         !
479         IF ( ln_dyninc ) THEN   
480            CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg )
481            CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg )
482            u_bkg(:,:,:) = u_bkg(:,:,:) * umask(:,:,:)
483            v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
484         ENDIF
485         !
486         IF ( ln_sshinc ) THEN
487            CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg )
488            ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1)
489         ENDIF
490         !
491         CALL iom_close( inum )
492         !
493      ENDIF
494      !
495      IF(lwp) WRITE(numout,*) '   ==>>>   Euler time step switch is ', neuler
496      !
497      IF( lk_asminc ) THEN                            !==  data assimilation  ==!
498         IF( ln_bkgwri )   CALL asm_bkg_wri( nit000 - 1, Kmm )      ! Output background fields
499         IF( ln_asmdin ) THEN                                  ! Direct initialization
500            IF( ln_trainc )   CALL tra_asm_inc( nit000 - 1, Kbb, Kmm, ts    , Krhs )      ! Tracers
501            IF( ln_dyninc )   CALL dyn_asm_inc( nit000 - 1, Kbb, Kmm, uu, vv, Krhs )      ! Dynamics
502            IF( ln_sshinc )   CALL ssh_asm_inc( nit000 - 1, Kbb, Kmm )                    ! SSH
503         ENDIF
504      ENDIF
505      !
506   END SUBROUTINE asm_inc_init
507   
508   
509   SUBROUTINE tra_asm_inc( kt, Kbb, Kmm, pts, Krhs )
510      !!----------------------------------------------------------------------
511      !!                    ***  ROUTINE tra_asm_inc  ***
512      !!         
513      !! ** Purpose : Apply the tracer (T and S) assimilation increments
514      !!
515      !! ** Method  : Direct initialization or Incremental Analysis Updating
516      !!
517      !! ** Action  :
518      !!----------------------------------------------------------------------
519      INTEGER                                  , INTENT(in   ) :: kt             ! Current time step
520      INTEGER                                  , INTENT(in   ) :: Kbb, Kmm, Krhs ! Time level indices
521      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers and RHS of tracer equation
522      !
523      INTEGER  :: ji, jj, jk
524      INTEGER  :: it
525      REAL(wp) :: zincwgt  ! IAU weight for current time step
526      REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values
527      !!----------------------------------------------------------------------
528      !
529      ! freezing point calculation taken from oc_fz_pt (but calculated for all depths)
530      ! used to prevent the applied increments taking the temperature below the local freezing point
531      DO jk = 1, jpkm1
532        CALL eos_fzp( pts(:,:,jk,jp_sal,Kmm), fzptnz(:,:,jk), gdept(:,:,jk,Kmm) )
533      END DO
534         !
535         !                             !--------------------------------------
536      IF ( ln_asmiau ) THEN            ! Incremental Analysis Updating
537         !                             !--------------------------------------
538         !
539         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
540            !
541            it = kt - nit000 + 1
542            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
543            !
544            IF(lwp) THEN
545               WRITE(numout,*) 
546               WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
547               WRITE(numout,*) '~~~~~~~~~~~~'
548            ENDIF
549            !
550            ! Update the tracer tendencies
551            DO jk = 1, jpkm1
552               IF (ln_temnofreeze) THEN
553                  ! Do not apply negative increments if the temperature will fall below freezing
554                  WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. &
555                     &   pts(:,:,jk,jp_tem,Kmm) + pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 
556                     pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 
557                  END WHERE
558               ELSE
559                  pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 
560               ENDIF
561               IF (ln_salfix) THEN
562                  ! Do not apply negative increments if the salinity will fall below a specified
563                  ! minimum value salfixmin
564                  WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. &
565                     &   pts(:,:,jk,jp_sal,Kmm) + pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 
566                     pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt
567                  END WHERE
568               ELSE
569                  pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt
570               ENDIF
571            END DO
572            !
573         ENDIF
574         !
575         IF ( kt == nitiaufin_r + 1  ) THEN   ! For bias crcn to work
576            DEALLOCATE( t_bkginc )
577            DEALLOCATE( s_bkginc )
578         ENDIF
579         !                             !--------------------------------------
580      ELSEIF ( ln_asmdin ) THEN        ! Direct Initialization
581         !                             !--------------------------------------
582         !           
583         IF ( kt == nitdin_r ) THEN
584            !
585            neuler = 0  ! Force Euler forward step
586            !
587            ! Initialize the now fields with the background + increment
588            IF (ln_temnofreeze) THEN
589               ! Do not apply negative increments if the temperature will fall below freezing
590               WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_tem,Kmm) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 
591                  pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
592               END WHERE
593            ELSE
594               pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
595            ENDIF
596            IF (ln_salfix) THEN
597               ! Do not apply negative increments if the salinity will fall below a specified
598               ! minimum value salfixmin
599               WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_sal,Kmm) + s_bkginc(:,:,:) > salfixmin ) 
600                  pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
601               END WHERE
602            ELSE
603               pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
604            ENDIF
605
606            pts(:,:,:,:,Kbb) = pts(:,:,:,:,Kmm)                 ! Update before fields
607
608            CALL eos( pts(:,:,:,:,Kbb), rhd, rhop, gdept_0(:,:,:) )  ! Before potential and in situ densities
609!!gm  fabien
610!            CALL eos( pts(:,:,:,:,Kbb), rhd, rhop )                ! Before potential and in situ densities
611!!gm
612
613            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)           &
614               &  CALL zps_hde    ( kt, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv,        &  ! Partial steps: before horizontal gradient
615               &                              rhd, gru , grv               )  ! of t, s, rd at the last ocean level
616            IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)                       &
617               &  CALL zps_hde_isf( nit000, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, gtui, gtvi,    &  ! Partial steps for top cell (ISF)
618               &                                  rhd, gru , grv , grui, grvi          )  ! of t, s, rd at the last ocean level
619
620            DEALLOCATE( t_bkginc )
621            DEALLOCATE( s_bkginc )
622            DEALLOCATE( t_bkg    )
623            DEALLOCATE( s_bkg    )
624         ENDIF
625         
626      ENDIF
627      ! Perhaps the following call should be in step
628      IF ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment
629      !
630   END SUBROUTINE tra_asm_inc
631
632
633   SUBROUTINE dyn_asm_inc( kt, Kbb, Kmm, puu, pvv, Krhs )
634      !!----------------------------------------------------------------------
635      !!                    ***  ROUTINE dyn_asm_inc  ***
636      !!         
637      !! ** Purpose : Apply the dynamics (u and v) assimilation increments.
638      !!
639      !! ** Method  : Direct initialization or Incremental Analysis Updating.
640      !!
641      !! ** Action  :
642      !!----------------------------------------------------------------------
643      INTEGER                             , INTENT( in )  ::  kt             ! ocean time-step index
644      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs ! ocean time level indices
645      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv       ! ocean velocities and RHS of momentum equation
646      !
647      INTEGER :: jk
648      INTEGER :: it
649      REAL(wp) :: zincwgt  ! IAU weight for current time step
650      !!----------------------------------------------------------------------
651      !
652      !                          !--------------------------------------------
653      IF ( ln_asmiau ) THEN      ! Incremental Analysis Updating
654         !                       !--------------------------------------------
655         !
656         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
657            !
658            it = kt - nit000 + 1
659            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
660            !
661            IF(lwp) THEN
662               WRITE(numout,*) 
663               WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
664               WRITE(numout,*) '~~~~~~~~~~~~'
665            ENDIF
666            !
667            ! Update the dynamic tendencies
668            DO jk = 1, jpkm1
669               puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + u_bkginc(:,:,jk) * zincwgt
670               pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + v_bkginc(:,:,jk) * zincwgt
671            END DO
672            !
673            IF ( kt == nitiaufin_r ) THEN
674               DEALLOCATE( u_bkginc )
675               DEALLOCATE( v_bkginc )
676            ENDIF
677            !
678         ENDIF
679         !                          !-----------------------------------------
680      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
681         !                          !-----------------------------------------
682         !         
683         IF ( kt == nitdin_r ) THEN
684            !
685            neuler = 0                    ! Force Euler forward step
686            !
687            ! Initialize the now fields with the background + increment
688            puu(:,:,:,Kmm) = u_bkg(:,:,:) + u_bkginc(:,:,:)
689            pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:) 
690            !
691            puu(:,:,:,Kbb) = puu(:,:,:,Kmm)         ! Update before fields
692            pvv(:,:,:,Kbb) = pvv(:,:,:,Kmm)
693            !
694            DEALLOCATE( u_bkg    )
695            DEALLOCATE( v_bkg    )
696            DEALLOCATE( u_bkginc )
697            DEALLOCATE( v_bkginc )
698         ENDIF
699         !
700      ENDIF
701      !
702   END SUBROUTINE dyn_asm_inc
703
704
705   SUBROUTINE ssh_asm_inc( kt, Kbb, Kmm )
706      !!----------------------------------------------------------------------
707      !!                    ***  ROUTINE ssh_asm_inc  ***
708      !!         
709      !! ** Purpose : Apply the sea surface height assimilation increment.
710      !!
711      !! ** Method  : Direct initialization or Incremental Analysis Updating.
712      !!
713      !! ** Action  :
714      !!----------------------------------------------------------------------
715      INTEGER, INTENT(IN) :: kt         ! Current time step
716      INTEGER, INTENT(IN) :: Kbb, Kmm   ! Current time step
717      !
718      INTEGER :: it
719      INTEGER :: jk
720      REAL(wp) :: zincwgt  ! IAU weight for current time step
721      !!----------------------------------------------------------------------
722      !
723      !                             !-----------------------------------------
724      IF ( ln_asmiau ) THEN         ! Incremental Analysis Updating
725         !                          !-----------------------------------------
726         !
727         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
728            !
729            it = kt - nit000 + 1
730            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
731            !
732            IF(lwp) THEN
733               WRITE(numout,*) 
734               WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
735                  &  kt,' with IAU weight = ', wgtiau(it)
736               WRITE(numout,*) '~~~~~~~~~~~~'
737            ENDIF
738            !
739            ! Save the tendency associated with the IAU weighted SSH increment
740            ! (applied in dynspg.*)
741#if defined key_asminc
742            ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
743#endif
744            !
745         ELSE IF( kt == nitiaufin_r+1 ) THEN
746            !
747            ! test on ssh_bkginc needed as ssh_asm_inc is called twice by time step
748            IF (ALLOCATED(ssh_bkginc)) DEALLOCATE( ssh_bkginc )
749            !
750#if defined key_asminc
751            ssh_iau(:,:) = 0._wp
752#endif
753            !
754         ENDIF
755         !                          !-----------------------------------------
756      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
757         !                          !-----------------------------------------
758         !
759         IF ( kt == nitdin_r ) THEN
760            !
761            neuler = 0                                   ! Force Euler forward step
762            !
763            ssh(:,:,Kmm) = ssh_bkg(:,:) + ssh_bkginc(:,:)   ! Initialize the now fields the background + increment
764            !
765            ssh(:,:,Kbb) = ssh(:,:,Kmm)                        ! Update before fields
766            e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
767!!gm why not e3u(:,:,:,Kbb), e3v(:,:,:,Kbb), gdept(:,:,:,Kbb) ????
768            !
769            DEALLOCATE( ssh_bkg    )
770            DEALLOCATE( ssh_bkginc )
771            !
772         ENDIF
773         !
774      ENDIF
775      !
776   END SUBROUTINE ssh_asm_inc
777
778
779   SUBROUTINE ssh_asm_div( kt, Kbb, Kmm, phdivn )
780      !!----------------------------------------------------------------------
781      !!                  ***  ROUTINE ssh_asm_div  ***
782      !!
783      !! ** Purpose :   ssh increment with z* is incorporated via a correction of the local divergence         
784      !!                across all the water column
785      !!
786      !! ** Method  :
787      !!                CAUTION : sshiau is positive (inflow) decreasing the
788      !!                          divergence and expressed in m/s
789      !!
790      !! ** Action  :   phdivn   decreased by the ssh increment
791      !!----------------------------------------------------------------------
792      INTEGER, INTENT(IN) :: kt                               ! ocean time-step index
793      INTEGER, INTENT(IN) :: Kbb, Kmm                         ! time level indices
794      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
795      !!
796      INTEGER  ::   jk                                        ! dummy loop index
797      REAL(wp), DIMENSION(:,:)  , POINTER       ::   ztim     ! local array
798      !!----------------------------------------------------------------------
799      !
800#if defined key_asminc
801      CALL ssh_asm_inc( kt, Kbb, Kmm ) !==   (calculate increments)
802      !
803      IF( ln_linssh ) THEN
804         phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t(:,:,1,Kmm) * tmask(:,:,1)
805      ELSE
806         ALLOCATE( ztim(jpi,jpj) )
807         ztim(:,:) = ssh_iau(:,:) / ( ht(:,:) + 1.0 - ssmask(:,:) )
808         DO jk = 1, jpkm1                                 
809            phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) 
810         END DO
811         !
812         DEALLOCATE(ztim)
813      ENDIF
814#endif
815      !
816   END SUBROUTINE ssh_asm_div
817
818
819   SUBROUTINE seaice_asm_inc( kt, kindic )
820      !!----------------------------------------------------------------------
821      !!                    ***  ROUTINE seaice_asm_inc  ***
822      !!         
823      !! ** Purpose : Apply the sea ice assimilation increment.
824      !!
825      !! ** Method  : Direct initialization or Incremental Analysis Updating.
826      !!
827      !! ** Action  :
828      !!
829      !!----------------------------------------------------------------------
830      INTEGER, INTENT(in)           ::   kt       ! Current time step
831      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation
832      !
833      INTEGER  ::   it
834      REAL(wp) ::   zincwgt   ! IAU weight for current time step
835#if defined key_si3
836      REAL(wp), DIMENSION(jpi,jpj) ::   zofrld, zohicif, zseaicendg, zhicifinc
837      REAL(wp) ::   zhicifmin = 0.5_wp      ! ice minimum depth in metres
838#endif
839      !!----------------------------------------------------------------------
840      !
841      !                             !-----------------------------------------
842      IF ( ln_asmiau ) THEN         ! Incremental Analysis Updating
843         !                          !-----------------------------------------
844         !
845         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
846            !
847            it = kt - nit000 + 1
848            zincwgt = wgtiau(it)      ! IAU weight for the current time step
849            ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments)
850            !
851            IF(lwp) THEN
852               WRITE(numout,*) 
853               WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
854               WRITE(numout,*) '~~~~~~~~~~~~'
855            ENDIF
856            !
857            ! Sea-ice : SI3 case
858            !
859#if defined key_si3
860            zofrld (:,:) = 1._wp - at_i(:,:)
861            zohicif(:,:) = hm_i(:,:)
862            !
863            at_i  (:,:) = 1. - MIN( MAX( 1.-at_i  (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
864            at_i_b(:,:) = 1. - MIN( MAX( 1.-at_i_b(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
865            fr_i(:,:) = at_i(:,:)        ! adjust ice fraction
866            !
867            zseaicendg(:,:) = zofrld(:,:) - (1. - at_i(:,:))   ! find out actual sea ice nudge applied
868            !
869            ! Nudge sea ice depth to bring it up to a required minimum depth
870            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 
871               zhicifinc(:,:) = (zhicifmin - hm_i(:,:)) * zincwgt   
872            ELSEWHERE
873               zhicifinc(:,:) = 0.0_wp
874            END WHERE
875            !
876            ! nudge ice depth
877            hm_i (:,:) = hm_i (:,:) + zhicifinc(:,:)
878            !
879            ! seaice salinity balancing (to add)
880#endif
881            !
882#if defined key_cice && defined key_asminc
883            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
884            ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt
885#endif
886            !
887            IF ( kt == nitiaufin_r ) THEN
888               DEALLOCATE( seaice_bkginc )
889            ENDIF
890            !
891         ELSE
892            !
893#if defined key_cice && defined key_asminc
894            ndaice_da(:,:) = 0._wp        ! Sea-ice : CICE case. Zero ice increment tendency into CICE
895#endif
896            !
897         ENDIF
898         !                          !-----------------------------------------
899      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
900         !                          !-----------------------------------------
901         !
902         IF ( kt == nitdin_r ) THEN
903            !
904            neuler = 0                    ! Force Euler forward step
905            !
906            ! Sea-ice : SI3 case
907            !
908#if defined key_si3
909            zofrld (:,:) = 1._wp - at_i(:,:)
910            zohicif(:,:) = hm_i(:,:)
911            !
912            ! Initialize the now fields the background + increment
913            at_i(:,:) = 1. - MIN( MAX( 1.-at_i(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp)
914            at_i_b(:,:) = at_i(:,:) 
915            fr_i(:,:) = at_i(:,:)        ! adjust ice fraction
916            !
917            zseaicendg(:,:) = zofrld(:,:) - (1. - at_i(:,:))   ! find out actual sea ice nudge applied
918            !
919            ! Nudge sea ice depth to bring it up to a required minimum depth
920            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 
921               zhicifinc(:,:) = zhicifmin - hm_i(:,:)
922            ELSEWHERE
923               zhicifinc(:,:) = 0.0_wp
924            END WHERE
925            !
926            ! nudge ice depth
927            hm_i (:,:) = hm_i (:,:) + zhicifinc(:,:)
928            !
929            ! seaice salinity balancing (to add)
930#endif
931            !
932#if defined key_cice && defined key_asminc
933            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
934           ndaice_da(:,:) = seaice_bkginc(:,:) / rdt
935#endif
936            IF ( .NOT. PRESENT(kindic) ) THEN
937               DEALLOCATE( seaice_bkginc )
938            END IF
939            !
940         ELSE
941            !
942#if defined key_cice && defined key_asminc
943            ndaice_da(:,:) = 0._wp     ! Sea-ice : CICE case. Zero ice increment tendency into CICE
944#endif
945            !
946         ENDIF
947
948!#if defined defined key_si3 || defined key_cice
949!
950!            IF (ln_seaicebal ) THEN       
951!             !! balancing salinity increments
952!             !! simple case from limflx.F90 (doesn't include a mass flux)
953!             !! assumption is that as ice concentration is reduced or increased
954!             !! the snow and ice depths remain constant
955!             !! note that snow is being created where ice concentration is being increased
956!             !! - could be more sophisticated and
957!             !! not do this (but would need to alter h_snow)
958!
959!             usave(:,:,:)=sb(:,:,:)   ! use array as a temporary store
960!
961!             DO jj = 1, jpj
962!               DO ji = 1, jpi
963!           ! calculate change in ice and snow mass per unit area
964!           ! positive values imply adding salt to the ocean (results from ice formation)
965!           ! fwf : ice formation and melting
966!
967!                 zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rdt
968!
969!           ! change salinity down to mixed layer depth
970!                 mld=hmld_kara(ji,jj)
971!
972!           ! prevent small mld
973!           ! less than 10m can cause salinity instability
974!                 IF (mld < 10) mld=10
975!
976!           ! set to bottom of a level
977!                 DO jk = jpk-1, 2, -1
978!                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN
979!                     mld=gdepw(ji,jj,jk+1)
980!                     jkmax=jk
981!                   ENDIF
982!                 ENDDO
983!
984!            ! avoid applying salinity balancing in shallow water or on land
985!            !
986!
987!            ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m)
988!
989!                 dsal_ocn=0.0_wp
990!                 sal_thresh=5.0_wp        ! minimum salinity threshold for salinity balancing
991!
992!                 if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) &
993!                              dsal_ocn = zfons / (rhop(ji,jj,1) * mld)
994!
995!           ! put increments in for levels in the mixed layer
996!           ! but prevent salinity below a threshold value
997!
998!                   DO jk = 1, jkmax             
999!
1000!                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN
1001!                           sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn
1002!                           sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn
1003!                     ENDIF
1004!
1005!                   ENDDO
1006!
1007!      !            !  salt exchanges at the ice/ocean interface
1008!      !            zpmess         = zfons / rdt_ice    ! rdt_ice is ice timestep
1009!      !
1010!      !! Adjust fsalt. A +ve fsalt means adding salt to ocean
1011!      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt 
1012!      !!               
1013!      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d)
1014!      !!                                                     ! E-P (kg m-2 s-2)
1015!      !            emp(ji,jj) = emp(ji,jj) + zpmess          ! E-P (kg m-2 s-2)
1016!               ENDDO !ji
1017!             ENDDO !jj!
1018!
1019!            ENDIF !ln_seaicebal
1020!
1021!#endif
1022         !
1023      ENDIF
1024      !
1025   END SUBROUTINE seaice_asm_inc
1026   
1027   !!======================================================================
1028END MODULE asminc
Note: See TracBrowser for help on using the repository browser.