source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 8031

Last change on this file since 8031 was 8031, checked in by lovato, 3 years ago

3.6 stable: bugfix for ticket #1528

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