source: branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 9331

Last change on this file since 9331 was 9331, checked in by dford, 2 years ago

Add balancing code.

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