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

Last change on this file since 13318 was 13318, checked in by dford, 5 months ago

Allow assimilation of PFT chlorophyll profiles. See Met Office utils ticket 346.

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