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_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 11953

Last change on this file since 11953 was 11953, checked in by dcarneir, 4 years ago

Compilation errors fixed in asminc.F90 and diaobs.F90

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