source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 9213

Last change on this file since 9213 was 9213, checked in by gm, 3 years ago

dev_merge_2017: nemogcm.F90 : updated in SAS & OFF + data assimilation initial calls (asm_bkg_wri , tra_asm_inc …) moved to asm_inc_init + closed sea : restructure namcfg & its control print + set ln_closea = false if domcfg file not read (ln_domcfg=F

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