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/branches/UKMO/NEMO_4.0.4_FOAM_IAU_SI3_SIC/src/OCE/ASM – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_FOAM_IAU_SI3_SIC/src/OCE/ASM/asminc.F90 @ 15412

Last change on this file since 15412 was 15412, checked in by dcarneir, 9 months ago

Final tested version of the SIC IAU code in GO8

File size: 53.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 phycst         ! physical constants
39   USE ice1D          ! sea-ice: thermodynamics variables
40   USE icetab         ! sea-ice: 3D <==> 2D, 2D <==> 1D
41   USE ice
42#endif
43   !
44   USE in_out_manager  ! I/O manager
45   USE iom             ! Library to read input files
46   USE lib_mpp         ! MPP library
47
48   IMPLICIT NONE
49   PRIVATE
50   
51   PUBLIC   asm_inc_init   !: Initialize the increment arrays and IAU weights
52   PUBLIC   tra_asm_inc    !: Apply the tracer (T and S) increments
53   PUBLIC   dyn_asm_inc    !: Apply the dynamic (u and v) increments
54   PUBLIC   ssh_asm_inc    !: Apply the SSH increment
55   PUBLIC   ssh_asm_div    !: Apply the SSH divergence
56   PUBLIC   seaice_asm_inc !: Apply the seaice increment
57
58#if defined key_asminc
59    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .TRUE.   !: Logical switch for assimilation increment interface
60#else
61    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE.  !: No assimilation increments
62#endif
63   LOGICAL, PUBLIC :: ln_bkgwri     !: No output of the background state fields
64   LOGICAL, PUBLIC :: ln_asmiau     !: No applying forcing with an assimilation increment
65   LOGICAL, PUBLIC :: ln_asmdin     !: No direct initialization
66   LOGICAL, PUBLIC :: ln_trainc     !: No tracer (T and S) assimilation increments
67   LOGICAL, PUBLIC :: ln_dyninc     !: No dynamics (u and v) assimilation increments
68   LOGICAL, PUBLIC :: ln_sshinc     !: No sea surface height assimilation increment
69   LOGICAL, PUBLIC :: ln_seaiceinc  !: No sea ice concentration increment
70   LOGICAL, PUBLIC :: ln_salfix     !: Apply minimum salinity check
71   LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
72   INTEGER, PUBLIC :: nn_divdmp     !: Apply divergence damping filter nn_divdmp times
73
74   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkg   , s_bkg      !: Background temperature and salinity
75   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkg   , v_bkg      !: Background u- & v- velocity components
76   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkginc, s_bkginc   !: Increment to the background T & S
77   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components
78   REAL(wp), PUBLIC, DIMENSION(:)    , ALLOCATABLE ::   wgtiau               !: IAU weights for each time step
79#if defined key_asminc
80   REAL(wp), PUBLIC, DIMENSION(:,:)  , ALLOCATABLE ::   ssh_iau              !: IAU-weighted sea surface height increment
81#endif
82   !                                !!! time steps relative to the cycle interval [0,nitend-nit000-1]
83   INTEGER , PUBLIC ::   nitbkg      !: Time step of the background state used in the Jb term
84   INTEGER , PUBLIC ::   nitdin      !: Time step of the background state for direct initialization
85   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval
86   INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval
87   !
88   INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting
89   !                                 !: = 1   Linear hat-like, centred in middle of IAU interval
90   REAL(wp), PUBLIC ::   salfixmin   !: Ensure that the salinity is larger than this  value if (ln_salfix)
91
92   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment
93   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc
94#if defined key_si3
95   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   a_i_bkginc          ! Increment to the background sea ice conc categories
96   LOGICAL, PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: incr_newice    ! Mask .TRUE.=DA positive ice increment to open water
97   REAL(wp) :: zhi_damin                                            ! Ice thickness for new sea ice from DA increment
98#endif
99#if defined key_cice && defined key_asminc
100   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ndaice_da             ! ice increment tendency into CICE
101#endif
102
103   !! * Substitutions
104#  include "vectopt_loop_substitute.h90"
105   !!----------------------------------------------------------------------
106   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
107   !! $Id$
108   !! Software governed by the CeCILL license (see ./LICENSE)
109   !!----------------------------------------------------------------------
110CONTAINS
111
112   SUBROUTINE asm_inc_init
113      !!----------------------------------------------------------------------
114      !!                    ***  ROUTINE asm_inc_init  ***
115      !!         
116      !! ** Purpose : Initialize the assimilation increment and IAU weights.
117      !!
118      !! ** Method  : Initialize the assimilation increment and IAU weights.
119      !!
120      !! ** Action  :
121      !!----------------------------------------------------------------------
122      INTEGER :: ji, jj, jk, jt, jl  ! dummy loop indices
123      INTEGER :: imid, inum          ! local integers
124      INTEGER :: ios                 ! Local integer output status for namelist read
125      INTEGER :: iiauper             ! Number of time steps in the IAU period
126      INTEGER :: icycper             ! Number of time steps in the cycle
127      REAL(KIND=dp) :: ditend_date     ! Date YYYYMMDD.HHMMSS of final time step
128      REAL(KIND=dp) :: ditbkg_date     ! Date YYYYMMDD.HHMMSS of background time step for Jb term
129      REAL(KIND=dp) :: ditdin_date     ! Date YYYYMMDD.HHMMSS of background time step for DI
130      REAL(KIND=dp) :: ditiaustr_date  ! Date YYYYMMDD.HHMMSS of IAU interval start time step
131      REAL(KIND=dp) :: ditiaufin_date  ! Date YYYYMMDD.HHMMSS of IAU interval final time step
132
133      REAL(wp) :: znorm        ! Normalization factor for IAU weights
134      REAL(wp) :: ztotwgt      ! Value of time-integrated IAU weights (should be equal to one)
135      REAL(wp) :: z_inc_dateb  ! Start date of interval on which increment is valid
136      REAL(wp) :: z_inc_datef  ! End date of interval on which increment is valid
137      REAL(wp) :: zdate_bkg    ! Date in background state file for DI
138      REAL(wp) :: zdate_inc    ! Time axis in increments file
139      !
140      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zhdiv   ! 2D workspace
141      REAL(wp) :: zremaining_increment
142
143      !!
144      NAMELIST/nam_asminc/ ln_bkgwri,                                      &
145         &                 ln_trainc, ln_dyninc, ln_sshinc,                &
146         &                 ln_seaiceinc, ln_asmdin, ln_asmiau,             &
147         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
148         &                 ln_salfix, salfixmin, nn_divdmp, zhi_damin
149      !!----------------------------------------------------------------------
150
151      !-----------------------------------------------------------------------
152      ! Read Namelist nam_asminc : assimilation increment interface
153      !-----------------------------------------------------------------------
154      ln_seaiceinc   = .FALSE.
155      ln_temnofreeze = .FALSE.
156
157      REWIND( numnam_ref )              ! Namelist nam_asminc in reference namelist : Assimilation increment
158      READ  ( numnam_ref, nam_asminc, IOSTAT = ios, ERR = 901)
159901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_asminc in reference namelist' )
160      REWIND( numnam_cfg )              ! Namelist nam_asminc in configuration namelist : Assimilation increment
161      READ  ( numnam_cfg, nam_asminc, IOSTAT = ios, ERR = 902 )
162902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_asminc in configuration namelist' )
163      IF(lwm) WRITE ( numond, nam_asminc )
164
165      ! Control print
166      IF(lwp) THEN
167         WRITE(numout,*)
168         WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :'
169         WRITE(numout,*) '~~~~~~~~~~~~'
170         WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters'
171         WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri
172         WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc
173         WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc
174         WRITE(numout,*) '      Logical switch for applying SSH increments               ln_sshinc = ', ln_sshinc
175         WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin
176         WRITE(numout,*) '      Logical switch for applying sea ice increments        ln_seaiceinc = ', ln_seaiceinc
177         WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau
178         WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg
179         WRITE(numout,*) '      Timestep of background for DI in [0,nitend-nit000-1]     nitdin    = ', nitdin
180         WRITE(numout,*) '      Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr
181         WRITE(numout,*) '      Timestep of end of IAU interval in [0,nitend-nit000-1]   nitiaufin = ', nitiaufin
182         WRITE(numout,*) '      Type of IAU weighting function                           niaufn    = ', niaufn
183         WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix
184         WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin
185         WRITE(numout,*) '      Minimum ice thickness for new ice from DA                zhi_damin = ', zhi_damin
186      ENDIF
187
188      nitbkg_r    = nitbkg    + nit000 - 1            ! Background time referenced to nit000
189      nitdin_r    = nitdin    + nit000 - 1            ! Background time for DI referenced to nit000
190      nitiaustr_r = nitiaustr + nit000 - 1            ! Start of IAU interval referenced to nit000
191      nitiaufin_r = nitiaufin + nit000 - 1            ! End of IAU interval referenced to nit000
192
193      iiauper     = nitiaufin_r - nitiaustr_r + 1     ! IAU interval length
194      icycper     = nitend      - nit000      + 1     ! Cycle interval length
195
196      CALL calc_date( nitend     , ditend_date    )   ! Date of final time step
197      CALL calc_date( nitbkg_r   , ditbkg_date    )   ! Background time for Jb referenced to ndate0
198      CALL calc_date( nitdin_r   , ditdin_date    )   ! Background time for DI referenced to ndate0
199      CALL calc_date( nitiaustr_r, ditiaustr_date )   ! IAU start time referenced to ndate0
200      CALL calc_date( nitiaufin_r, ditiaufin_date )   ! IAU end time referenced to ndate0
201
202      IF(lwp) THEN
203         WRITE(numout,*)
204         WRITE(numout,*) '   Time steps referenced to current cycle:'
205         WRITE(numout,*) '       iitrst      = ', nit000 - 1
206         WRITE(numout,*) '       nit000      = ', nit000
207         WRITE(numout,*) '       nitend      = ', nitend
208         WRITE(numout,*) '       nitbkg_r    = ', nitbkg_r
209         WRITE(numout,*) '       nitdin_r    = ', nitdin_r
210         WRITE(numout,*) '       nitiaustr_r = ', nitiaustr_r
211         WRITE(numout,*) '       nitiaufin_r = ', nitiaufin_r
212         WRITE(numout,*)
213         WRITE(numout,*) '   Dates referenced to current cycle:'
214         WRITE(numout,*) '       ndastp         = ', ndastp
215         WRITE(numout,*) '       ndate0         = ', ndate0
216         WRITE(numout,*) '       nn_time0       = ', nn_time0
217         WRITE(numout,*) '       ditend_date    = ', ditend_date
218         WRITE(numout,*) '       ditbkg_date    = ', ditbkg_date
219         WRITE(numout,*) '       ditdin_date    = ', ditdin_date
220         WRITE(numout,*) '       ditiaustr_date = ', ditiaustr_date
221         WRITE(numout,*) '       ditiaufin_date = ', ditiaufin_date
222      ENDIF
223
224
225      IF ( ( ln_asmdin ).AND.( ln_asmiau ) )   &
226         & CALL ctl_stop( ' ln_asmdin and ln_asmiau :', &
227         &                ' Choose Direct Initialization OR Incremental Analysis Updating')
228
229      IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
230           .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) &
231         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', &
232         &                ' but ln_asmdin and ln_asmiau are both set to .false. :', &
233         &                ' Inconsistent options')
234
235      IF ( ( niaufn /= 0 ).AND.( niaufn /= 1 ) ) &
236         & CALL ctl_stop( ' niaufn /= 0 or niaufn /=1 :',  &
237         &                ' Type IAU weighting function is invalid')
238
239      IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) &
240         &                     )  &
241         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', &
242         &                ' The assimilation increments are not applied')
243
244      IF ( ( ln_asmiau ).AND.( nitiaustr == nitiaufin ) ) &
245         & CALL ctl_stop( ' nitiaustr = nitiaufin :',  &
246         &                ' IAU interval is of zero length')
247
248      IF ( ( ln_asmiau ).AND.( ( nitiaustr_r < nit000 ).OR.( nitiaufin_r > nitend ) ) ) &
249         & CALL ctl_stop( ' nitiaustr or nitiaufin :',  &
250         &                ' IAU starting or final time step is outside the cycle interval', &
251         &                 ' Valid range nit000 to nitend')
252
253      IF ( ( nitbkg_r < nit000 - 1 ).OR.( nitbkg_r > nitend ) ) &
254         & CALL ctl_stop( ' nitbkg :',  &
255         &                ' Background time step is outside the cycle interval')
256
257      IF ( ( nitdin_r < nit000 - 1 ).OR.( nitdin_r > nitend ) ) &
258         & CALL ctl_stop( ' nitdin :',  &
259         &                ' Background time step for Direct Initialization is outside', &
260         &                ' the cycle interval')
261
262      IF ( nstop > 0 ) RETURN       ! if there are any errors then go no further
263
264      !--------------------------------------------------------------------
265      ! Initialize the Incremental Analysis Updating weighting function
266      !--------------------------------------------------------------------
267
268      IF( ln_asmiau ) THEN
269         !
270         ALLOCATE( wgtiau( icycper ) )
271         !
272         wgtiau(:) = 0._wp
273         !
274         !                                !---------------------------------------------------------
275         IF( niaufn == 0 ) THEN           ! Constant IAU forcing
276            !                             !---------------------------------------------------------
277            DO jt = 1, iiauper
278               wgtiau(jt+nitiaustr-1) = 1.0 / REAL( iiauper )
279            END DO
280            !                             !---------------------------------------------------------
281         ELSEIF ( niaufn == 1 ) THEN      ! Linear hat-like, centred in middle of IAU interval
282            !                             !---------------------------------------------------------
283            ! Compute the normalization factor
284            znorm = 0._wp
285            IF( MOD( iiauper, 2 ) == 0 ) THEN   ! Even number of time steps in IAU interval
286               imid = iiauper / 2 
287               DO jt = 1, imid
288                  znorm = znorm + REAL( jt )
289               END DO
290               znorm = 2.0 * znorm
291            ELSE                                ! Odd number of time steps in IAU interval
292               imid = ( iiauper + 1 ) / 2       
293               DO jt = 1, imid - 1
294                  znorm = znorm + REAL( jt )
295               END DO
296               znorm = 2.0 * znorm + REAL( imid )
297            ENDIF
298            znorm = 1.0 / znorm
299            !
300            DO jt = 1, imid - 1
301               wgtiau(jt+nitiaustr-1) = REAL( jt ) * znorm
302            END DO
303            DO jt = imid, iiauper
304               wgtiau(jt+nitiaustr-1) = REAL( iiauper - jt + 1 ) * znorm
305            END DO
306            !
307         ENDIF
308
309         ! Test that the integral of the weights over the weighting interval equals 1
310          IF(lwp) THEN
311             WRITE(numout,*)
312             WRITE(numout,*) 'asm_inc_init : IAU weights'
313             WRITE(numout,*) '~~~~~~~~~~~~'
314             WRITE(numout,*) '             time step         IAU  weight'
315             WRITE(numout,*) '             =========     ====================='
316             ztotwgt = 0.0
317             DO jt = 1, icycper
318                ztotwgt = ztotwgt + wgtiau(jt)
319                WRITE(numout,*) '         ', jt, '       ', wgtiau(jt) 
320             END DO   
321             WRITE(numout,*) '         ==================================='
322             WRITE(numout,*) '         Time-integrated weight = ', ztotwgt
323             WRITE(numout,*) '         ==================================='
324          ENDIF
325         
326      ENDIF
327
328      !--------------------------------------------------------------------
329      ! Allocate and initialize the increment arrays
330      !--------------------------------------------------------------------
331
332      ALLOCATE( t_bkginc     (jpi,jpj,jpk) )   ;   t_bkginc     (:,:,:) = 0._wp
333      ALLOCATE( s_bkginc     (jpi,jpj,jpk) )   ;   s_bkginc     (:,:,:) = 0._wp
334      ALLOCATE( u_bkginc     (jpi,jpj,jpk) )   ;   u_bkginc     (:,:,:) = 0._wp
335      ALLOCATE( v_bkginc     (jpi,jpj,jpk) )   ;   v_bkginc     (:,:,:) = 0._wp
336      ALLOCATE( ssh_bkginc   (jpi,jpj)     )   ;   ssh_bkginc   (:,:)   = 0._wp
337      ALLOCATE( seaice_bkginc(jpi,jpj)     )   ;   seaice_bkginc(:,:)   = 0._wp
338#if defined key_si3
339      ALLOCATE( a_i_bkginc   (jpi,jpj,jpl) )   ;   a_i_bkginc   (:,:,:) = 0._wp
340      ALLOCATE( incr_newice(jpi,jpj,jpl) )     ;   incr_newice(:,:,:) = .FALSE.
341#endif
342#if defined key_asminc
343      ALLOCATE( ssh_iau      (jpi,jpj)     )   ;   ssh_iau      (:,:)   = 0._wp
344#endif
345#if defined key_cice && defined key_asminc
346      ALLOCATE( ndaice_da    (jpi,jpj)     )   ;   ndaice_da    (:,:)   = 0._wp
347#endif
348      !
349      IF ( ln_trainc .OR. ln_dyninc .OR.   &       !--------------------------------------
350         & ln_sshinc .OR. ln_seaiceinc   ) THEN    ! Read the increments from file
351         !                                         !--------------------------------------
352         CALL iom_open( c_asminc, inum )
353         !
354         CALL iom_get( inum, 'time'       , zdate_inc   ) 
355         CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb )
356         CALL iom_get( inum, 'z_inc_datef', z_inc_datef )
357         z_inc_dateb = zdate_inc
358         z_inc_datef = zdate_inc
359         !
360         IF(lwp) THEN
361            WRITE(numout,*) 
362            WRITE(numout,*) 'asm_inc_init : Assimilation increments valid between dates ', z_inc_dateb,' and ', z_inc_datef
363            WRITE(numout,*) '~~~~~~~~~~~~'
364         ENDIF
365         !
366         IF ( ( z_inc_dateb < ndastp + nn_time0*0.0001_wp ) .OR.   &
367            & ( z_inc_datef > ditend_date ) ) &
368            &    CALL ctl_warn( ' Validity time of assimilation increments is ', &
369            &                   ' outside the assimilation interval' )
370
371         IF ( ( ln_asmdin ).AND.( zdate_inc /= ditdin_date ) ) &
372            & CALL ctl_warn( ' Validity time of assimilation increments does ', &
373            &                ' not agree with Direct Initialization time' )
374
375         IF ( ln_trainc ) THEN   
376            CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 )
377            CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 )
378            ! Apply the masks
379            t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:)
380            s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:)
381            ! Set missing increments to 0.0 rather than 1e+20
382            ! to allow for differences in masks
383            WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0
384            WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0
385         ENDIF
386
387         IF ( ln_dyninc ) THEN   
388            CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )             
389            CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )             
390            ! Apply the masks
391            u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:)
392            v_bkginc(:,:,:) = v_bkginc(:,:,:) * vmask(:,:,:)
393            ! Set missing increments to 0.0 rather than 1e+20
394            ! to allow for differences in masks
395            WHERE( ABS( u_bkginc(:,:,:) ) > 1.0e+10 ) u_bkginc(:,:,:) = 0.0
396            WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0
397         ENDIF
398       
399         IF ( ln_sshinc ) THEN
400            CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 )
401            ! Apply the masks
402            ssh_bkginc(:,:) = ssh_bkginc(:,:) * tmask(:,:,1)
403            ! Set missing increments to 0.0 rather than 1e+20
404            ! to allow for differences in masks
405            WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0
406         ENDIF
407
408         IF ( ln_seaiceinc ) THEN
409            CALL iom_get( inum, jpdom_autoglo, 'bckinseaice', seaice_bkginc, 1 )
410            ! Apply the masks
411            seaice_bkginc(:,:) = seaice_bkginc(:,:) * tmask(:,:,1)
412            ! Set missing increments to 0.0 rather than 1e+20
413            ! to allow for differences in masks
414            ! very small increments are also set to zero
415            WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 .OR. &
416              &    ABS( seaice_bkginc(:,:) ) < epsi03 ) seaice_bkginc(:,:) = 0.0_wp
417           
418            IF (lwp) THEN
419               WRITE(numout,*)
420               WRITE(numout,*) 'asm_inc_init : Converting seaice_bkginc to a_i_bkginc'
421               WRITE(numout,*) '~~~~~~~~~~~~'
422            END IF
423
424            !single category increment for sea ice conc
425            !convert single category increment to multi category
426            a_i_bkginc = 0.0_wp
427            at_i = SUM( a_i(:,:,:), DIM=3 )
428
429            ! ensure zhi_damin lies within 1st category
430            IF ( zhi_damin > hi_max(1) ) zhi_damin = hi_max(1) - epsi02
431           
432            DO jj = 1, jpj
433               DO ji = 1, jpi
434                  IF ( seaice_bkginc(ji,jj) > 0.0_wp) THEN
435                     !positive ice concentration increments are always
436                     !added to the thinnest category of ice
437                     a_i_bkginc(ji,jj,1) = seaice_bkginc(ji,jj)
438                  ELSE
439                     !negative increments are first removed from the thinnest
440                     !available category until it reaches zero concentration
441                     !and then progressively removed from thicker categories
442                     zremaining_increment = seaice_bkginc(ji,jj)
443                     DO jl = 1, jpl
444                        ! assign as much increment as possible to current category
445                        a_i_bkginc(ji,jj,jl) = -MIN( a_i(ji,jj,jl), -zremaining_increment )
446                        ! update remaining amount of increment
447                        zremaining_increment = zremaining_increment - a_i_bkginc(ji,jj,jl)
448                     END DO
449                  END IF
450               END DO
451            END DO
452            ! find model points where DA new ice should be added to open water
453            DO jl = 1,jpl
454               WHERE ( at_i(:,:) < epsi02 .AND. seaice_bkginc(:,:) > 0.0_wp )
455                  incr_newice(:,:,jl) = .TRUE.
456               END WHERE
457            END DO
458         ENDIF
459         !
460         CALL iom_close( inum )
461         !
462      ENDIF
463      !
464      !                                            !--------------------------------------
465      IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN    ! Apply divergence damping filter
466         !                                         !--------------------------------------
467         ALLOCATE( zhdiv(jpi,jpj) ) 
468         !
469         DO jt = 1, nn_divdmp
470            !
471            DO jk = 1, jpkm1           ! zhdiv = e1e1 * div
472               zhdiv(:,:) = 0._wp
473               DO jj = 2, jpjm1
474                  DO ji = fs_2, fs_jpim1   ! vector opt.
475                     zhdiv(ji,jj) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * u_bkginc(ji  ,jj,jk)    &
476                        &            - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * u_bkginc(ji-1,jj,jk)    &
477                        &            + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * v_bkginc(ji,jj  ,jk)    &
478                        &            - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * v_bkginc(ji,jj-1,jk)  ) / e3t_n(ji,jj,jk)
479                  END DO
480               END DO
481               CALL lbc_lnk( 'asminc', zhdiv, 'T', 1. )   ! lateral boundary cond. (no sign change)
482               !
483               DO jj = 2, jpjm1
484                  DO ji = fs_2, fs_jpim1   ! vector opt.
485                     u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk)                         &
486                        &               + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji  ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
487                     v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk)                         &
488                        &               + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj  ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 
489                  END DO
490               END DO
491            END DO
492            !
493         END DO
494         !
495         DEALLOCATE( zhdiv ) 
496         !
497      ENDIF
498      !
499      !                             !-----------------------------------------------------
500      IF ( ln_asmdin ) THEN         ! Allocate and initialize the background state arrays
501         !                          !-----------------------------------------------------
502         !
503         ALLOCATE( t_bkg  (jpi,jpj,jpk) )   ;   t_bkg  (:,:,:) = 0._wp
504         ALLOCATE( s_bkg  (jpi,jpj,jpk) )   ;   s_bkg  (:,:,:) = 0._wp
505         ALLOCATE( u_bkg  (jpi,jpj,jpk) )   ;   u_bkg  (:,:,:) = 0._wp
506         ALLOCATE( v_bkg  (jpi,jpj,jpk) )   ;   v_bkg  (:,:,:) = 0._wp
507         ALLOCATE( ssh_bkg(jpi,jpj)     )   ;   ssh_bkg(:,:)   = 0._wp
508         !
509         !
510         !--------------------------------------------------------------------
511         ! Read from file the background state at analysis time
512         !--------------------------------------------------------------------
513         !
514         CALL iom_open( c_asmdin, inum )
515         !
516         CALL iom_get( inum, 'rdastp', zdate_bkg ) 
517         !
518         IF(lwp) THEN
519            WRITE(numout,*) 
520            WRITE(numout,*) '   ==>>>  Assimilation background state valid at : ', zdate_bkg
521            WRITE(numout,*)
522         ENDIF
523         !
524         IF ( zdate_bkg /= ditdin_date )   &
525            & CALL ctl_warn( ' Validity time of assimilation background state does', &
526            &                ' not agree with Direct Initialization time' )
527         !
528         IF ( ln_trainc ) THEN   
529            CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg )
530            CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg )
531            t_bkg(:,:,:) = t_bkg(:,:,:) * tmask(:,:,:)
532            s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:)
533         ENDIF
534         !
535         IF ( ln_dyninc ) THEN   
536            CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg )
537            CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg )
538            u_bkg(:,:,:) = u_bkg(:,:,:) * umask(:,:,:)
539            v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
540         ENDIF
541         !
542         IF ( ln_sshinc ) THEN
543            CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg )
544            ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1)
545         ENDIF
546         !
547         CALL iom_close( inum )
548         !
549      ENDIF
550      !
551      IF(lwp) WRITE(numout,*) '   ==>>>   Euler time step switch is ', neuler
552      !
553      IF( lk_asminc ) THEN                            !==  data assimilation  ==!
554         IF( ln_bkgwri )   CALL asm_bkg_wri( nit000 - 1 )      ! Output background fields
555         IF( ln_asmdin ) THEN                                  ! Direct initialization
556            IF( ln_trainc )   CALL tra_asm_inc( nit000 - 1 )      ! Tracers
557            IF( ln_dyninc )   CALL dyn_asm_inc( nit000 - 1 )      ! Dynamics
558            IF( ln_sshinc )   CALL ssh_asm_inc( nit000 - 1 )      ! SSH
559         ENDIF
560      ENDIF
561      !
562   END SUBROUTINE asm_inc_init
563   
564   
565   SUBROUTINE tra_asm_inc( kt )
566      !!----------------------------------------------------------------------
567      !!                    ***  ROUTINE tra_asm_inc  ***
568      !!         
569      !! ** Purpose : Apply the tracer (T and S) assimilation increments
570      !!
571      !! ** Method  : Direct initialization or Incremental Analysis Updating
572      !!
573      !! ** Action  :
574      !!----------------------------------------------------------------------
575      INTEGER, INTENT(IN) ::   kt   ! Current time step
576      !
577      INTEGER  :: ji, jj, jk
578      INTEGER  :: it
579      REAL(wp) :: zincwgt  ! IAU weight for current time step
580      REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values
581      !!----------------------------------------------------------------------
582      !
583      ! freezing point calculation taken from oc_fz_pt (but calculated for all depths)
584      ! used to prevent the applied increments taking the temperature below the local freezing point
585      DO jk = 1, jpkm1
586        CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), gdept_n(:,:,jk) )
587      END DO
588         !
589         !                             !--------------------------------------
590      IF ( ln_asmiau ) THEN            ! Incremental Analysis Updating
591         !                             !--------------------------------------
592         !
593         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
594            !
595            it = kt - nit000 + 1
596            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
597            !
598            IF(lwp) THEN
599               WRITE(numout,*) 
600               WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
601               WRITE(numout,*) '~~~~~~~~~~~~'
602            ENDIF
603            !
604            ! Update the tracer tendencies
605            DO jk = 1, jpkm1
606               IF (ln_temnofreeze) THEN
607                  ! Do not apply negative increments if the temperature will fall below freezing
608                  WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. &
609                     &   tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 
610                     tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 
611                  END WHERE
612               ELSE
613                  tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 
614               ENDIF
615               IF (ln_salfix) THEN
616                  ! Do not apply negative increments if the salinity will fall below a specified
617                  ! minimum value salfixmin
618                  WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. &
619                     &   tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 
620                     tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
621                  END WHERE
622               ELSE
623                  tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
624               ENDIF
625            END DO
626            !
627         ENDIF
628         !
629         IF ( kt == nitiaufin_r + 1  ) THEN   ! For bias crcn to work
630            DEALLOCATE( t_bkginc )
631            DEALLOCATE( s_bkginc )
632         ENDIF
633         !                             !--------------------------------------
634      ELSEIF ( ln_asmdin ) THEN        ! Direct Initialization
635         !                             !--------------------------------------
636         !           
637         IF ( kt == nitdin_r ) THEN
638            !
639            neuler = 0  ! Force Euler forward step
640            !
641            ! Initialize the now fields with the background + increment
642            IF (ln_temnofreeze) THEN
643               ! Do not apply negative increments if the temperature will fall below freezing
644               WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 
645                  tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
646               END WHERE
647            ELSE
648               tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
649            ENDIF
650            IF (ln_salfix) THEN
651               ! Do not apply negative increments if the salinity will fall below a specified
652               ! minimum value salfixmin
653               WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 
654                  tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
655               END WHERE
656            ELSE
657               tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
658            ENDIF
659
660            tsb(:,:,:,:) = tsn(:,:,:,:)                 ! Update before fields
661
662            CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )  ! Before potential and in situ densities
663!!gm  fabien
664!            CALL eos( tsb, rhd, rhop )                ! Before potential and in situ densities
665!!gm
666
667            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)      &
668               &  CALL zps_hde    ( kt, jpts, tsb, gtsu, gtsv,        &  ! Partial steps: before horizontal gradient
669               &                              rhd, gru , grv          )  ! of t, s, rd at the last ocean level
670            IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)      &
671               &  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi,    &    ! Partial steps for top cell (ISF)
672               &                                  rhd, gru , grv , grui, grvi     )    ! of t, s, rd at the last ocean level
673
674            DEALLOCATE( t_bkginc )
675            DEALLOCATE( s_bkginc )
676            DEALLOCATE( t_bkg    )
677            DEALLOCATE( s_bkg    )
678         ENDIF
679         
680      ENDIF
681      ! Perhaps the following call should be in step
682      IF ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment
683      !
684   END SUBROUTINE tra_asm_inc
685
686
687   SUBROUTINE dyn_asm_inc( kt )
688      !!----------------------------------------------------------------------
689      !!                    ***  ROUTINE dyn_asm_inc  ***
690      !!         
691      !! ** Purpose : Apply the dynamics (u and v) assimilation increments.
692      !!
693      !! ** Method  : Direct initialization or Incremental Analysis Updating.
694      !!
695      !! ** Action  :
696      !!----------------------------------------------------------------------
697      INTEGER, INTENT(IN) :: kt   ! Current time step
698      !
699      INTEGER :: jk
700      INTEGER :: it
701      REAL(wp) :: zincwgt  ! IAU weight for current time step
702      !!----------------------------------------------------------------------
703      !
704      !                          !--------------------------------------------
705      IF ( ln_asmiau ) THEN      ! Incremental Analysis Updating
706         !                       !--------------------------------------------
707         !
708         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
709            !
710            it = kt - nit000 + 1
711            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
712            !
713            IF(lwp) THEN
714               WRITE(numout,*) 
715               WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
716               WRITE(numout,*) '~~~~~~~~~~~~'
717            ENDIF
718            !
719            ! Update the dynamic tendencies
720            DO jk = 1, jpkm1
721               ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt
722               va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt
723            END DO
724            !
725            IF ( kt == nitiaufin_r ) THEN
726               DEALLOCATE( u_bkginc )
727               DEALLOCATE( v_bkginc )
728            ENDIF
729            !
730         ENDIF
731         !                          !-----------------------------------------
732      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
733         !                          !-----------------------------------------
734         !         
735         IF ( kt == nitdin_r ) THEN
736            !
737            neuler = 0                    ! Force Euler forward step
738            !
739            ! Initialize the now fields with the background + increment
740            un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:)
741            vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:) 
742            !
743            ub(:,:,:) = un(:,:,:)         ! Update before fields
744            vb(:,:,:) = vn(:,:,:)
745            !
746            DEALLOCATE( u_bkg    )
747            DEALLOCATE( v_bkg    )
748            DEALLOCATE( u_bkginc )
749            DEALLOCATE( v_bkginc )
750         ENDIF
751         !
752      ENDIF
753      !
754   END SUBROUTINE dyn_asm_inc
755
756
757   SUBROUTINE ssh_asm_inc( kt )
758      !!----------------------------------------------------------------------
759      !!                    ***  ROUTINE ssh_asm_inc  ***
760      !!         
761      !! ** Purpose : Apply the sea surface height assimilation increment.
762      !!
763      !! ** Method  : Direct initialization or Incremental Analysis Updating.
764      !!
765      !! ** Action  :
766      !!----------------------------------------------------------------------
767      INTEGER, INTENT(IN) :: kt   ! Current time step
768      !
769      INTEGER :: it
770      INTEGER :: jk
771      REAL(wp) :: zincwgt  ! IAU weight for current time step
772      !!----------------------------------------------------------------------
773      !
774      !                             !-----------------------------------------
775      IF ( ln_asmiau ) THEN         ! Incremental Analysis Updating
776         !                          !-----------------------------------------
777         !
778         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
779            !
780            it = kt - nit000 + 1
781            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
782            !
783            IF(lwp) THEN
784               WRITE(numout,*) 
785               WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
786                  &  kt,' with IAU weight = ', wgtiau(it)
787               WRITE(numout,*) '~~~~~~~~~~~~'
788            ENDIF
789            !
790            ! Save the tendency associated with the IAU weighted SSH increment
791            ! (applied in dynspg.*)
792#if defined key_asminc
793            ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
794#endif
795            !
796         ELSE IF( kt == nitiaufin_r+1 ) THEN
797            !
798            ! test on ssh_bkginc needed as ssh_asm_inc is called twice by time step
799            IF (ALLOCATED(ssh_bkginc)) DEALLOCATE( ssh_bkginc )
800            !
801#if defined key_asminc
802            ssh_iau(:,:) = 0._wp
803#endif
804            !
805         ENDIF
806         !                          !-----------------------------------------
807      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
808         !                          !-----------------------------------------
809         !
810         IF ( kt == nitdin_r ) THEN
811            !
812            neuler = 0                                   ! Force Euler forward step
813            !
814            sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)   ! Initialize the now fields the background + increment
815            !
816            sshb(:,:) = sshn(:,:)                        ! Update before fields
817            e3t_b(:,:,:) = e3t_n(:,:,:)
818!!gm why not e3u_b, e3v_b, gdept_b ????
819            !
820            DEALLOCATE( ssh_bkg    )
821            DEALLOCATE( ssh_bkginc )
822            !
823         ENDIF
824         !
825      ENDIF
826      !
827   END SUBROUTINE ssh_asm_inc
828
829
830   SUBROUTINE ssh_asm_div( kt, phdivn )
831      !!----------------------------------------------------------------------
832      !!                  ***  ROUTINE ssh_asm_div  ***
833      !!
834      !! ** Purpose :   ssh increment with z* is incorporated via a correction of the local divergence         
835      !!                across all the water column
836      !!
837      !! ** Method  :
838      !!                CAUTION : sshiau is positive (inflow) decreasing the
839      !!                          divergence and expressed in m/s
840      !!
841      !! ** Action  :   phdivn   decreased by the ssh increment
842      !!----------------------------------------------------------------------
843      INTEGER, INTENT(IN) :: kt                               ! ocean time-step index
844      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
845      !!
846      INTEGER  ::   jk                                        ! dummy loop index
847      REAL(wp), DIMENSION(:,:)  , POINTER       ::   ztim     ! local array
848      !!----------------------------------------------------------------------
849      !
850#if defined key_asminc
851      CALL ssh_asm_inc( kt ) !==   (calculate increments)
852      !
853      IF( ln_linssh ) THEN
854         phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t_n(:,:,1) * tmask(:,:,1)
855      ELSE
856         ALLOCATE( ztim(jpi,jpj) )
857         ztim(:,:) = ssh_iau(:,:) / ( ht_n(:,:) + 1.0 - ssmask(:,:) )
858         DO jk = 1, jpkm1                                 
859            phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) 
860         END DO
861         !
862         DEALLOCATE(ztim)
863      ENDIF
864#endif
865      !
866   END SUBROUTINE ssh_asm_div
867
868
869   SUBROUTINE seaice_asm_inc( kt, kindic )
870      !!----------------------------------------------------------------------
871      !!                    ***  ROUTINE seaice_asm_inc  ***
872      !!         
873      !! ** Purpose : Apply the sea ice assimilation increment.
874      !!
875      !! ** Method  : Direct initialization or Incremental Analysis Updating.
876      !!
877      !! ** Action  :
878      !!
879      !!----------------------------------------------------------------------
880      INTEGER, INTENT(in)           ::   kt       ! Current time step
881      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation
882      !
883      INTEGER  ::   it, jl
884      REAL(wp) ::   zincwgt                            ! IAU weight for current time step
885#if defined key_si3
886      LOGICAL, SAVE :: logical_newice = .TRUE.         ! Flag to add new ice from DA to open water
887      REAL(wp), DIMENSION(jpi,jpj,jpl) :: da_i         ! change in ice concentration
888      REAL(wp), DIMENSION(jpi,jpj,jpl) :: dv_i         ! change in ice volume
889      REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i       ! inverse of old ice concentration
890      REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_v_i       ! inverse of old ice volume
891      REAL(wp), DIMENSION(jpi,jpj)     :: zhi_damin_2D ! array with DA thickness for incr_newice
892#endif
893      !!----------------------------------------------------------------------
894      !
895      !                             !-----------------------------------------
896      IF ( ln_asmiau ) THEN         ! Incremental Analysis Updating
897         !                          !-----------------------------------------
898         !
899         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
900            !
901            it = kt - nit000 + 1
902            zincwgt = wgtiau(it)      ! IAU weight for the current time step
903            ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments)
904            !
905            IF(lwp) THEN
906               WRITE(numout,*) 
907               WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
908               WRITE(numout,*) '~~~~~~~~~~~~'
909            ENDIF
910            !
911            ! Sea-ice : SI3 case
912            !
913#if defined key_si3
914            ! ensure zhi_damin lies within 1st category
915            IF ( zhi_damin > hi_max(1) ) zhi_damin = hi_max(1) - epsi02
916
917            ! compute the inverse of key sea ice variables
918            ! to be used later in the code
919            WHERE( a_i(:,:,:) > epsi10 )
920               z1_a_i(:,:,:) = 1.0_wp / a_i(:,:,:)
921               z1_v_i(:,:,:) = 1.0_wp / v_i(:,:,:)
922            ELSEWHERE
923               z1_a_i(:,:,:) = 0.0_wp
924               z1_v_i(:,:,:) = 0.0_wp
925            END WHERE
926
927            ! add positive concentration increments to regions where ice
928            ! is already present and bound them to 1
929            ! ice volume is added based on zhi_damin
930            WHERE ( .NOT. incr_newice .AND. a_i_bkginc(:,:,:) > 0.0_wp )
931               a_i(:,:,:) = a_i(:,:,:) + MIN( 1.0_wp - a_i(:,:,:), a_i_bkginc(:,:,:) * zincwgt )
932               v_i(:,:,:) = v_i(:,:,:) + MIN( 1.0_wp - a_i(:,:,:), a_i_bkginc(:,:,:) * zincwgt ) * zhi_damin
933            END WHERE
934
935            ! add negative concentration increments to regions where ice
936            ! is already present and bound them to 0
937            ! in this case ice volume is changed based on the current thickness
938            WHERE ( .NOT. incr_newice .AND. a_i_bkginc(:,:,:) < 0.0_wp )
939               a_i(:,:,:) = MAX( a_i(:,:,:) + a_i_bkginc(:,:,:) * zincwgt, 0.0_wp )
940               v_i(:,:,:) = a_i(:,:,:) * h_i(:,:,:)
941            END WHERE
942           
943            ! compute changes in ice concentration and volume
944            WHERE ( incr_newice )
945               da_i(:,:,:) = 1.0_wp
946               dv_i(:,:,:) = 1.0_wp
947            ELSEWHERE
948               da_i(:,:,:) = a_i(:,:,:) * z1_a_i(:,:,:)
949               dv_i(:,:,:) = v_i(:,:,:) * z1_v_i(:,:,:)
950            END WHERE
951
952            ! initialise thermodynamics of new ice being added to open water
953            ! just do it once since next IAU steps assume that new ice has
954            ! already been added in
955            IF ( kt == nitiaustr_r .AND. logical_newice ) THEN
956
957               ! assign zhi_damin to ice forming in open water
958               WHERE ( ANY( incr_newice, DIM=3 ) )
959                  zhi_damin_2D(:,:) = zhi_damin
960               ELSEWHERE
961                  zhi_damin_2D(:,:) = 0.0_wp
962               END WHERE
963
964               ! add ice concentration and volume
965               ! ensure the other prognostic variables are set to zero
966               WHERE ( incr_newice )
967                  a_i(:,:,:) = MIN( 1.0_wp, a_i_bkginc(:,:,:) * zincwgt )
968                  v_i(:,:,:) = MIN( 1.0_wp, a_i_bkginc(:,:,:) * zincwgt ) * zhi_damin
969                  v_s (:,:,:) = 0.0_wp
970                  a_ip(:,:,:) = 0.0_wp
971                  v_ip(:,:,:) = 0.0_wp
972                  sv_i(:,:,:) = 0.0_wp
973               END WHERE
974               DO jl = 1, nlay_i
975                  WHERE ( incr_newice )
976                     e_i(:,:,jl,:) = 0.0_wp
977                  END WHERE
978               END DO
979               DO jl = 1, nlay_s
980                  WHERE ( incr_newice )
981                     e_s(:,:,jl,:) = 0.0_wp
982                  END WHERE
983               END DO
984
985               ! Initialisation of the salt content and ice enthalpy
986               ! set flag of new ice to false after this
987               call init_new_ice_thd( zhi_damin_2D )
988               incr_newice(:,:,:) = .FALSE.
989            END IF
990
991            ! maintain equivalent values for key prognostic variables
992            v_s(:,:,:) = v_s(:,:,:) * da_i(:,:,:)
993            DO jl = 1, nlay_s
994               e_s(:,:,jl,:) = e_s(:,:,jl,:) * da_i(:,:,:)
995            END DO
996            a_ip (:,:,:) = a_ip(:,:,:) * da_i(:,:,:)
997            v_ip (:,:,:) = v_ip(:,:,:) * da_i(:,:,:)
998 
999            ! ice volume dependent variables
1000            sv_i (:,:,:) = sv_i(:,:,:) * dv_i(:,:,:)
1001            DO jl = 1, nlay_i
1002               e_i(:,:,jl,:) = e_i(:,:,jl,:) * dv_i(:,:,:)
1003            END DO
1004#endif
1005
1006#if defined key_cice && defined key_asminc
1007            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
1008            ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt
1009#endif
1010            !
1011            IF ( kt == nitiaufin_r ) THEN
1012               DEALLOCATE( seaice_bkginc )
1013#if defined key_si3
1014               DEALLOCATE( incr_newice )
1015               DEALLOCATE( a_i_bkginc )
1016#endif
1017            ENDIF
1018            !
1019         ELSE
1020            !
1021#if defined key_cice && defined key_asminc
1022            ndaice_da(:,:) = 0._wp        ! Sea-ice : CICE case. Zero ice increment tendency into CICE
1023#endif
1024            !
1025         ENDIF
1026         !                          !-----------------------------------------
1027      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
1028         !                          !-----------------------------------------
1029         !
1030         IF ( kt == nitdin_r ) THEN
1031            !
1032            neuler = 0                    ! Force Euler forward step
1033
1034            IF(lwp) THEN
1035               WRITE(numout,*)
1036               WRITE(numout,*) 'seaice_asm_inc : sea ice direct initialization at time step = ', kt
1037               WRITE(numout,*) '~~~~~~~~~~~~'
1038            ENDIF
1039            !
1040#if defined key_cice && defined key_asminc
1041            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
1042           ndaice_da(:,:) = seaice_bkginc(:,:) / rdt
1043#endif
1044            IF ( .NOT. PRESENT(kindic) ) THEN
1045               DEALLOCATE( seaice_bkginc )
1046            END IF
1047            !
1048         ELSE
1049            !
1050#if defined key_cice && defined key_asminc
1051            ndaice_da(:,:) = 0._wp     ! Sea-ice : CICE case. Zero ice increment tendency into CICE
1052#endif
1053            !
1054         ENDIF
1055         !
1056      ENDIF
1057      !
1058   END SUBROUTINE seaice_asm_inc
1059
1060
1061   SUBROUTINE init_new_ice_thd( hi_new )
1062      !!----------------------------------------------------------------------
1063      !!                  ***  ROUTINE init_new_ice_thd  ***
1064      !!
1065      !! ** Purpose :   Initialise thermodynamics of new ice
1066      !!                forming at 1st category with thickness hi_new
1067      !!
1068      !! ** Method  :   Apply SI3 thermodynamic equations to initialise
1069      !!                thermodynamic of new ice
1070      !!
1071      !! ** Action  :   update thermodynamic variables
1072      !!----------------------------------------------------------------------
1073      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    ::   hi_new  ! total thickness of new ice
1074
1075      INTEGER :: jj, ji, jk, jl
1076      REAL(wp) ::   ztmelts   ! melting point
1077
1078      REAL(wp), DIMENSION(jpij) ::   zh_newice   ! 1d version of hi_new
1079      REAL(wp), DIMENSION(jpij) ::   zs_newice   ! salinity of new ice
1080      !!----------------------------------------------------------------------
1081
1082      ! Identify grid points where new ice forms
1083      npti = 0   ;   nptidx(:) = 0
1084      DO jj = 1, jpj
1085         DO ji = 1, jpi
1086            IF ( hi_new(ji,jj) > 0._wp ) THEN
1087               npti = npti + 1
1088               nptidx( npti ) = (jj - 1) * jpi + ji
1089            ENDIF
1090         END DO
1091      END DO
1092
1093      ! Move from 2-D to 1-D vectors
1094      IF ( npti > 0 ) THEN
1095         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
1096         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
1097         CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , hi_new        )
1098         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d    (1:npti) , sss_m         )
1099         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d   (1:npti) , t_bo          )
1100         DO jk = 1, nlay_i
1101            CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,1) )
1102         END DO
1103
1104         ! --- Salinity of new ice --- !
1105         SELECT CASE ( nn_icesal )
1106         CASE ( 1 )                    ! Sice = constant
1107            zs_newice(1:npti) = rn_icesal
1108         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
1109            DO ji = 1, npti
1110               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_1d(ji) )
1111            END DO
1112         CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
1113            zs_newice(1:npti) =   2.3
1114         END SELECT
1115
1116         ! --- Update ice salt content --- !
1117         DO ji = 1, npti
1118            sv_i_2d(ji,1) = sv_i_2d(ji,1) + zs_newice(ji) * ( v_i_2d(ji,1) )
1119         END DO
1120
1121         ! --- Heat content of new ice --- !
1122         ! We assume that new ice is formed at the seawater freezing point
1123         DO ji = 1, npti
1124               ztmelts       = - rTmlt * zs_newice(ji)                  ! Melting point (C)
1125               e_i_1d(ji,:)  =   rhoi * (  rcpi  * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     &
1126                  &                      + rLfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
1127                  &                      - rcp   *         ztmelts )
1128         END DO
1129
1130         ! Change units for e_i
1131         DO jk = 1, nlay_i
1132            e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) * v_i_2d(1:npti,1) * r1_nlay_i 
1133         END DO
1134
1135         ! Reforming full thermodynamic variables
1136         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
1137         DO jk = 1, nlay_i
1138               CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,1) )
1139         END DO
1140      END IF
1141
1142   END SUBROUTINE init_new_ice_thd
1143
1144
1145   !!======================================================================
1146END MODULE asminc
Note: See TracBrowser for help on using the repository browser.