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/trunk/src/SWE – NEMO

source: NEMO/trunk/src/SWE/asminc.F90 @ 13295

Last change on this file since 13295 was 13295, checked in by acc, 4 years ago

Replace do-loop macros in the trunk with alternative forms with greater flexibility for extra halo applications. This alters a lot of routines but does not change any behaviour or results. do_loop_substitute.h90 is greatly simplified by this change. SETTE results are identical to those with the previous revision

File size: 48.8 KB
Line 
1MODULE asminc
2   !!======================================================================
3   !!                       ***  MODULE asminc  ***
4   !! Assimilation increment : Apply an increment generated by data
5   !!                          assimilation
6   !!======================================================================
7   !! History :       ! 2007-03  (M. Martin)  Met Office version
8   !!                 ! 2007-04  (A. Weaver)  calc_date original code
9   !!                 ! 2007-04  (A. Weaver)  Merge with OPAVAR/NEMOVAR
10   !!   NEMO     3.3  ! 2010-05  (D. Lea)  Update to work with NEMO v3.2
11   !!             -   ! 2010-05  (D. Lea)  add calc_month_len routine based on day_init
12   !!            3.4  ! 2012-10  (A. Weaver and K. Mogensen) Fix for direct initialization
13   !!                 ! 2014-09  (D. Lea)  Local calc_date removed use routine from OBS
14   !!                 ! 2015-11  (D. Lea)  Handle non-zero initial time of day
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!   asm_inc_init   : Initialize the increment arrays and IAU weights
19   !!   tra_asm_inc    : Apply the tracer (T and S) increments
20   !!   dyn_asm_inc    : Apply the dynamic (u and v) increments
21   !!   ssh_asm_inc    : Apply the SSH increment
22   !!   ssh_asm_div    : Apply divergence associated with SSH increment
23   !!   seaice_asm_inc : Apply the seaice increment
24   !!----------------------------------------------------------------------
25   USE oce             ! Dynamics and active tracers defined in memory
26   USE par_oce         ! Ocean space and time domain variables
27   USE dom_oce         ! Ocean space and time domain
28   USE domvvl          ! domain: variable volume level
29   USE ldfdyn          ! lateral diffusion: eddy viscosity coefficients
30   USE eosbn2          ! Equation of state - in situ and potential density
31   USE zpshde          ! Partial step : Horizontal Derivative
32   USE asmpar          ! Parameters for the assmilation interface
33   USE asmbkg          !
34   USE c1d             ! 1D initialization
35   USE sbc_oce         ! Surface boundary condition variables.
36   USE diaobs   , ONLY : calc_date     ! Compute the calendar date on a given step
37#if defined key_si3
38   USE ice      , ONLY : hm_i, at_i, at_i_b
39#endif
40   !
41   USE in_out_manager  ! I/O manager
42   USE iom             ! Library to read input files
43   USE lib_mpp         ! MPP library
44
45   IMPLICIT NONE
46   PRIVATE
47   
48   PUBLIC   asm_inc_init   !: Initialize the increment arrays and IAU weights
49   PUBLIC   tra_asm_inc    !: Apply the tracer (T and S) increments
50   PUBLIC   dyn_asm_inc    !: Apply the dynamic (u and v) increments
51   PUBLIC   ssh_asm_inc    !: Apply the SSH increment
52   PUBLIC   ssh_asm_div    !: Apply the SSH divergence
53   PUBLIC   seaice_asm_inc !: Apply the seaice increment
54
55#if defined key_asminc
56    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .TRUE.   !: Logical switch for assimilation increment interface
57#else
58    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE.  !: No assimilation increments
59#endif
60   LOGICAL, PUBLIC :: ln_bkgwri     !: No output of the background state fields
61   LOGICAL, PUBLIC :: ln_asmiau     !: No applying forcing with an assimilation increment
62   LOGICAL, PUBLIC :: ln_asmdin     !: No direct initialization
63   LOGICAL, PUBLIC :: ln_trainc     !: No tracer (T and S) assimilation increments
64   LOGICAL, PUBLIC :: ln_dyninc     !: No dynamics (u and v) assimilation increments
65   LOGICAL, PUBLIC :: ln_sshinc     !: No sea surface height assimilation increment
66   LOGICAL, PUBLIC :: ln_seaiceinc  !: No sea ice concentration increment
67   LOGICAL, PUBLIC :: ln_salfix     !: Apply minimum salinity check
68   LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
69   INTEGER, PUBLIC :: nn_divdmp     !: Apply divergence damping filter nn_divdmp times
70
71   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkg   , s_bkg      !: Background temperature and salinity
72   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkg   , v_bkg      !: Background u- & v- velocity components
73   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkginc, s_bkginc   !: Increment to the background T & S
74   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components
75   REAL(wp), PUBLIC, DIMENSION(:)    , ALLOCATABLE ::   wgtiau               !: IAU weights for each time step
76#if defined key_asminc
77   REAL(wp), PUBLIC, DIMENSION(:,:)  , ALLOCATABLE ::   ssh_iau              !: IAU-weighted sea surface height increment
78#endif
79   !                                !!! time steps relative to the cycle interval [0,nitend-nit000-1]
80   INTEGER , PUBLIC ::   nitbkg      !: Time step of the background state used in the Jb term
81   INTEGER , PUBLIC ::   nitdin      !: Time step of the background state for direct initialization
82   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval
83   INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval
84   !
85   INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting
86   !                                 !: = 1   Linear hat-like, centred in middle of IAU interval
87   REAL(wp), PUBLIC ::   salfixmin   !: Ensure that the salinity is larger than this  value if (ln_salfix)
88
89   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment
90   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc
91#if defined key_cice && defined key_asminc
92   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ndaice_da             ! ice increment tendency into CICE
93#endif
94
95   !! * Substitutions
96#  include "do_loop_substitute.h90"
97!!st10
98#  include "domzgr_substitute.h90"
99!!st10
100   !!----------------------------------------------------------------------
101   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
102   !! $Id: asminc.F90 12614 2020-03-26 14:59:52Z gm $
103   !! Software governed by the CeCILL license (see ./LICENSE)
104   !!----------------------------------------------------------------------
105CONTAINS
106
107   SUBROUTINE asm_inc_init( Kbb, Kmm, Krhs )
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, INTENT(in) ::  Kbb, Kmm, Krhs  ! time level indices
118      !
119      INTEGER :: ji, jj, jk, jt  ! dummy loop indices
120      INTEGER :: imid, inum      ! local integers
121      INTEGER :: ios             ! Local integer output status for namelist read
122      INTEGER :: iiauper         ! Number of time steps in the IAU period
123      INTEGER :: icycper         ! Number of time steps in the cycle
124      REAL(KIND=dp) :: ditend_date     ! Date YYYYMMDD.HHMMSS of final time step
125      REAL(KIND=dp) :: ditbkg_date     ! Date YYYYMMDD.HHMMSS of background time step for Jb term
126      REAL(KIND=dp) :: ditdin_date     ! Date YYYYMMDD.HHMMSS of background time step for DI
127      REAL(KIND=dp) :: ditiaustr_date  ! Date YYYYMMDD.HHMMSS of IAU interval start time step
128      REAL(KIND=dp) :: ditiaufin_date  ! Date YYYYMMDD.HHMMSS of IAU interval final time step
129
130      REAL(wp) :: znorm        ! Normalization factor for IAU weights
131      REAL(wp) :: ztotwgt      ! Value of time-integrated IAU weights (should be equal to one)
132      REAL(wp) :: z_inc_dateb  ! Start date of interval on which increment is valid
133      REAL(wp) :: z_inc_datef  ! End date of interval on which increment is valid
134      REAL(wp) :: zdate_bkg    ! Date in background state file for DI
135      REAL(wp) :: zdate_inc    ! Time axis in increments file
136      !
137      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zhdiv   ! 2D workspace
138      !!
139      NAMELIST/nam_asminc/ ln_bkgwri,                                      &
140         &                 ln_trainc, ln_dyninc, ln_sshinc,                &
141         &                 ln_asmdin, ln_asmiau,                           &
142         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
143         &                 ln_salfix, salfixmin, nn_divdmp
144      !!----------------------------------------------------------------------
145
146      !-----------------------------------------------------------------------
147      ! Read Namelist nam_asminc : assimilation increment interface
148      !-----------------------------------------------------------------------
149      ln_seaiceinc   = .FALSE.
150      ln_temnofreeze = .FALSE.
151
152      READ  ( numnam_ref, nam_asminc, IOSTAT = ios, ERR = 901)
153901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_asminc in reference namelist' )
154      READ  ( numnam_cfg, nam_asminc, IOSTAT = ios, ERR = 902 )
155902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_asminc in configuration namelist' )
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_2D( 0, 0, 0, 0 )
419                  zhdiv(ji,jj) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * u_bkginc(ji  ,jj,jk)    &
420                     &            - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * u_bkginc(ji-1,jj,jk)    &
421                     &            + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * v_bkginc(ji,jj  ,jk)    &
422                     &            - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk)  ) &
423                     &            / e3t(ji,jj,jk,Kmm)
424               END_2D
425               CALL lbc_lnk( 'asminc', zhdiv, 'T', 1. )   ! lateral boundary cond. (no sign change)
426               !
427               DO_2D( 0, 0, 0, 0 )
428                  u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk)                         &
429                     &               + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji  ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
430                  v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk)                         &
431                     &               + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj  ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 
432               END_2D
433            END DO
434            !
435         END DO
436         !
437         DEALLOCATE( zhdiv ) 
438         !
439      ENDIF
440      !
441      !                             !-----------------------------------------------------
442      IF ( ln_asmdin ) THEN         ! Allocate and initialize the background state arrays
443         !                          !-----------------------------------------------------
444         !
445         ALLOCATE( t_bkg  (jpi,jpj,jpk) )   ;   t_bkg  (:,:,:) = 0._wp
446         ALLOCATE( s_bkg  (jpi,jpj,jpk) )   ;   s_bkg  (:,:,:) = 0._wp
447         ALLOCATE( u_bkg  (jpi,jpj,jpk) )   ;   u_bkg  (:,:,:) = 0._wp
448         ALLOCATE( v_bkg  (jpi,jpj,jpk) )   ;   v_bkg  (:,:,:) = 0._wp
449         ALLOCATE( ssh_bkg(jpi,jpj)     )   ;   ssh_bkg(:,:)   = 0._wp
450         !
451         !
452         !--------------------------------------------------------------------
453         ! Read from file the background state at analysis time
454         !--------------------------------------------------------------------
455         !
456         CALL iom_open( c_asmdin, inum )
457         !
458         CALL iom_get( inum, 'rdastp', zdate_bkg ) 
459         !
460         IF(lwp) THEN
461            WRITE(numout,*) 
462            WRITE(numout,*) '   ==>>>  Assimilation background state valid at : ', zdate_bkg
463            WRITE(numout,*)
464         ENDIF
465         !
466         IF ( zdate_bkg /= ditdin_date )   &
467            & CALL ctl_warn( ' Validity time of assimilation background state does', &
468            &                ' not agree with Direct Initialization time' )
469         !
470         IF ( ln_trainc ) THEN   
471            CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg )
472            CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg )
473            t_bkg(:,:,:) = t_bkg(:,:,:) * tmask(:,:,:)
474            s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:)
475         ENDIF
476         !
477         IF ( ln_dyninc ) THEN   
478            CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg )
479            CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg )
480            u_bkg(:,:,:) = u_bkg(:,:,:) * umask(:,:,:)
481            v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
482         ENDIF
483         !
484         IF ( ln_sshinc ) THEN
485            CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg )
486            ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1)
487         ENDIF
488         !
489         CALL iom_close( inum )
490         !
491      ENDIF
492      !
493      IF(lwp) WRITE(numout,*) '   ==>>>   Euler time step switch is ', l_1st_euler
494      !
495      IF( lk_asminc ) THEN                            !==  data assimilation  ==!
496         IF( ln_bkgwri )   CALL asm_bkg_wri( nit000 - 1, Kmm )      ! Output background fields
497         IF( ln_asmdin ) THEN                                  ! Direct initialization
498            IF( ln_trainc )   CALL tra_asm_inc( nit000 - 1, Kbb, Kmm, ts    , Krhs )      ! Tracers
499            IF( ln_dyninc )   CALL dyn_asm_inc( nit000 - 1, Kbb, Kmm, uu, vv, Krhs )      ! Dynamics
500            IF( ln_sshinc )   CALL ssh_asm_inc( nit000 - 1, Kbb, Kmm )                    ! SSH
501         ENDIF
502      ENDIF
503      !
504   END SUBROUTINE asm_inc_init
505   
506   
507   SUBROUTINE tra_asm_inc( kt, Kbb, Kmm, pts, Krhs )
508      !!----------------------------------------------------------------------
509      !!                    ***  ROUTINE tra_asm_inc  ***
510      !!         
511      !! ** Purpose : Apply the tracer (T and S) assimilation increments
512      !!
513      !! ** Method  : Direct initialization or Incremental Analysis Updating
514      !!
515      !! ** Action  :
516      !!----------------------------------------------------------------------
517      INTEGER                                  , INTENT(in   ) :: kt             ! Current time step
518      INTEGER                                  , INTENT(in   ) :: Kbb, Kmm, Krhs ! Time level indices
519      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers and RHS of tracer equation
520      !
521      INTEGER  :: ji, jj, jk
522      INTEGER  :: it
523      REAL(wp) :: zincwgt  ! IAU weight for current time step
524      REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values
525      !!----------------------------------------------------------------------
526      !
527      ! freezing point calculation taken from oc_fz_pt (but calculated for all depths)
528      ! used to prevent the applied increments taking the temperature below the local freezing point
529      DO jk = 1, jpkm1
530        CALL eos_fzp( pts(:,:,jk,jp_sal,Kmm), fzptnz(:,:,jk), gdept(:,:,jk,Kmm) )
531      END DO
532         !
533         !                             !--------------------------------------
534      IF ( ln_asmiau ) THEN            ! Incremental Analysis Updating
535         !                             !--------------------------------------
536         !
537         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
538            !
539            it = kt - nit000 + 1
540            zincwgt = wgtiau(it) / rn_Dt   ! IAU weight for the current time step
541            !
542            IF(lwp) THEN
543               WRITE(numout,*) 
544               WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
545               WRITE(numout,*) '~~~~~~~~~~~~'
546            ENDIF
547            !
548            ! Update the tracer tendencies
549            DO jk = 1, jpkm1
550               IF (ln_temnofreeze) THEN
551                  ! Do not apply negative increments if the temperature will fall below freezing
552                  WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. &
553                     &   pts(:,:,jk,jp_tem,Kmm) + pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 
554                     pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 
555                  END WHERE
556               ELSE
557                  pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 
558               ENDIF
559               IF (ln_salfix) THEN
560                  ! Do not apply negative increments if the salinity will fall below a specified
561                  ! minimum value salfixmin
562                  WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. &
563                     &   pts(:,:,jk,jp_sal,Kmm) + pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 
564                     pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt
565                  END WHERE
566               ELSE
567                  pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt
568               ENDIF
569            END DO
570            !
571         ENDIF
572         !
573         IF ( kt == nitiaufin_r + 1  ) THEN   ! For bias crcn to work
574            DEALLOCATE( t_bkginc )
575            DEALLOCATE( s_bkginc )
576         ENDIF
577         !                             !--------------------------------------
578      ELSEIF ( ln_asmdin ) THEN        ! Direct Initialization
579         !                             !--------------------------------------
580         !           
581         IF ( kt == nitdin_r ) THEN
582            !
583            l_1st_euler = .TRUE.  ! Force Euler forward step
584            !
585            ! Initialize the now fields with the background + increment
586            IF (ln_temnofreeze) THEN
587               ! Do not apply negative increments if the temperature will fall below freezing
588               WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_tem,Kmm) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 
589                  pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
590               END WHERE
591            ELSE
592               pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
593            ENDIF
594            IF (ln_salfix) THEN
595               ! Do not apply negative increments if the salinity will fall below a specified
596               ! minimum value salfixmin
597               WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_sal,Kmm) + s_bkginc(:,:,:) > salfixmin ) 
598                  pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
599               END WHERE
600            ELSE
601               pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
602            ENDIF
603
604            pts(:,:,:,:,Kbb) = pts(:,:,:,:,Kmm)                 ! Update before fields
605
606            CALL eos( pts(:,:,:,:,Kbb), rhd, rhop, gdept_0(:,:,:) )  ! Before potential and in situ densities
607!!gm  fabien
608!            CALL eos( pts(:,:,:,:,Kbb), rhd, rhop )                ! Before potential and in situ densities
609!!gm
610
611            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)           &
612               &  CALL zps_hde    ( kt, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv,        &  ! Partial steps: before horizontal gradient
613               &                              rhd, gru , grv               )  ! of t, s, rd at the last ocean level
614            IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)                       &
615               &  CALL zps_hde_isf( nit000, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, gtui, gtvi,    &  ! Partial steps for top cell (ISF)
616               &                                  rhd, gru , grv , grui, grvi          )  ! of t, s, rd at the last ocean level
617
618            DEALLOCATE( t_bkginc )
619            DEALLOCATE( s_bkginc )
620            DEALLOCATE( t_bkg    )
621            DEALLOCATE( s_bkg    )
622         ENDIF
623         
624      ENDIF
625      ! Perhaps the following call should be in step
626      IF ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment
627      !
628   END SUBROUTINE tra_asm_inc
629
630
631   SUBROUTINE dyn_asm_inc( kt, Kbb, Kmm, puu, pvv, Krhs )
632      !!----------------------------------------------------------------------
633      !!                    ***  ROUTINE dyn_asm_inc  ***
634      !!         
635      !! ** Purpose : Apply the dynamics (u and v) assimilation increments.
636      !!
637      !! ** Method  : Direct initialization or Incremental Analysis Updating.
638      !!
639      !! ** Action  :
640      !!----------------------------------------------------------------------
641      INTEGER                             , INTENT( in )  ::  kt             ! ocean time-step index
642      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs ! ocean time level indices
643      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv       ! ocean velocities and RHS of momentum equation
644      !
645      INTEGER :: jk
646      INTEGER :: it
647      REAL(wp) :: zincwgt  ! IAU weight for current time step
648      !!----------------------------------------------------------------------
649      !
650      !                          !--------------------------------------------
651      IF ( ln_asmiau ) THEN      ! Incremental Analysis Updating
652         !                       !--------------------------------------------
653         !
654         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
655            !
656            it = kt - nit000 + 1
657            zincwgt = wgtiau(it) / rn_Dt   ! IAU weight for the current time step
658            !
659            IF(lwp) THEN
660               WRITE(numout,*) 
661               WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
662               WRITE(numout,*) '~~~~~~~~~~~~'
663            ENDIF
664            !
665            ! Update the dynamic tendencies
666            DO jk = 1, jpkm1
667               puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + u_bkginc(:,:,jk) * zincwgt
668               pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + v_bkginc(:,:,jk) * zincwgt
669            END DO
670            !
671            IF ( kt == nitiaufin_r ) THEN
672               DEALLOCATE( u_bkginc )
673               DEALLOCATE( v_bkginc )
674            ENDIF
675            !
676         ENDIF
677         !                          !-----------------------------------------
678      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
679         !                          !-----------------------------------------
680         !         
681         IF ( kt == nitdin_r ) THEN
682            !
683            l_1st_euler = .TRUE.                    ! Force Euler forward step
684            !
685            ! Initialize the now fields with the background + increment
686            puu(:,:,:,Kmm) = u_bkg(:,:,:) + u_bkginc(:,:,:)
687            pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:) 
688            !
689            puu(:,:,:,Kbb) = puu(:,:,:,Kmm)         ! Update before fields
690            pvv(:,:,:,Kbb) = pvv(:,:,:,Kmm)
691            !
692            DEALLOCATE( u_bkg    )
693            DEALLOCATE( v_bkg    )
694            DEALLOCATE( u_bkginc )
695            DEALLOCATE( v_bkginc )
696         ENDIF
697         !
698      ENDIF
699      !
700   END SUBROUTINE dyn_asm_inc
701
702
703   SUBROUTINE ssh_asm_inc( kt, Kbb, Kmm )
704      !!----------------------------------------------------------------------
705      !!                    ***  ROUTINE ssh_asm_inc  ***
706      !!         
707      !! ** Purpose : Apply the sea surface height assimilation increment.
708      !!
709      !! ** Method  : Direct initialization or Incremental Analysis Updating.
710      !!
711      !! ** Action  :
712      !!----------------------------------------------------------------------
713      INTEGER, INTENT(IN) :: kt         ! Current time step
714      INTEGER, INTENT(IN) :: Kbb, Kmm   ! Current time step
715      !
716      INTEGER :: it
717      INTEGER :: jk
718      REAL(wp) :: zincwgt  ! IAU weight for current time step
719      !!----------------------------------------------------------------------
720      !
721      !                             !-----------------------------------------
722      IF ( ln_asmiau ) THEN         ! Incremental Analysis Updating
723         !                          !-----------------------------------------
724         !
725         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
726            !
727            it = kt - nit000 + 1
728            zincwgt = wgtiau(it) / rn_Dt   ! IAU weight for the current time step
729            !
730            IF(lwp) THEN
731               WRITE(numout,*) 
732               WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
733                  &  kt,' with IAU weight = ', wgtiau(it)
734               WRITE(numout,*) '~~~~~~~~~~~~'
735            ENDIF
736            !
737            ! Save the tendency associated with the IAU weighted SSH increment
738            ! (applied in dynspg.*)
739#if defined key_asminc
740            ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
741#endif
742            !
743         ELSE IF( kt == nitiaufin_r+1 ) THEN
744            !
745            ! test on ssh_bkginc needed as ssh_asm_inc is called twice by time step
746            IF (ALLOCATED(ssh_bkginc)) DEALLOCATE( ssh_bkginc )
747            !
748#if defined key_asminc
749            ssh_iau(:,:) = 0._wp
750#endif
751            !
752         ENDIF
753         !                          !-----------------------------------------
754      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
755         !                          !-----------------------------------------
756         !
757         IF ( kt == nitdin_r ) THEN
758            !
759            l_1st_euler = .TRUE.                            ! Force Euler forward step
760            !
761            ssh(:,:,Kmm) = ssh_bkg(:,:) + ssh_bkginc(:,:)   ! Initialize the now fields the background + increment
762            !
763            ssh(:,:,Kbb) = ssh(:,:,Kmm)                        ! Update before fields
764!!st11
765#if ! defined key_qco
766            e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
767#endif
768!!st11
769!!gm why not e3u(:,:,:,Kbb), e3v(:,:,:,Kbb), gdept(:,:,jk,Kbb) ????
770            !
771            DEALLOCATE( ssh_bkg    )
772            DEALLOCATE( ssh_bkginc )
773            !
774         ENDIF
775         !
776      ENDIF
777      !
778   END SUBROUTINE ssh_asm_inc
779
780
781   SUBROUTINE ssh_asm_div( kt, Kbb, Kmm, phdivn )
782      !!----------------------------------------------------------------------
783      !!                  ***  ROUTINE ssh_asm_div  ***
784      !!
785      !! ** Purpose :   ssh increment with z* is incorporated via a correction of the local divergence         
786      !!                across all the water column
787      !!
788      !! ** Method  :
789      !!                CAUTION : sshiau is positive (inflow) decreasing the
790      !!                          divergence and expressed in m/s
791      !!
792      !! ** Action  :   phdivn   decreased by the ssh increment
793      !!----------------------------------------------------------------------
794      INTEGER, INTENT(IN) :: kt                               ! ocean time-step index
795      INTEGER, INTENT(IN) :: Kbb, Kmm                         ! time level indices
796      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
797      !!
798      INTEGER  ::   jk                                        ! dummy loop index
799      REAL(wp), DIMENSION(:,:)  , POINTER       ::   ztim     ! local array
800      !!----------------------------------------------------------------------
801      !
802#if defined key_asminc
803      CALL ssh_asm_inc( kt, Kbb, Kmm ) !==   (calculate increments)
804      !
805      IF( ln_linssh ) THEN
806         phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t(:,:,1,Kmm) * tmask(:,:,1)
807      ELSE
808         ALLOCATE( ztim(jpi,jpj) )
809         ztim(:,:) = ssh_iau(:,:) / ( ht(:,:) + 1.0 - ssmask(:,:) )
810         DO jk = 1, jpkm1                                 
811            phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) 
812         END DO
813         !
814         DEALLOCATE(ztim)
815      ENDIF
816#endif
817      !
818   END SUBROUTINE ssh_asm_div
819
820
821   SUBROUTINE seaice_asm_inc( kt, kindic )
822      !!----------------------------------------------------------------------
823      !!                    ***  ROUTINE seaice_asm_inc  ***
824      !!         
825      !! ** Purpose : Apply the sea ice assimilation increment.
826      !!
827      !! ** Method  : Direct initialization or Incremental Analysis Updating.
828      !!
829      !! ** Action  :
830      !!
831      !!----------------------------------------------------------------------
832      INTEGER, INTENT(in)           ::   kt       ! Current time step
833      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation
834      !
835      INTEGER  ::   it
836      REAL(wp) ::   zincwgt   ! IAU weight for current time step
837#if defined key_si3
838      REAL(wp), DIMENSION(jpi,jpj) ::   zofrld, zohicif, zseaicendg, zhicifinc
839      REAL(wp) ::   zhicifmin = 0.5_wp      ! ice minimum depth in metres
840#endif
841      !!----------------------------------------------------------------------
842      !
843      !                             !-----------------------------------------
844      IF ( ln_asmiau ) THEN         ! Incremental Analysis Updating
845         !                          !-----------------------------------------
846         !
847         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
848            !
849            it = kt - nit000 + 1
850            zincwgt = wgtiau(it)      ! IAU weight for the current time step
851            ! note this is not a tendency so should not be divided by rn_Dt (as with the tracer and other increments)
852            !
853            IF(lwp) THEN
854               WRITE(numout,*) 
855               WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
856               WRITE(numout,*) '~~~~~~~~~~~~'
857            ENDIF
858            !
859            ! Sea-ice : SI3 case
860            !
861#if defined key_si3
862            zofrld (:,:) = 1._wp - at_i(:,:)
863            zohicif(:,:) = hm_i(:,:)
864            !
865            at_i  (:,:) = 1. - MIN( MAX( 1.-at_i  (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
866            at_i_b(:,:) = 1. - MIN( MAX( 1.-at_i_b(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
867            fr_i(:,:) = at_i(:,:)        ! adjust ice fraction
868            !
869            zseaicendg(:,:) = zofrld(:,:) - (1. - at_i(:,:))   ! find out actual sea ice nudge applied
870            !
871            ! Nudge sea ice depth to bring it up to a required minimum depth
872            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 
873               zhicifinc(:,:) = (zhicifmin - hm_i(:,:)) * zincwgt   
874            ELSEWHERE
875               zhicifinc(:,:) = 0.0_wp
876            END WHERE
877            !
878            ! nudge ice depth
879            hm_i (:,:) = hm_i (:,:) + zhicifinc(:,:)
880            !
881            ! seaice salinity balancing (to add)
882#endif
883            !
884#if defined key_cice && defined key_asminc
885            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
886            ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rn_Dt
887#endif
888            !
889            IF ( kt == nitiaufin_r ) THEN
890               DEALLOCATE( seaice_bkginc )
891            ENDIF
892            !
893         ELSE
894            !
895#if defined key_cice && defined key_asminc
896            ndaice_da(:,:) = 0._wp        ! Sea-ice : CICE case. Zero ice increment tendency into CICE
897#endif
898            !
899         ENDIF
900         !                          !-----------------------------------------
901      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
902         !                          !-----------------------------------------
903         !
904         IF ( kt == nitdin_r ) THEN
905            !
906            l_1st_euler = .TRUE.              ! Force Euler forward step
907            !
908            ! Sea-ice : SI3 case
909            !
910#if defined key_si3
911            zofrld (:,:) = 1._wp - at_i(:,:)
912            zohicif(:,:) = hm_i(:,:)
913            !
914            ! Initialize the now fields the background + increment
915            at_i(:,:) = 1. - MIN( MAX( 1.-at_i(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp)
916            at_i_b(:,:) = at_i(:,:) 
917            fr_i(:,:) = at_i(:,:)        ! adjust ice fraction
918            !
919            zseaicendg(:,:) = zofrld(:,:) - (1. - at_i(:,:))   ! find out actual sea ice nudge applied
920            !
921            ! Nudge sea ice depth to bring it up to a required minimum depth
922            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 
923               zhicifinc(:,:) = zhicifmin - hm_i(:,:)
924            ELSEWHERE
925               zhicifinc(:,:) = 0.0_wp
926            END WHERE
927            !
928            ! nudge ice depth
929            hm_i (:,:) = hm_i (:,:) + zhicifinc(:,:)
930            !
931            ! seaice salinity balancing (to add)
932#endif
933            !
934#if defined key_cice && defined key_asminc
935            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
936           ndaice_da(:,:) = seaice_bkginc(:,:) / rn_Dt
937#endif
938            IF ( .NOT. PRESENT(kindic) ) THEN
939               DEALLOCATE( seaice_bkginc )
940            END IF
941            !
942         ELSE
943            !
944#if defined key_cice && defined key_asminc
945            ndaice_da(:,:) = 0._wp     ! Sea-ice : CICE case. Zero ice increment tendency into CICE
946#endif
947            !
948         ENDIF
949
950!#if defined defined key_si3 || defined key_cice
951!
952!            IF (ln_seaicebal ) THEN       
953!             !! balancing salinity increments
954!             !! simple case from limflx.F90 (doesn't include a mass flux)
955!             !! assumption is that as ice concentration is reduced or increased
956!             !! the snow and ice depths remain constant
957!             !! note that snow is being created where ice concentration is being increased
958!             !! - could be more sophisticated and
959!             !! not do this (but would need to alter h_snow)
960!
961!             usave(:,:,:)=sb(:,:,:)   ! use array as a temporary store
962!
963!             DO jj = 1, jpj
964!               DO ji = 1, jpi
965!           ! calculate change in ice and snow mass per unit area
966!           ! positive values imply adding salt to the ocean (results from ice formation)
967!           ! fwf : ice formation and melting
968!
969!                 zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rn_Dt
970!
971!           ! change salinity down to mixed layer depth
972!                 mld=hmld_kara(ji,jj)
973!
974!           ! prevent small mld
975!           ! less than 10m can cause salinity instability
976!                 IF (mld < 10) mld=10
977!
978!           ! set to bottom of a level
979!                 DO jk = jpk-1, 2, -1
980!                   IF ((mld > gdepw(ji,jj,jk,Kmm)) .and. (mld < gdepw(ji,jj,jk+1,Kmm))) THEN
981!                     mld=gdepw(ji,jj,jk+1,Kmm)
982!                     jkmax=jk
983!                   ENDIF
984!                 ENDDO
985!
986!            ! avoid applying salinity balancing in shallow water or on land
987!            !
988!
989!            ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m)
990!
991!                 dsal_ocn=0.0_wp
992!                 sal_thresh=5.0_wp        ! minimum salinity threshold for salinity balancing
993!
994!                 if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) &
995!                              dsal_ocn = zfons / (rhop(ji,jj,1) * mld)
996!
997!           ! put increments in for levels in the mixed layer
998!           ! but prevent salinity below a threshold value
999!
1000!                   DO jk = 1, jkmax             
1001!
1002!                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN
1003!                           sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn
1004!                           sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn
1005!                     ENDIF
1006!
1007!                   ENDDO
1008!
1009!      !            !  salt exchanges at the ice/ocean interface
1010!      !            zpmess         = zfons / rDt_ice    ! rDt_ice is ice timestep
1011!      !
1012!      !! Adjust fsalt. A +ve fsalt means adding salt to ocean
1013!      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt 
1014!      !!               
1015!      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d)
1016!      !!                                                     ! E-P (kg m-2 s-2)
1017!      !            emp(ji,jj) = emp(ji,jj) + zpmess          ! E-P (kg m-2 s-2)
1018!               ENDDO !ji
1019!             ENDDO !jj!
1020!
1021!            ENDIF !ln_seaicebal
1022!
1023!#endif
1024         !
1025      ENDIF
1026      !
1027   END SUBROUTINE seaice_asm_inc
1028   
1029   !!======================================================================
1030END MODULE asminc
Note: See TracBrowser for help on using the repository browser.