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/UKMO/dev_r5518_GO6_package_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/dev_r5518_GO6_package_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 10005

Last change on this file since 10005 was 10005, checked in by timgraham, 4 years ago

Included all of the functional changes from STARTHOUR branch in branches/2014

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