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

source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 5168

Last change on this file since 5168 was 5168, checked in by pabouttier, 7 years ago

Fixed several bugs described in Ticket #1360

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