source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 8400

Last change on this file since 8400 was 8400, checked in by timgraham, 3 years ago

GMED ticket 335:

  • Merge dev_r5518_GO6_package_inc_asm into package branch to make everything easier for data assimilation
  • No effect on configs without data assimilation
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#if defined key_cice && defined key_asminc
42   USE sbc_ice, ONLY : & ! CICE Ice model variables
43   & ndaice_da, nfresh_da, nfsalt_da
44#endif
45   USE sbc_oce          ! Surface boundary condition variables.
46
47   IMPLICIT NONE
48   PRIVATE
49   
50   PUBLIC   asm_inc_init   !: Initialize the increment arrays and IAU weights
51   PUBLIC   calc_date      !: Compute the calendar date YYYYMMDD on a given step
52   PUBLIC   tra_asm_inc    !: Apply the tracer (T and S) increments
53   PUBLIC   dyn_asm_inc    !: Apply the dynamic (u and v) increments
54   PUBLIC   ssh_asm_inc    !: Apply the SSH increment
55   PUBLIC   seaice_asm_inc !: Apply the seaice increment
56
57#if defined key_asminc
58    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .TRUE.   !: Logical switch for assimilation increment interface
59#else
60    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE.  !: No assimilation increments
61#endif
62   LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.      !: No output of the background state fields
63   LOGICAL, PUBLIC :: ln_asmiau = .FALSE.      !: No applying forcing with an assimilation increment
64   LOGICAL, PUBLIC :: ln_asmdin = .FALSE.      !: No direct initialization
65   LOGICAL, PUBLIC :: ln_trainc = .FALSE.      !: No tracer (T and S) assimilation increments
66   LOGICAL, PUBLIC :: ln_dyninc = .FALSE.      !: No dynamics (u and v) assimilation increments
67   LOGICAL, PUBLIC :: ln_sshinc = .FALSE.      !: No sea surface height assimilation increment
68   LOGICAL, PUBLIC :: ln_seaiceinc             !: No sea ice concentration increment
69   LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check
70   LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
71   INTEGER, PUBLIC :: nn_divdmp                !: Apply divergence damping filter nn_divdmp times
72
73   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkg   , s_bkg      !: Background temperature and salinity
74   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkg   , v_bkg      !: Background u- & v- velocity components
75   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkginc, s_bkginc   !: Increment to the background T & S
76   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components
77   REAL(wp), PUBLIC, DIMENSION(:)    , ALLOCATABLE ::   wgtiau               !: IAU weights for each time step
78#if defined key_asminc
79   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssh_iau           !: IAU-weighted sea surface height increment
80#endif
81   !                                !!! time steps relative to the cycle interval [0,nitend-nit000-1]
82   INTEGER , PUBLIC ::   nitbkg      !: Time step of the background state used in the Jb term
83   INTEGER , PUBLIC ::   nitdin      !: Time step of the background state for direct initialization
84   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval
85   INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval
86   !
87   INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting
88   !                                 !: = 1   Linear hat-like, centred in middle of IAU interval
89   REAL(wp), PUBLIC ::   salfixmin   !: Ensure that the salinity is larger than this  value if (ln_salfix)
90
91   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment
92   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc
93
94   !! * Substitutions
95#  include "domzgr_substitute.h90"
96#  include "ldfdyn_substitute.h90"
97#  include "vectopt_loop_substitute.h90"
98   !!----------------------------------------------------------------------
99   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
100   !! $Id$
101   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
102   !!----------------------------------------------------------------------
103CONTAINS
104
105   SUBROUTINE asm_inc_init
106      !!----------------------------------------------------------------------
107      !!                    ***  ROUTINE asm_inc_init  ***
108      !!         
109      !! ** Purpose : Initialize the assimilation increment and IAU weights.
110      !!
111      !! ** Method  : Initialize the assimilation increment and IAU weights.
112      !!
113      !! ** Action  :
114      !!----------------------------------------------------------------------
115      INTEGER :: ji, jj, jk, jt  ! dummy loop indices
116      INTEGER :: imid, inum      ! local integers
117      INTEGER :: ios             ! Local integer output status for namelist read
118      INTEGER :: iiauper         ! Number of time steps in the IAU period
119      INTEGER :: icycper         ! Number of time steps in the cycle
120      INTEGER :: iitend_date     ! Date YYYYMMDD of final time step
121      INTEGER :: iitbkg_date     ! Date YYYYMMDD of background time step for Jb term
122      INTEGER :: iitdin_date     ! Date YYYYMMDD of background time step for DI
123      INTEGER :: iitiaustr_date  ! Date YYYYMMDD of IAU interval start time step
124      INTEGER :: iitiaufin_date  ! Date YYYYMMDD of IAU interval final time step
125      !
126      REAL(wp) :: znorm        ! Normalization factor for IAU weights
127      REAL(wp) :: ztotwgt      ! Value of time-integrated IAU weights (should be equal to one)
128      REAL(wp) :: z_inc_dateb  ! Start date of interval on which increment is valid
129      REAL(wp) :: z_inc_datef  ! End date of interval on which increment is valid
130      REAL(wp) :: zdate_bkg    ! Date in background state file for DI
131      REAL(wp) :: zdate_inc    ! Time axis in increments file
132      !
133      REAL(wp), POINTER, DIMENSION(:,:) ::   hdiv   ! 2D workspace
134      !!
135      NAMELIST/nam_asminc/ ln_bkgwri,                                      &
136         &                 ln_trainc, ln_dyninc, ln_sshinc,                &
137         &                 ln_asmdin, ln_asmiau,                           &
138         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
139         &                 ln_salfix, salfixmin, nn_divdmp,                &
140         &                 ln_seaiceinc, ln_temnofreeze
141      !!----------------------------------------------------------------------
142
143      !-----------------------------------------------------------------------
144      ! Read Namelist nam_asminc : assimilation increment interface
145      !-----------------------------------------------------------------------
146      ln_seaiceinc = .FALSE.
147      ln_temnofreeze = .FALSE.
148
149      REWIND( numnam_ref )              ! Namelist nam_asminc in reference namelist : Assimilation increment
150      READ  ( numnam_ref, nam_asminc, IOSTAT = ios, ERR = 901)
151901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_asminc in reference namelist', lwp )
152
153      REWIND( numnam_cfg )              ! Namelist nam_asminc in configuration namelist : Assimilation increment
154      READ  ( numnam_cfg, nam_asminc, IOSTAT = ios, ERR = 902 )
155902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_asminc in configuration namelist', lwp )
156      IF(lwm) WRITE ( numond, nam_asminc )
157
158      ! Control print
159      IF(lwp) THEN
160         WRITE(numout,*)
161         WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :'
162         WRITE(numout,*) '~~~~~~~~~~~~'
163         WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters'
164         WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri
165         WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc
166         WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc
167         WRITE(numout,*) '      Logical switch for applying SSH increments               ln_sshinc = ', ln_sshinc
168         WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin
169         WRITE(numout,*) '      Logical switch for applying sea ice increments        ln_seaiceinc = ', ln_seaiceinc
170         WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau
171         WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg
172         WRITE(numout,*) '      Timestep of background for DI in [0,nitend-nit000-1]     nitdin    = ', nitdin
173         WRITE(numout,*) '      Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr
174         WRITE(numout,*) '      Timestep of end of IAU interval in [0,nitend-nit000-1]   nitiaufin = ', nitiaufin
175         WRITE(numout,*) '      Type of IAU weighting function                           niaufn    = ', niaufn
176         WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix
177         WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin
178      ENDIF
179
180      nitbkg_r    = nitbkg    + nit000 - 1  ! Background time referenced to nit000
181      nitdin_r    = nitdin    + nit000 - 1  ! Background time for DI referenced to nit000
182      nitiaustr_r = nitiaustr + nit000 - 1  ! Start of IAU interval referenced to nit000
183      nitiaufin_r = nitiaufin + nit000 - 1  ! End of IAU interval referenced to nit000
184
185      iiauper = nitiaufin_r - nitiaustr_r + 1  ! IAU interval length
186      icycper = nitend      - nit000      + 1  ! Cycle interval length
187
188      CALL calc_date( nit000, nitend     , ndate0, iitend_date    )     ! Date of final time step
189      CALL calc_date( nit000, nitbkg_r   , ndate0, iitbkg_date    )     ! Background time for Jb referenced to ndate0
190      CALL calc_date( nit000, nitdin_r   , ndate0, iitdin_date    )     ! Background time for DI referenced to ndate0
191      CALL calc_date( nit000, nitiaustr_r, ndate0, iitiaustr_date )     ! IAU start time referenced to ndate0
192      CALL calc_date( nit000, nitiaufin_r, ndate0, iitiaufin_date )     ! IAU end time referenced to ndate0
193      !
194      IF(lwp) THEN
195         WRITE(numout,*)
196         WRITE(numout,*) '   Time steps referenced to current cycle:'
197         WRITE(numout,*) '       iitrst      = ', nit000 - 1
198         WRITE(numout,*) '       nit000      = ', nit000
199         WRITE(numout,*) '       nitend      = ', nitend
200         WRITE(numout,*) '       nitbkg_r    = ', nitbkg_r
201         WRITE(numout,*) '       nitdin_r    = ', nitdin_r
202         WRITE(numout,*) '       nitiaustr_r = ', nitiaustr_r
203         WRITE(numout,*) '       nitiaufin_r = ', nitiaufin_r
204         WRITE(numout,*)
205         WRITE(numout,*) '   Dates referenced to current cycle:'
206         WRITE(numout,*) '       ndastp         = ', ndastp
207         WRITE(numout,*) '       ndate0         = ', ndate0
208         WRITE(numout,*) '       iitend_date    = ', iitend_date
209         WRITE(numout,*) '       iitbkg_date    = ', iitbkg_date
210         WRITE(numout,*) '       iitdin_date    = ', iitdin_date
211         WRITE(numout,*) '       iitiaustr_date = ', iitiaustr_date
212         WRITE(numout,*) '       iitiaufin_date = ', iitiaufin_date
213      ENDIF
214
215      IF ( nacc /= 0 ) &
216         & CALL ctl_stop( ' nacc /= 0 and key_asminc :',  &
217         &                ' Assimilation increments have only been implemented', &
218         &                ' for synchronous time stepping' )
219
220      IF ( ( ln_asmdin ).AND.( ln_asmiau ) )   &
221         & CALL ctl_stop( ' ln_asmdin and ln_asmiau :', &
222         &                ' Choose Direct Initialization OR Incremental Analysis Updating')
223
224      IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
225           .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) &
226         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', &
227         &                ' but ln_asmdin and ln_asmiau are both set to .false. :', &
228         &                ' Inconsistent options')
229
230      IF ( ( niaufn /= 0 ).AND.( niaufn /= 1 ) ) &
231         & CALL ctl_stop( ' niaufn /= 0 or niaufn /=1 :',  &
232         &                ' Type IAU weighting function is invalid')
233
234      IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) &
235         &                     )  &
236         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', &
237         &                ' The assimilation increments are not applied')
238
239      IF ( ( ln_asmiau ).AND.( nitiaustr == nitiaufin ) ) &
240         & CALL ctl_stop( ' nitiaustr = nitiaufin :',  &
241         &                ' IAU interval is of zero length')
242
243      IF ( ( ln_asmiau ).AND.( ( nitiaustr_r < nit000 ).OR.( nitiaufin_r > nitend ) ) ) &
244         & CALL ctl_stop( ' nitiaustr or nitiaufin :',  &
245         &                ' IAU starting or final time step is outside the cycle interval', &
246         &                 ' Valid range nit000 to nitend')
247
248      IF ( ( nitbkg_r < nit000 - 1 ).OR.( nitbkg_r > nitend ) ) &
249         & CALL ctl_stop( ' nitbkg :',  &
250         &                ' Background time step is outside the cycle interval')
251
252      IF ( ( nitdin_r < nit000 - 1 ).OR.( nitdin_r > nitend ) ) &
253         & CALL ctl_stop( ' nitdin :',  &
254         &                ' Background time step for Direct Initialization is outside', &
255         &                ' the cycle interval')
256
257      IF ( nstop > 0 ) RETURN       ! if there are any errors then go no further
258
259      !--------------------------------------------------------------------
260      ! Initialize the Incremental Analysis Updating weighting function
261      !--------------------------------------------------------------------
262
263      IF ( ln_asmiau ) THEN
264
265         ALLOCATE( wgtiau( icycper ) )
266
267         wgtiau(:) = 0.0
268
269         IF ( niaufn == 0 ) THEN
270
271            !---------------------------------------------------------
272            ! Constant IAU forcing
273            !---------------------------------------------------------
274
275            DO jt = 1, iiauper
276               wgtiau(jt+nitiaustr-1) = 1.0 / REAL( iiauper )
277            END DO
278
279         ELSEIF ( niaufn == 1 ) THEN
280
281            !---------------------------------------------------------
282            ! Linear hat-like, centred in middle of IAU interval
283            !---------------------------------------------------------
284
285            ! Compute the normalization factor
286            znorm = 0.0
287            IF ( MOD( iiauper, 2 ) == 0 ) THEN  ! Even number of time steps in IAU interval
288               imid = iiauper / 2 
289               DO jt = 1, imid
290                  znorm = znorm + REAL( jt )
291               END DO
292               znorm = 2.0 * znorm
293            ELSE                               ! Odd number of time steps in IAU interval
294               imid = ( iiauper + 1 ) / 2       
295               DO jt = 1, imid - 1
296                  znorm = znorm + REAL( jt )
297               END DO
298               znorm = 2.0 * znorm + REAL( imid )
299            ENDIF
300            znorm = 1.0 / znorm
301
302            DO jt = 1, imid - 1
303               wgtiau(jt+nitiaustr-1) = REAL( jt ) * znorm
304            END DO
305            DO jt = imid, iiauper
306               wgtiau(jt+nitiaustr-1) = REAL( iiauper - jt + 1 ) * znorm
307            END DO
308
309         ENDIF
310
311         ! Test that the integral of the weights over the weighting interval equals 1
312          IF(lwp) THEN
313             WRITE(numout,*)
314             WRITE(numout,*) 'asm_inc_init : IAU weights'
315             WRITE(numout,*) '~~~~~~~~~~~~'
316             WRITE(numout,*) '             time step         IAU  weight'
317             WRITE(numout,*) '             =========     ====================='
318             ztotwgt = 0.0
319             DO jt = 1, icycper
320                ztotwgt = ztotwgt + wgtiau(jt)
321                WRITE(numout,*) '         ', jt, '       ', wgtiau(jt) 
322             END DO   
323             WRITE(numout,*) '         ==================================='
324             WRITE(numout,*) '         Time-integrated weight = ', ztotwgt
325             WRITE(numout,*) '         ==================================='
326          ENDIF
327         
328      ENDIF
329
330      !--------------------------------------------------------------------
331      ! Allocate and initialize the increment arrays
332      !--------------------------------------------------------------------
333
334      ALLOCATE( t_bkginc(jpi,jpj,jpk) )
335      ALLOCATE( s_bkginc(jpi,jpj,jpk) )
336      ALLOCATE( u_bkginc(jpi,jpj,jpk) )
337      ALLOCATE( v_bkginc(jpi,jpj,jpk) )
338      ALLOCATE( ssh_bkginc(jpi,jpj)   )
339      ALLOCATE( seaice_bkginc(jpi,jpj))
340#if defined key_asminc
341      ALLOCATE( ssh_iau(jpi,jpj)      )
342#endif
343      t_bkginc(:,:,:) = 0.0
344      s_bkginc(:,:,:) = 0.0
345      u_bkginc(:,:,:) = 0.0
346      v_bkginc(:,:,:) = 0.0
347      ssh_bkginc(:,:) = 0.0
348      seaice_bkginc(:,:) = 0.0
349#if defined key_asminc
350      ssh_iau(:,:)    = 0.0
351#endif
352      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN
353
354         !--------------------------------------------------------------------
355         ! Read the increments from file
356         !--------------------------------------------------------------------
357
358         CALL iom_open( c_asminc, inum )
359
360         CALL iom_get( inum, 'time', zdate_inc ) 
361
362         CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb )
363         CALL iom_get( inum, 'z_inc_datef', z_inc_datef )
364         z_inc_dateb = zdate_inc
365         z_inc_datef = zdate_inc
366
367         IF(lwp) THEN
368            WRITE(numout,*) 
369            WRITE(numout,*) 'asm_inc_init : Assimilation increments valid ', &
370               &            ' between dates ', NINT( z_inc_dateb ),' and ',  &
371               &            NINT( z_inc_datef )
372            WRITE(numout,*) '~~~~~~~~~~~~'
373         ENDIF
374
375         IF (     ( NINT( z_inc_dateb ) < ndastp      ) &
376            & .OR.( NINT( z_inc_datef ) > iitend_date ) ) &
377            & CALL ctl_warn( ' Validity time of assimilation increments is ', &
378            &                ' outside the assimilation interval' )
379
380         IF ( ( ln_asmdin ).AND.( NINT( zdate_inc ) /= iitdin_date ) ) &
381            & CALL ctl_warn( ' Validity time of assimilation increments does ', &
382            &                ' not agree with Direct Initialization time' )
383
384         IF ( ln_trainc ) THEN   
385            CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 )
386            CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 )
387            ! Apply the masks
388            t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:)
389            s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:)
390            ! Set missing increments to 0.0 rather than 1e+20
391            ! to allow for differences in masks
392            WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0
393            WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0
394         ENDIF
395
396         IF ( ln_dyninc ) THEN   
397            CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )             
398            CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )             
399            ! Apply the masks
400            u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:)
401            v_bkginc(:,:,:) = v_bkginc(:,:,:) * vmask(:,:,:)
402            ! Set missing increments to 0.0 rather than 1e+20
403            ! to allow for differences in masks
404            WHERE( ABS( u_bkginc(:,:,:) ) > 1.0e+10 ) u_bkginc(:,:,:) = 0.0
405            WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0
406         ENDIF
407       
408         IF ( ln_sshinc ) THEN
409            CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 )
410            ! Apply the masks
411            ssh_bkginc(:,:) = ssh_bkginc(:,:) * tmask(:,:,1)
412            ! Set missing increments to 0.0 rather than 1e+20
413            ! to allow for differences in masks
414            WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0
415         ENDIF
416
417         IF ( ln_seaiceinc ) THEN
418            CALL iom_get( inum, jpdom_autoglo, 'bckinseaice', seaice_bkginc, 1 )
419            ! Apply the masks
420            seaice_bkginc(:,:) = seaice_bkginc(:,:) * tmask(:,:,1)
421            ! Set missing increments to 0.0 rather than 1e+20
422            ! to allow for differences in masks
423            WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0
424         ENDIF
425
426         CALL iom_close( inum )
427 
428      ENDIF
429
430      !-----------------------------------------------------------------------
431      ! Apply divergence damping filter
432      !-----------------------------------------------------------------------
433
434      IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN
435
436         CALL wrk_alloc(jpi,jpj,hdiv) 
437
438         DO  jt = 1, nn_divdmp
439
440            DO jk = 1, jpkm1
441
442               hdiv(:,:) = 0._wp
443
444               DO jj = 2, jpjm1
445                  DO ji = fs_2, fs_jpim1   ! vector opt.
446                     hdiv(ji,jj) =   &
447                        (  e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) * u_bkginc(ji  ,jj  ,jk)     &
448                         - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) * u_bkginc(ji-1,jj  ,jk)     &
449                         + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) * v_bkginc(ji  ,jj  ,jk)     &
450                         - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) * v_bkginc(ji  ,jj-1,jk)  )  &
451                         / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
452                  END DO
453               END DO
454
455               CALL lbc_lnk( hdiv, 'T', 1. )   ! lateral boundary cond. (no sign change)
456
457               DO jj = 2, jpjm1
458                  DO ji = fs_2, fs_jpim1   ! vector opt.
459                     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)   &
460                                                                        - e1t(ji  ,jj)*e2t(ji  ,jj) * hdiv(ji  ,jj) ) &
461                                                                      / e1u(ji,jj) * umask(ji,jj,jk) 
462                     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)   &
463                                                                        - e1t(ji,jj  )*e2t(ji,jj  ) * hdiv(ji,jj  ) ) &
464                                                                      / e2v(ji,jj) * vmask(ji,jj,jk) 
465                  END DO
466               END DO
467
468            END DO
469
470         END DO
471
472         CALL wrk_dealloc(jpi,jpj,hdiv) 
473
474      ENDIF
475
476
477
478      !-----------------------------------------------------------------------
479      ! Allocate and initialize the background state arrays
480      !-----------------------------------------------------------------------
481
482      IF ( ln_asmdin ) THEN
483
484         ALLOCATE( t_bkg(jpi,jpj,jpk) )
485         ALLOCATE( s_bkg(jpi,jpj,jpk) )
486         ALLOCATE( u_bkg(jpi,jpj,jpk) )
487         ALLOCATE( v_bkg(jpi,jpj,jpk) )
488         ALLOCATE( ssh_bkg(jpi,jpj)   )
489
490         t_bkg(:,:,:) = 0.0
491         s_bkg(:,:,:) = 0.0
492         u_bkg(:,:,:) = 0.0
493         v_bkg(:,:,:) = 0.0
494         ssh_bkg(:,:) = 0.0
495
496         !--------------------------------------------------------------------
497         ! Read from file the background state at analysis time
498         !--------------------------------------------------------------------
499
500         CALL iom_open( c_asmdin, inum )
501
502         CALL iom_get( inum, 'rdastp', zdate_bkg ) 
503       
504         IF(lwp) THEN
505            WRITE(numout,*) 
506            WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', &
507               &  NINT( zdate_bkg )
508            WRITE(numout,*) '~~~~~~~~~~~~'
509         ENDIF
510
511         IF ( NINT( zdate_bkg ) /= iitdin_date ) &
512            & CALL ctl_warn( ' Validity time of assimilation background state does', &
513            &                ' not agree with Direct Initialization time' )
514
515         IF ( ln_trainc ) THEN   
516            CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg )
517            CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg )
518            t_bkg(:,:,:) = t_bkg(:,:,:) * tmask(:,:,:)
519            s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:)
520         ENDIF
521
522         IF ( ln_dyninc ) THEN   
523            CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg )
524            CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg )
525            u_bkg(:,:,:) = u_bkg(:,:,:) * umask(:,:,:)
526            v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
527         ENDIF
528       
529         IF ( ln_sshinc ) THEN
530            CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg )
531            ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1)
532         ENDIF
533
534         CALL iom_close( inum )
535
536      ENDIF
537      !
538   END SUBROUTINE asm_inc_init
539
540
541   SUBROUTINE calc_date( kit000, kt, kdate0, kdate )
542      !!----------------------------------------------------------------------
543      !!                    ***  ROUTINE calc_date  ***
544      !!         
545      !! ** Purpose : Compute the calendar date YYYYMMDD at a given time step.
546      !!
547      !! ** Method  : Compute the calendar date YYYYMMDD at a given time step.
548      !!
549      !! ** Action  :
550      !!----------------------------------------------------------------------
551      INTEGER, INTENT(IN) :: kit000  ! Initial time step
552      INTEGER, INTENT(IN) :: kt      ! Current time step referenced to kit000
553      INTEGER, INTENT(IN) :: kdate0  ! Initial date
554      INTEGER, INTENT(OUT) :: kdate  ! Current date reference to kdate0
555      !
556      INTEGER :: iyea0    ! Initial year
557      INTEGER :: imon0    ! Initial month
558      INTEGER :: iday0    ! Initial day
559      INTEGER :: iyea     ! Current year
560      INTEGER :: imon     ! Current month
561      INTEGER :: iday     ! Current day
562      INTEGER :: idaystp  ! Number of days between initial and current date
563      INTEGER :: idaycnt  ! Day counter
564
565      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
566
567      !-----------------------------------------------------------------------
568      ! Compute the calendar date YYYYMMDD
569      !-----------------------------------------------------------------------
570
571      ! Initial date
572      iyea0 =   kdate0 / 10000
573      imon0 = ( kdate0 - ( iyea0 * 10000 ) ) / 100
574      iday0 =   kdate0 - ( iyea0 * 10000 ) - ( imon0 * 100 ) 
575
576      ! Check that kt >= kit000 - 1
577      IF ( kt < kit000 - 1 ) CALL ctl_stop( ' kt must be >= kit000 - 1')
578
579      ! If kt = kit000 - 1 then set the date to the restart date
580      IF ( kt == kit000 - 1 ) THEN
581
582         kdate = ndastp
583         RETURN
584
585      ENDIF
586
587      ! Compute the number of days from the initial date
588      idaystp = INT( REAL( kt - kit000 ) * rdt / 86400. )
589   
590      iday    = iday0
591      imon    = imon0
592      iyea    = iyea0
593      idaycnt = 0
594
595      CALL calc_month_len( iyea, imonth_len )
596
597      DO WHILE ( idaycnt < idaystp )
598         iday = iday + 1
599         IF ( iday > imonth_len(imon) )  THEN
600            iday = 1
601            imon = imon + 1
602         ENDIF
603         IF ( imon > 12 ) THEN
604            imon = 1
605            iyea = iyea + 1
606            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
607         ENDIF                 
608         idaycnt = idaycnt + 1
609      END DO
610      !
611      kdate = iyea * 10000 + imon * 100 + iday
612      !
613   END SUBROUTINE
614
615
616   SUBROUTINE calc_month_len( iyear, imonth_len )
617      !!----------------------------------------------------------------------
618      !!                    ***  ROUTINE calc_month_len  ***
619      !!         
620      !! ** Purpose : Compute the number of days in a months given a year.
621      !!
622      !! ** Method  :
623      !!----------------------------------------------------------------------
624      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
625      INTEGER :: iyear         !: year
626      !!----------------------------------------------------------------------
627      !
628      ! length of the month of the current year (from nleapy, read in namelist)
629      IF ( nleapy < 2 ) THEN
630         imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
631         IF ( nleapy == 1 ) THEN   ! we are using calendar with leap years
632            IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
633               imonth_len(2) = 29
634            ENDIF
635         ENDIF
636      ELSE
637         imonth_len(:) = nleapy   ! all months with nleapy days per year
638      ENDIF
639      !
640   END SUBROUTINE
641
642
643   SUBROUTINE tra_asm_inc( kt )
644      !!----------------------------------------------------------------------
645      !!                    ***  ROUTINE tra_asm_inc  ***
646      !!         
647      !! ** Purpose : Apply the tracer (T and S) assimilation increments
648      !!
649      !! ** Method  : Direct initialization or Incremental Analysis Updating
650      !!
651      !! ** Action  :
652      !!----------------------------------------------------------------------
653      INTEGER, INTENT(IN) :: kt               ! Current time step
654      !
655      INTEGER :: ji,jj,jk
656      INTEGER :: it
657      REAL(wp) :: zincwgt  ! IAU weight for current time step
658      REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values
659      !!----------------------------------------------------------------------
660
661      ! freezing point calculation taken from oc_fz_pt (but calculated for all depths)
662      ! used to prevent the applied increments taking the temperature below the local freezing point
663
664      DO jk = 1, jpkm1
665        CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), fsdept(:,:,jk) )
666      END DO
667
668      IF ( ln_asmiau ) THEN
669
670         !--------------------------------------------------------------------
671         ! Incremental Analysis Updating
672         !--------------------------------------------------------------------
673
674         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
675
676            it = kt - nit000 + 1
677            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
678
679            IF(lwp) THEN
680               WRITE(numout,*) 
681               WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
682               WRITE(numout,*) '~~~~~~~~~~~~'
683            ENDIF
684
685            ! Update the tracer tendencies
686            DO jk = 1, jpkm1
687               IF (ln_temnofreeze) THEN
688                  ! Do not apply negative increments if the temperature will fall below freezing
689                  WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. &
690                     &   tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 
691                     tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 
692                  END WHERE
693               ELSE
694                  tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 
695               ENDIF
696               IF (ln_salfix) THEN
697                  ! Do not apply negative increments if the salinity will fall below a specified
698                  ! minimum value salfixmin
699                  WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. &
700                     &   tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 
701                     tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
702                  END WHERE
703               ELSE
704                  tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
705               ENDIF
706            END DO
707
708         ENDIF
709
710         IF ( kt == nitiaufin_r + 1  ) THEN   ! For bias crcn to work
711            DEALLOCATE( t_bkginc )
712            DEALLOCATE( s_bkginc )
713         ENDIF
714
715
716      ELSEIF ( ln_asmdin ) THEN
717
718         !--------------------------------------------------------------------
719         ! Direct Initialization
720         !--------------------------------------------------------------------
721           
722         IF ( kt == nitdin_r ) THEN
723
724            neuler = 0  ! Force Euler forward step
725
726            ! Initialize the now fields with the background + increment
727            IF (ln_temnofreeze) THEN
728               ! Do not apply negative increments if the temperature will fall below freezing
729               WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 
730                  tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
731               END WHERE
732            ELSE
733               tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
734            ENDIF
735            IF (ln_salfix) THEN
736               ! Do not apply negative increments if the salinity will fall below a specified
737               ! minimum value salfixmin
738               WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 
739                  tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
740               END WHERE
741            ELSE
742               tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
743            ENDIF
744
745            tsb(:,:,:,:) = tsn(:,:,:,:)                 ! Update before fields
746
747            CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )  ! Before potential and in situ densities
748!!gm  fabien
749!            CALL eos( tsb, rhd, rhop )                ! Before potential and in situ densities
750!!gm
751
752
753            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)      &
754               &  CALL zps_hde    ( kt, jpts, tsb, gtsu, gtsv,        &  ! Partial steps: before horizontal gradient
755               &                              rhd, gru , grv          )  ! of t, s, rd at the last ocean level
756            IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)      &
757               &  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,    &    ! Partial steps for top cell (ISF)
758               &                                  rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
759               &                           gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
760
761#if defined key_zdfkpp
762            CALL eos( tsn, rhd, fsdept_n(:,:,:) )                      ! Compute rhd
763!!gm fabien            CALL eos( tsn, rhd )                      ! Compute rhd
764#endif
765
766            DEALLOCATE( t_bkginc )
767            DEALLOCATE( s_bkginc )
768            DEALLOCATE( t_bkg    )
769            DEALLOCATE( s_bkg    )
770         ENDIF
771         
772      ENDIF
773      ! Perhaps the following call should be in step
774      IF   ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment
775      !
776   END SUBROUTINE tra_asm_inc
777
778
779   SUBROUTINE dyn_asm_inc( kt )
780      !!----------------------------------------------------------------------
781      !!                    ***  ROUTINE dyn_asm_inc  ***
782      !!         
783      !! ** Purpose : Apply the dynamics (u and v) assimilation increments.
784      !!
785      !! ** Method  : Direct initialization or Incremental Analysis Updating.
786      !!
787      !! ** Action  :
788      !!----------------------------------------------------------------------
789      INTEGER, INTENT(IN) :: kt   ! Current time step
790      !
791      INTEGER :: jk
792      INTEGER :: it
793      REAL(wp) :: zincwgt  ! IAU weight for current time step
794      !!----------------------------------------------------------------------
795
796      IF ( ln_asmiau ) THEN
797
798         !--------------------------------------------------------------------
799         ! Incremental Analysis Updating
800         !--------------------------------------------------------------------
801
802         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
803
804            it = kt - nit000 + 1
805            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
806
807            IF(lwp) THEN
808               WRITE(numout,*) 
809               WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', &
810                  &  kt,' with IAU weight = ', wgtiau(it)
811               WRITE(numout,*) '~~~~~~~~~~~~'
812            ENDIF
813
814            ! Update the dynamic tendencies
815            DO jk = 1, jpkm1
816               ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt
817               va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt
818            END DO
819           
820            IF ( kt == nitiaufin_r ) THEN
821               DEALLOCATE( u_bkginc )
822               DEALLOCATE( v_bkginc )
823            ENDIF
824
825         ENDIF
826
827      ELSEIF ( ln_asmdin ) THEN 
828
829         !--------------------------------------------------------------------
830         ! Direct Initialization
831         !--------------------------------------------------------------------
832         
833         IF ( kt == nitdin_r ) THEN
834
835            neuler = 0                    ! Force Euler forward step
836
837            ! Initialize the now fields with the background + increment
838            un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:)
839            vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:) 
840
841            ub(:,:,:) = un(:,:,:)         ! Update before fields
842            vb(:,:,:) = vn(:,:,:)
843 
844            DEALLOCATE( u_bkg    )
845            DEALLOCATE( v_bkg    )
846            DEALLOCATE( u_bkginc )
847            DEALLOCATE( v_bkginc )
848         ENDIF
849         !
850      ENDIF
851      !
852   END SUBROUTINE dyn_asm_inc
853
854
855   SUBROUTINE ssh_asm_inc( kt )
856      !!----------------------------------------------------------------------
857      !!                    ***  ROUTINE ssh_asm_inc  ***
858      !!         
859      !! ** Purpose : Apply the sea surface height assimilation increment.
860      !!
861      !! ** Method  : Direct initialization or Incremental Analysis Updating.
862      !!
863      !! ** Action  :
864      !!----------------------------------------------------------------------
865      INTEGER, INTENT(IN) :: kt   ! Current time step
866      !
867      INTEGER :: it
868      INTEGER :: jk
869      REAL(wp) :: zincwgt  ! IAU weight for current time step
870      !!----------------------------------------------------------------------
871
872      IF ( ln_asmiau ) THEN
873
874         !--------------------------------------------------------------------
875         ! Incremental Analysis Updating
876         !--------------------------------------------------------------------
877
878         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
879
880            it = kt - nit000 + 1
881            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
882
883            IF(lwp) THEN
884               WRITE(numout,*) 
885               WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
886                  &  kt,' with IAU weight = ', wgtiau(it)
887               WRITE(numout,*) '~~~~~~~~~~~~'
888            ENDIF
889
890            ! Save the tendency associated with the IAU weighted SSH increment
891            ! (applied in dynspg.*)
892#if defined key_asminc
893            ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
894#endif
895            IF ( kt == nitiaufin_r ) THEN
896               DEALLOCATE( ssh_bkginc )
897            ENDIF
898
899         ELSE
900#if defined key_asminc
901            ssh_iau(:,:) = 0.0
902#endif
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.