[[BR]] == '''Conservative Leap-Frog + Asselin filter'' == Last edited [[Timestamp]] [[PageOutline]] '''Purpose''': [[BR]] 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 : [[BR]] - Flow chart : general structure must be changed in a similar way to what have been done for vvl, see [wiki:2009WP/2009Stream3/VVL] [[BR]] - Forcing evaluation : the momentum and tracers forcing are evaluated at kt+1/2, no-more at kt. [[BR]] - Forcing time discretisation : 0.5 * ( Q(k+1/2) + Q(kt-1/2) ), instead of simply Q(kt). [[BR]] - 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 [[BR]] [[Image(NEW_LF_RA.jpg, 50%)]][[BR]] Illustration of forcing integration methods. Old formulation (top) and new formulation (bottom). '''Status''': [[BR]] ongoing work (NEMO-Paris). [[BR]] '''Module involved''': [[BR]] ''fldread.F90'' [[BR]] ''sbcmod.F90'' ; trasbc.F90 , traqsr.F90 , dynzdf_exp.F90 , dynzdf_imp.F90 [[BR]] ''sshnxt.F90'' , tranxt.F90 , dynnxt.F90[[BR]] ''step.F90'' ==> the version of step introduced for vvl after revision v3.1 should be used. [[BR]] Header associated with the modification (to be included in all modified modules) :[[BR]] {{{ !! 3.2 ! 2009-06 (M. Leclair, G. Madec) Forcing averaged over 2 time steps }}} '''Pending issues''': [[BR]] - 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... [[BR]] - Check the consistency with the change planned in the ''emps'' field (switch to ''fsalt'', a true salt flux). [[BR]] - modification required in limsbc, limsbc_2 and CICE interface : conversion of emps into a salt flux (new name 'fsalt' ?) [[BR]] - the modifications in trasbc interfere with the new runoff scheme (see [wiki:2009WP/2009Stream3#Substream3.2:fromNEMOteamshort-term]) [[BR]] ---- == 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. [[BR]] [[Image(NEW_LF_RA_FlowCart.jpg, 50%)]][[BR]] Step routine flow chart, form Leclair and Madec (OM 2009).[[BR]] ---- == 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.[[BR]] [[BR]] ---- == 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.[[BR]] === • 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 }}} [[BR]] === • 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 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 }}} [[BR]] === • 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 }}} [[BR]] ---- === 3. Modified filter === Blah blah blah... [[BR]]