S3.1-d Conservative Leap-Frog + Asselin filter

Last edited Timestamp?

Purpose:
Introduce in the trunk the modified Leap-Frog + Robert-Asselin filter time stepping scheme (thereafter LF-RA) described in (Leclair and Madec Ocean Modelling 2009, in press). The scheme allows a exact conservation of the heat and salt contents when using variable volume thickness (key_vvl). The non-conservation of the standard LF-RA due to the time diffusion of the forcing through the filter. In the modified scheme the exact conservation is achieved by removing the forcing from the filter. This requires 4 modifications regarding the forcing evaluation, its time discretisation, the filter, and the general flowchart :

  • Flow chart : general structure must be changed in a similar way to what have been done for vvl, see 2009WP/2009Stream3/VVL
  • Forcing evaluation : the momentum and tracers forcing are evaluated at kt+½, no-more at kt.
  • Forcing time discretisation : 0.5 * ( Q(k+½) + Q(kt-½) ), instead of simply Q(kt).
  • Filter : the RA filter no-more act on the forcing, i.e. -afp * rdt * ( Q(kt+½) - Q(kt-½) ) is removed from the filter

NB: While introducing the updated time stepping strategy we found that the current version of enhanced vertical diffusion (evd) was a source small divergence of two consecutive time step. To remove this source, the static instability should be tested at both before and now time step. See ticket #401.

New versus old time discretisation of the forcing
Illustration of forcing integration methods. Old formulation (top) and new formulation (bottom).

Status:
ongoing work (NEMO-Paris).

Module involved: starting from the trunk after release v3.1, including vvl update (June 14th 2009)
step.F90 =⇒ the version introduced for vvl after revision v3.1 should be used.
fldread.F90
sbcmod.F90 ; trasbc.F90 , traqsr.F90 , dynzdf_exp.F90 , dynzdf_imp.F90
sshnxt.F90 , tranxt.F90 , dynnxt.F90 =⇒ the versions introduced for vvl after revision v3.1 should be used.
Header associated with the modification (to be included in all modified modules) :

   !!            3.2  ! 2009-06  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps

Pending issues:

  • The associated changes in TOP are still unclear. A priori, in PISCES p4zopt.F90 should be modified. Nevertheless, the TOP time stepping strategy seems to be a forward in time scheme. If yes, then the evaluation of the forcing at kt+½ is an improvement (if not to say a bug fix!). If not, action is required…
  • Check the consistency with the change planned in the emps field (switch to fsalt, a true salt flux).
  • modification required in limsbc, limsbc_2 and CICE interface : conversion of emps into a salt flux (new name 'fsalt' ?)
  • the modifications in trasbc interfere with the new runoff scheme (see 2009WP/2009Stream3)

0. Flow Chart

The flow chart of step used must be the one introduce in the trunk after the v3.1, to restore the key_vvl option. Its main characteristic is that the ssh is time-stepped first, followed immediately by the computation of the now vertical velocity.

Flow chart of step with the new Leap-Frog + Robert-Asselin filter
Step routine flow chart, form Leclair and Madec (OM 2009).


1. Forcing evaluation

In order to evaluate the forcing at kt+½, all the input data in the surface module should be interpolated in time at kt+½. This have been introduced in fldread.F90 module since release 3.0. There is no required additional action.


2. Forcing : time discretisation

Each time a forcing term appear in the 3 time evolution equation solved in NEMO/OPA (ssh, momentum, tracer), it is the mean of kt+½ and kt+½ which is used. Therefore, the surface module (including the sea-ice!) must now provide both kt+½ and kt-½ forcing fields, and all the module involving a forcing term have to be modified.

surface module (SBC)

- sbc_oce.F90 ====

add before field (denoted by a '_b' sufix) for all ocean forcing (utau_b, vtau_b, qns_b, qsr_b, emp_b, emps_b)

the following lines

   !!----------------------------------------------------------------------
   !!              Ocean Surface Boundary Condition fields
   !!----------------------------------------------------------------------
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   utau      !: sea surface i-stress (ocean referential)     [N/m2]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   vtau      !: sea surface j-stress (ocean referential)     [N/m2]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   wndm      !: wind speed module at T-point (=|U10m-Uoce|)  [m/s]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr       !: sea heat flux:     solar                     [W/m2]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns       !: sea heat flux: non solar                     [W/m2]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr_tot   !: total     solar heat flux (over sea and ice) [W/m2]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns_tot   !: total non solar heat flux (over sea and ice) [W/m2]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp       !: freshwater budget: volume flux               [Kg/m2/s]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emps      !: freshwater budget: concentration/dillution   [Kg/m2/s]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp_tot   !: total evaporation - (liquid + solid) precpitation over oce and ice
   ...

becomes

   !!----------------------------------------------------------------------
   !!              Ocean Surface Boundary Condition fields
   !!----------------------------------------------------------------------
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   utau   , utau_b   !: sea surface i-stress (ocean referential)     [N/m2]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   vtau   , vtau_b   !: sea surface j-stress (ocean referential)     [N/m2]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   wndm              !: wind speed module at T-point (=|U10m-Uoce|)  [m/s]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr    , qsr_b    !: sea heat flux:     solar                     [W/m2]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns    , qns_b    !: sea heat flux: non solar                     [W/m2]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr_tot           !: total     solar heat flux (over sea and ice) [W/m2]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns_tot           !: total non solar heat flux (over sea and ice) [W/m2]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp    , emp_b    !: freshwater budget: volume flux               [Kg/m2/s]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emps   , emps_b   !: freshwater budget: concentration/dillution   [Kg/m2/s]
   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp_tot           !: total E-P-R over ocean and ice               [Kg/m2/s]
   ...

- sbcmod.F90

add the swap of surface fields and the read/write in the restart file

Beginning of sbs routine, introduce the swap, but only for kt ≥ nit000+1. If no restart, the '_b' fields is known only at the end of the sbc calculation. Therefore we have choosen to regroup the setting of '_b' fields at nit000 at the end of sbc routine.

CAUTION : here the swap must be done at each kt, what ever the nn_fsbc value is. Otherwise the ocean forcing will be wrong.

A the very begining of sbc following lines

      IF( kt == nit000 )   CALL sbc_init         ! Read namsbc namelist : surface module

becomes

      !                                                ! ---------------------------------------- !
      !                                                !          Initialisation / Swap           !
      !                                                ! ---------------------------------------- !
      IF( kt == nit000 ) THEN         
         CALL sbc_init                                 ! Read namsbc namelist : surface module
      ELSE                             
         utau_b(:,:) = utau(:,:)                       ! Swap the ocean forcing fields
         utau_b(:,:) = utau(:,:)
         qns_b (:,:) = qns (:,:)
         qsr_b (:,:) = qsr (:,:)
         emp_b (:,:) = emp (:,:)
         emps_b(:,:) = emps(:,:)
      ENDIF

A the very end of sbc, we add the setting of the forcing field at nit000-1 (either from nit000 forcing or read in restart file) and the write in the restart file at nitend. This corresponds to the addition of the following lines

      !                                                ! ---------------------------------------- !
      IF( kt == nit000 ) THEN                          !    set before forcing field at nit000    !
         !                                             ! ---------------------------------------- !
         IF( ln_rstart .AND.    &                                !* Restart: read in restart file
            & iom_varid( numror, 'utau_b', ldstop = .FALSE. ) > 0 ) THEN 
            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file'
            CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b )   ! before i-stress  (U-point)
            CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b )   ! before j-stress  (V-point)
            CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b  )   ! before non solar heat flux (T-point)
            CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b  )   ! before     solar heat flux (T-point)
            CALL iom_get( numror, jpdom_autoglo, 'emp_b' , emp_b  )   ! before     freshwater flux (T-point)
            CALL iom_get( numror, jpdom_autoglo, 'emps_b', emp_b  )   ! before C/D freshwater flux (T-point)
            !
         ELSE                                                    !* no restart: set from nit000 values
            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields set to nit000'
            utau_b(:,:) = utau(:,:) 
            utau_b(:,:) = utau(:,:)
            qns_b (:,:) = qns (:,:)
            qsr_b (:,:) = qsr (:,:)
            emp_b (:,:) = emp (:,:)
            emps_b(:,:) = emps(:,:)
         ENDIF
      ENDIF

      !                                                ! ---------------------------------------- !
      IF( lrst_oce ) THEN                              !      Write in the ocean restart file     !
         !                                             ! ---------------------------------------- !
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ',   &
            &                    'at it= ', kt,' date= ', ndastp
         IF(lwp) WRITE(numout,*) '~~~~'
         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau )    ! 
         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , vtau )
         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns  )
         CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  )
         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  )
         CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emp  )
         !
      ENDIF


tracer forcing term (TRA)

- trasbc.F90

applied heat and freshwater/salt fluxes (qns, emp, emps) as the mean of 2 consecutive time-steps (kt+½ and kt-½)

The following lines

      ! Concentration dillution effect on (t,s)
      DO jj = 2, jpj
         DO ji = fs_2, fs_jpim1   ! vector opt.
#if ! defined key_zco
            zse3t = 1. / fse3t(ji,jj,1)
#endif
            IF( lk_vvl) THEN
               zta = ro0cpr * qns(ji,jj) * zse3t &                   ! temperature : heat flux
                &    - emp(ji,jj) * zsrau * tn(ji,jj,1)  * zse3t     ! & cooling/heating effet of EMP flux
               zsa = 0.e0                                            ! No salinity concent./dilut. effect
            ELSE
               zta = ro0cpr * qns(ji,jj) * zse3t     ! temperature : heat flux
               zsa = emps(ji,jj) * zsrau * sn(ji,jj,1)   * zse3t     ! salinity :  concent./dilut. effect
            ENDIF
            ta(ji,jj,1) = ta(ji,jj,1) + zta                          ! add the trend to the general tracer trend
            sa(ji,jj,1) = sa(ji,jj,1) + zsa
         END DO
      END DO

have been modified to : (1) incorporate the time mean of the forcing ; (2) put the IF-ELSE-ENDIF outside the DO-loop for vector optimisation ; (3) remove the key_zco optimisation (increase readability) as zco is less and less used ; (4) correct a bug in vvl case, where salt flux associated with freezing/melting were not taken into account. The resulting code the becomes

      !                             ! ---------------------- !
      IF( lk_vvl ) THEN             !  Variable Volume case  !
         !                          ! ---------------------- !
         DO jj = 2, jpj
            DO ji = fs_2, fs_jpim1   ! vector opt.
               zse3t = 0.5 / fse3t(ji,jj,1)
               !                           ! temperature: heat flux + cooling/heating effet of EMP flux
               ta(ji,jj,1) = ta(ji,jj,1) + (  ro0cpr * ( qns(ji,jj) - qns_b(ji,jj) )                & 
                  &                         - zsrau  * ( emp(ji,jj) - emp_b(ji,jj) ) * tn(ji,jj,1)  ) * zse3t
               !                           ! salinity: salt flux
               sa(ji,jj,1) = sa(ji,jj,1) + ( emps(ji,jj) + emps_b(ji,jj) ) * zse3t

!!gm BUG : in key_vvl emps must be modified to only include the salt flux due to with sea-ice freezing/melting
!!gm       otherwise this flux will be missing  ==> modification required in limsbc,  limsbc_2 and CICE interface.

            END DO
         END DO
         !                          ! ---------------------- !
      ELSE                          !  Constant Volume case  !
         !                          ! ---------------------- !
         DO jj = 2, jpj
            DO ji = fs_2, fs_jpim1   ! vector opt.
               zse3t = 0.5 / fse3t(ji,jj,1)
               !                           ! temperature: heat flux 
               ta(ji,jj,1) = ta(ji,jj,1) + ro0cpr * ( qns (ji,jj) + qns_b (ji,jj) )                * zse3t 
               !                           ! salinity: salt flux + concent./dilut. effect (both in emps)
               sa(ji,jj,1) = sa(ji,jj,1) + zsrau  * ( emps(ji,jj) + emps_b(ji,jj) )  * sn(ji,jj,1) * zse3t
            END DO
         END DO
         !
      ENDIF

CAUTION 1 : bug in the salt flux in key_vvl ==⇒ modif in limsbc_2.F90 (emps becomes a salt flux in vvl,limsbc, limsbc_2 and CICE interface must be modified)

CAUTION 2 : bug in with no solar light penetration (ln_traqsr=.false.) =⇒ modif in trasbc : qns = qns + qsr. Therefore in trasbc, the following lines:

     IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration

becomes

      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
         qsr(:,:) = 0.e0                     ! qsr set to zero
         IF( kt == nit000 ) THEN             ! idem on before field at nit000
            qns_b(:,:) = qns_b(:,:) + qsr_b(:,:)  
            qsr_b(:,:) = 0.e0                     
         ENDIF
      ENDIF

- traqsr.F90

applied the solar heat flux (qsr) as the mean of 2 consecutive time-steps (kt+½ and kt-½)

            IF( nn_chldta ==1 ) THEN                             !*  Variable Chlorophyll
...
               zsi0r = 1.e0 / rn_si0
               zcoef  = ( 1. - rn_abs ) / 3.e0                        ! equi-partition in R-G-B
               ze0(:,:,1) = rn_abs  * qsr(:,:)
               ze1(:,:,1) = zcoef * qsr(:,:)
               ze2(:,:,1) = zcoef * qsr(:,:)
               ze3(:,:,1) = zcoef * qsr(:,:)
               zea(:,:,1) =         qsr(:,:)
...
                  ta(:,:,jk) = ta(:,:,jk) + etot3(:,:,jk) * qsr(:,:)

         ENDIF
         !                                                ! ------------------------- !
         IF( ln_qsr_2bd ) THEN                            !  2 band light penetration !
            !                                             ! ------------------------- !
            !
            DO jk = 1, nksr
               DO jj = 2, jpjm1
                  DO ji = fs_2, fs_jpim1   ! vector opt.
                     ta(ji,jj,jk) = ta(ji,jj,jk) + etot3(ji,jj,jk) * qsr(ji,jj)
                  END DO
               END DO
            END DO
            !
         ENDIF

becomes

            IF( nn_chldta ==1 ) THEN                             !*  Variable Chlorophyll
...
               zsi0r = 1.e0 / rn_si0
               zcoef  = ( 1. - rn_abs ) / 3.e0                        ! equi-partition in R-G-B
               zea(:,:,1) =    0.5 * ( qsr_b(:,:) + qsr(:,:) )
               ze0(:,:,1) = rn_abs * zea(:,:)
               ze1(:,:,1) = zcoef  * zea(:,:)
               ze2(:,:,1) = zcoef  * zea(:,:)
               ze3(:,:,1) = zcoef  * zea(:,:)
...
            ELSE                                                 !*  Constant Chlorophyll
               DO jk = 1, nksr
                  ta(:,:,jk) = ta(:,:,jk) + etot3(:,:,jk) * 0.5 * ( qsr_b(:,:) + qsr(:,:) )
               END DO
            ENDIF

         ENDIF
         !                                                ! ------------------------- !
         IF( ln_qsr_2bd ) THEN                            !  2 band light penetration !
            !                                             ! ------------------------- !
            !
            DO jk = 1, nksr
               ta(:,:,jk) = ta(:,:,jk) + etot3(:,:,jk) * 0.5 * ( qsr_b(:,:) + qsr(:,:) )
            END DO
            !
         ENDIF

CAUTION 1 : the case where the light penetration is provided by the bio-model ( i.e. lk_qsr_bio = ln_qsr_bio =.true. is not yet treated. It probably requires changes in p4zopt.F90 and trcopt.F90. Not completely obvious since this depends on the TOP time-stepping strategy.


momentum forcing term (DYN)

- dynzdf_exp.F90

applied the stress (utau, vtau) as the mean of 2 consecutive time-steps (kt+½ and kt-½)

      zrau0r = 1. / rau0
...
         ! Surface boundary condition
         DO ji = 2, jpim1
            zwy(ji,1) = utau(ji,jj) * zrau0r
            zww(ji,1) = vtau(ji,jj) * zrau0r
         END DO  

becomes

      zrau0r = 0.5 / rau0
...
         DO ji = 2, jpim1            ! Surface boundary condition
            zwy(ji,1) = ( utau_b(ji,jj) + utau(ji,jj) ) * zrau0r
            zww(ji,1) = ( utau_b(ji,jj) + vtau(ji,jj) ) * zrau0r
         END DO  

- dynzdf_imp.F90

applied the stress as the mean of 2 consecutive time-steps (kt+½ and kt-½)

      REAL(wp) ::   zrau0r, z2dtf, zcoef, zzws, zrhs   ! temporary scalars
...
      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
      DO jj = 2, jpjm1   
         DO ji = fs_2, fs_jpim1   ! vector opt.
!!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
!!!         ua(ji,jj,1) = ub(ji,jj,1)  &
!!!                      + p2dt * ( ua(ji,jj,1) + utau(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) )
            z2dtf = p2dt / ( fse3u(ji,jj,1)*rau0 )
            ua(ji,jj,1) = ub(ji,jj,1)  &
                         + p2dt *  ua(ji,jj,1) + z2dtf * utau(ji,jj)
         END DO
      END DO
...
      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
      DO jj = 2, jpjm1
         DO ji = fs_2, fs_jpim1   ! vector opt.
!!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
!!!         va(ji,jj,1) = vb(ji,jj,1)  &
!!!                      + p2dt * ( va(ji,jj,1) + vtau(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) )
            z2dtf = p2dt / ( fse3v(ji,jj,1)*rau0 )
            va(ji,jj,1) = vb(ji,jj,1)  &
                         + p2dt * va(ji,jj,1) + z2dtf * vtau(ji,jj)
         END DO
      END DO

becomes

      REAL(wp) ::   zrau0r, zcoef, zzws, zrhs     ! temporary scalars
...
      !                     !* second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
      DO jj = 2, jpjm1
         DO ji = fs_2, fs_jpim1   ! vector opt.
            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5 * ( utau_b(ji,jj) + utau(ji,jj)   &
               &                                                    / ( fse3u(ji,jj,1) * rau0 )     )
         END DO
      END DO
...
      !                    !* second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
      DO jj = 2, jpjm1
         DO ji = fs_2, fs_jpim1   ! vector opt.
            va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
               &                                                    / ( fse3v(ji,jj,1) * rau0 )       )
         END DO
      END DO



3. Modified filter

- sshwzv.F90 (currently named wzvmod)

to be explained

- tranxt.F90

to be explained

- dynnxt.F90

to be explained


Last modified 11 years ago Last modified on 2009-06-25T07:21:10+02:00

Attachments (2)

Download all attachments as: .zip