source: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 8648

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

Change source of MLD for direct initialisation of logchl and pCO2 to work.

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