New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
asminc.F90 in branches/UKMO/dev_r5518_assim_inertial_osci/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/dev_r5518_assim_inertial_osci/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 15804

Last change on this file since 15804 was 15804, checked in by jenniewaters, 2 years ago

initial testing of inertial oscillation assimilation.

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