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

Last change on this file since 12489 was 12489, checked in by davestorkey, 6 months ago

Preparation for new timestepping scheme #2390.
Main changes:

  1. Initial euler timestep now handled in stp and not in TRA/DYN routines.
  2. Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
  3. Renaming of a few miscellaneous parameters, eg. atfp → rn_atfp (namelist parameter used everywhere) and rau0 → rho0.

This version gives bit-comparable results to the previous version of the trunk.

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