[[BR]] == '''S3.1-d 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]] 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. [[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''': starting from the trunk after release v3.1, including vvl update (June 14th 2009) [[BR]] ''step.F90'' ==> the version introduced for vvl after revision v3.1 should be used. [[BR]] ''fldread.F90'' [[BR]] ''sbcmod.F90'' ; ''trasbc.F90'' , ''traqsr.F90'' , ''dynzdf_exp.F90'' , ''dynzdf_imp.F90'' [[BR]] ''sshnxt.F90'' , ''tranxt.F90'' , ''dynnxt.F90'' ==> the versions 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) ==== ===== - 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''' (TRA) ==== ==== - trasbc.F90 ==== applied heat and freshwater/salt fluxes (qns, 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 }}} ==== - traqsr.F90 ==== applied the solar heat flux (qsr) as the mean of 2 consecutive time-steps (kt+1/2 and kt-1/2) {{{ 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. [[BR]] ==== • '''momentum forcing term''' (DYN) ==== ==== - dynzdf_exp.F90 ==== applied the stress (utau, vtau) 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 === ==== - sshwzv.F90 (currently named wzvmod) ==== '''to be explained''' ==== - tranxt.F90 ==== '''to be explained''' ==== - dynnxt.F90 ==== '''to be explained''' [[BR]]