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_package/src/OCE/ASM – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_FOAM_package/src/OCE/ASM/asminc.F90 @ 15486

Last change on this file since 15486 was 15486, checked in by dcarneir, 8 months ago

SIC assimilation interface (IAU) for SI3

File size: 54.3 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 && defined key_asminc
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   REAL(wp) :: zhi_damin                                            ! Ice thickness for new sea ice from DA increment
95#if defined key_si3 && defined key_asminc
96   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   a_i_bkginc          ! Increment to the background sea ice conc categories
97   LOGICAL, PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: incr_newice    ! Mask .TRUE.=DA positive ice increment to open water
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 && defined key_asminc
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(:,:) ) < 1.0e-03_wp ) seaice_bkginc(:,:) = 0.0_wp
417           
418#if defined key_si3 && defined key_asminc
419            IF (lwp) THEN
420               WRITE(numout,*)
421               WRITE(numout,*) 'asm_inc_init : Converting single category increment to multi-category'
422               WRITE(numout,*) '~~~~~~~~~~~~'
423            END IF
424
425            !single category increment for sea ice conc
426            !convert single category increment to multi category
427            at_i = SUM( a_i(:,:,:), DIM=3 )
428
429            ! ensure zhi_damin lies within 1st category
430            IF ( zhi_damin > hi_max(1) ) THEN
431               IF (lwp) THEN
432                  WRITE(numout,*)
433                  WRITE(numout,*) 'minimum DA thickness > hmax(1): setting minimum DA thickness to hmax(1)'
434                  WRITE(numout,*) '~~~~~~~~~~~~'
435               END IF
436               zhi_damin = hi_max(1) - epsi02
437            END IF
438           
439            DO jj = 1, jpj
440               DO ji = 1, jpi
441                  IF ( seaice_bkginc(ji,jj) > 0.0_wp ) THEN
442                     !positive ice concentration increments are always
443                     !added to the thinnest category of ice
444                     a_i_bkginc(ji,jj,1) = seaice_bkginc(ji,jj)
445                  ELSE
446                     !negative increments are first removed from the thinnest
447                     !available category until it reaches zero concentration
448                     !and then progressively removed from thicker categories
449
450                     !NOTE: might consider the remote possibility that zremaining_increment
451                     !left over is not zero, particulary if SI3 is being run
452                     !as a single category
453                     zremaining_increment = seaice_bkginc(ji,jj)
454                     DO jl = 1, jpl
455                        ! assign as much increment as possible to current category
456                        a_i_bkginc(ji,jj,jl) = -MIN( a_i(ji,jj,jl), -zremaining_increment )
457                        ! update remaining amount of increment
458                        zremaining_increment = zremaining_increment - a_i_bkginc(ji,jj,jl)
459                     END DO
460                  END IF
461               END DO
462            END DO
463            ! find model points where DA new ice should be added to open water
464            ! in any ice category
465            DO jl = 1,jpl
466               WHERE ( at_i(:,:) < epsi02 .AND. seaice_bkginc(:,:) > 0.0_wp )
467                  incr_newice(:,:,jl) = .TRUE.
468               END WHERE
469            END DO
470         ENDIF
471#endif
472         !
473         CALL iom_close( inum )
474         !
475      ENDIF
476      !
477      !                                            !--------------------------------------
478      IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN    ! Apply divergence damping filter
479         !                                         !--------------------------------------
480         ALLOCATE( zhdiv(jpi,jpj) ) 
481         !
482         DO jt = 1, nn_divdmp
483            !
484            DO jk = 1, jpkm1           ! zhdiv = e1e1 * div
485               zhdiv(:,:) = 0._wp
486               DO jj = 2, jpjm1
487                  DO ji = fs_2, fs_jpim1   ! vector opt.
488                     zhdiv(ji,jj) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * u_bkginc(ji  ,jj,jk)    &
489                        &            - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * u_bkginc(ji-1,jj,jk)    &
490                        &            + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * v_bkginc(ji,jj  ,jk)    &
491                        &            - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * v_bkginc(ji,jj-1,jk)  ) / e3t_n(ji,jj,jk)
492                  END DO
493               END DO
494               CALL lbc_lnk( 'asminc', zhdiv, 'T', 1. )   ! lateral boundary cond. (no sign change)
495               !
496               DO jj = 2, jpjm1
497                  DO ji = fs_2, fs_jpim1   ! vector opt.
498                     u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk)                         &
499                        &               + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji  ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
500                     v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk)                         &
501                        &               + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj  ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 
502                  END DO
503               END DO
504            END DO
505            !
506         END DO
507         !
508         DEALLOCATE( zhdiv ) 
509         !
510      ENDIF
511      !
512      !                             !-----------------------------------------------------
513      IF ( ln_asmdin ) THEN         ! Allocate and initialize the background state arrays
514         !                          !-----------------------------------------------------
515         !
516         ALLOCATE( t_bkg  (jpi,jpj,jpk) )   ;   t_bkg  (:,:,:) = 0._wp
517         ALLOCATE( s_bkg  (jpi,jpj,jpk) )   ;   s_bkg  (:,:,:) = 0._wp
518         ALLOCATE( u_bkg  (jpi,jpj,jpk) )   ;   u_bkg  (:,:,:) = 0._wp
519         ALLOCATE( v_bkg  (jpi,jpj,jpk) )   ;   v_bkg  (:,:,:) = 0._wp
520         ALLOCATE( ssh_bkg(jpi,jpj)     )   ;   ssh_bkg(:,:)   = 0._wp
521         !
522         !
523         !--------------------------------------------------------------------
524         ! Read from file the background state at analysis time
525         !--------------------------------------------------------------------
526         !
527         CALL iom_open( c_asmdin, inum )
528         !
529         CALL iom_get( inum, 'rdastp', zdate_bkg ) 
530         !
531         IF(lwp) THEN
532            WRITE(numout,*) 
533            WRITE(numout,*) '   ==>>>  Assimilation background state valid at : ', zdate_bkg
534            WRITE(numout,*)
535         ENDIF
536         !
537         IF ( zdate_bkg /= ditdin_date )   &
538            & CALL ctl_warn( ' Validity time of assimilation background state does', &
539            &                ' not agree with Direct Initialization time' )
540         !
541         IF ( ln_trainc ) THEN   
542            CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg )
543            CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg )
544            t_bkg(:,:,:) = t_bkg(:,:,:) * tmask(:,:,:)
545            s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:)
546         ENDIF
547         !
548         IF ( ln_dyninc ) THEN   
549            CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg )
550            CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg )
551            u_bkg(:,:,:) = u_bkg(:,:,:) * umask(:,:,:)
552            v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
553         ENDIF
554         !
555         IF ( ln_sshinc ) THEN
556            CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg )
557            ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1)
558         ENDIF
559         !
560         CALL iom_close( inum )
561         !
562      ENDIF
563      !
564      IF(lwp) WRITE(numout,*) '   ==>>>   Euler time step switch is ', neuler
565      !
566      IF( lk_asminc ) THEN                            !==  data assimilation  ==!
567         IF( ln_bkgwri )   CALL asm_bkg_wri( nit000 - 1 )      ! Output background fields
568         IF( ln_asmdin ) THEN                                  ! Direct initialization
569            IF( ln_trainc )   CALL tra_asm_inc( nit000 - 1 )      ! Tracers
570            IF( ln_dyninc )   CALL dyn_asm_inc( nit000 - 1 )      ! Dynamics
571            IF( ln_sshinc )   CALL ssh_asm_inc( nit000 - 1 )      ! SSH
572         ENDIF
573      ENDIF
574      !
575   END SUBROUTINE asm_inc_init
576   
577   
578   SUBROUTINE tra_asm_inc( kt )
579      !!----------------------------------------------------------------------
580      !!                    ***  ROUTINE tra_asm_inc  ***
581      !!         
582      !! ** Purpose : Apply the tracer (T and S) assimilation increments
583      !!
584      !! ** Method  : Direct initialization or Incremental Analysis Updating
585      !!
586      !! ** Action  :
587      !!----------------------------------------------------------------------
588      INTEGER, INTENT(IN) ::   kt   ! Current time step
589      !
590      INTEGER  :: ji, jj, jk
591      INTEGER  :: it
592      REAL(wp) :: zincwgt  ! IAU weight for current time step
593      REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values
594      !!----------------------------------------------------------------------
595      !
596      ! freezing point calculation taken from oc_fz_pt (but calculated for all depths)
597      ! used to prevent the applied increments taking the temperature below the local freezing point
598      DO jk = 1, jpkm1
599        CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), gdept_n(:,:,jk) )
600      END DO
601         !
602         !                             !--------------------------------------
603      IF ( ln_asmiau ) THEN            ! Incremental Analysis Updating
604         !                             !--------------------------------------
605         !
606         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
607            !
608            it = kt - nit000 + 1
609            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
610            !
611            IF(lwp) THEN
612               WRITE(numout,*) 
613               WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
614               WRITE(numout,*) '~~~~~~~~~~~~'
615            ENDIF
616            !
617            ! Update the tracer tendencies
618            DO jk = 1, jpkm1
619               IF (ln_temnofreeze) THEN
620                  ! Do not apply negative increments if the temperature will fall below freezing
621                  WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. &
622                     &   tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 
623                     tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 
624                  END WHERE
625               ELSE
626                  tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 
627               ENDIF
628               IF (ln_salfix) THEN
629                  ! Do not apply negative increments if the salinity will fall below a specified
630                  ! minimum value salfixmin
631                  WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. &
632                     &   tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 
633                     tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
634                  END WHERE
635               ELSE
636                  tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
637               ENDIF
638            END DO
639            !
640         ENDIF
641         !
642         IF ( kt == nitiaufin_r + 1  ) THEN   ! For bias crcn to work
643            DEALLOCATE( t_bkginc )
644            DEALLOCATE( s_bkginc )
645         ENDIF
646         !                             !--------------------------------------
647      ELSEIF ( ln_asmdin ) THEN        ! Direct Initialization
648         !                             !--------------------------------------
649         !           
650         IF ( kt == nitdin_r ) THEN
651            !
652            neuler = 0  ! Force Euler forward step
653            !
654            ! Initialize the now fields with the background + increment
655            IF (ln_temnofreeze) THEN
656               ! Do not apply negative increments if the temperature will fall below freezing
657               WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 
658                  tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
659               END WHERE
660            ELSE
661               tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
662            ENDIF
663            IF (ln_salfix) THEN
664               ! Do not apply negative increments if the salinity will fall below a specified
665               ! minimum value salfixmin
666               WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 
667                  tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
668               END WHERE
669            ELSE
670               tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
671            ENDIF
672
673            tsb(:,:,:,:) = tsn(:,:,:,:)                 ! Update before fields
674
675            CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )  ! Before potential and in situ densities
676!!gm  fabien
677!            CALL eos( tsb, rhd, rhop )                ! Before potential and in situ densities
678!!gm
679
680            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)      &
681               &  CALL zps_hde    ( kt, jpts, tsb, gtsu, gtsv,        &  ! Partial steps: before horizontal gradient
682               &                              rhd, gru , grv          )  ! of t, s, rd at the last ocean level
683            IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)      &
684               &  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi,    &    ! Partial steps for top cell (ISF)
685               &                                  rhd, gru , grv , grui, grvi     )    ! of t, s, rd at the last ocean level
686
687            DEALLOCATE( t_bkginc )
688            DEALLOCATE( s_bkginc )
689            DEALLOCATE( t_bkg    )
690            DEALLOCATE( s_bkg    )
691         ENDIF
692         
693      ENDIF
694      ! Perhaps the following call should be in step
695      IF ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment
696      !
697   END SUBROUTINE tra_asm_inc
698
699
700   SUBROUTINE dyn_asm_inc( kt )
701      !!----------------------------------------------------------------------
702      !!                    ***  ROUTINE dyn_asm_inc  ***
703      !!         
704      !! ** Purpose : Apply the dynamics (u and v) assimilation increments.
705      !!
706      !! ** Method  : Direct initialization or Incremental Analysis Updating.
707      !!
708      !! ** Action  :
709      !!----------------------------------------------------------------------
710      INTEGER, INTENT(IN) :: kt   ! Current time step
711      !
712      INTEGER :: jk
713      INTEGER :: it
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) / rdt   ! IAU weight for the current time step
725            !
726            IF(lwp) THEN
727               WRITE(numout,*) 
728               WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
729               WRITE(numout,*) '~~~~~~~~~~~~'
730            ENDIF
731            !
732            ! Update the dynamic tendencies
733            DO jk = 1, jpkm1
734               ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt
735               va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt
736            END DO
737            !
738            IF ( kt == nitiaufin_r ) THEN
739               DEALLOCATE( u_bkginc )
740               DEALLOCATE( v_bkginc )
741            ENDIF
742            !
743         ENDIF
744         !                          !-----------------------------------------
745      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
746         !                          !-----------------------------------------
747         !         
748         IF ( kt == nitdin_r ) THEN
749            !
750            neuler = 0                    ! Force Euler forward step
751            !
752            ! Initialize the now fields with the background + increment
753            un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:)
754            vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:) 
755            !
756            ub(:,:,:) = un(:,:,:)         ! Update before fields
757            vb(:,:,:) = vn(:,:,:)
758            !
759            DEALLOCATE( u_bkg    )
760            DEALLOCATE( v_bkg    )
761            DEALLOCATE( u_bkginc )
762            DEALLOCATE( v_bkginc )
763         ENDIF
764         !
765      ENDIF
766      !
767   END SUBROUTINE dyn_asm_inc
768
769
770   SUBROUTINE ssh_asm_inc( kt )
771      !!----------------------------------------------------------------------
772      !!                    ***  ROUTINE ssh_asm_inc  ***
773      !!         
774      !! ** Purpose : Apply the sea surface height assimilation increment.
775      !!
776      !! ** Method  : Direct initialization or Incremental Analysis Updating.
777      !!
778      !! ** Action  :
779      !!----------------------------------------------------------------------
780      INTEGER, INTENT(IN) :: kt   ! Current time step
781      !
782      INTEGER :: it
783      INTEGER :: jk
784      REAL(wp) :: zincwgt  ! IAU weight for current time step
785      !!----------------------------------------------------------------------
786      !
787      !                             !-----------------------------------------
788      IF ( ln_asmiau ) THEN         ! Incremental Analysis Updating
789         !                          !-----------------------------------------
790         !
791         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
792            !
793            it = kt - nit000 + 1
794            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
795            !
796            IF(lwp) THEN
797               WRITE(numout,*) 
798               WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
799                  &  kt,' with IAU weight = ', wgtiau(it)
800               WRITE(numout,*) '~~~~~~~~~~~~'
801            ENDIF
802            !
803            ! Save the tendency associated with the IAU weighted SSH increment
804            ! (applied in dynspg.*)
805#if defined key_asminc
806            ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
807#endif
808            !
809         ELSE IF( kt == nitiaufin_r+1 ) THEN
810            !
811            ! test on ssh_bkginc needed as ssh_asm_inc is called twice by time step
812            IF (ALLOCATED(ssh_bkginc)) DEALLOCATE( ssh_bkginc )
813            !
814#if defined key_asminc
815            ssh_iau(:,:) = 0._wp
816#endif
817            !
818         ENDIF
819         !                          !-----------------------------------------
820      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
821         !                          !-----------------------------------------
822         !
823         IF ( kt == nitdin_r ) THEN
824            !
825            neuler = 0                                   ! Force Euler forward step
826            !
827            sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)   ! Initialize the now fields the background + increment
828            !
829            sshb(:,:) = sshn(:,:)                        ! Update before fields
830            e3t_b(:,:,:) = e3t_n(:,:,:)
831!!gm why not e3u_b, e3v_b, gdept_b ????
832            !
833            DEALLOCATE( ssh_bkg    )
834            DEALLOCATE( ssh_bkginc )
835            !
836         ENDIF
837         !
838      ENDIF
839      !
840   END SUBROUTINE ssh_asm_inc
841
842
843   SUBROUTINE ssh_asm_div( kt, phdivn )
844      !!----------------------------------------------------------------------
845      !!                  ***  ROUTINE ssh_asm_div  ***
846      !!
847      !! ** Purpose :   ssh increment with z* is incorporated via a correction of the local divergence         
848      !!                across all the water column
849      !!
850      !! ** Method  :
851      !!                CAUTION : sshiau is positive (inflow) decreasing the
852      !!                          divergence and expressed in m/s
853      !!
854      !! ** Action  :   phdivn   decreased by the ssh increment
855      !!----------------------------------------------------------------------
856      INTEGER, INTENT(IN) :: kt                               ! ocean time-step index
857      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
858      !!
859      INTEGER  ::   jk                                        ! dummy loop index
860      REAL(wp), DIMENSION(:,:)  , POINTER       ::   ztim     ! local array
861      !!----------------------------------------------------------------------
862      !
863#if defined key_asminc
864      CALL ssh_asm_inc( kt ) !==   (calculate increments)
865      !
866      IF( ln_linssh ) THEN
867         phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t_n(:,:,1) * tmask(:,:,1)
868      ELSE
869         ALLOCATE( ztim(jpi,jpj) )
870         ztim(:,:) = ssh_iau(:,:) / ( ht_n(:,:) + 1.0 - ssmask(:,:) )
871         DO jk = 1, jpkm1                                 
872            phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) 
873         END DO
874         !
875         DEALLOCATE(ztim)
876      ENDIF
877#endif
878      !
879   END SUBROUTINE ssh_asm_div
880
881
882   SUBROUTINE seaice_asm_inc( kt, kindic )
883      !!----------------------------------------------------------------------
884      !!                    ***  ROUTINE seaice_asm_inc  ***
885      !!         
886      !! ** Purpose : Apply the sea ice assimilation increment.
887      !!
888      !! ** Method  : Direct initialization or Incremental Analysis Updating.
889      !!
890      !! ** Action  :
891      !!
892      !!----------------------------------------------------------------------
893      INTEGER, INTENT(in)           ::   kt       ! Current time step
894      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation
895      !
896      INTEGER  ::   it, jk
897      REAL(wp) ::   zincwgt                            ! IAU weight for current time step
898#if defined key_si3 && defined key_asminc
899      REAL(wp), DIMENSION(jpi,jpj,jpl) :: da_i         ! change in ice concentration
900      REAL(wp), DIMENSION(jpi,jpj,jpl) :: dv_i         ! change in ice volume
901      REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i       ! inverse of ice concentration before current IAU step
902      REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_v_i       ! inverse of ice volume before current IAU step
903      REAL(wp), DIMENSION(jpi,jpj)     :: zhi_damin_2D ! array with DA thickness for incr_newice
904#endif
905      !!----------------------------------------------------------------------
906      !
907      !                             !-----------------------------------------
908      IF ( ln_asmiau ) THEN         ! Incremental Analysis Updating
909         !                          !-----------------------------------------
910         !
911         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
912            !
913            it = kt - nit000 + 1
914            zincwgt = wgtiau(it)      ! IAU weight for the current time step
915            ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments)
916            !
917            IF(lwp) THEN
918               WRITE(numout,*) 
919               WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
920               WRITE(numout,*) '~~~~~~~~~~~~'
921            ENDIF
922            !
923            ! Sea-ice : SI3 case
924            !
925#if defined key_si3 && defined key_asminc
926            ! compute the inverse of key sea ice variables
927            ! to be used later in the code
928            WHERE( a_i(:,:,:) > epsi10 )
929               z1_a_i(:,:,:) = 1.0_wp / a_i(:,:,:)
930               z1_v_i(:,:,:) = 1.0_wp / v_i(:,:,:)
931            ELSEWHERE
932               z1_a_i(:,:,:) = 0.0_wp
933               z1_v_i(:,:,:) = 0.0_wp
934            END WHERE
935
936            ! add positive concentration increments to regions where ice
937            ! is already present and bound them to 1
938            ! ice volume is added based on zhi_damin
939            WHERE ( .NOT. incr_newice .AND. a_i_bkginc(:,:,:) > 0.0_wp )
940               a_i(:,:,:) = a_i(:,:,:) + MIN( 1.0_wp - a_i(:,:,:), a_i_bkginc(:,:,:) * zincwgt )
941               v_i(:,:,:) = v_i(:,:,:) + MIN( 1.0_wp - a_i(:,:,:), a_i_bkginc(:,:,:) * zincwgt ) * zhi_damin
942            END WHERE
943
944            ! add negative concentration increments to regions where ice
945            ! is already present and bound them to 0
946            ! in this case ice volume is changed based on the current thickness
947            WHERE ( .NOT. incr_newice .AND. a_i_bkginc(:,:,:) < 0.0_wp )
948               a_i(:,:,:) = MAX( a_i(:,:,:) + a_i_bkginc(:,:,:) * zincwgt, 0.0_wp )
949               v_i(:,:,:) = a_i(:,:,:) * h_i(:,:,:)
950            END WHERE
951           
952            ! compute changes in ice concentration and volume
953            WHERE ( incr_newice )
954               da_i(:,:,:) = 1.0_wp
955               dv_i(:,:,:) = 1.0_wp
956            ELSEWHERE
957               da_i(:,:,:) = a_i(:,:,:) * z1_a_i(:,:,:)
958               dv_i(:,:,:) = v_i(:,:,:) * z1_v_i(:,:,:)
959            END WHERE
960
961            ! initialise thermodynamics of new ice being added to open water
962            ! just do it once since next IAU steps assume that new ice has
963            ! already been added in
964            IF ( kt == nitiaustr_r ) THEN
965
966               ! assign zhi_damin to ice forming in open water
967               WHERE ( ANY( incr_newice, DIM=3 ) )
968                  zhi_damin_2D(:,:) = zhi_damin
969               ELSEWHERE
970                  zhi_damin_2D(:,:) = 0.0_wp
971               END WHERE
972
973               ! add ice concentration and volume
974               ! ensure the other prognostic variables are set to zero
975               WHERE ( incr_newice )
976                  a_i(:,:,:) = MIN( 1.0_wp, a_i_bkginc(:,:,:) * zincwgt )
977                  v_i(:,:,:) = MIN( 1.0_wp, a_i_bkginc(:,:,:) * zincwgt ) * zhi_damin
978                  v_s (:,:,:) = 0.0_wp
979                  a_ip(:,:,:) = 0.0_wp
980                  v_ip(:,:,:) = 0.0_wp
981                  sv_i(:,:,:) = 0.0_wp
982               END WHERE
983               DO jk = 1, nlay_i
984                  WHERE ( incr_newice )
985                     e_i(:,:,jk,:) = 0.0_wp
986                  END WHERE
987               END DO
988               DO jk = 1, nlay_s
989                  WHERE ( incr_newice )
990                     e_s(:,:,jk,:) = 0.0_wp
991                  END WHERE
992               END DO
993
994               ! Initialisation of the salt content and ice enthalpy
995               ! set flag of new ice to false after this
996               CALL init_new_ice_thd( zhi_damin_2D )
997               incr_newice(:,:,:) = .FALSE.
998            END IF
999
1000            ! maintain equivalent values for key prognostic variables
1001            v_s(:,:,:) = v_s(:,:,:) * da_i(:,:,:)
1002            DO jk = 1, nlay_s
1003               e_s(:,:,jk,:) = e_s(:,:,jk,:) * da_i(:,:,:)
1004            END DO
1005            a_ip (:,:,:) = a_ip(:,:,:) * da_i(:,:,:)
1006            v_ip (:,:,:) = v_ip(:,:,:) * da_i(:,:,:)
1007 
1008            ! ice volume dependent variables
1009            sv_i (:,:,:) = sv_i(:,:,:) * dv_i(:,:,:)
1010            DO jk = 1, nlay_i
1011               e_i(:,:,jk,:) = e_i(:,:,jk,:) * dv_i(:,:,:)
1012            END DO
1013#endif
1014
1015#if defined key_cice && defined key_asminc
1016            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
1017            ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt
1018#endif
1019            !
1020            IF ( kt == nitiaufin_r ) THEN
1021               DEALLOCATE( seaice_bkginc )
1022#if defined key_si3 && defined key_asminc
1023               DEALLOCATE( incr_newice )
1024               DEALLOCATE( a_i_bkginc )
1025#endif
1026            ENDIF
1027            !
1028         ELSE
1029            !
1030#if defined key_cice && defined key_asminc
1031            ndaice_da(:,:) = 0._wp        ! Sea-ice : CICE case. Zero ice increment tendency into CICE
1032#endif
1033            !
1034         ENDIF
1035         !                          !-----------------------------------------
1036      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
1037         !                          !-----------------------------------------
1038         !
1039         IF ( kt == nitdin_r ) THEN
1040            !
1041            neuler = 0                    ! Force Euler forward step
1042
1043            IF(lwp) THEN
1044               WRITE(numout,*)
1045               WRITE(numout,*) 'seaice_asm_inc : sea ice direct initialization at time step = ', kt
1046               WRITE(numout,*) '~~~~~~~~~~~~'
1047            ENDIF
1048            !
1049#if defined key_cice && defined key_asminc
1050            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
1051           ndaice_da(:,:) = seaice_bkginc(:,:) / rdt
1052#endif
1053            IF ( .NOT. PRESENT(kindic) ) THEN
1054               DEALLOCATE( seaice_bkginc )
1055            END IF
1056            !
1057         ELSE
1058            !
1059#if defined key_cice && defined key_asminc
1060            ndaice_da(:,:) = 0._wp     ! Sea-ice : CICE case. Zero ice increment tendency into CICE
1061#endif
1062            !
1063         ENDIF
1064         !
1065      ENDIF
1066      !
1067   END SUBROUTINE seaice_asm_inc
1068
1069
1070   SUBROUTINE init_new_ice_thd( hi_new )
1071      !!----------------------------------------------------------------------
1072      !!                  ***  ROUTINE init_new_ice_thd  ***
1073      !!
1074      !! ** Purpose :   Initialise thermodynamics of new ice
1075      !!                forming at 1st category with thickness hi_new
1076      !!
1077      !! ** Method  :   Apply SI3  equations to initialise
1078      !!                thermodynamics of new ice
1079      !!
1080      !! ** Action  :   update sea ice thermodynamics
1081      !!----------------------------------------------------------------------
1082      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    ::   hi_new  ! total thickness of new ice
1083
1084      INTEGER :: jj, ji, jk
1085      REAL(wp) ::   ztmelts         ! melting point
1086      REAL(wp) ::   Sice_Fz=2.3_wp  ! Salinity of the ice = F(z) [multiyear ice]
1087
1088      REAL(wp), DIMENSION(jpij) ::   zh_newice   ! 1d version of hi_new
1089      REAL(wp), DIMENSION(jpij) ::   zs_newice   ! salinity of new ice
1090      !!----------------------------------------------------------------------
1091
1092      ! Identify grid points where new ice forms
1093      npti = 0   ;   nptidx(:) = 0
1094      DO jj = 1, jpj
1095         DO ji = 1, jpi
1096            IF ( hi_new(ji,jj) > 0._wp ) THEN
1097               npti = npti + 1
1098               nptidx( npti ) = (jj - 1) * jpi + ji
1099            ENDIF
1100         END DO
1101      END DO
1102
1103      ! Move from 2-D to 1-D vectors
1104      IF ( npti > 0 ) THEN
1105         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
1106         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
1107         CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , hi_new        )
1108         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d    (1:npti) , sss_m         )
1109         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d   (1:npti) , t_bo          )
1110         DO jk = 1, nlay_i
1111            CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,1) )
1112         END DO
1113
1114         ! --- Salinity of new ice --- !
1115         SELECT CASE ( nn_icesal )
1116         CASE ( 1 )                    ! Sice = constant
1117            zs_newice(1:npti) = rn_icesal
1118         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
1119            DO ji = 1, npti
1120               zs_newice(ji) = MIN(  4.606_wp + 0.91_wp / zh_newice(ji) , rn_simax , 0.5_wp * sss_1d(ji) )
1121            END DO
1122         CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
1123            zs_newice(1:npti) = Sice_Fz
1124         END SELECT
1125
1126         ! --- Update ice salt content --- !
1127         DO ji = 1, npti
1128            sv_i_2d(ji,1) = sv_i_2d(ji,1) + zs_newice(ji) * ( v_i_2d(ji,1) )
1129         END DO
1130
1131         ! --- Heat content of new ice --- !
1132         ! We assume that new ice is formed at the seawater freezing point
1133         DO ji = 1, npti
1134               ztmelts       = - rTmlt * zs_newice(ji)                  ! Melting point (C)
1135               e_i_1d(ji,:)  =   rhoi * (  rcpi  * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                        &
1136                  &                      + rLfus * ( 1.0_wp - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
1137                  &                      - rcp   *         ztmelts )
1138         END DO
1139
1140         ! Change units for e_i
1141         DO jk = 1, nlay_i
1142            e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) * v_i_2d(1:npti,1) * r1_nlay_i 
1143         END DO
1144
1145         ! Reforming full thermodynamic variables
1146         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
1147         DO jk = 1, nlay_i
1148               CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,1) )
1149         END DO
1150      END IF
1151
1152   END SUBROUTINE init_new_ice_thd
1153
1154
1155   !!======================================================================
1156END MODULE asminc
Note: See TracBrowser for help on using the repository browser.