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

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

Fix for direct initialization case (although note that this is untested in the vvl case and may still contain bugs).

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