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

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

source: NEMO/trunk/src/OCE/ASM/asminc.F90 @ 13295

Last change on this file since 13295 was 13295, checked in by acc, 4 years ago

Replace do-loop macros in the trunk with alternative forms with greater flexibility for extra halo applications. This alters a lot of routines but does not change any behaviour or results. do_loop_substitute.h90 is greatly simplified by this change. SETTE results are identical to those with the previous revision

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