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

Last change on this file since 10574 was 10574, checked in by dford, 19 months ago

Merge in functionality added to GO6 at r10149.

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