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

Last change on this file since 6140 was 6140, checked in by timgraham, 5 years ago

Merge of branches/2015/dev_merge_2015 back into trunk. Merge excludes NEMOGCM/TOOLS/OBSTOOLS/ for now due to issues with the change of file type. Will sort these manually with further commits.

Branch merged as follows:
In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
Small conflicts due to bug fixes applied to trunk since the dev_merge_2015 was copied. Bug fixes were applied to the branch as well so these were easy to resolve.
Branch committed at this stage

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge —reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2015/dev_merge_2015
to merge the branch into the trunk and then commit - no conflicts at this stage.

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