source: branches/2017/dev_r8789_sbc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 8832

Last change on this file since 8832 was 8832, checked in by charris, 3 years ago

Fix for compilation problem without key_asminc.

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