source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 10622

Last change on this file since 10622 was 10622, checked in by dford, 18 months ago

Implement biogeochemistry assimilation for FABM-ERSEM.

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