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

source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ASM/asminc.F90 @ 10001

Last change on this file since 10001 was 10001, checked in by gm, 6 years ago

#1911 (ENHANCE-04): RK3 branch - step I.1 and I.2 (see wiki page)

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