source: branches/UKMO/dev_r5518_v3.6_asm_nemovar_community/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 7803

Last change on this file since 7803 was 7803, checked in by dford, 3 years ago

Merge in changes to allow a maximum chlorophyll assimilation increment to be specified. See internal UKMO ticket 690.

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