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

source: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 9537

Last change on this file since 9537 was 9537, checked in by jcastill, 6 years ago

Bug fix for key_asminc
WARNING - this fix is applied to revision 9196 of the branch, so the fix at revision 9473 is not included - it will be included in the next revision: this is done so we can use the same version as the operational suite, which at the moment does not contain the fix included in r9473

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