source: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 8428

Last change on this file since 8428 was 8428, checked in by dford, 3 years ago

Add logchl assimilation, with nitrogen balancing for HadOCC.

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