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

Last change on this file since 8436 was 8436, checked in by dford, 4 years ago

Implement initial version of surface chlorophyll assimilation for MEDUSA.

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