source: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 8440

Last change on this file since 8440 was 8440, checked in by dford, 4 years ago

Add more variables to assimilation background, so it includes all required elements of the background state.

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