source: branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 12611

Last change on this file since 12611 was 12611, checked in by dcarneir, 8 months ago

Including sea ice thickness variable

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