New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
asminc.F90 in branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 5 years ago

GMED 450 add flush after prints

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