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

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

source: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 9180

Last change on this file since 9180 was 9180, checked in by kingr, 6 years ago

Adding changes required to write out time-averaged assimilation background.

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