source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_sshiau/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 14161

Last change on this file since 14161 was 14161, checked in by kingr, 12 months ago

Allow heat/salt conservation when assimilating SSH to be namelist controlled.

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