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

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

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

Last change on this file since 9023 was 9023, checked in by timgraham, 6 years ago

Merged METO_MERCATOR branch and resolved all conflicts in OPA_SRC

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