source: trunk/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 5836

Last change on this file since 5836 was 5836, checked in by cetlod, 6 years ago

merge the simplification branch onto the trunk, see ticket #1612

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