New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
asminc.F90 in branches/UKMO/dev_r5518_v3.4_asm_nemovar_community_bgc_ersem/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/dev_r5518_v3.4_asm_nemovar_community_bgc_ersem/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 7716

Last change on this file since 7716 was 7716, checked in by dford, 6 years ago

Tidy up and simplify use of arrays for logchl assimilation balancing.

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