source: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 12613

Last change on this file since 12613 was 12613, checked in by dcarneir, 19 months ago

Including sea ice thickness in GO6 branch

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