source: trunk/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 8547

Last change on this file since 8547 was 8547, checked in by timgraham, 3 years ago

Correct typo in last commit for #1849

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