source: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis3/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 11070

Last change on this file since 11070 was 11070, checked in by rrenshaw, 18 months ago

code fixes

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