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

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

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

Last change on this file since 8456 was 8456, checked in by dford, 7 years ago

Add pCO2/fCO2 assimilation.

File size: 78.6 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    = 5   !: 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         ENDIF
753
754         CALL iom_close( inum )
755         
756         DO jt = 1, jptra
757            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
758         END DO
759     
760      ELSE IF ( ln_pco2inc .OR. ln_fco2inc ) THEN
761
762         ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
763         tracer_bkg(:,:,:,:) = 0.0
764
765         CALL iom_open( c_asmbkg, inum )
766         
767#if defined key_hadocc
768         CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
769         CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
770#elif defined key_medusa
771         CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
772         CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
773#endif
774
775         CALL iom_close( inum )
776         
777         DO jt = 1, jptra
778            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
779         END DO
780     
781      ENDIF
782#endif
783      !
784   END SUBROUTINE asm_inc_init
785
786
787   SUBROUTINE calc_date( kit000, kt, kdate0, kdate )
788      !!----------------------------------------------------------------------
789      !!                    ***  ROUTINE calc_date  ***
790      !!         
791      !! ** Purpose : Compute the calendar date YYYYMMDD at a given time step.
792      !!
793      !! ** Method  : Compute the calendar date YYYYMMDD at a given time step.
794      !!
795      !! ** Action  :
796      !!----------------------------------------------------------------------
797      INTEGER, INTENT(IN) :: kit000  ! Initial time step
798      INTEGER, INTENT(IN) :: kt      ! Current time step referenced to kit000
799      INTEGER, INTENT(IN) :: kdate0  ! Initial date
800      INTEGER, INTENT(OUT) :: kdate  ! Current date reference to kdate0
801      !
802      INTEGER :: iyea0    ! Initial year
803      INTEGER :: imon0    ! Initial month
804      INTEGER :: iday0    ! Initial day
805      INTEGER :: iyea     ! Current year
806      INTEGER :: imon     ! Current month
807      INTEGER :: iday     ! Current day
808      INTEGER :: idaystp  ! Number of days between initial and current date
809      INTEGER :: idaycnt  ! Day counter
810
811      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
812
813      !-----------------------------------------------------------------------
814      ! Compute the calendar date YYYYMMDD
815      !-----------------------------------------------------------------------
816
817      ! Initial date
818      iyea0 =   kdate0 / 10000
819      imon0 = ( kdate0 - ( iyea0 * 10000 ) ) / 100
820      iday0 =   kdate0 - ( iyea0 * 10000 ) - ( imon0 * 100 ) 
821
822      ! Check that kt >= kit000 - 1
823      IF ( kt < kit000 - 1 ) CALL ctl_stop( ' kt must be >= kit000 - 1')
824
825      ! If kt = kit000 - 1 then set the date to the restart date
826      IF ( kt == kit000 - 1 ) THEN
827
828         kdate = ndastp
829         RETURN
830
831      ENDIF
832
833      ! Compute the number of days from the initial date
834      idaystp = INT( REAL( kt - kit000 ) * rdt / 86400. )
835   
836      iday    = iday0
837      imon    = imon0
838      iyea    = iyea0
839      idaycnt = 0
840
841      CALL calc_month_len( iyea, imonth_len )
842
843      DO WHILE ( idaycnt < idaystp )
844         iday = iday + 1
845         IF ( iday > imonth_len(imon) )  THEN
846            iday = 1
847            imon = imon + 1
848         ENDIF
849         IF ( imon > 12 ) THEN
850            imon = 1
851            iyea = iyea + 1
852            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
853         ENDIF                 
854         idaycnt = idaycnt + 1
855      END DO
856      !
857      kdate = iyea * 10000 + imon * 100 + iday
858      !
859   END SUBROUTINE
860
861
862   SUBROUTINE calc_month_len( iyear, imonth_len )
863      !!----------------------------------------------------------------------
864      !!                    ***  ROUTINE calc_month_len  ***
865      !!         
866      !! ** Purpose : Compute the number of days in a months given a year.
867      !!
868      !! ** Method  :
869      !!----------------------------------------------------------------------
870      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
871      INTEGER :: iyear         !: year
872      !!----------------------------------------------------------------------
873      !
874      ! length of the month of the current year (from nleapy, read in namelist)
875      IF ( nleapy < 2 ) THEN
876         imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
877         IF ( nleapy == 1 ) THEN   ! we are using calendar with leap years
878            IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
879               imonth_len(2) = 29
880            ENDIF
881         ENDIF
882      ELSE
883         imonth_len(:) = nleapy   ! all months with nleapy days per year
884      ENDIF
885      !
886   END SUBROUTINE
887
888
889   SUBROUTINE tra_asm_inc( kt )
890      !!----------------------------------------------------------------------
891      !!                    ***  ROUTINE tra_asm_inc  ***
892      !!         
893      !! ** Purpose : Apply the tracer (T and S) assimilation increments
894      !!
895      !! ** Method  : Direct initialization or Incremental Analysis Updating
896      !!
897      !! ** Action  :
898      !!----------------------------------------------------------------------
899      INTEGER, INTENT(IN) :: kt               ! Current time step
900      !
901      INTEGER :: ji,jj,jk
902      INTEGER :: it
903      REAL(wp) :: zincwgt  ! IAU weight for current time step
904      REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values
905      !!----------------------------------------------------------------------
906
907      ! freezing point calculation taken from oc_fz_pt (but calculated for all depths)
908      ! used to prevent the applied increments taking the temperature below the local freezing point
909
910      DO jk = 1, jpkm1
911        CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), fsdept(:,:,jk) )
912      END DO
913
914      IF ( ln_asmiau ) THEN
915
916         !--------------------------------------------------------------------
917         ! Incremental Analysis Updating
918         !--------------------------------------------------------------------
919
920         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
921
922            it = kt - nit000 + 1
923            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
924
925            IF(lwp) THEN
926               WRITE(numout,*) 
927               WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
928               WRITE(numout,*) '~~~~~~~~~~~~'
929            ENDIF
930
931            ! Update the tracer tendencies
932            DO jk = 1, jpkm1
933               IF (ln_temnofreeze) THEN
934                  ! Do not apply negative increments if the temperature will fall below freezing
935                  WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. &
936                     &   tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 
937                     tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 
938                  END WHERE
939               ELSE
940                  tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 
941               ENDIF
942               IF (ln_salfix) THEN
943                  ! Do not apply negative increments if the salinity will fall below a specified
944                  ! minimum value salfixmin
945                  WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. &
946                     &   tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 
947                     tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
948                  END WHERE
949               ELSE
950                  tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
951               ENDIF
952            END DO
953
954         ENDIF
955
956         IF ( kt == nitiaufin_r + 1  ) THEN   ! For bias crcn to work
957            DEALLOCATE( t_bkginc )
958            DEALLOCATE( s_bkginc )
959         ENDIF
960
961
962      ELSEIF ( ln_asmdin ) THEN
963
964         !--------------------------------------------------------------------
965         ! Direct Initialization
966         !--------------------------------------------------------------------
967           
968         IF ( kt == nitdin_r ) THEN
969
970            neuler = 0  ! Force Euler forward step
971
972            ! Initialize the now fields with the background + increment
973            IF (ln_temnofreeze) THEN
974               ! Do not apply negative increments if the temperature will fall below freezing
975               WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 
976                  tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
977               END WHERE
978            ELSE
979               tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
980            ENDIF
981            IF (ln_salfix) THEN
982               ! Do not apply negative increments if the salinity will fall below a specified
983               ! minimum value salfixmin
984               WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 
985                  tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
986               END WHERE
987            ELSE
988               tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
989            ENDIF
990
991            tsb(:,:,:,:) = tsn(:,:,:,:)                 ! Update before fields
992
993            CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )  ! Before potential and in situ densities
994!!gm  fabien
995!            CALL eos( tsb, rhd, rhop )                ! Before potential and in situ densities
996!!gm
997
998
999            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)      &
1000               &  CALL zps_hde    ( kt, jpts, tsb, gtsu, gtsv,        &  ! Partial steps: before horizontal gradient
1001               &                              rhd, gru , grv          )  ! of t, s, rd at the last ocean level
1002            IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)      &
1003               &  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,    &    ! Partial steps for top cell (ISF)
1004               &                                  rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
1005               &                           gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
1006
1007#if defined key_zdfkpp
1008            CALL eos( tsn, rhd, fsdept_n(:,:,:) )                      ! Compute rhd
1009!!gm fabien            CALL eos( tsn, rhd )                      ! Compute rhd
1010#endif
1011
1012            DEALLOCATE( t_bkginc )
1013            DEALLOCATE( s_bkginc )
1014            DEALLOCATE( t_bkg    )
1015            DEALLOCATE( s_bkg    )
1016         ENDIF
1017         
1018      ENDIF
1019      ! Perhaps the following call should be in step
1020      IF   ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment
1021      !
1022   END SUBROUTINE tra_asm_inc
1023
1024
1025   SUBROUTINE dyn_asm_inc( kt )
1026      !!----------------------------------------------------------------------
1027      !!                    ***  ROUTINE dyn_asm_inc  ***
1028      !!         
1029      !! ** Purpose : Apply the dynamics (u and v) assimilation increments.
1030      !!
1031      !! ** Method  : Direct initialization or Incremental Analysis Updating.
1032      !!
1033      !! ** Action  :
1034      !!----------------------------------------------------------------------
1035      INTEGER, INTENT(IN) :: kt   ! Current time step
1036      !
1037      INTEGER :: jk
1038      INTEGER :: it
1039      REAL(wp) :: zincwgt  ! IAU weight for current time step
1040      !!----------------------------------------------------------------------
1041
1042      IF ( ln_asmiau ) THEN
1043
1044         !--------------------------------------------------------------------
1045         ! Incremental Analysis Updating
1046         !--------------------------------------------------------------------
1047
1048         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1049
1050            it = kt - nit000 + 1
1051            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
1052
1053            IF(lwp) THEN
1054               WRITE(numout,*) 
1055               WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', &
1056                  &  kt,' with IAU weight = ', wgtiau(it)
1057               WRITE(numout,*) '~~~~~~~~~~~~'
1058            ENDIF
1059
1060            ! Update the dynamic tendencies
1061            DO jk = 1, jpkm1
1062               ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt
1063               va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt
1064            END DO
1065           
1066            IF ( kt == nitiaufin_r ) THEN
1067               DEALLOCATE( u_bkginc )
1068               DEALLOCATE( v_bkginc )
1069            ENDIF
1070
1071         ENDIF
1072
1073      ELSEIF ( ln_asmdin ) THEN 
1074
1075         !--------------------------------------------------------------------
1076         ! Direct Initialization
1077         !--------------------------------------------------------------------
1078         
1079         IF ( kt == nitdin_r ) THEN
1080
1081            neuler = 0                    ! Force Euler forward step
1082
1083            ! Initialize the now fields with the background + increment
1084            un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:)
1085            vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:) 
1086
1087            ub(:,:,:) = un(:,:,:)         ! Update before fields
1088            vb(:,:,:) = vn(:,:,:)
1089 
1090            DEALLOCATE( u_bkg    )
1091            DEALLOCATE( v_bkg    )
1092            DEALLOCATE( u_bkginc )
1093            DEALLOCATE( v_bkginc )
1094         ENDIF
1095         !
1096      ENDIF
1097      !
1098   END SUBROUTINE dyn_asm_inc
1099
1100
1101   SUBROUTINE ssh_asm_inc( kt )
1102      !!----------------------------------------------------------------------
1103      !!                    ***  ROUTINE ssh_asm_inc  ***
1104      !!         
1105      !! ** Purpose : Apply the sea surface height assimilation increment.
1106      !!
1107      !! ** Method  : Direct initialization or Incremental Analysis Updating.
1108      !!
1109      !! ** Action  :
1110      !!----------------------------------------------------------------------
1111      INTEGER, INTENT(IN) :: kt   ! Current time step
1112      !
1113      INTEGER :: it
1114      INTEGER :: jk
1115      REAL(wp) :: zincwgt  ! IAU weight for current time step
1116      !!----------------------------------------------------------------------
1117
1118      IF ( ln_asmiau ) THEN
1119
1120         !--------------------------------------------------------------------
1121         ! Incremental Analysis Updating
1122         !--------------------------------------------------------------------
1123
1124         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1125
1126            it = kt - nit000 + 1
1127            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
1128
1129            IF(lwp) THEN
1130               WRITE(numout,*) 
1131               WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
1132                  &  kt,' with IAU weight = ', wgtiau(it)
1133               WRITE(numout,*) '~~~~~~~~~~~~'
1134            ENDIF
1135
1136            ! Save the tendency associated with the IAU weighted SSH increment
1137            ! (applied in dynspg.*)
1138#if defined key_asminc
1139            ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
1140#endif
1141            IF ( kt == nitiaufin_r ) THEN
1142               DEALLOCATE( ssh_bkginc )
1143            ENDIF
1144
1145         ELSE
1146#if defined key_asminc
1147            ssh_iau(:,:) = 0.0
1148#endif
1149         ENDIF
1150
1151      ELSEIF ( ln_asmdin ) THEN
1152
1153         !--------------------------------------------------------------------
1154         ! Direct Initialization
1155         !--------------------------------------------------------------------
1156
1157         IF ( kt == nitdin_r ) THEN
1158
1159            neuler = 0                    ! Force Euler forward step
1160
1161            ! Initialize the now fields the background + increment
1162            sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) 
1163
1164            ! Update before fields
1165            sshb(:,:) = sshn(:,:)         
1166
1167            IF( lk_vvl ) THEN
1168               DO jk = 1, jpk
1169                  fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
1170               END DO
1171            ENDIF
1172
1173            DEALLOCATE( ssh_bkg    )
1174            DEALLOCATE( ssh_bkginc )
1175
1176         ENDIF
1177         !
1178      ENDIF
1179      !
1180   END SUBROUTINE ssh_asm_inc
1181
1182
1183   SUBROUTINE seaice_asm_inc( kt, kindic )
1184      !!----------------------------------------------------------------------
1185      !!                    ***  ROUTINE seaice_asm_inc  ***
1186      !!         
1187      !! ** Purpose : Apply the sea ice assimilation increment.
1188      !!
1189      !! ** Method  : Direct initialization or Incremental Analysis Updating.
1190      !!
1191      !! ** Action  :
1192      !!
1193      !!----------------------------------------------------------------------
1194      IMPLICIT NONE
1195      !
1196      INTEGER, INTENT(in)           ::   kt   ! Current time step
1197      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation
1198      !
1199      INTEGER  ::   it
1200      REAL(wp) ::   zincwgt   ! IAU weight for current time step
1201#if defined key_lim2
1202      REAL(wp), DIMENSION(jpi,jpj) ::   zofrld, zohicif, zseaicendg, zhicifinc  ! LIM
1203      REAL(wp) ::   zhicifmin = 0.5_wp      ! ice minimum depth in metres
1204#endif
1205      !!----------------------------------------------------------------------
1206
1207      IF ( ln_asmiau ) THEN
1208
1209         !--------------------------------------------------------------------
1210         ! Incremental Analysis Updating
1211         !--------------------------------------------------------------------
1212
1213         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1214
1215            it = kt - nit000 + 1
1216            zincwgt = wgtiau(it)      ! IAU weight for the current time step
1217            ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments)
1218
1219            IF(lwp) THEN
1220               WRITE(numout,*) 
1221               WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', &
1222                  &  kt,' with IAU weight = ', wgtiau(it)
1223               WRITE(numout,*) '~~~~~~~~~~~~'
1224            ENDIF
1225
1226            ! Sea-ice : LIM-3 case (to add)
1227
1228#if defined key_lim2
1229            ! Sea-ice : LIM-2 case
1230            zofrld (:,:) = frld(:,:)
1231            zohicif(:,:) = hicif(:,:)
1232            !
1233            frld  = MIN( MAX( frld (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
1234            pfrld = MIN( MAX( pfrld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
1235            fr_i(:,:) = 1.0_wp - frld(:,:)        ! adjust ice fraction
1236            !
1237            zseaicendg(:,:) = zofrld(:,:) - frld(:,:)   ! find out actual sea ice nudge applied
1238            !
1239            ! Nudge sea ice depth to bring it up to a required minimum depth
1240            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 
1241               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt   
1242            ELSEWHERE
1243               zhicifinc(:,:) = 0.0_wp
1244            END WHERE
1245            !
1246            ! nudge ice depth
1247            hicif (:,:) = hicif (:,:) + zhicifinc(:,:)
1248            phicif(:,:) = phicif(:,:) + zhicifinc(:,:)       
1249            !
1250            ! seaice salinity balancing (to add)
1251#endif
1252
1253#if defined key_cice && defined key_asminc
1254            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
1255            ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt
1256#endif
1257
1258            IF ( kt == nitiaufin_r ) THEN
1259               DEALLOCATE( seaice_bkginc )
1260            ENDIF
1261
1262         ELSE
1263
1264#if defined key_cice && defined key_asminc
1265            ! Sea-ice : CICE case. Zero ice increment tendency into CICE
1266            ndaice_da(:,:) = 0.0_wp
1267#endif
1268
1269         ENDIF
1270
1271      ELSEIF ( ln_asmdin ) THEN
1272
1273         !--------------------------------------------------------------------
1274         ! Direct Initialization
1275         !--------------------------------------------------------------------
1276
1277         IF ( kt == nitdin_r ) THEN
1278
1279            neuler = 0                    ! Force Euler forward step
1280
1281            ! Sea-ice : LIM-3 case (to add)
1282
1283#if defined key_lim2
1284            ! Sea-ice : LIM-2 case.
1285            zofrld(:,:)=frld(:,:)
1286            zohicif(:,:)=hicif(:,:)
1287            !
1288            ! Initialize the now fields the background + increment
1289            frld (:,:) = MIN( MAX( frld(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp)
1290            pfrld(:,:) = frld(:,:) 
1291            fr_i (:,:) = 1.0_wp - frld(:,:)                ! adjust ice fraction
1292            zseaicendg(:,:) = zofrld(:,:) - frld(:,:)      ! find out actual sea ice nudge applied
1293            !
1294            ! Nudge sea ice depth to bring it up to a required minimum depth
1295            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 
1296               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt   
1297            ELSEWHERE
1298               zhicifinc(:,:) = 0.0_wp
1299            END WHERE
1300            !
1301            ! nudge ice depth
1302            hicif (:,:) = hicif (:,:) + zhicifinc(:,:)
1303            phicif(:,:) = phicif(:,:)       
1304            !
1305            ! seaice salinity balancing (to add)
1306#endif
1307 
1308#if defined key_cice && defined key_asminc
1309            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
1310           ndaice_da(:,:) = seaice_bkginc(:,:) / rdt
1311#endif
1312           IF ( .NOT. PRESENT(kindic) ) THEN
1313              DEALLOCATE( seaice_bkginc )
1314           END IF
1315
1316         ELSE
1317
1318#if defined key_cice && defined key_asminc
1319            ! Sea-ice : CICE case. Zero ice increment tendency into CICE
1320            ndaice_da(:,:) = 0.0_wp
1321#endif
1322         
1323         ENDIF
1324
1325!#if defined defined key_lim2 || defined key_cice
1326!
1327!            IF (ln_seaicebal ) THEN       
1328!             !! balancing salinity increments
1329!             !! simple case from limflx.F90 (doesn't include a mass flux)
1330!             !! assumption is that as ice concentration is reduced or increased
1331!             !! the snow and ice depths remain constant
1332!             !! note that snow is being created where ice concentration is being increased
1333!             !! - could be more sophisticated and
1334!             !! not do this (but would need to alter h_snow)
1335!
1336!             usave(:,:,:)=sb(:,:,:)   ! use array as a temporary store
1337!
1338!             DO jj = 1, jpj
1339!               DO ji = 1, jpi
1340!           ! calculate change in ice and snow mass per unit area
1341!           ! positive values imply adding salt to the ocean (results from ice formation)
1342!           ! fwf : ice formation and melting
1343!
1344!                 zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rdt
1345!
1346!           ! change salinity down to mixed layer depth
1347!                 mld=hmld_kara(ji,jj)
1348!
1349!           ! prevent small mld
1350!           ! less than 10m can cause salinity instability
1351!                 IF (mld < 10) mld=10
1352!
1353!           ! set to bottom of a level
1354!                 DO jk = jpk-1, 2, -1
1355!                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN
1356!                     mld=gdepw(ji,jj,jk+1)
1357!                     jkmax=jk
1358!                   ENDIF
1359!                 ENDDO
1360!
1361!            ! avoid applying salinity balancing in shallow water or on land
1362!            !
1363!
1364!            ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m)
1365!
1366!                 dsal_ocn=0.0_wp
1367!                 sal_thresh=5.0_wp        ! minimum salinity threshold for salinity balancing
1368!
1369!                 if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) &
1370!                              dsal_ocn = zfons / (rhop(ji,jj,1) * mld)
1371!
1372!           ! put increments in for levels in the mixed layer
1373!           ! but prevent salinity below a threshold value
1374!
1375!                   DO jk = 1, jkmax             
1376!
1377!                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN
1378!                           sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn
1379!                           sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn
1380!                     ENDIF
1381!
1382!                   ENDDO
1383!
1384!      !            !  salt exchanges at the ice/ocean interface
1385!      !            zpmess         = zfons / rdt_ice    ! rdt_ice is ice timestep
1386!      !
1387!      !! Adjust fsalt. A +ve fsalt means adding salt to ocean
1388!      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt 
1389!      !!               
1390!      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d)
1391!      !!                                                     ! E-P (kg m-2 s-2)
1392!      !            emp(ji,jj) = emp(ji,jj) + zpmess          ! E-P (kg m-2 s-2)
1393!               ENDDO !ji
1394!             ENDDO !jj!
1395!
1396!            ENDIF !ln_seaicebal
1397!
1398!#endif
1399
1400      ENDIF
1401
1402   END SUBROUTINE seaice_asm_inc
1403
1404   SUBROUTINE logchl_asm_inc( kt )
1405      !!----------------------------------------------------------------------
1406      !!                    ***  ROUTINE logchl_asm_inc  ***
1407      !!         
1408      !! ** Purpose : Apply the chlorophyll assimilation increments.
1409      !!
1410      !! ** Method  : Calculate increments to state variables using nitrogen
1411      !!              balancing.
1412      !!              Direct initialization or Incremental Analysis Updating.
1413      !!
1414      !! ** Action  :
1415      !!----------------------------------------------------------------------
1416      INTEGER, INTENT(IN) :: kt   ! Current time step
1417      !
1418      INTEGER  :: jk              ! Loop counter
1419      INTEGER  :: it              ! Index
1420      REAL(wp) :: zincwgt         ! IAU weight for current time step
1421      REAL(wp) :: zincper         ! IAU interval in seconds
1422      !!----------------------------------------------------------------------
1423
1424      IF ( kt <= nit000 ) THEN
1425
1426         zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
1427
1428#if defined key_medusa && defined key_foam_medusa
1429         CALL asm_logchl_bal_medusa( logchl_bkginc, zincper, mld_choice_bgc, &
1430            &                        rn_maxchlinc, ln_logchlbal,             &
1431            &                        pgrow_avg_bkg, ploss_avg_bkg,           &
1432            &                        phyt_avg_bkg, mld_max_bkg,              &
1433            &                        tracer_bkg, logchl_balinc )
1434#elif defined key_hadocc
1435         CALL asm_logchl_bal_hadocc( logchl_bkginc, zincper, mld_choice_bgc, &
1436            &                        rn_maxchlinc, ln_logchlbal,             &
1437            &                        pgrow_avg_bkg, ploss_avg_bkg,           &
1438            &                        phyt_avg_bkg, mld_max_bkg,              &
1439            &                        chl_bkg, cchl_p_bkg,                    &
1440            &                        tracer_bkg, logchl_balinc )
1441#else
1442         CALL ctl_stop( 'Attempting to assimilate logchl, ', &
1443            &           'but not defined a biogeochemical model' )
1444#endif
1445
1446      ENDIF
1447
1448      IF ( ln_asmiau ) THEN
1449
1450         !--------------------------------------------------------------------
1451         ! Incremental Analysis Updating
1452         !--------------------------------------------------------------------
1453
1454         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1455
1456            it = kt - nit000 + 1
1457            zincwgt = wgtiau(it)      ! IAU weight for the current time step
1458            ! note this is not a tendency so should not be divided by rdt
1459
1460            IF(lwp) THEN
1461               WRITE(numout,*) 
1462               WRITE(numout,*) 'logchl_asm_inc : logchl IAU at time step = ', &
1463                  &  kt,' with IAU weight = ', wgtiau(it)
1464               WRITE(numout,*) '~~~~~~~~~~~~'
1465            ENDIF
1466
1467            ! Update the biogeochemical variables
1468            ! Add directly to trn and trb, rather than to tra, as not a tendency
1469#if defined key_medusa && defined key_foam_medusa
1470            DO jk = 1, jpkm1
1471               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
1472                  &                          logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
1473               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
1474                  &                          logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
1475            END DO
1476#elif defined key_hadocc
1477            DO jk = 1, jpkm1
1478               trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
1479                  &                          logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
1480               trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
1481                  &                          logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
1482            END DO
1483#endif
1484
1485            ! Do not deallocate arrays - needed by asm_bal_wri
1486            ! which is called at end of model run
1487         ENDIF
1488
1489      ELSEIF ( ln_asmdin ) THEN 
1490
1491         !--------------------------------------------------------------------
1492         ! Direct Initialization
1493         !--------------------------------------------------------------------
1494         
1495         IF ( kt == nitdin_r ) THEN
1496
1497            neuler = 0                    ! Force Euler forward step
1498
1499#if defined key_medusa && defined key_foam_medusa
1500            ! Initialize the now fields with the background + increment
1501            ! Background currently is what the model is initialised with
1502            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with chlorophyll assimilation', &
1503               &           ' Background state is taken from model rather than background file' )
1504            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1505               &                         logchl_balinc(:,:,:,jp_msa0:jp_msa1)
1506            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1507#elif defined key_hadocc
1508            ! Initialize the now fields with the background + increment
1509            ! Background currently is what the model is initialised with
1510            CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', &
1511               &           ' Background state is taken from model rather than background file' )
1512            trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1513               &                         logchl_balinc(:,:,:,jp_had0:jp_had1)
1514            trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1515#endif
1516 
1517            ! Do not deallocate arrays - needed by asm_bal_wri
1518            ! which is called at end of model run
1519         ENDIF
1520         !
1521      ENDIF
1522      !
1523   END SUBROUTINE logchl_asm_inc
1524
1525
1526   SUBROUTINE pco2_asm_inc( kt )
1527      !!----------------------------------------------------------------------
1528      !!                    ***  ROUTINE pco2_asm_inc  ***
1529      !!         
1530      !! ** Purpose : Apply the pco2/fco2 assimilation increments.
1531      !!
1532      !! ** Method  : Calculate increments to state variables using carbon
1533      !!              balancing.
1534      !!              Direct initialization or Incremental Analysis Updating.
1535      !!
1536      !! ** Action  :
1537      !!----------------------------------------------------------------------
1538      INTEGER, INTENT(IN)                   :: kt           ! Current time step
1539      !
1540      INTEGER                               :: ji, jj, jk   ! Loop counters
1541      INTEGER                               :: it           ! Index
1542      INTEGER                               :: jkmax        ! Index of mixed layer depth
1543      REAL(wp)                              :: zincwgt      ! IAU weight for current time step
1544      REAL(wp), DIMENSION(jpi,jpj)          :: zmld         ! Mixed layer depth for balancing
1545      REAL(wp), DIMENSION(jpi,jpj)          :: pco2_bkginc_temp ! For fCO2 conversion
1546      REAL(wp), DIMENSION(jpi,jpj)          :: dic_bkg_temp ! DIC background
1547      REAL(wp), DIMENSION(jpi,jpj)          :: alk_bkg_temp ! Alkalinity background
1548      REAL(wp), DIMENSION(jpi,jpj)          :: tem_bkg_temp ! Temperature background
1549      REAL(wp), DIMENSION(jpi,jpj)          :: sal_bkg_temp ! Salinity background
1550      REAL(wp), DIMENSION(1)                :: patm         ! Atmospheric pressure
1551      REAL(wp), DIMENSION(1)                :: phyd         ! Hydrostatic pressure
1552
1553      ! Coefficients for fCO2 to pCO2 conversion
1554      REAL(wp)                              :: zcoef_fco2_1 = -1636.75
1555      REAL(wp)                              :: zcoef_fco2_2 = 12.0408
1556      REAL(wp)                              :: zcoef_fco2_3 = 0.0327957
1557      REAL(wp)                              :: zcoef_fco2_4 = 0.0000316528
1558      REAL(wp)                              :: zcoef_fco2_5 = 2.0
1559      REAL(wp)                              :: zcoef_fco2_6 = 57.7
1560      REAL(wp)                              :: zcoef_fco2_7 = 0.118
1561      REAL(wp)                              :: zcoef_fco2_8 = 82.0578
1562      !!----------------------------------------------------------------------
1563
1564      IF ( kt <= nit000 ) THEN
1565
1566         !--------------------------------------------------------------------
1567         ! DIC and alkalinity balancing
1568         !--------------------------------------------------------------------
1569
1570         IF ( ln_fco2inc ) THEN
1571#if defined key_medusa && defined key_roam
1572            ! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine
1573            patm(1) = 1.0
1574            phyd(1) = 0.0
1575            DO jj = 1, jpj
1576               DO ji = 1, jpi
1577                  CALL f2pCO2( pco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) )
1578               END DO
1579            END DO
1580#else
1581            ! If assimilating fCO2, then convert to pCO2 using temperature
1582            ! See flux_gas.F90 within HadOCC for details of calculation
1583            pco2_bkginc_temp(:,:) = pco2_bkginc(:,:) /                                                             &
1584               &                    EXP((zcoef_fco2_1                                                            + &
1585               &                         zcoef_fco2_2 * (tsn(:,:,1,1)+rt0)                                       - &
1586               &                         zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)                    + &
1587               &                         zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + &
1588               &                         zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0)))      / &
1589               &                        (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0)))
1590#endif
1591         ELSE
1592            pco2_bkginc_temp(:,:) = pco2_bkginc(:,:)
1593         ENDIF
1594
1595         ! Account for physics assimilation if required
1596         IF ( ln_trainc ) THEN
1597            tem_bkg_temp(:,:) = tsn(:,:,1,1) + t_bkginc(:,:,1)
1598            sal_bkg_temp(:,:) = tsn(:,:,1,2) + s_bkginc(:,:,1)
1599         ELSE
1600            tem_bkg_temp(:,:) = tsn(:,:,1,1)
1601            sal_bkg_temp(:,:) = tsn(:,:,1,2)
1602         ENDIF
1603
1604#if defined key_medusa
1605         ! Account for logchl balancing if required
1606         IF ( ln_logchlinc .AND. ln_logchlbal ) THEN
1607            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + logchl_balinc(:,:,1,jpdic)
1608            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + logchl_balinc(:,:,1,jpalk)
1609         ELSE
1610            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
1611            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk)
1612         ENDIF
1613
1614         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1615            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
1616            &               pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) )
1617
1618#elif defined key_hadocc
1619         ! Account for logchl balancing if required
1620         IF ( ln_logchlinc .AND. ln_logchlbal ) THEN
1621            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + logchl_balinc(:,:,1,jp_had_dic)
1622            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + logchl_balinc(:,:,1,jp_had_alk)
1623         ELSE
1624            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
1625            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk)
1626         ENDIF
1627
1628         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1629            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
1630            &               pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) )
1631
1632#else
1633         CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
1634            &           'but not defined a biogeochemical model' )
1635#endif
1636
1637         ! Select mixed layer
1638         SELECT CASE( mld_choice_bgc )
1639         CASE ( 1 )                   ! Turbocline/mixing depth [W points]
1640            zmld(:,:) = hmld(:,:)
1641         CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
1642            zmld(:,:) = hmlp(:,:)
1643         CASE ( 3 )                   ! Kara MLD [Interpolated]
1644#if defined key_karaml
1645            IF ( ln_kara ) THEN
1646               zmld(:,:) = hmld_kara(:,:)
1647            ELSE
1648               CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
1649                  &           ' but ln_kara=.false.' )
1650            ENDIF
1651#else
1652            CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
1653               &           ' but is not defined' )
1654#endif
1655         CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
1656            !zmld(:,:) = hmld_tref(:,:)
1657            CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', &
1658               &           ' but is not available in this version' )
1659         CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
1660            zmld(:,:) = hmlpt(:,:)
1661         END SELECT
1662
1663         ! Propagate through mixed layer
1664         DO jj = 1, jpj
1665            DO ji = 1, jpi
1666               !
1667               jkmax = jpk-1
1668               DO jk = jpk-1, 1, -1
1669                  IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
1670                     & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
1671                     zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
1672                     jkmax = jk
1673                  ENDIF
1674               END DO
1675               !
1676               DO jk = 2, jkmax
1677                  pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:)
1678               END DO
1679               !
1680            END DO
1681         END DO
1682
1683      ENDIF
1684
1685      IF ( ln_asmiau ) THEN
1686
1687         !--------------------------------------------------------------------
1688         ! Incremental Analysis Updating
1689         !--------------------------------------------------------------------
1690
1691         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1692
1693            it = kt - nit000 + 1
1694            zincwgt = wgtiau(it)      ! IAU weight for the current time step
1695            ! note this is not a tendency so should not be divided by rdt
1696
1697            IF(lwp) THEN
1698               WRITE(numout,*) 
1699               IF ( ln_pco2inc ) THEN
1700                  WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
1701                     &  kt,' with IAU weight = ', wgtiau(it)
1702               ELSE IF ( ln_fco2inc ) THEN
1703                  WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', &
1704                     &  kt,' with IAU weight = ', wgtiau(it)
1705               ENDIF
1706               WRITE(numout,*) '~~~~~~~~~~~~'
1707            ENDIF
1708
1709            ! Update the biogeochemical variables
1710            ! Add directly to trn and trb, rather than to tra, as not a tendency
1711#if defined key_medusa && defined key_foam_medusa
1712            DO jk = 1, jpkm1
1713               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
1714                  &                          pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
1715               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
1716                  &                          pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
1717            END DO
1718#elif defined key_hadocc
1719            DO jk = 1, jpkm1
1720               trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
1721                  &                          pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
1722               trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
1723                  &                          pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
1724            END DO
1725#endif
1726
1727            ! Do not deallocate arrays - needed by asm_bal_wri
1728            ! which is called at end of model run
1729
1730         ENDIF
1731
1732      ELSEIF ( ln_asmdin ) THEN 
1733
1734         !--------------------------------------------------------------------
1735         ! Direct Initialization
1736         !--------------------------------------------------------------------
1737         
1738         IF ( kt == nitdin_r ) THEN
1739
1740            neuler = 0                    ! Force Euler forward step
1741
1742#if defined key_medusa && defined key_foam_medusa
1743            ! Initialize the now fields with the background + increment
1744            ! Background currently is what the model is initialised with
1745            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pCO2 assimilation', &
1746               &           ' Background state is taken from model rather than background file' )
1747            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1748               &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1)
1749            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1750#elif defined key_hadocc
1751            ! Initialize the now fields with the background + increment
1752            ! Background currently is what the model is initialised with
1753            CALL ctl_warn( ' Doing direct initialisation of HadOCC with pCO2 assimilation', &
1754               &           ' Background state is taken from model rather than background file' )
1755            trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1756               &                         pco2_balinc(:,:,:,jp_had0:jp_had1)
1757            trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1758#endif
1759 
1760            ! Do not deallocate arrays - needed by asm_bal_wri
1761            ! which is called at end of model run
1762         ENDIF
1763         !
1764      ENDIF
1765      !
1766   END SUBROUTINE pco2_asm_inc
1767   
1768   !!======================================================================
1769END MODULE asminc
Note: See TracBrowser for help on using the repository browser.