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

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

source: branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 12609

Last change on this file since 12609 was 12609, checked in by dcarneir, 2 years ago

Making SIT changes in global configs to be compatible with build in shelf configs

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