Changeset 10267
- Timestamp:
- 2018-10-31T19:09:43+01:00 (6 years ago)
- Location:
- NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE
- Files:
-
- 2 edited
-
icedyn.F90 (modified) (2 diffs)
-
icedyn_adv_umx.F90 (modified) (27 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn.F90
r9604 r10267 71 71 INTEGER, INTENT(in) :: kt ! ice time step 72 72 !! 73 INTEGER :: ji, jj, jl ! dummy loop indices 74 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhmax 73 INTEGER :: ji, jj, jl ! dummy loop indices 74 REAL(wp) :: zcoefu, zcoefv 75 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhmax 76 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdivu_i 75 77 !!-------------------------------------------------------------------- 76 78 ! … … 118 120 119 121 CASE ( np_dynADV ) !== pure advection ==! (prescribed velocities) 122 ALLOCATE( zdivu_i(jpi,jpj) ) 123 120 124 u_ice(:,:) = rn_uice * umask(:,:,1) 121 125 v_ice(:,:) = rn_vice * vmask(:,:,1) 122 !!CALL RANDOM_NUMBER(u_ice(:,:)) 123 !!CALL RANDOM_NUMBER(v_ice(:,:)) 126 !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1) 127 !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 128 ! --- monotonicity test from Lipscomb et al 2004 --- ! 129 ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length 130 ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s 131 !DO jj = 1, jpj 132 ! DO ji = 1, jpi 133 ! zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. ) 134 ! zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. ) 135 ! u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1., zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1) 136 ! v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1., zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1) 137 ! END DO 138 !END DO 139 ! --- 124 140 CALL ice_dyn_adv ( kt ) ! -- advection of ice 141 142 ! diagnostics: divergence at T points 143 DO jj = 2, jpjm1 144 DO ji = 2, jpim1 145 zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 146 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 147 END DO 148 END DO 149 CALL lbc_lnk( zdivu_i, 'T', 1. ) 150 IF( iom_use('icediv') ) CALL iom_put( "icediv" , zdivu_i(:,:) ) 151 152 DEALLOCATE( zdivu_i ) 125 153 126 154 END SELECT -
NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv_umx.F90
r9929 r10267 14 14 !! ultimate_x(_y) : compute a tracer value at velocity points using ULTIMATE scheme at various orders 15 15 !! macho : ??? 16 !! nonosc _2d: compute monotonic tracer fluxes by a non-oscillatory algorithm16 !! nonosc : compute monotonic tracer fluxes by a non-oscillatory algorithm 17 17 !!---------------------------------------------------------------------- 18 18 USE phycst ! physical constant … … 20 20 USE sbc_oce , ONLY : nn_fsbc ! update frequency of surface boundary condition 21 21 USE ice ! sea-ice variables 22 USE icevar ! sea-ice: operations 22 23 ! 23 24 USE in_out_manager ! I/O manager … … 43 44 CONTAINS 44 45 45 SUBROUTINE ice_dyn_adv_umx( k _order, kt, pu_ice, pv_ice, &46 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )46 SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, & 47 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 47 48 !!---------------------------------------------------------------------- 48 49 !! *** ROUTINE ice_dyn_adv_umx *** … … 54 55 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 55 56 !!---------------------------------------------------------------------- 56 INTEGER , INTENT(in ) :: k _order ! order of the scheme (1-5 or 20)57 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 57 58 INTEGER , INTENT(in ) :: kt ! time step 58 59 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity … … 70 71 ! 71 72 INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices 72 INTEGER :: initad ! number of sub-timestep for the advection 73 REAL(wp) :: zcfl , zusnit, zdt ! - - 74 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zudy, zvdx, zcu_box, zcv_box 73 INTEGER :: icycle ! number of sub-timestep for the advection 74 REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers 75 REAL(wp) :: zcfl , zdt 76 REAL(wp) :: zeps = 0.1_wp ! shift in concentration to avoid division by 0 77 ! ! must be >= 0.01 and the best seems to be 0.1 78 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box, zfu, zfv 79 REAL(wp), DIMENSION(jpi,jpj) :: z1_a, zh_i, zh_s, zs_i, zo_i, ze_i, ze_s, z1_ap, zh_ip 75 80 !!---------------------------------------------------------------------- 76 81 ! 77 82 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 78 !79 ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) )80 83 ! 81 84 ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! … … 84 87 IF( lk_mpp ) CALL mpp_max( zcfl ) 85 88 86 IF( zcfl > 0.5 ) THEN ; i nitad = 2 ; zusnit = 0.5_wp87 ELSE ; i nitad = 1 ; zusnit = 1.0_wp89 IF( zcfl > 0.5 ) THEN ; icycle = 2 90 ELSE ; icycle = 1 88 91 ENDIF 89 92 90 zdt = rdt_ice / REAL(i nitad)93 zdt = rdt_ice / REAL(icycle) 91 94 92 95 ! --- transport --- ! … … 112 115 !== advection ==! 113 116 !---------------! 114 DO jt = 1, initad 115 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) ) ! Open water area 117 DO jt = 1, icycle 118 119 zamsk = 1._wp 120 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) ) ! Open water area 121 116 122 DO jl = 1, jpl 117 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) ) ! Ice area 118 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) ) ! Ice volume 119 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, psv_i(:,:,jl) ) ! Salt content 120 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i(:,:,jl) ) ! Age content 123 ! to avoid a problem with the limiter nonosc when A gets close to 0 124 pa_i(:,:,jl) = pa_i(:,:,jl) + zeps * tmask(:,:,1) 125 ! 126 WHERE( pa_i(:,:,jl) > epsi20 ) ; z1_a(:,:) = 1._wp / pa_i(:,:,jl) 127 ELSEWHERE ; z1_a(:,:) = 0. !!; pa_i(:,:,jl) = 0. ; pv_i(:,:,jl) = 0. 128 END WHERE 129 ! 130 zamsk = 1._wp 131 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 132 !!zfu = zudy ; zfv = zvdx 133 !!CALL adv_umx( kn_umx, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, pv_i(:,:,jl), pv_i(:,:,jl) ) 134 zamsk = 0._wp ; zh_i(:,:) = pv_i (:,:,jl) * z1_a(:,:) 135 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_i(:,:), pv_i (:,:,jl) ) ! Ice volume 136 zamsk = 0._wp ; zh_s(:,:) = pv_s (:,:,jl) * z1_a(:,:) 137 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_s(:,:), pv_s (:,:,jl) ) ! Snw volume 138 zamsk = 0._wp ; zs_i(:,:) = psv_i(:,:,jl) * z1_a(:,:) 139 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zs_i(:,:), psv_i(:,:,jl) ) ! Salt content 140 zamsk = 0._wp ; zo_i(:,:) = poa_i(:,:,jl) * z1_a(:,:) 141 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zo_i(:,:), poa_i(:,:,jl) ) ! Age content 121 142 DO jk = 1, nlay_i 122 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) ) ! Ice heat content123 END DO124 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) ) ! Snow volume143 zamsk = 0._wp ; ze_i(:,:) = pe_i(:,:,jk,jl) * z1_a(:,:) 144 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 145 END DO 125 146 DO jk = 1, nlay_s 126 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,jk,jl) ) ! Snow heat content 127 END DO 128 IF ( ln_pnd_H12 ) THEN 129 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) ) ! Melt pond fraction 130 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) ) ! Melt pond volume 147 zamsk = 0._wp ; ze_s(:,:) = pe_s(:,:,jk,jl) * z1_a(:,:) 148 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 149 END DO 150 ! 151 IF ( ln_pnd_H12 ) THEN ! melt ponds 152 ! to avoid a problem with the limiter nonosc when A gets close to 0 153 pa_ip(:,:,jl) = pa_ip(:,:,jl) + zeps * tmask(:,:,1) 154 ! 155 WHERE( pa_ip(:,:,jl) > epsi20 ) ; z1_ap(:,:) = 1._wp / pa_ip(:,:,jl) 156 ELSEWHERE ; z1_ap(:,:) = 0. 157 END WHERE 158 ! 159 zamsk = 1._wp 160 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 161 zamsk = 0._wp ; zh_ip(:,:) = pv_ip(:,:,jl) * z1_a(:,:) 162 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zfu , zfv , zcu_box, zcv_box, zh_ip(:,:), pv_ip(:,:,jl) ) ! mp volume 131 163 ENDIF 132 END DO 133 END DO 134 ! 135 DEALLOCATE( zudy, zvdx, zcu_box, zcv_box ) 164 ! 165 ! to avoid a problem with the limiter nonosc when A gets close to 0 166 DO jj = 2, jpjm1 167 DO ji = fs_2, fs_jpim1 168 !pa_i(ji,jj,jl) = ( pa_i(ji,jj,jl) - zeps ) * tmask(ji,jj,1) 169 pa_i(ji,jj,jl) = ( pa_i(ji,jj,jl) - zeps & 170 & + zeps * ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt & 171 & ) * tmask(ji,jj,1) 172 IF ( ln_pnd_H12 ) THEN ! melt ponds 173 pa_ip(ji,jj,jl) = ( pa_ip(ji,jj,jl) - zeps & 174 & + zeps * ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt & 175 & ) * tmask(ji,jj,1) 176 ENDIF 177 END DO 178 END DO 179 CALL lbc_lnk( pa_i(:,:,jl), 'T', 1. ) 180 ! 181 END DO 182 183 END DO 136 184 ! 137 185 END SUBROUTINE ice_dyn_adv_umx 138 186 139 187 140 SUBROUTINE adv_umx( k_order, kt, pdt, puc, pvc, pubox, pvbox, ptc)188 SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pfu, pfv ) 141 189 !!---------------------------------------------------------------------- 142 190 !! *** ROUTINE adv_umx *** … … 151 199 !! ** Action : - pt the after advective tracer 152 200 !!---------------------------------------------------------------------- 153 INTEGER , INTENT(in ) :: k_order ! order of the ULTIMATE scheme 154 INTEGER , INTENT(in ) :: kt ! number of iteration 155 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 156 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 157 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity 158 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: ptc ! tracer content field 201 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 202 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 203 INTEGER , INTENT(in ) :: jt ! number of sub-iteration 204 INTEGER , INTENT(in ) :: kt ! number of iteration 205 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 206 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2 207 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 or u*a*e2u 208 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity 209 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer field 210 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: ptc ! tracer content field 211 REAL(wp), DIMENSION(jpi,jpj), INTENT( out), OPTIONAL :: pfu, pfv ! high order fluxes 159 212 ! 160 213 INTEGER :: ji, jj ! dummy loop indices 161 214 REAL(wp) :: ztra ! local scalar 162 REAL(wp), DIMENSION(jpi,jpj) :: zfu_ups, zfu_ho, zt_u, zt_ups 163 REAL(wp), DIMENSION(jpi,jpj) :: zfv_ups, zfv_ho, zt_v, ztrd 164 !!---------------------------------------------------------------------- 165 ! 166 ! upstream advection with initial mass fluxes & intermediate update 167 ! -------------------------------------------------------------------- 168 DO jj = 1, jpjm1 ! upstream tracer flux in the i and j direction 169 DO ji = 1, fs_jpim1 ! vector opt. 170 zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * ptc(ji,jj) + MIN( puc(ji,jj), 0._wp ) * ptc(ji+1,jj) 171 zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * ptc(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * ptc(ji,jj+1) 172 END DO 173 END DO 174 175 DO jj = 2, jpjm1 ! total intermediate advective trends 176 DO ji = fs_2, fs_jpim1 ! vector opt. 177 ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj ) & 178 & + zfv_ups(ji,jj) - zfv_ups(ji ,jj-1) ) * r1_e1e2t(ji,jj) 179 ! 180 ztrd(ji,jj) = ztra ! upstream trend [ -div(uh) or -div(uhT) ] 181 zt_ups (ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) ! guess after content field with monotonic scheme 182 END DO 183 END DO 184 CALL lbc_lnk( zt_ups, 'T', 1. ) ! Lateral boundary conditions (unchanged sign) 215 !!clem REAL(wp) :: zeps = 1.e-02 ! local scalar 216 INTEGER :: kn_limiter = 1 ! 1=nonosc ; 2=superbee ; 3=h3(rachid) 217 REAL(wp), DIMENSION(jpi,jpj) :: zfu_ho , zfv_ho , zt_u, zt_v, zpt 218 REAL(wp), DIMENSION(jpi,jpj) :: zfu_ups, zfv_ups, zt_ups ! only for nonosc 219 !!---------------------------------------------------------------------- 220 ! 221 ! add a constant value to avoid problems with zeros 222 DO jj = 1, jpj 223 DO ji = 1, jpi 224 zpt(ji,jj) = pt(ji,jj) !!clem + zeps * tmask(ji,jj,1) 225 END DO 226 END DO 227 228 ! upstream (_ups) advection with initial mass fluxes 229 ! --------------------------------------------------- 230 ! fluxes 231 DO jj = 1, jpjm1 232 DO ji = 1, fs_jpim1 233 zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj) 234 zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1) 235 END DO 236 END DO 237 ! guess after content field with upstream scheme 238 DO jj = 2, jpjm1 239 DO ji = fs_2, fs_jpim1 240 ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj ) & 241 & + zfv_ups(ji,jj) - zfv_ups(ji ,jj-1) ) * r1_e1e2t(ji,jj) 242 zt_ups(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 243 END DO 244 END DO 245 CALL lbc_lnk( zt_ups, 'T', 1. ) 185 246 186 247 ! High order (_ho) fluxes 187 248 ! ----------------------- 188 SELECT CASE( k_order ) 189 CASE ( 20 ) ! centered second order 190 DO jj = 1, jpjm1 191 DO ji = 1, fs_jpim1 ! vector opt. 192 zfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( ptc(ji,jj) + ptc(ji+1,jj) ) 193 zfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj) + ptc(ji,jj+1) ) 194 END DO 195 END DO 196 ! 197 CASE ( 1:5 ) ! 1st to 5th order ULTIMATE-MACHO scheme 198 CALL macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v ) 199 ! 200 DO jj = 1, jpjm1 201 DO ji = 1, fs_jpim1 ! vector opt. 202 zfu_ho(ji,jj) = puc(ji,jj) * zt_u(ji,jj) 203 zfv_ho(ji,jj) = pvc(ji,jj) * zt_v(ji,jj) 204 END DO 205 END DO 249 SELECT CASE( kn_umx ) 250 ! 251 CASE ( 20 ) !== centered second order ==! 252 ! 253 CALL cen2( kn_limiter, jt, kt, pdt, zpt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho, & 254 & zt_ups, zfu_ups, zfv_ups ) 255 ! 256 CASE ( 1:5 ) !== 1st to 5th order ULTIMATE-MACHO scheme ==! 257 ! 258 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, & 259 & zt_ups, zfu_ups, zfv_ups ) 206 260 ! 207 261 END SELECT 208 262 209 ! antidiffusive flux : high order minus low order210 ! ------------------------ --------------------------211 DO jj = 1, jpjm1212 DO j i = 1, fs_jpim1 ! vector opt.213 zfu_ho(ji,jj) = zfu_ho(ji,jj) - zfu_ups(ji,jj)214 zfv_ho(ji,jj) = zfv_ho(ji,jj) - zfv_ups(ji,jj)215 END DO216 END DO217 218 ! monotonicity algorithm219 ! -------------------------220 CALL nonosc_2d( ptc, zfu_ho, zfv_ho, zt_ups, pdt )263 ! output high order fluxes 264 ! ------------------------ 265 IF( PRESENT(pfu) ) THEN 266 DO jj = 1, jpjm1 267 DO ji = 1, fs_jpim1 268 pfu(ji,jj) = zfu_ho(ji,jj) !!clem - zeps * puc(ji,jj) 269 pfv(ji,jj) = zfv_ho(ji,jj) !!clem - zeps * pvc(ji,jj) 270 END DO 271 END DO 272 !!CALL lbc_lnk( pfu, 'U', -1. ) ! clem: not needed I think 273 !!CALL lbc_lnk( pfv, 'V', -1. ) 274 ENDIF 221 275 222 276 ! final trend with corrected fluxes 223 277 ! ------------------------------------ 224 278 DO jj = 2, jpjm1 225 DO ji = fs_2, fs_jpim1 ! vector opt. 226 ztra = ztrd(ji,jj) - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj ) & 227 & + zfv_ho(ji,jj) - zfv_ho(ji ,jj-1) ) * r1_e1e2t(ji,jj) 279 DO ji = fs_2, fs_jpim1 280 ztra = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) & ! Div(uaH) or Div(ua) 281 !!clem & + ( puc (ji,jj) - puc (ji-1,jj) + pvc (ji,jj) - pvc (ji,jj-1) ) * zeps & ! epsi * Div(ua) or epsi * Div(u) 282 & ) * r1_e1e2t(ji,jj) 228 283 ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 229 284 END DO … … 233 288 END SUBROUTINE adv_umx 234 289 235 236 SUBROUTINE macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v ) 290 SUBROUTINE cen2( kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, & 291 & pt_ups, pfu_ups, pfv_ups ) 292 !!--------------------------------------------------------------------- 293 !! *** ROUTINE macho *** 294 !! 295 !! ** Purpose : compute 296 !! 297 !! ** Method : ... ??? 298 !! TIM = transient interpolation Modeling 299 !! 300 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 301 !!---------------------------------------------------------------------- 302 INTEGER , INTENT(in ) :: kn_limiter ! limiter 303 INTEGER , INTENT(in ) :: jt ! number of sub-iteration 304 INTEGER , INTENT(in ) :: kt ! number of iteration 305 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 306 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields 307 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu, pv ! 2 ice velocity components 308 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc, pvc ! 2 ice velocity * A components 309 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptc ! tracer content at before time step 310 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes 311 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt_ups ! upstream guess of tracer content 312 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pfu_ups, pfv_ups ! upstream fluxes 313 ! 314 INTEGER :: ji, jj ! dummy loop indices 315 LOGICAL :: ll_xy = .TRUE. 316 REAL(wp), DIMENSION(jpi,jpj) :: zzt 317 !!---------------------------------------------------------------------- 318 ! 319 IF( .NOT.ll_xy ) THEN !-- no alternate directions --! 320 ! 321 DO jj = 1, jpjm1 322 DO ji = 1, fs_jpim1 323 pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) ) 324 pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) ) 325 END DO 326 END DO 327 IF ( kn_limiter == 1 ) THEN 328 CALL nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 329 ELSEIF( kn_limiter == 2 ) THEN 330 CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 331 CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 332 ELSEIF( kn_limiter == 3 ) THEN 333 CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 334 CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 335 ENDIF 336 ! 337 ELSE !-- alternate directions --! 338 ! 339 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 340 ! 341 ! flux in x-direction 342 DO jj = 1, jpjm1 343 DO ji = 1, fs_jpim1 344 pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) ) 345 END DO 346 END DO 347 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 348 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 349 350 ! first guess of tracer content from u-flux 351 DO jj = 2, jpjm1 352 DO ji = fs_2, fs_jpim1 ! vector opt. 353 zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 354 END DO 355 END DO 356 CALL lbc_lnk( zzt, 'T', 1. ) 357 358 ! flux in y-direction 359 DO jj = 1, jpjm1 360 DO ji = 1, fs_jpim1 361 pfv_ho(ji,jj) = 0.5 * pv(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) ) 362 END DO 363 END DO 364 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 365 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 366 367 ELSE !== even ice time step: adv_y then adv_x ==! 368 ! 369 ! flux in y-direction 370 DO jj = 1, jpjm1 371 DO ji = 1, fs_jpim1 372 pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) ) 373 END DO 374 END DO 375 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 376 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 377 ! 378 ! first guess of tracer content from v-flux 379 DO jj = 2, jpjm1 380 DO ji = fs_2, fs_jpim1 ! vector opt. 381 zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 382 END DO 383 END DO 384 CALL lbc_lnk( zzt, 'T', 1. ) 385 ! 386 ! flux in x-direction 387 DO jj = 1, jpjm1 388 DO ji = 1, fs_jpim1 389 pfu_ho(ji,jj) = 0.5 * pu(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) ) 390 END DO 391 END DO 392 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 393 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 394 395 ENDIF 396 IF( kn_limiter == 1 ) CALL nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 397 398 ENDIF 399 400 END SUBROUTINE cen2 401 402 403 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, & 404 & pt_ups, pfu_ups, pfv_ups ) 405 !!--------------------------------------------------------------------- 406 !! *** ROUTINE macho *** 407 !! 408 !! ** Purpose : compute 409 !! 410 !! ** Method : ... ??? 411 !! TIM = transient interpolation Modeling 412 !! 413 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 414 !!---------------------------------------------------------------------- 415 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 416 INTEGER , INTENT(in ) :: kn_limiter ! limiter 417 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 418 INTEGER , INTENT(in ) :: jt ! number of sub-iteration 419 INTEGER , INTENT(in ) :: kt ! number of iteration 420 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 421 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields 422 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu, pv ! 2 ice velocity components 423 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc, pvc ! 2 ice velocity * A components 424 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity 425 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptc ! tracer content at before time step 426 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_u, pt_v ! tracer at u- and v-points 427 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes 428 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt_ups ! upstream guess of tracer content 429 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pfu_ups, pfv_ups ! upstream fluxes 430 ! 431 INTEGER :: ji, jj ! dummy loop indices 432 REAL(wp), DIMENSION(jpi,jpj) :: zzt 433 !!---------------------------------------------------------------------- 434 ! 435 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 436 ! 437 ! !-- ultimate interpolation of pt at u-point --! 438 CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 439 ! !-- limiter in x --! 440 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 441 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 442 ! !-- advective form update in zzt --! 443 DO jj = 2, jpjm1 444 DO ji = fs_2, fs_jpim1 ! vector opt. 445 zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj) & 446 & - pt (ji,jj) * pdt * ( pu (ji,jj) - pu (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk 447 zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 448 END DO 449 END DO 450 CALL lbc_lnk( zzt, 'T', 1. ) 451 ! !-- ultimate interpolation of pt at v-point --! 452 CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho ) 453 ! !-- limiter in y --! 454 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 455 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 456 ! 457 ELSE !== even ice time step: adv_y then adv_x ==! 458 ! 459 ! !-- ultimate interpolation of pt at v-point --! 460 CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 461 ! !-- limiter in y --! 462 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 463 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 464 ! !-- advective form update in zzt --! 465 DO jj = 2, jpjm1 466 DO ji = fs_2, fs_jpim1 467 zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj) & 468 & - pt (ji,jj) * pdt * ( pv (ji,jj) - pv (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk 469 zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 470 END DO 471 END DO 472 CALL lbc_lnk( zzt, 'T', 1. ) 473 ! !-- ultimate interpolation of pt at u-point --! 474 CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho ) 475 ! !-- limiter in x --! 476 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 477 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 478 ! 479 ENDIF 480 IF( kn_limiter == 1 ) CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 481 ! 482 END SUBROUTINE macho 483 484 485 SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 237 486 !!--------------------------------------------------------------------- 238 487 !! *** ROUTINE ultimate_x *** … … 245 494 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 246 495 !!---------------------------------------------------------------------- 247 INTEGER , INTENT(in ) :: k_order ! order of the ULTIMATE scheme 248 INTEGER , INTENT(in ) :: kt ! number of iteration 249 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 250 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptc ! tracer fields 251 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc, pvc ! 2 ice velocity components 252 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity 253 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_u, pt_v ! tracer at u- and v-points 254 ! 255 INTEGER :: ji, jj ! dummy loop indices 256 REAL(wp) :: zc_box ! - - 257 REAL(wp), DIMENSION(jpi,jpj) :: zzt 258 !!---------------------------------------------------------------------- 259 ! 260 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 261 ! 262 ! !-- ultimate interpolation of pt at u-point --! 263 CALL ultimate_x( k_order, pdt, ptc, puc, pt_u ) 264 ! 265 ! !-- advective form update in zzt --! 266 DO jj = 2, jpjm1 267 DO ji = fs_2, fs_jpim1 ! vector opt. 268 zzt(ji,jj) = ptc(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj) & 269 & - ptc (ji,jj) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e1e2t(ji,jj) 270 zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 271 END DO 272 END DO 273 CALL lbc_lnk( zzt, 'T', 1. ) 274 ! 275 ! !-- ultimate interpolation of pt at v-point --! 276 CALL ultimate_y( k_order, pdt, zzt, pvc, pt_v ) 277 ! 278 ELSE !== even ice time step: adv_y then adv_x ==! 279 ! 280 ! !-- ultimate interpolation of pt at v-point --! 281 CALL ultimate_y( k_order, pdt, ptc, pvc, pt_v ) 282 ! 283 ! !-- advective form update in zzt --! 284 DO jj = 2, jpjm1 285 DO ji = fs_2, fs_jpim1 286 zzt(ji,jj) = ptc(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj) & 287 & - ptc (ji,jj) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e1e2t(ji,jj) 288 zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 289 END DO 290 END DO 291 CALL lbc_lnk( zzt, 'T', 1. ) 292 ! 293 ! !-- ultimate interpolation of pt at u-point --! 294 CALL ultimate_x( k_order, pdt, zzt, puc, pt_u ) 295 ! 296 ENDIF 297 ! 298 END SUBROUTINE macho 299 300 301 SUBROUTINE ultimate_x( k_order, pdt, pt, puc, pt_u ) 302 !!--------------------------------------------------------------------- 303 !! *** ROUTINE ultimate_x *** 304 !! 305 !! ** Purpose : compute 306 !! 307 !! ** Method : ... ??? 308 !! TIM = transient interpolation Modeling 309 !! 310 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 311 !!---------------------------------------------------------------------- 312 INTEGER , INTENT(in ) :: k_order ! ocean time-step index 496 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 313 497 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 314 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc ! ice i-velocity component 498 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu ! ice i-velocity component 499 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc ! ice i-velocity * A component 315 500 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields 316 501 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_u ! tracer at u-point 317 ! 318 INTEGER :: ji, jj ! dummy loop indices 502 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pfu_ho ! high order flux 503 ! 504 INTEGER :: ji, jj ! dummy loop indices 319 505 REAL(wp) :: zcu, zdx2, zdx4 ! - - 320 REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu2, ztu3, ztu4506 REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu2, ztu3, ztu4 321 507 !!---------------------------------------------------------------------- 322 508 ! … … 346 532 ! 347 533 ! 348 SELECT CASE (k _order)534 SELECT CASE (kn_umx ) 349 535 ! 350 536 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) … … 352 538 DO jj = 2, jpjm1 353 539 DO ji = 1, fs_jpim1 ! vector opt. 354 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( 355 & - SIGN( 1._wp, pu c(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) )540 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) & 541 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) ) 356 542 END DO 357 543 END DO … … 361 547 DO jj = 2, jpjm1 362 548 DO ji = 1, fs_jpim1 ! vector opt. 363 zcu = pu c(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)549 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 364 550 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) & 365 551 & - zcu * ( pt(ji+1,jj) - pt(ji,jj) ) ) … … 371 557 DO jj = 2, jpjm1 372 558 DO ji = 1, fs_jpim1 ! vector opt. 373 zcu = pu c(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)559 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 374 560 zdx2 = e1u(ji,jj) * e1u(ji,jj) 375 561 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) … … 385 571 DO jj = 2, jpjm1 386 572 DO ji = 1, fs_jpim1 ! vector opt. 387 zcu = pu c(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)573 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 388 574 zdx2 = e1u(ji,jj) * e1u(ji,jj) 389 575 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) … … 399 585 DO jj = 2, jpjm1 400 586 DO ji = 1, fs_jpim1 ! vector opt. 401 zcu = pu c(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)587 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 402 588 zdx2 = e1u(ji,jj) * e1u(ji,jj) 403 589 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) … … 413 599 ! 414 600 END SELECT 601 ! !-- High order flux in i-direction --! 602 DO jj = 1, jpjm1 603 DO ji = 1, fs_jpim1 ! vector opt. 604 pfu_ho(ji,jj) = puc(ji,jj) * pt_u(ji,jj) 605 END DO 606 END DO 415 607 ! 416 608 END SUBROUTINE ultimate_x 417 609 418 610 419 SUBROUTINE ultimate_y( k _order, pdt, pt, pvc, pt_v)611 SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 420 612 !!--------------------------------------------------------------------- 421 613 !! *** ROUTINE ultimate_y *** … … 428 620 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 429 621 !!---------------------------------------------------------------------- 430 INTEGER , INTENT(in ) :: k _order ! ocean time-step index622 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 431 623 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 432 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pvc ! ice j-velocity component 624 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pv ! ice j-velocity component 625 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pvc ! ice j-velocity*A component 433 626 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields 434 627 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_v ! tracer at v-point 628 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pfv_ho ! high order flux 435 629 ! 436 630 INTEGER :: ji, jj ! dummy loop indices 437 631 REAL(wp) :: zcv, zdy2, zdy4 ! - - 438 REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv2, ztv3, ztv4632 REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv2, ztv3, ztv4 439 633 !!---------------------------------------------------------------------- 440 634 ! … … 466 660 ! 467 661 ! 468 SELECT CASE (k _order)662 SELECT CASE (kn_umx ) 469 663 ! 470 664 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 471 665 DO jj = 1, jpjm1 472 666 DO ji = fs_2, fs_jpim1 473 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( 474 & - SIGN( 1._wp, pv c(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )667 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) & 668 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) 475 669 END DO 476 670 END DO … … 479 673 DO jj = 1, jpjm1 480 674 DO ji = fs_2, fs_jpim1 481 zcv = pv c(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)675 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 482 676 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) & 483 677 & - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) ) … … 489 683 DO jj = 1, jpjm1 490 684 DO ji = fs_2, fs_jpim1 491 zcv = pv c(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)685 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 492 686 zdy2 = e2v(ji,jj) * e2v(ji,jj) 493 687 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) … … 502 696 DO jj = 1, jpjm1 503 697 DO ji = fs_2, fs_jpim1 504 zcv = pv c(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)698 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 505 699 zdy2 = e2v(ji,jj) * e2v(ji,jj) 506 700 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) … … 515 709 DO jj = 1, jpjm1 516 710 DO ji = fs_2, fs_jpim1 517 zcv = pv c(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)711 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 518 712 zdy2 = e2v(ji,jj) * e2v(ji,jj) 519 713 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) … … 529 723 ! 530 724 END SELECT 725 ! !-- High order flux in j-direction --! 726 DO jj = 1, jpjm1 727 DO ji = 1, fs_jpim1 ! vector opt. 728 pfv_ho(ji,jj) = pvc(ji,jj) * pt_v(ji,jj) 729 END DO 730 END DO 531 731 ! 532 732 END SUBROUTINE ultimate_y 533 534 535 SUBROUTINE nonosc_2d( p bef, paa, pbb, paft, pdt)733 734 735 SUBROUTINE nonosc_2d( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 536 736 !!--------------------------------------------------------------------- 537 737 !! *** ROUTINE nonosc *** … … 541 741 !! 542 742 !! ** Method : ... ??? 543 !! warning : p bef and paftmust be masked, but the boundaries743 !! warning : ptc and pt_ups must be masked, but the boundaries 544 744 !! conditions on the fluxes are not necessary zalezak (1979) 545 745 !! drange (1995) multi-dimensional forward-in-time and upstream- 546 746 !! in-space based differencing for fluid 547 747 !!---------------------------------------------------------------------- 548 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 549 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pbef, paft ! before & after field 550 REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: paa, pbb ! monotonic fluxes in the 2 directions 748 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 749 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: ptc, pt_ups ! before field & upstream guess of after field 750 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pfv_ups, pfu_ups ! upstream flux 751 REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: pfv_ho, pfu_ho ! monotonic flux 551 752 ! 552 753 INTEGER :: ji, jj ! dummy loop indices 553 INTEGER :: ikm1 ! local integer 554 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zsml, z1_dt ! local scalars 555 REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - 556 REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zmsk, zdiv 557 !!---------------------------------------------------------------------- 558 ! 754 REAL(wp) :: zpos, zneg, zbig, zsml, z1_dt ! local scalars 755 REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - 756 REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zdiv 757 !!---------------------------------------------------------------------- 559 758 zbig = 1.e+40_wp 560 759 zsml = 1.e-15_wp … … 563 762 DO jj = 2, jpjm1 564 763 DO ji = fs_2, fs_jpim1 ! vector opt. 565 zdiv(ji,jj) = - ( paa(ji,jj) - paa(ji-1,jj ) & 566 & + pbb(ji,jj) - pbb(ji ,jj-1) ) 764 zdiv(ji,jj) = - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) + pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) 567 765 END DO 568 766 END DO 569 767 CALL lbc_lnk( zdiv, 'T', 1. ) ! Lateral boundary conditions (unchanged sign) 570 768 571 ! Determine ice masks for before and after tracers 572 WHERE( pbef(:,:) == 0._wp .AND. paft(:,:) == 0._wp .AND. zdiv(:,:) == 0._wp ) ; zmsk(:,:) = 0._wp 573 ELSEWHERE ; zmsk(:,:) = 1._wp * tmask(:,:,1) 574 END WHERE 769 ! antidiffusive flux : high order minus low order 770 ! -------------------------------------------------- 771 DO jj = 1, jpjm1 772 DO ji = 1, fs_jpim1 ! vector opt. 773 pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_ups(ji,jj) 774 pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_ups(ji,jj) 775 END DO 776 END DO 575 777 576 778 ! Search local extrema 577 779 ! -------------------- 578 ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 579 ! zbup(:,:) = MAX( pbef(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ), & 580 ! & paft(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ) ) 581 ! zbdo(:,:) = MIN( pbef(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ), & 582 ! & paft(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ) ) 583 zbup(:,:) = MAX( pbef(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ), & 584 & paft(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ) ) 585 zbdo(:,:) = MIN( pbef(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ), & 586 & paft(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ) ) 780 ! max/min of ptc & pt_ups with large negative/positive value (-/+zbig) outside ice cover 781 DO jj = 1, jpj 782 DO ji = fs_2, fs_jpim1 783 IF( ptc(ji,jj) == 0._wp .AND. pt_ups(ji,jj) == 0._wp .AND. zdiv(ji,jj) == 0._wp ) THEN 784 zbup(ji,jj) = -zbig 785 zbdo(ji,jj) = zbig 786 ELSE 787 zbup(ji,jj) = MAX( ptc(ji,jj) , pt_ups(ji,jj) ) 788 zbdo(ji,jj) = MIN( ptc(ji,jj) , pt_ups(ji,jj) ) 789 ENDIF 790 END DO 791 END DO 587 792 588 793 z1_dt = 1._wp / pdt … … 590 795 DO ji = fs_2, fs_jpim1 ! vector opt. 591 796 ! 592 zup = MAX( zbup(ji,jj), zbup(ji-1,jj ), zbup(ji+1,jj ), & ! search max/min in neighbourhood 593 & zbup(ji ,jj-1), zbup(ji ,jj+1) ) 594 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), & 595 & zbdo(ji ,jj-1), zbdo(ji ,jj+1) ) 596 ! 597 zpos = MAX( 0., paa(ji-1,jj ) ) - MIN( 0., paa(ji ,jj ) ) & ! positive/negative part of the flux 598 & + MAX( 0., pbb(ji ,jj-1) ) - MIN( 0., pbb(ji ,jj ) ) 599 zneg = MAX( 0., paa(ji ,jj ) ) - MIN( 0., paa(ji-1,jj ) ) & 600 & + MAX( 0., pbb(ji ,jj ) ) - MIN( 0., pbb(ji ,jj-1) ) 601 ! 602 zbt = e1e2t(ji,jj) * z1_dt ! up & down beta terms 603 zbetup(ji,jj) = ( zup - paft(ji,jj) ) / ( zpos + zsml ) * zbt 604 zbetdo(ji,jj) = ( paft(ji,jj) - zdo ) / ( zneg + zsml ) * zbt 797 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 798 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) ) 799 ! 800 zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji ,jj) ) & 801 & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj ) ) ! positive/negative part of the flux 802 zneg = MAX( 0., pfu_ho(ji ,jj) ) - MIN( 0., pfu_ho(ji-1,jj) ) & 803 & + MAX( 0., pfv_ho(ji,jj ) ) - MIN( 0., pfv_ho(ji,jj-1) ) 804 ! 805 ! ! up & down beta terms 806 !!clem zbetup(ji,jj) = ( zup - pt_ups(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt 807 !!clem zbetdo(ji,jj) = ( pt_ups(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt 808 IF( zpos >= epsi20 ) THEN 809 zbetup(ji,jj) = ( zup - pt_ups(ji,jj) ) / zpos * e1e2t(ji,jj) * z1_dt 810 ELSE 811 zbetup(ji,jj) = zbig 812 ENDIF 813 ! 814 IF( zneg >= epsi20 ) THEN 815 zbetdo(ji,jj) = ( pt_ups(ji,jj) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 816 ELSE 817 zbetdo(ji,jj) = zbig 818 ENDIF 819 ! 605 820 END DO 606 821 END DO 607 822 CALL lbc_lnk_multi( zbetup, 'T', 1., zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign) 608 823 609 ! monotonic flux in the i & j direction (paa & pbb)610 ! --------------------------------- ----611 DO jj = 2, jpjm1824 ! monotonic flux in the y direction 825 ! --------------------------------- 826 DO jj = 1, jpjm1 612 827 DO ji = 1, fs_jpim1 ! vector opt. 613 828 zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) ) 614 829 zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) ) 615 zcu = 0.5 + SIGN( 0.5 , p aa(ji,jj) )616 ! 617 p aa(ji,jj) = paa(ji,jj) * ( zcu * zau + ( 1._wp - zcu) * zbu)618 END DO 619 END DO 620 ! 830 zcu = 0.5 + SIGN( 0.5 , pfu_ho(ji,jj) ) 831 ! 832 pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_ups(ji,jj) 833 END DO 834 END DO 835 621 836 DO jj = 1, jpjm1 622 DO ji = fs_2, fs_jpim1 ! vector opt.837 DO ji = 1, fs_jpim1 ! vector opt. 623 838 zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) ) 624 839 zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) ) 625 zcv = 0.5 + SIGN( 0.5 , p bb(ji,jj) )626 ! 627 p bb(ji,jj) = pbb(ji,jj) * ( zcv * zav + ( 1._wp - zcv) * zbv)840 zcv = 0.5 + SIGN( 0.5 , pfv_ho(ji,jj) ) 841 ! 842 pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_ups(ji,jj) 628 843 END DO 629 844 END DO 630 845 ! 631 846 END SUBROUTINE nonosc_2d 847 848 SUBROUTINE limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 849 !!--------------------------------------------------------------------- 850 !! *** ROUTINE limiter_x *** 851 !! 852 !! ** Purpose : compute flux limiter 853 !!---------------------------------------------------------------------- 854 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 855 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pu ! ice i-velocity => u*e2 856 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: puc ! ice i-velocity *A => u*e2*a 857 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pt ! ice tracer 858 REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: pfu_ho ! high order flux 859 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ), OPTIONAL :: pfu_ups ! upstream flux 860 ! 861 REAL(wp) :: Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr 862 INTEGER :: ji, jj ! dummy loop indices 863 REAL(wp), DIMENSION (jpi,jpj) :: zslpx ! tracer slopes 864 !!---------------------------------------------------------------------- 865 ! 866 DO jj = 2, jpjm1 867 DO ji = fs_2, fs_jpim1 ! vector opt. 868 zslpx(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * umask(ji,jj,1) 869 END DO 870 END DO 871 CALL lbc_lnk( zslpx, 'U', -1.) ! lateral boundary cond. 872 873 DO jj = 2, jpjm1 874 DO ji = fs_2, fs_jpim1 ! vector opt. 875 uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj) 876 877 Rjm = zslpx(ji-1,jj) 878 Rj = zslpx(ji ,jj) 879 Rjp = zslpx(ji+1,jj) 880 881 IF( PRESENT(pfu_ups) ) THEN 882 883 IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm 884 ELSE ; Rr = Rjp 885 ENDIF 886 887 zh3 = pfu_ho(ji,jj) - pfu_ups(ji,jj) 888 IF( Rj > 0. ) THEN 889 zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(puc(ji,jj)), & 890 & MIN( 2. * Rr * 0.5 * ABS(puc(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) ) 891 ELSE 892 zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(puc(ji,jj)), & 893 & MIN(-2. * Rr * 0.5 * ABS(puc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) ) 894 ENDIF 895 pfu_ho(ji,jj) = pfu_ups(ji,jj) + zlimiter 896 897 ELSE 898 IF( Rj /= 0. ) THEN 899 IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj 900 ELSE ; Cr = Rjp / Rj 901 ENDIF 902 ELSE 903 Cr = 0. 904 !IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm * 1.e20 905 !ELSE ; Cr = Rjp * 1.e20 906 !ENDIF 907 ENDIF 908 909 ! -- superbee -- 910 zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 911 ! -- van albada 2 -- 912 !!zpsi = 2.*Cr / (Cr*Cr+1.) 913 914 ! -- sweby (with beta=1) -- 915 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 916 ! -- van Leer -- 917 !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 918 ! -- ospre -- 919 !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 920 ! -- koren -- 921 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 922 ! -- charm -- 923 !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 924 !ELSE ; zpsi = 0. 925 !ENDIF 926 ! -- van albada 1 -- 927 !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 928 ! -- smart -- 929 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 930 ! -- umist -- 931 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 932 933 ! high order flux corrected by the limiter 934 pfu_ho(ji,jj) = pfu_ho(ji,jj) - ABS( puc(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5 935 936 ENDIF 937 END DO 938 END DO 939 CALL lbc_lnk( pfu_ho, 'U', -1.) ! lateral boundary cond. 940 ! 941 END SUBROUTINE limiter_x 942 943 SUBROUTINE limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 944 !!--------------------------------------------------------------------- 945 !! *** ROUTINE limiter_y *** 946 !! 947 !! ** Purpose : compute flux limiter 948 !!---------------------------------------------------------------------- 949 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 950 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pv ! ice i-velocity => u*e2 951 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pvc ! ice i-velocity *A => u*e2*a 952 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pt ! ice tracer 953 REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: pfv_ho ! high order flux 954 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ), OPTIONAL :: pfv_ups ! upstream flux 955 ! 956 REAL(wp) :: Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr 957 INTEGER :: ji, jj ! dummy loop indices 958 REAL(wp), DIMENSION (jpi,jpj) :: zslpy ! tracer slopes 959 !!---------------------------------------------------------------------- 960 ! 961 DO jj = 2, jpjm1 962 DO ji = fs_2, fs_jpim1 ! vector opt. 963 zslpy(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * vmask(ji,jj,1) 964 END DO 965 END DO 966 CALL lbc_lnk( zslpy, 'V', -1.) ! lateral boundary cond. 967 968 DO jj = 2, jpjm1 969 DO ji = fs_2, fs_jpim1 ! vector opt. 970 vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj) 971 972 Rjm = zslpy(ji,jj-1) 973 Rj = zslpy(ji,jj ) 974 Rjp = zslpy(ji,jj+1) 975 976 IF( PRESENT(pfv_ups) ) THEN 977 978 IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm 979 ELSE ; Rr = Rjp 980 ENDIF 981 982 zh3 = pfv_ho(ji,jj) - pfv_ups(ji,jj) 983 IF( Rj > 0. ) THEN 984 zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pvc(ji,jj)), & 985 & MIN( 2. * Rr * 0.5 * ABS(pvc(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) ) 986 ELSE 987 zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pvc(ji,jj)), & 988 & MIN(-2. * Rr * 0.5 * ABS(pvc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) ) 989 ENDIF 990 pfv_ho(ji,jj) = pfv_ups(ji,jj) + zlimiter 991 992 ELSE 993 994 IF( Rj /= 0. ) THEN 995 IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj 996 ELSE ; Cr = Rjp / Rj 997 ENDIF 998 ELSE 999 Cr = 0. 1000 !IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm * 1.e20 1001 !ELSE ; Cr = Rjp * 1.e20 1002 !ENDIF 1003 ENDIF 1004 1005 ! -- superbee -- 1006 zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 1007 ! -- van albada 2 -- 1008 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1009 1010 ! -- sweby (with beta=1) -- 1011 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 1012 ! -- van Leer -- 1013 !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 1014 ! -- ospre -- 1015 !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 1016 ! -- koren -- 1017 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 1018 ! -- charm -- 1019 !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 1020 !ELSE ; zpsi = 0. 1021 !ENDIF 1022 ! -- van albada 1 -- 1023 !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 1024 ! -- smart -- 1025 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 1026 ! -- umist -- 1027 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 1028 1029 ! high order flux corrected by the limiter 1030 pfv_ho(ji,jj) = pfv_ho(ji,jj) - ABS( pvc(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5 1031 1032 ENDIF 1033 END DO 1034 END DO 1035 CALL lbc_lnk( pfv_ho, 'V', -1.) ! lateral boundary cond. 1036 ! 1037 END SUBROUTINE limiter_y 632 1038 633 1039 #else
Note: See TracChangeset
for help on using the changeset viewer.
