source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_sit/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 11952

Last change on this file since 11952 was 11952, checked in by dcarneir, 13 months ago

Compilation errors fixed in asminc.F90

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