Version 5 (modified by gm, 15 years ago) (diff) |
---|
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+1/2, no-more at kt.
- Forcing time discretisation : 0.5 * ( Q(k+1/2) + Q(kt-1/2) ), instead of simply Q(kt).
- Filter : the RA filter no-more act on the forcing, i.e. -afp * rdt * ( Q(kt+1/2) - Q(kt-1/2) ) is removed from the filter
Illustration of forcing integration methods. Old formulation (top) and new formulation (bottom).
Status:
ongoing work (NEMO-Paris).
Module involved:
fldread.F90
sbcmod.F90 ; trasbc.F90 , traqsr.F90 , dynzdf_exp.F90 , dynzdf_imp.F90
sshnxt.F90 , tranxt.F90 , dynnxt.F90
step.F90 ==> the version of step 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+1/2 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.
Step routine flow chart, form Leclair and Madec (OM 2009).
1. Forcing evaluation
In order to evaluate the forcing at kt+1/2, all the input data in the surface module should be interpolated in time at kt+1/2. 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+1/2 and kt+1/2 which is used. Therefore, the surface module (including the sea-ice!) must now provide both kt+1/2 and kt-1/2 forcing fields, and all the module involving a forcing term have to be modified.
• surface module (sbc_oce.F90 and sbcmod.F90)
- 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
- trasbc.F90
applied the qns, qsr, emp, emps as the mean of 2 consecutive time-steps (kt+1/2 and kt-1/2)
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 DOhave 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
• momentum forcing term
- dynzdf_exp.F90
applied the stress as the mean of 2 consecutive time-steps (kt+1/2 and kt-1/2)
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+1/2 and kt-1/2)
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
Blah blah blah...
Attachments (2)
-
NEW_LF_RA.jpg
(53.6 KB) -
added by gm 15 years ago.
New versus old time discretisation of the forcing
-
NEW_LF_RA_FlowCart.jpg
(79.1 KB) -
added by gm 15 years ago.
Flow chart of step with the new Leap-Frog + Robert-Asselin filter
Download all attachments as: .zip