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 @ 12594

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

Removing commented code in sit_asm_inc

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