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 @ 15242

Last change on this file since 15242 was 15242, checked in by dcarneir, 3 years ago

Leaving only the essential bits in the code

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