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/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 3784

Last change on this file since 3784 was 3784, checked in by gm, 9 years ago

dev_v3_4_STABLE_2012: #1046 : bug correction to compile LIM3 with ASM

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