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/2015/dev_r5056_CMCC4_simplification/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/2015/dev_r5056_CMCC4_simplification/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 5282

Last change on this file since 5282 was 5282, checked in by diovino, 9 years ago

Dev. branch CMCC4_simplification ticket #1456

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