MODULE icedyn_adv_umx !!============================================================================== !! *** MODULE icedyn_adv_umx *** !! sea-ice : advection using the ULTIMATE-MACHO scheme !!============================================================================== !! History : 3.6 ! 2014-11 (C. Rousset, G. Madec) Original code !! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube] !!---------------------------------------------------------------------- #if defined key_si3 !!---------------------------------------------------------------------- !! 'key_si3' SI3 sea-ice model !!---------------------------------------------------------------------- !! ice_dyn_adv_umx : update the tracer trend with the 3D advection trends using a TVD scheme !! ultimate_x(_y) : compute a tracer value at velocity points using ULTIMATE scheme at various orders !! macho : ??? !! nonosc : compute monotonic tracer fluxes by a non-oscillatory algorithm !!---------------------------------------------------------------------- USE phycst ! physical constant USE dom_oce ! ocean domain USE sbc_oce , ONLY : nn_fsbc ! update frequency of surface boundary condition USE ice ! sea-ice variables USE icevar ! sea-ice: operations ! USE in_out_manager ! I/O manager USE lib_mpp ! MPP library USE lib_fortran ! fortran utilities (glob_sum + no signed zero) USE lbclnk ! lateral boundary conditions (or mpp links) IMPLICIT NONE PRIVATE PUBLIC ice_dyn_adv_umx ! called by icedyn_adv.F90 REAL(wp) :: z1_6 = 1._wp / 6._wp ! =1/6 REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120 !! * Substitutions # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/ICE 4.0 , NEMO Consortium (2018) !! $Id: icedyn_adv_umx.F90 4499 2014-02-18 15:14:31Z timgraham $ !! Software governed by the CeCILL licence (./LICENSE) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, & & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) !!---------------------------------------------------------------------- !! *** ROUTINE ice_dyn_adv_umx *** !! !! ** Purpose : Compute the now trend due to total advection of !! tracers and add it to the general trend of tracer equations !! using an "Ultimate-Macho" scheme !! !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) INTEGER , INTENT(in ) :: kt ! time step REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_i ! salt content REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content ! INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices INTEGER :: icycle ! number of sub-timestep for the advection REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers REAL(wp) :: zcfl , zdt REAL(wp) :: zeps = 0.1_wp ! shift in concentration to avoid division by 0 ! ! must be >= 0.01 and the best seems to be 0.1 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box, zfu, zfv REAL(wp), DIMENSION(jpi,jpj) :: z1_a, zh_i, zh_s, zs_i, zo_i, ze_i, ze_s, z1_ap, zh_ip !!---------------------------------------------------------------------- ! IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' ! ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! zcfl = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) zcfl = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) IF( lk_mpp ) CALL mpp_max( zcfl ) IF( zcfl > 0.5 ) THEN ; icycle = 2 ELSE ; icycle = 1 ENDIF zdt = rdt_ice / REAL(icycle) ! --- transport --- ! zudy(:,:) = pu_ice(:,:) * e2u(:,:) zvdx(:,:) = pv_ice(:,:) * e1v(:,:) ! --- define velocity for advection: u*grad(H) --- ! DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 IF ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN ; zcu_box(ji,jj) = 0._wp ELSEIF( pu_ice(ji,jj) > 0._wp ) THEN ; zcu_box(ji,jj) = pu_ice(ji-1,jj) ELSE ; zcu_box(ji,jj) = pu_ice(ji ,jj) ENDIF IF ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN ; zcv_box(ji,jj) = 0._wp ELSEIF( pv_ice(ji,jj) > 0._wp ) THEN ; zcv_box(ji,jj) = pv_ice(ji,jj-1) ELSE ; zcv_box(ji,jj) = pv_ice(ji,jj ) ENDIF END DO END DO !---------------! !== advection ==! !---------------! DO jt = 1, icycle zamsk = 1._wp CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) ) ! Open water area DO jl = 1, jpl ! to avoid a problem with the limiter nonosc when A gets close to 0 pa_i(:,:,jl) = pa_i(:,:,jl) + zeps * tmask(:,:,1) ! WHERE( pa_i(:,:,jl) > epsi20 ) ; z1_a(:,:) = 1._wp / pa_i(:,:,jl) ELSEWHERE ; z1_a(:,:) = 0. !!; pa_i(:,:,jl) = 0. ; pv_i(:,:,jl) = 0. END WHERE ! zamsk = 1._wp CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), zfu, zfv ) ! Ice area !!zfu = zudy ; zfv = zvdx !!CALL adv_umx( kn_umx, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, pv_i(:,:,jl), pv_i(:,:,jl) ) zamsk = 0._wp ; zh_i(:,:) = pv_i (:,:,jl) * z1_a(:,:) CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_i(:,:), pv_i (:,:,jl) ) ! Ice volume zamsk = 0._wp ; zh_s(:,:) = pv_s (:,:,jl) * z1_a(:,:) CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_s(:,:), pv_s (:,:,jl) ) ! Snw volume zamsk = 0._wp ; zs_i(:,:) = psv_i(:,:,jl) * z1_a(:,:) CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zs_i(:,:), psv_i(:,:,jl) ) ! Salt content zamsk = 0._wp ; zo_i(:,:) = poa_i(:,:,jl) * z1_a(:,:) CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zo_i(:,:), poa_i(:,:,jl) ) ! Age content DO jk = 1, nlay_i zamsk = 0._wp ; ze_i(:,:) = pe_i(:,:,jk,jl) * z1_a(:,:) CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu, zfv, zcu_box, zcv_box, ze_i(:,:), pe_i(:,:,jk,jl) ) ! Ice heat content END DO DO jk = 1, nlay_s zamsk = 0._wp ; ze_s(:,:) = pe_s(:,:,jk,jl) * z1_a(:,:) CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu, zfv, zcu_box, zcv_box, ze_s(:,:), pe_s(:,:,jk,jl) ) ! Snw heat content END DO ! IF ( ln_pnd_H12 ) THEN ! melt ponds ! to avoid a problem with the limiter nonosc when A gets close to 0 pa_ip(:,:,jl) = pa_ip(:,:,jl) + zeps * tmask(:,:,1) ! WHERE( pa_ip(:,:,jl) > epsi20 ) ; z1_ap(:,:) = 1._wp / pa_ip(:,:,jl) ELSEWHERE ; z1_ap(:,:) = 0. END WHERE ! zamsk = 1._wp CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), zfu, zfv ) ! mp fraction zamsk = 0._wp ; zh_ip(:,:) = pv_ip(:,:,jl) * z1_a(:,:) CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_ip(:,:), pv_ip(:,:,jl) ) ! mp volume ENDIF ! ! to avoid a problem with the limiter nonosc when A gets close to 0 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 !pa_i(ji,jj,jl) = ( pa_i(ji,jj,jl) - zeps ) * tmask(ji,jj,1) pa_i(ji,jj,jl) = ( pa_i(ji,jj,jl) - zeps & & + zeps * ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt & & ) * tmask(ji,jj,1) IF ( ln_pnd_H12 ) THEN ! melt ponds pa_ip(ji,jj,jl) = ( pa_ip(ji,jj,jl) - zeps & & + zeps * ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt & & ) * tmask(ji,jj,1) ENDIF END DO END DO CALL lbc_lnk( pa_i(:,:,jl), 'T', 1. ) ! END DO END DO ! END SUBROUTINE ice_dyn_adv_umx SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pfu, pfv ) !!---------------------------------------------------------------------- !! *** ROUTINE adv_umx *** !! !! ** Purpose : Compute the now trend due to total advection of !! tracers and add it to the general trend of tracer equations !! !! ** Method : TVD scheme, i.e. 2nd order centered scheme with !! corrected flux (monotonic correction) !! note: - this advection scheme needs a leap-frog time scheme !! !! ** Action : - pt the after advective tracer !!---------------------------------------------------------------------- REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) INTEGER , INTENT(in ) :: jt ! number of sub-iteration INTEGER , INTENT(in ) :: kt ! number of iteration REAL(wp) , INTENT(in ) :: pdt ! tracer time-step REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 or u*a*e2u REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer field REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: ptc ! tracer content field REAL(wp), DIMENSION(jpi,jpj), INTENT( out), OPTIONAL :: pfu, pfv ! high order fluxes ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: ztra ! local scalar !!clem REAL(wp) :: zeps = 1.e-02 ! local scalar INTEGER :: kn_limiter = 1 ! 1=nonosc ; 2=superbee ; 3=h3(rachid) REAL(wp), DIMENSION(jpi,jpj) :: zfu_ho , zfv_ho , zt_u, zt_v, zpt REAL(wp), DIMENSION(jpi,jpj) :: zfu_ups, zfv_ups, zt_ups ! only for nonosc !!---------------------------------------------------------------------- ! ! add a constant value to avoid problems with zeros DO jj = 1, jpj DO ji = 1, jpi zpt(ji,jj) = pt(ji,jj) !!clem + zeps * tmask(ji,jj,1) END DO END DO ! upstream (_ups) advection with initial mass fluxes ! --------------------------------------------------- ! fluxes DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj) zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1) END DO END DO ! guess after content field with upstream scheme DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj ) & & + zfv_ups(ji,jj) - zfv_ups(ji ,jj-1) ) * r1_e1e2t(ji,jj) zt_ups(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) END DO END DO CALL lbc_lnk( zt_ups, 'T', 1. ) ! High order (_ho) fluxes ! ----------------------- SELECT CASE( kn_umx ) ! CASE ( 20 ) !== centered second order ==! ! CALL cen2( kn_limiter, jt, kt, pdt, zpt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho, & & zt_ups, zfu_ups, zfv_ups ) ! CASE ( 1:5 ) !== 1st to 5th order ULTIMATE-MACHO scheme ==! ! CALL macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, zpt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho, & & zt_ups, zfu_ups, zfv_ups ) ! END SELECT ! output high order fluxes ! ------------------------ IF( PRESENT(pfu) ) THEN DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 pfu(ji,jj) = zfu_ho(ji,jj) !!clem - zeps * puc(ji,jj) pfv(ji,jj) = zfv_ho(ji,jj) !!clem - zeps * pvc(ji,jj) END DO END DO !!CALL lbc_lnk( pfu, 'U', -1. ) ! clem: not needed I think !!CALL lbc_lnk( pfv, 'V', -1. ) ENDIF ! final trend with corrected fluxes ! ------------------------------------ DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ztra = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) & ! Div(uaH) or Div(ua) !!clem & + ( puc (ji,jj) - puc (ji-1,jj) + pvc (ji,jj) - pvc (ji,jj-1) ) * zeps & ! epsi * Div(ua) or epsi * Div(u) & ) * r1_e1e2t(ji,jj) ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) END DO END DO CALL lbc_lnk( ptc, 'T', 1. ) ! END SUBROUTINE adv_umx SUBROUTINE cen2( kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, & & pt_ups, pfu_ups, pfv_ups ) !!--------------------------------------------------------------------- !! *** ROUTINE macho *** !! !! ** Purpose : compute !! !! ** Method : ... ??? !! TIM = transient interpolation Modeling !! !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kn_limiter ! limiter INTEGER , INTENT(in ) :: jt ! number of sub-iteration INTEGER , INTENT(in ) :: kt ! number of iteration REAL(wp) , INTENT(in ) :: pdt ! tracer time-step REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu, pv ! 2 ice velocity components REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc, pvc ! 2 ice velocity * A components REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptc ! tracer content at before time step REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt_ups ! upstream guess of tracer content REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pfu_ups, pfv_ups ! upstream fluxes ! INTEGER :: ji, jj ! dummy loop indices LOGICAL :: ll_xy = .TRUE. REAL(wp), DIMENSION(jpi,jpj) :: zzt !!---------------------------------------------------------------------- ! IF( .NOT.ll_xy ) THEN !-- no alternate directions --! ! DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) ) pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) ) END DO END DO IF ( kn_limiter == 1 ) THEN CALL nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) ELSEIF( kn_limiter == 2 ) THEN CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) ELSEIF( kn_limiter == 3 ) THEN CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) ENDIF ! ELSE !-- alternate directions --! ! IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! ! ! flux in x-direction DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) ) END DO END DO IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) ! first guess of tracer content from u-flux DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) END DO END DO CALL lbc_lnk( zzt, 'T', 1. ) ! flux in y-direction DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 pfv_ho(ji,jj) = 0.5 * pv(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) ) END DO END DO IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) ELSE !== even ice time step: adv_y then adv_x ==! ! ! flux in y-direction DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) ) END DO END DO IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) ! ! first guess of tracer content from v-flux DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) END DO END DO CALL lbc_lnk( zzt, 'T', 1. ) ! ! flux in x-direction DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 pfu_ho(ji,jj) = 0.5 * pu(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) ) END DO END DO IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) ENDIF IF( kn_limiter == 1 ) CALL nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) ENDIF END SUBROUTINE cen2 SUBROUTINE macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, pt_u, pt_v, pfu_ho, pfv_ho, & & pt_ups, pfu_ups, pfv_ups ) !!--------------------------------------------------------------------- !! *** ROUTINE macho *** !! !! ** Purpose : compute !! !! ** Method : ... ??? !! TIM = transient interpolation Modeling !! !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. !!---------------------------------------------------------------------- REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) INTEGER , INTENT(in ) :: kn_limiter ! limiter INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) INTEGER , INTENT(in ) :: jt ! number of sub-iteration INTEGER , INTENT(in ) :: kt ! number of iteration REAL(wp) , INTENT(in ) :: pdt ! tracer time-step REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu, pv ! 2 ice velocity components REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc, pvc ! 2 ice velocity * A components REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptc ! tracer content at before time step REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_u, pt_v ! tracer at u- and v-points REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt_ups ! upstream guess of tracer content REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pfu_ups, pfv_ups ! upstream fluxes ! INTEGER :: ji, jj ! dummy loop indices REAL(wp), DIMENSION(jpi,jpj) :: zzt !!---------------------------------------------------------------------- ! IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! ! ! !-- ultimate interpolation of pt at u-point --! CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) ! !-- limiter in x --! IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) ! !-- advective form update in zzt --! DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj) & & - pt (ji,jj) * pdt * ( pu (ji,jj) - pu (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) END DO END DO CALL lbc_lnk( zzt, 'T', 1. ) ! !-- ultimate interpolation of pt at v-point --! CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho ) ! !-- limiter in y --! IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) ! ELSE !== even ice time step: adv_y then adv_x ==! ! ! !-- ultimate interpolation of pt at v-point --! CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) ! !-- limiter in y --! IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) ! !-- advective form update in zzt --! DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj) & & - pt (ji,jj) * pdt * ( pv (ji,jj) - pv (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) END DO END DO CALL lbc_lnk( zzt, 'T', 1. ) ! !-- ultimate interpolation of pt at u-point --! CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho ) ! !-- limiter in x --! IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) ! ENDIF IF( kn_limiter == 1 ) CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) ! END SUBROUTINE macho SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) !!--------------------------------------------------------------------- !! *** ROUTINE ultimate_x *** !! !! ** Purpose : compute !! !! ** Method : ... ??? !! TIM = transient interpolation Modeling !! !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) REAL(wp) , INTENT(in ) :: pdt ! tracer time-step REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu ! ice i-velocity component REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc ! ice i-velocity * A component REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_u ! tracer at u-point REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pfu_ho ! high order flux ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: zcu, zdx2, zdx4 ! - - REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu2, ztu3, ztu4 !!---------------------------------------------------------------------- ! ! !-- Laplacian in i-direction --! DO jj = 2, jpjm1 ! First derivative (gradient) DO ji = 1, fs_jpim1 ztu1(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1) END DO ! ! Second derivative (Laplacian) DO ji = fs_2, fs_jpim1 ztu2(ji,jj) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj) END DO END DO CALL lbc_lnk( ztu2, 'T', 1. ) ! ! !-- BiLaplacian in i-direction --! DO jj = 2, jpjm1 ! Third derivative DO ji = 1, fs_jpim1 ztu3(ji,jj) = ( ztu2(ji+1,jj) - ztu2(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1) END DO ! ! Fourth derivative DO ji = fs_2, fs_jpim1 ztu4(ji,jj) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj) END DO END DO CALL lbc_lnk( ztu4, 'T', 1. ) ! ! SELECT CASE (kn_umx ) ! CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) ! DO jj = 2, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) & & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) ) END DO END DO ! CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) ! DO jj = 2, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) & & - zcu * ( pt(ji+1,jj) - pt(ji,jj) ) ) END DO END DO ! CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) ! DO jj = 2, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) zdx2 = e1u(ji,jj) * e1u(ji,jj) !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj) + pt (ji,jj) & & - zcu * ( pt (ji+1,jj) - pt (ji,jj) ) ) & & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj) + ztu2(ji,jj) & & - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) ) ) END DO END DO ! CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) ! DO jj = 2, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) zdx2 = e1u(ji,jj) * e1u(ji,jj) !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj) + pt (ji,jj) & & - zcu * ( pt (ji+1,jj) - pt (ji,jj) ) ) & & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj) + ztu2(ji,jj) & & - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) ) ) END DO END DO ! CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) ! DO jj = 2, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) zdx2 = e1u(ji,jj) * e1u(ji,jj) !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) zdx4 = zdx2 * zdx2 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj) + pt (ji,jj) & & - zcu * ( pt (ji+1,jj) - pt (ji,jj) ) ) & & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj) + ztu2(ji,jj) & & - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) ) & & + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj) + ztu4(ji,jj) & & - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj) - ztu4(ji,jj) ) ) ) END DO END DO ! END SELECT ! !-- High order flux in i-direction --! DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. pfu_ho(ji,jj) = puc(ji,jj) * pt_u(ji,jj) END DO END DO ! END SUBROUTINE ultimate_x SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) !!--------------------------------------------------------------------- !! *** ROUTINE ultimate_y *** !! !! ** Purpose : compute !! !! ** Method : ... ??? !! TIM = transient interpolation Modeling !! !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) REAL(wp) , INTENT(in ) :: pdt ! tracer time-step REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pv ! ice j-velocity component REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pvc ! ice j-velocity*A component REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_v ! tracer at v-point REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pfv_ho ! high order flux ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: zcv, zdy2, zdy4 ! - - REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv2, ztv3, ztv4 !!---------------------------------------------------------------------- ! ! !-- Laplacian in j-direction --! DO jj = 1, jpjm1 ! First derivative (gradient) DO ji = fs_2, fs_jpim1 ztv1(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) END DO END DO DO jj = 2, jpjm1 ! Second derivative (Laplacian) DO ji = fs_2, fs_jpim1 ztv2(ji,jj) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj) END DO END DO CALL lbc_lnk( ztv2, 'T', 1. ) ! ! !-- BiLaplacian in j-direction --! DO jj = 1, jpjm1 ! First derivative DO ji = fs_2, fs_jpim1 ztv3(ji,jj) = ( ztv2(ji,jj+1) - ztv2(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) END DO END DO DO jj = 2, jpjm1 ! Second derivative DO ji = fs_2, fs_jpim1 ztv4(ji,jj) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj) END DO END DO CALL lbc_lnk( ztv4, 'T', 1. ) ! ! SELECT CASE (kn_umx ) ! CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) DO jj = 1, jpjm1 DO ji = fs_2, fs_jpim1 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) & & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) END DO END DO ! CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) DO jj = 1, jpjm1 DO ji = fs_2, fs_jpim1 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) & & - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) ) END DO END DO CALL lbc_lnk( pt_v, 'V', 1. ) ! CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) DO jj = 1, jpjm1 DO ji = fs_2, fs_jpim1 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) zdy2 = e2v(ji,jj) * e2v(ji,jj) !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) & & - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) & & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1) + ztv2(ji,jj) & & - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) ) END DO END DO ! CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) DO jj = 1, jpjm1 DO ji = fs_2, fs_jpim1 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) zdy2 = e2v(ji,jj) * e2v(ji,jj) !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) & & - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) & & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1) + ztv2(ji,jj) & & - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) ) END DO END DO ! CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) DO jj = 1, jpjm1 DO ji = fs_2, fs_jpim1 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) zdy2 = e2v(ji,jj) * e2v(ji,jj) !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) zdy4 = zdy2 * zdy2 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) & & - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) & & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1) + ztv2(ji,jj) & & - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) & & + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1) + ztv4(ji,jj) & & - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1) - ztv4(ji,jj) ) ) ) END DO END DO ! END SELECT ! !-- High order flux in j-direction --! DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. pfv_ho(ji,jj) = pvc(ji,jj) * pt_v(ji,jj) END DO END DO ! END SUBROUTINE ultimate_y SUBROUTINE nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) !!--------------------------------------------------------------------- !! *** ROUTINE nonosc *** !! !! ** Purpose : compute monotonic tracer fluxes from the upstream !! scheme and the before field by a nonoscillatory algorithm !! !! ** Method : ... ??? !! warning : ptc and pt_ups must be masked, but the boundaries !! conditions on the fluxes are not necessary zalezak (1979) !! drange (1995) multi-dimensional forward-in-time and upstream- !! in-space based differencing for fluid !!---------------------------------------------------------------------- REAL(wp) , INTENT(in ) :: pdt ! tracer time-step REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: ptc, pt_ups ! before field & upstream guess of after field REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pfv_ups, pfu_ups ! upstream flux REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: pfv_ho, pfu_ho ! monotonic flux ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: zpos, zneg, zbig, zsml, z1_dt ! local scalars REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zdiv !!---------------------------------------------------------------------- zbig = 1.e+40_wp zsml = 1.e-15_wp ! test on divergence DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zdiv(ji,jj) = - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) + pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) END DO END DO CALL lbc_lnk( zdiv, 'T', 1. ) ! Lateral boundary conditions (unchanged sign) ! antidiffusive flux : high order minus low order ! -------------------------------------------------- DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_ups(ji,jj) pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_ups(ji,jj) END DO END DO ! Search local extrema ! -------------------- ! max/min of ptc & pt_ups with large negative/positive value (-/+zbig) outside ice cover DO jj = 1, jpj DO ji = fs_2, fs_jpim1 IF( ptc(ji,jj) == 0._wp .AND. pt_ups(ji,jj) == 0._wp .AND. zdiv(ji,jj) == 0._wp ) THEN zbup(ji,jj) = -zbig zbdo(ji,jj) = zbig ELSE zbup(ji,jj) = MAX( ptc(ji,jj) , pt_ups(ji,jj) ) zbdo(ji,jj) = MIN( ptc(ji,jj) , pt_ups(ji,jj) ) ENDIF END DO END DO z1_dt = 1._wp / pdt DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! zup = MAX( zbup(ji,jj), zbup(ji-1,jj), zbup(ji+1,jj), zbup(ji,jj-1), zbup(ji,jj+1) ) ! search max/min in neighbourhood zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) ) ! zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji ,jj) ) & & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj ) ) ! positive/negative part of the flux zneg = MAX( 0., pfu_ho(ji ,jj) ) - MIN( 0., pfu_ho(ji-1,jj) ) & & + MAX( 0., pfv_ho(ji,jj ) ) - MIN( 0., pfv_ho(ji,jj-1) ) ! ! ! up & down beta terms !!clem zbetup(ji,jj) = ( zup - pt_ups(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt !!clem zbetdo(ji,jj) = ( pt_ups(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt IF( zpos >= epsi20 ) THEN zbetup(ji,jj) = ( zup - pt_ups(ji,jj) ) / zpos * e1e2t(ji,jj) * z1_dt ELSE zbetup(ji,jj) = zbig ENDIF ! IF( zneg >= epsi20 ) THEN zbetdo(ji,jj) = ( pt_ups(ji,jj) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt ELSE zbetdo(ji,jj) = zbig ENDIF ! END DO END DO CALL lbc_lnk_multi( zbetup, 'T', 1., zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign) ! monotonic flux in the y direction ! --------------------------------- DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) ) zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) ) zcu = 0.5 + SIGN( 0.5 , pfu_ho(ji,jj) ) ! pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_ups(ji,jj) END DO END DO DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) ) zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) ) zcv = 0.5 + SIGN( 0.5 , pfv_ho(ji,jj) ) ! pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_ups(ji,jj) END DO END DO ! END SUBROUTINE nonosc_2d SUBROUTINE limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) !!--------------------------------------------------------------------- !! *** ROUTINE limiter_x *** !! !! ** Purpose : compute flux limiter !!---------------------------------------------------------------------- REAL(wp) , INTENT(in ) :: pdt ! tracer time-step REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pu ! ice i-velocity => u*e2 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: puc ! ice i-velocity *A => u*e2*a REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pt ! ice tracer REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: pfu_ho ! high order flux REAL(wp), DIMENSION (jpi,jpj), INTENT(in ), OPTIONAL :: pfu_ups ! upstream flux ! REAL(wp) :: Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr INTEGER :: ji, jj ! dummy loop indices REAL(wp), DIMENSION (jpi,jpj) :: zslpx ! tracer slopes !!---------------------------------------------------------------------- ! DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zslpx(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * umask(ji,jj,1) END DO END DO CALL lbc_lnk( zslpx, 'U', -1.) ! lateral boundary cond. DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj) Rjm = zslpx(ji-1,jj) Rj = zslpx(ji ,jj) Rjp = zslpx(ji+1,jj) IF( PRESENT(pfu_ups) ) THEN IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm ELSE ; Rr = Rjp ENDIF zh3 = pfu_ho(ji,jj) - pfu_ups(ji,jj) IF( Rj > 0. ) THEN zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(puc(ji,jj)), & & MIN( 2. * Rr * 0.5 * ABS(puc(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) ) ELSE zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(puc(ji,jj)), & & MIN(-2. * Rr * 0.5 * ABS(puc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) ) ENDIF pfu_ho(ji,jj) = pfu_ups(ji,jj) + zlimiter ELSE IF( Rj /= 0. ) THEN IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj ELSE ; Cr = Rjp / Rj ENDIF ELSE Cr = 0. !IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm * 1.e20 !ELSE ; Cr = Rjp * 1.e20 !ENDIF ENDIF ! -- superbee -- zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) ! -- van albada 2 -- !!zpsi = 2.*Cr / (Cr*Cr+1.) ! -- sweby (with beta=1) -- !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) ! -- van Leer -- !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) ! -- ospre -- !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) ! -- koren -- !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) ! -- charm -- !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) !ELSE ; zpsi = 0. !ENDIF ! -- van albada 1 -- !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) ! -- smart -- !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) ! -- umist -- !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) ! high order flux corrected by the limiter pfu_ho(ji,jj) = pfu_ho(ji,jj) - ABS( puc(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5 ENDIF END DO END DO CALL lbc_lnk( pfu_ho, 'U', -1.) ! lateral boundary cond. ! END SUBROUTINE limiter_x SUBROUTINE limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) !!--------------------------------------------------------------------- !! *** ROUTINE limiter_y *** !! !! ** Purpose : compute flux limiter !!---------------------------------------------------------------------- REAL(wp) , INTENT(in ) :: pdt ! tracer time-step REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pv ! ice i-velocity => u*e2 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pvc ! ice i-velocity *A => u*e2*a REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pt ! ice tracer REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: pfv_ho ! high order flux REAL(wp), DIMENSION (jpi,jpj), INTENT(in ), OPTIONAL :: pfv_ups ! upstream flux ! REAL(wp) :: Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr INTEGER :: ji, jj ! dummy loop indices REAL(wp), DIMENSION (jpi,jpj) :: zslpy ! tracer slopes !!---------------------------------------------------------------------- ! DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zslpy(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * vmask(ji,jj,1) END DO END DO CALL lbc_lnk( zslpy, 'V', -1.) ! lateral boundary cond. DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj) Rjm = zslpy(ji,jj-1) Rj = zslpy(ji,jj ) Rjp = zslpy(ji,jj+1) IF( PRESENT(pfv_ups) ) THEN IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm ELSE ; Rr = Rjp ENDIF zh3 = pfv_ho(ji,jj) - pfv_ups(ji,jj) IF( Rj > 0. ) THEN zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pvc(ji,jj)), & & MIN( 2. * Rr * 0.5 * ABS(pvc(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) ) ELSE zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pvc(ji,jj)), & & MIN(-2. * Rr * 0.5 * ABS(pvc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) ) ENDIF pfv_ho(ji,jj) = pfv_ups(ji,jj) + zlimiter ELSE IF( Rj /= 0. ) THEN IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj ELSE ; Cr = Rjp / Rj ENDIF ELSE Cr = 0. !IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm * 1.e20 !ELSE ; Cr = Rjp * 1.e20 !ENDIF ENDIF ! -- superbee -- zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) ! -- van albada 2 -- !!zpsi = 2.*Cr / (Cr*Cr+1.) ! -- sweby (with beta=1) -- !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) ! -- van Leer -- !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) ! -- ospre -- !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) ! -- koren -- !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) ! -- charm -- !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) !ELSE ; zpsi = 0. !ENDIF ! -- van albada 1 -- !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) ! -- smart -- !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) ! -- umist -- !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) ! high order flux corrected by the limiter pfv_ho(ji,jj) = pfv_ho(ji,jj) - ABS( pvc(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5 ENDIF END DO END DO CALL lbc_lnk( pfv_ho, 'V', -1.) ! lateral boundary cond. ! END SUBROUTINE limiter_y #else !!---------------------------------------------------------------------- !! Default option Dummy module NO SI3 sea-ice model !!---------------------------------------------------------------------- #endif !!====================================================================== END MODULE icedyn_adv_umx