New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
asminc.F90 in branches/UKMO/dev_r5518_GO6_package_inc_asm/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/dev_r5518_GO6_package_inc_asm/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 8043

Last change on this file since 8043 was 8043, checked in by jwhile, 7 years ago

Further merge of r5518_sshiau_off

File size: 49.1 KB
Line 
1MODULE asminc
2   !!======================================================================
3   !!                       ***  MODULE asminc  ***
4   !! Assimilation increment : Apply an increment generated by data
5   !!                          assimilation
6   !!======================================================================
7   !! History :       ! 2007-03  (M. Martin)  Met Office version
8   !!                 ! 2007-04  (A. Weaver)  calc_date original code
9   !!                 ! 2007-04  (A. Weaver)  Merge with OPAVAR/NEMOVAR
10   !!   NEMO     3.3  ! 2010-05  (D. Lea)  Update to work with NEMO v3.2
11   !!             -   ! 2010-05  (D. Lea)  add calc_month_len routine based on day_init
12   !!            3.4  ! 2012-10  (A. Weaver and K. Mogensen) Fix for direct initialization
13   !!----------------------------------------------------------------------
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      !
774   END SUBROUTINE tra_asm_inc
775
776
777   SUBROUTINE dyn_asm_inc( kt )
778      !!----------------------------------------------------------------------
779      !!                    ***  ROUTINE dyn_asm_inc  ***
780      !!         
781      !! ** Purpose : Apply the dynamics (u and v) assimilation increments.
782      !!
783      !! ** Method  : Direct initialization or Incremental Analysis Updating.
784      !!
785      !! ** Action  :
786      !!----------------------------------------------------------------------
787      INTEGER, INTENT(IN) :: kt   ! Current time step
788      !
789      INTEGER :: jk
790      INTEGER :: it
791      REAL(wp) :: zincwgt  ! IAU weight for current time step
792      !!----------------------------------------------------------------------
793
794      IF ( ln_asmiau ) THEN
795
796         !--------------------------------------------------------------------
797         ! Incremental Analysis Updating
798         !--------------------------------------------------------------------
799
800         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
801
802            it = kt - nit000 + 1
803            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
804
805            IF(lwp) THEN
806               WRITE(numout,*) 
807               WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', &
808                  &  kt,' with IAU weight = ', wgtiau(it)
809               WRITE(numout,*) '~~~~~~~~~~~~'
810            ENDIF
811
812            ! Update the dynamic tendencies
813            DO jk = 1, jpkm1
814               ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt
815               va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt
816            END DO
817           
818            IF ( kt == nitiaufin_r ) THEN
819               DEALLOCATE( u_bkginc )
820               DEALLOCATE( v_bkginc )
821            ENDIF
822
823         ENDIF
824
825      ELSEIF ( ln_asmdin ) THEN 
826
827         !--------------------------------------------------------------------
828         ! Direct Initialization
829         !--------------------------------------------------------------------
830         
831         IF ( kt == nitdin_r ) THEN
832
833            neuler = 0                    ! Force Euler forward step
834
835            ! Initialize the now fields with the background + increment
836            un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:)
837            vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:) 
838
839            ub(:,:,:) = un(:,:,:)         ! Update before fields
840            vb(:,:,:) = vn(:,:,:)
841 
842            DEALLOCATE( u_bkg    )
843            DEALLOCATE( v_bkg    )
844            DEALLOCATE( u_bkginc )
845            DEALLOCATE( v_bkginc )
846         ENDIF
847         !
848      ENDIF
849      !
850   END SUBROUTINE dyn_asm_inc
851
852
853   SUBROUTINE ssh_asm_inc( kt )
854      !!----------------------------------------------------------------------
855      !!                    ***  ROUTINE ssh_asm_inc  ***
856      !!         
857      !! ** Purpose : Apply the sea surface height assimilation increment.
858      !!
859      !! ** Method  : Direct initialization or Incremental Analysis Updating.
860      !!
861      !! ** Action  :
862      !!----------------------------------------------------------------------
863      INTEGER, INTENT(IN) :: kt   ! Current time step
864      !
865      INTEGER :: it
866      INTEGER :: jk
867      REAL(wp) :: zincwgt  ! IAU weight for current time step
868      !!----------------------------------------------------------------------
869
870      IF ( ln_asmiau ) THEN
871
872         !--------------------------------------------------------------------
873         ! Incremental Analysis Updating
874         !--------------------------------------------------------------------
875
876         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
877
878            it = kt - nit000 + 1
879            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
880
881            IF(lwp) THEN
882               WRITE(numout,*) 
883               WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
884                  &  kt,' with IAU weight = ', wgtiau(it)
885               WRITE(numout,*) '~~~~~~~~~~~~'
886            ENDIF
887
888            ! Save the tendency associated with the IAU weighted SSH increment
889            ! (applied in dynspg.*)
890#if defined key_asminc
891            ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
892#endif
893            IF ( kt == nitiaufin_r ) THEN
894               DEALLOCATE( ssh_bkginc )
895            ENDIF
896
897         ELSE
898
899            ssh_iau(:,:) = 0.0
900
901         ENDIF
902
903      ELSEIF ( ln_asmdin ) THEN
904
905         !--------------------------------------------------------------------
906         ! Direct Initialization
907         !--------------------------------------------------------------------
908
909         IF ( kt == nitdin_r ) THEN
910
911            neuler = 0                    ! Force Euler forward step
912
913            ! Initialize the now fields the background + increment
914            sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) 
915
916            ! Update before fields
917            sshb(:,:) = sshn(:,:)         
918
919            IF( lk_vvl ) THEN
920               DO jk = 1, jpk
921                  fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
922               END DO
923            ENDIF
924
925            DEALLOCATE( ssh_bkg    )
926            DEALLOCATE( ssh_bkginc )
927
928         ENDIF
929         !
930      ENDIF
931      !
932   END SUBROUTINE ssh_asm_inc
933
934
935   SUBROUTINE seaice_asm_inc( kt, kindic )
936      !!----------------------------------------------------------------------
937      !!                    ***  ROUTINE seaice_asm_inc  ***
938      !!         
939      !! ** Purpose : Apply the sea ice assimilation increment.
940      !!
941      !! ** Method  : Direct initialization or Incremental Analysis Updating.
942      !!
943      !! ** Action  :
944      !!
945      !!----------------------------------------------------------------------
946      IMPLICIT NONE
947      !
948      INTEGER, INTENT(in)           ::   kt   ! Current time step
949      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation
950      !
951      INTEGER  ::   it
952      REAL(wp) ::   zincwgt   ! IAU weight for current time step
953#if defined key_lim2
954      REAL(wp), DIMENSION(jpi,jpj) ::   zofrld, zohicif, zseaicendg, zhicifinc  ! LIM
955      REAL(wp) ::   zhicifmin = 0.5_wp      ! ice minimum depth in metres
956#endif
957      !!----------------------------------------------------------------------
958
959      IF ( ln_asmiau ) THEN
960
961         !--------------------------------------------------------------------
962         ! Incremental Analysis Updating
963         !--------------------------------------------------------------------
964
965         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
966
967            it = kt - nit000 + 1
968            zincwgt = wgtiau(it)      ! IAU weight for the current time step
969            ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments)
970
971            IF(lwp) THEN
972               WRITE(numout,*) 
973               WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', &
974                  &  kt,' with IAU weight = ', wgtiau(it)
975               WRITE(numout,*) '~~~~~~~~~~~~'
976            ENDIF
977
978            ! Sea-ice : LIM-3 case (to add)
979
980#if defined key_lim2
981            ! Sea-ice : LIM-2 case
982            zofrld (:,:) = frld(:,:)
983            zohicif(:,:) = hicif(:,:)
984            !
985            frld  = MIN( MAX( frld (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
986            pfrld = MIN( MAX( pfrld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
987            fr_i(:,:) = 1.0_wp - frld(:,:)        ! adjust ice fraction
988            !
989            zseaicendg(:,:) = zofrld(:,:) - frld(:,:)   ! find out actual sea ice nudge applied
990            !
991            ! Nudge sea ice depth to bring it up to a required minimum depth
992            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 
993               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt   
994            ELSEWHERE
995               zhicifinc(:,:) = 0.0_wp
996            END WHERE
997            !
998            ! nudge ice depth
999            hicif (:,:) = hicif (:,:) + zhicifinc(:,:)
1000            phicif(:,:) = phicif(:,:) + zhicifinc(:,:)       
1001            !
1002            ! seaice salinity balancing (to add)
1003#endif
1004
1005#if defined key_cice && defined key_asminc
1006            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
1007            ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt
1008#endif
1009
1010            IF ( kt == nitiaufin_r ) THEN
1011               DEALLOCATE( seaice_bkginc )
1012            ENDIF
1013
1014         ELSE
1015
1016#if defined key_cice && defined key_asminc
1017            ! Sea-ice : CICE case. Zero ice increment tendency into CICE
1018            ndaice_da(:,:) = 0.0_wp
1019#endif
1020
1021         ENDIF
1022
1023      ELSEIF ( ln_asmdin ) THEN
1024
1025         !--------------------------------------------------------------------
1026         ! Direct Initialization
1027         !--------------------------------------------------------------------
1028
1029         IF ( kt == nitdin_r ) THEN
1030
1031            neuler = 0                    ! Force Euler forward step
1032
1033            ! Sea-ice : LIM-3 case (to add)
1034
1035#if defined key_lim2
1036            ! Sea-ice : LIM-2 case.
1037            zofrld(:,:)=frld(:,:)
1038            zohicif(:,:)=hicif(:,:)
1039            !
1040            ! Initialize the now fields the background + increment
1041            frld (:,:) = MIN( MAX( frld(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp)
1042            pfrld(:,:) = frld(:,:) 
1043            fr_i (:,:) = 1.0_wp - frld(:,:)                ! adjust ice fraction
1044            zseaicendg(:,:) = zofrld(:,:) - frld(:,:)      ! find out actual sea ice nudge applied
1045            !
1046            ! Nudge sea ice depth to bring it up to a required minimum depth
1047            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 
1048               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt   
1049            ELSEWHERE
1050               zhicifinc(:,:) = 0.0_wp
1051            END WHERE
1052            !
1053            ! nudge ice depth
1054            hicif (:,:) = hicif (:,:) + zhicifinc(:,:)
1055            phicif(:,:) = phicif(:,:)       
1056            !
1057            ! seaice salinity balancing (to add)
1058#endif
1059 
1060#if defined key_cice && defined key_asminc
1061            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
1062           ndaice_da(:,:) = seaice_bkginc(:,:) / rdt
1063#endif
1064           IF ( .NOT. PRESENT(kindic) ) THEN
1065              DEALLOCATE( seaice_bkginc )
1066           END IF
1067
1068         ELSE
1069
1070#if defined key_cice && defined key_asminc
1071            ! Sea-ice : CICE case. Zero ice increment tendency into CICE
1072            ndaice_da(:,:) = 0.0_wp
1073#endif
1074         
1075         ENDIF
1076
1077!#if defined defined key_lim2 || defined key_cice
1078!
1079!            IF (ln_seaicebal ) THEN       
1080!             !! balancing salinity increments
1081!             !! simple case from limflx.F90 (doesn't include a mass flux)
1082!             !! assumption is that as ice concentration is reduced or increased
1083!             !! the snow and ice depths remain constant
1084!             !! note that snow is being created where ice concentration is being increased
1085!             !! - could be more sophisticated and
1086!             !! not do this (but would need to alter h_snow)
1087!
1088!             usave(:,:,:)=sb(:,:,:)   ! use array as a temporary store
1089!
1090!             DO jj = 1, jpj
1091!               DO ji = 1, jpi
1092!           ! calculate change in ice and snow mass per unit area
1093!           ! positive values imply adding salt to the ocean (results from ice formation)
1094!           ! fwf : ice formation and melting
1095!
1096!                 zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rdt
1097!
1098!           ! change salinity down to mixed layer depth
1099!                 mld=hmld_kara(ji,jj)
1100!
1101!           ! prevent small mld
1102!           ! less than 10m can cause salinity instability
1103!                 IF (mld < 10) mld=10
1104!
1105!           ! set to bottom of a level
1106!                 DO jk = jpk-1, 2, -1
1107!                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN
1108!                     mld=gdepw(ji,jj,jk+1)
1109!                     jkmax=jk
1110!                   ENDIF
1111!                 ENDDO
1112!
1113!            ! avoid applying salinity balancing in shallow water or on land
1114!            !
1115!
1116!            ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m)
1117!
1118!                 dsal_ocn=0.0_wp
1119!                 sal_thresh=5.0_wp        ! minimum salinity threshold for salinity balancing
1120!
1121!                 if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) &
1122!                              dsal_ocn = zfons / (rhop(ji,jj,1) * mld)
1123!
1124!           ! put increments in for levels in the mixed layer
1125!           ! but prevent salinity below a threshold value
1126!
1127!                   DO jk = 1, jkmax             
1128!
1129!                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN
1130!                           sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn
1131!                           sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn
1132!                     ENDIF
1133!
1134!                   ENDDO
1135!
1136!      !            !  salt exchanges at the ice/ocean interface
1137!      !            zpmess         = zfons / rdt_ice    ! rdt_ice is ice timestep
1138!      !
1139!      !! Adjust fsalt. A +ve fsalt means adding salt to ocean
1140!      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt 
1141!      !!               
1142!      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d)
1143!      !!                                                     ! E-P (kg m-2 s-2)
1144!      !            emp(ji,jj) = emp(ji,jj) + zpmess          ! E-P (kg m-2 s-2)
1145!               ENDDO !ji
1146!             ENDDO !jj!
1147!
1148!            ENDIF !ln_seaicebal
1149!
1150!#endif
1151
1152      ENDIF
1153
1154   END SUBROUTINE seaice_asm_inc
1155   
1156   !!======================================================================
1157END MODULE asminc
Note: See TracBrowser for help on using the repository browser.