Changeset 10413 for NEMO/trunk/src/ICE/icedyn_adv_umx.F90
- Timestamp:
- 2018-12-18T18:59:59+01:00 (6 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/ICE/icedyn_adv_umx.F90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icedyn_adv_umx.F90
r10344 r10413 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 … … 34 35 REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120 35 36 37 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z1_ai, amaxu, amaxv 38 39 LOGICAL ll_dens 40 41 ! advect H all the way (and get V=H*A at the end) 42 LOGICAL :: ll_thickness = .FALSE. 43 44 ! look for 9 points around in nonosc limiter 45 LOGICAL :: ll_9points = .FALSE. ! false=better h? 46 47 ! use HgradU in nonosc limiter 48 LOGICAL :: ll_HgradU = .TRUE. ! no effect? 49 50 ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme 51 LOGICAL :: ll_neg = .TRUE. ! keep TRUE 52 53 ! limit the fluxes 54 LOGICAL :: ll_zeroup1 = .FALSE. ! false ok if Hbig otherwise needed for 2D sinon on a des valeurs de H trop fortes !! 55 LOGICAL :: ll_zeroup2 = .FALSE. ! false ok for 1D, 2D, 3D 56 LOGICAL :: ll_zeroup4 = .FALSE. ! false ok for 1D, 2D, 3D 57 LOGICAL :: ll_zeroup5 = .FALSE. ! false ok for 1D, 2D 58 59 ! fluxes that are limited are u*H, then (u*H)*(ua/u) is used for V (only for nonosc) 60 LOGICAL :: ll_clem = .TRUE. ! simpler than rachid and works 61 62 ! First advect H as H*=Hdiv(u), then use H* for H(n+1)=H(n)-div(uH*) 63 LOGICAL :: ll_gurvan = .FALSE. ! must be false for 1D case !! 64 65 ! First guess as div(uH) (-true-) or Hdiv(u)+ugradH (-false-) 66 LOGICAL :: ll_1stguess_clem = .FALSE. ! better negative values but less good h 67 68 ! advect (or not) open water. If not, retrieve it from advection of A 69 LOGICAL :: ll_ADVopw = .FALSE. ! 70 71 ! alternate directions for upstream 72 LOGICAL :: ll_upsxy = .TRUE. 73 74 ! alternate directions for high order 75 LOGICAL :: ll_hoxy = .TRUE. 76 77 ! prelimiter: use it to avoid overshoot in H 78 LOGICAL :: ll_prelimiter_zalesak = .TRUE. ! from: Zalesak(1979) eq. 14 => true is better for 1D but false is better in 3D (for h and negative values) => pb in x-y? 79 LOGICAL :: ll_prelimiter_devore = .FALSE. ! from: Devore eq. 11 (worth than zalesak) 80 81 ! iterate on the limiter (only nonosc) 82 LOGICAL :: ll_limiter_it2 = .FALSE. 83 84 36 85 !! * Substitutions 37 86 # include "vectopt_loop_substitute.h90" … … 39 88 !! NEMO/ICE 4.0 , NEMO Consortium (2018) 40 89 !! $Id$ 41 !! Software governed by the CeCILL licen se (see./LICENSE)90 !! Software governed by the CeCILL licence (./LICENSE) 42 91 !!---------------------------------------------------------------------- 43 92 CONTAINS 44 93 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 )94 SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, & 95 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 47 96 !!---------------------------------------------------------------------- 48 97 !! *** ROUTINE ice_dyn_adv_umx *** … … 54 103 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 55 104 !!---------------------------------------------------------------------- 56 INTEGER , INTENT(in ) :: k _order ! order of the scheme (1-5 or 20)105 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 57 106 INTEGER , INTENT(in ) :: kt ! time step 58 107 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity … … 70 119 ! 71 120 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 121 INTEGER :: icycle ! number of sub-timestep for the advection 122 REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers 123 REAL(wp) :: zcfl , zdt 124 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box, zua_ho, zva_ho 125 REAL(wp), DIMENSION(jpi,jpj) :: zhvar 126 REAL(wp), DIMENSION(jpi,jpj) :: zai_b, zai_a, z1_ai_b 75 127 !!---------------------------------------------------------------------- 76 128 ! 77 129 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 78 130 ! 79 ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) )80 131 ! 81 132 ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! … … 84 135 IF( lk_mpp ) CALL mpp_max( zcfl ) 85 136 86 IF( zcfl > 0.5 ) THEN ; i nitad = 2 ; zusnit = 0.5_wp87 ELSE ; i nitad = 1 ; zusnit = 1.0_wp88 ENDIF 89 90 zdt = rdt_ice / REAL(i nitad)137 IF( zcfl > 0.5 ) THEN ; icycle = 2 138 ELSE ; icycle = 1 139 ENDIF 140 141 zdt = rdt_ice / REAL(icycle) 91 142 92 143 ! --- transport --- ! … … 109 160 END DO 110 161 162 IF(.NOT. ALLOCATED(z1_ai)) ALLOCATE(z1_ai(jpi,jpj)) 163 IF( ll_zeroup2 ) THEN 164 IF(.NOT. ALLOCATED(amaxu)) ALLOCATE(amaxu (jpi,jpj)) 165 IF(.NOT. ALLOCATED(amaxv)) ALLOCATE(amaxv (jpi,jpj)) 166 ENDIF 111 167 !---------------! 112 168 !== advection ==! 113 169 !---------------! 114 DO jt = 1, initad 115 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) ) ! Open water area 170 DO jt = 1, icycle 171 172 IF( ll_ADVopw ) THEN 173 ll_dens=.FALSE. 174 zamsk = 1._wp 175 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) ) ! Open water area 176 ELSE 177 zai_b(:,:) = SUM( pa_i(:,:,:), dim=3 ) 178 ENDIF 179 116 180 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 181 ! 182 WHERE( pa_i(:,:,jl) >= epsi20 ) ; z1_ai_b(:,:) = 1._wp / pa_i(:,:,jl) 183 ELSEWHERE ; z1_ai_b(:,:) = 0. 184 END WHERE 185 ! 186 IF( ll_zeroup2 ) THEN 187 DO jj = 2, jpjm1 188 DO ji = fs_2, fs_jpim1 189 amaxu(ji,jj)=MAX( pa_i(ji,jj,jl), pa_i(ji,jj-1,jl), pa_i(ji,jj+1,jl), & 190 & pa_i(ji+1,jj,jl), pa_i(ji+1,jj-1,jl), pa_i(ji+1,jj+1,jl) ) 191 amaxv(ji,jj)=MAX( pa_i(ji,jj,jl), pa_i(ji-1,jj,jl), pa_i(ji+1,jj,jl), & 192 & pa_i(ji,jj+1,jl), pa_i(ji-1,jj+1,jl), pa_i(ji+1,jj+1,jl) ) 193 END DO 194 END DO 195 CALL lbc_lnk_multi(amaxu, 'T', 1., amaxv, 'T', 1.) 196 ENDIF 197 ! 198 zamsk = 1._wp 199 ll_dens=.TRUE. 200 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), zua_ho, zva_ho ) ! Ice area 201 ll_dens=.FALSE. 202 203 WHERE( pa_i(:,:,jl) >= epsi20 ) ; z1_ai(:,:) = 1._wp / pa_i(:,:,jl) 204 ELSEWHERE ; z1_ai(:,:) = 0. 205 END WHERE 206 207 IF( ll_thickness ) THEN 208 zua_ho(:,:) = zudy(:,:) 209 zva_ho(:,:) = zvdx(:,:) 210 ENDIF 211 212 zamsk = 0._wp ; zhvar(:,:) = pv_i (:,:,jl) * z1_ai_b(:,:) 213 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_i (:,:,jl) ) ! Ice volume 214 IF( ll_thickness ) pv_i(:,:,jl) = zhvar(:,:) * pa_i(:,:,jl) 215 216 zamsk = 0._wp ; zhvar(:,:) = pv_s (:,:,jl) * z1_ai_b(:,:) 217 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_s (:,:,jl) ) ! Snw volume 218 IF( ll_thickness ) pv_s(:,:,jl) = zhvar(:,:) * pa_i(:,:,jl) 219 220 zamsk = 0._wp ; zhvar(:,:) = psv_i(:,:,jl) * z1_ai_b(:,:) 221 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), psv_i(:,:,jl) ) ! Salt content 222 223 zamsk = 0._wp ; zhvar(:,:) = poa_i(:,:,jl) * z1_ai_b(:,:) 224 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), poa_i(:,:,jl) ) ! Age content 225 121 226 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 content 123 END DO 124 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) ) ! Snow volume 227 zamsk = 0._wp ; zhvar(:,:) = pe_i(:,:,jk,jl) * z1_ai_b(:,:) 228 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar(:,:), pe_i(:,:,jk,jl) ) ! Ice heat content 229 END DO 230 125 231 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 131 ENDIF 132 END DO 133 END DO 134 ! 135 DEALLOCATE( zudy, zvdx, zcu_box, zcv_box ) 232 zamsk = 0._wp ; zhvar(:,:) = pe_s(:,:,jk,jl) * z1_ai_b(:,:) 233 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar(:,:), pe_s(:,:,jk,jl) ) ! Snw heat content 234 END DO 235 ! 236 IF ( ln_pnd_H12 ) THEN ! melt ponds (must be the last ones to be advected because of z1_ai_b...) 237 ! 238 WHERE( pa_ip(:,:,jl) >= epsi20 ) ; z1_ai_b(:,:) = 1._wp / pa_ip(:,:,jl) 239 ELSEWHERE ; z1_ai_b(:,:) = 0. 240 END WHERE 241 ! 242 zamsk = 1._wp 243 ll_dens=.TRUE. 244 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), zua_ho, zva_ho ) ! mp fractio!n 245 ll_dens=.FALSE. 246 247 WHERE( pa_ip(:,:,jl) >= epsi20 ) ; z1_ai(:,:) = 1._wp / pa_ip(:,:,jl) 248 ELSEWHERE ; z1_ai(:,:) = 0. 249 END WHERE 250 251 zamsk = 0._wp ; zhvar(:,:) = pv_ip(:,:,jl) * z1_ai_b(:,:) 252 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar(:,:), pv_ip(:,:,jl) ) ! mp volume 253 ENDIF 254 ! 255 ! 256 END DO 257 258 IF( .NOT. ll_ADVopw ) THEN 259 zai_a(:,:) = SUM( pa_i(:,:,:), dim=3 ) 260 DO jj = 2, jpjm1 261 DO ji = fs_2, fs_jpim1 262 pato_i(ji,jj) = pato_i(ji,jj) - ( zai_a(ji,jj) - zai_b(ji,jj) ) & 263 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 264 END DO 265 END DO 266 CALL lbc_lnk( pato_i(:,:), 'T', 1. ) 267 ENDIF 268 269 END DO 136 270 ! 137 271 END SUBROUTINE ice_dyn_adv_umx 138 272 139 273 140 SUBROUTINE adv_umx( k_order, kt, pdt, puc, pvc, pubox, pvbox, ptc)274 SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho ) 141 275 !!---------------------------------------------------------------------- 142 276 !! *** ROUTINE adv_umx *** … … 151 285 !! ** Action : - pt the after advective tracer 152 286 !!---------------------------------------------------------------------- 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 287 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 288 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 289 INTEGER , INTENT(in ) :: jt ! number of sub-iteration 290 INTEGER , INTENT(in ) :: kt ! number of iteration 291 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 292 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2 293 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 or u*a*e2u 294 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity 295 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pt ! tracer field 296 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: ptc ! tracer content field 297 REAL(wp), DIMENSION(jpi,jpj), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes 159 298 ! 160 299 INTEGER :: ji, jj ! dummy loop indices 161 300 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 301 INTEGER :: kn_limiter = 1 ! 1=nonosc ; 2=superbee ; 3=h3(rachid) 302 REAL(wp), DIMENSION(jpi,jpj) :: zfu_ho , zfv_ho , zt_u, zt_v, zpt 303 REAL(wp), DIMENSION(jpi,jpj) :: zfu_ups, zfv_ups, zt_ups ! only for nonosc 304 !!---------------------------------------------------------------------- 305 ! 306 ! upstream (_ups) advection with initial mass fluxes 307 ! --------------------------------------------------- 308 IF( ll_clem ) zfu_ups=0.; zfv_ups=0. 309 310 IF( ll_gurvan .AND. pamsk==0. ) THEN 311 DO jj = 2, jpjm1 312 DO ji = fs_2, fs_jpim1 313 pt(ji,jj) = ( pt (ji,jj) + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) & 314 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 315 END DO 316 END DO 317 CALL lbc_lnk( pt, 'T', 1. ) 318 ENDIF 319 174 320 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) 185 321 IF( .NOT. ll_upsxy ) THEN 322 323 ! fluxes 324 DO jj = 1, jpjm1 325 DO ji = 1, fs_jpim1 326 IF( ll_clem ) THEN 327 zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj) 328 zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1) 329 ELSE 330 zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj) 331 zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1) 332 ENDIF 333 END DO 334 END DO 335 336 ELSE 337 ! 1 if advection of A 338 ! z1_ai already defined IF advection of other tracers 339 IF( pamsk == 1. ) z1_ai(:,:) = 1._wp 340 ! 341 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 342 ! flux in x-direction 343 DO jj = 1, jpjm1 344 DO ji = 1, fs_jpim1 345 IF( ll_clem ) THEN 346 zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj) 347 ELSE 348 zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * pt(ji+1,jj) 349 ENDIF 350 END DO 351 END DO 352 353 ! first guess of tracer content from u-flux 354 DO jj = 2, jpjm1 355 DO ji = fs_2, fs_jpim1 ! vector opt. 356 IF( ll_clem ) THEN 357 IF( ll_gurvan ) THEN 358 zpt(ji,jj) = ( pt(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 359 ELSE 360 zpt(ji,jj) = ( pt(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 361 & + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 362 & * tmask(ji,jj,1) 363 ENDIF 364 ELSE 365 zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) & 366 & * tmask(ji,jj,1) 367 ENDIF 368 !! IF( ji==26 .AND. jj==86) THEN 369 !! WRITE(numout,*) '************************' 370 !! WRITE(numout,*) 'zpt upstream',zpt(ji,jj) 371 !! ENDIF 372 END DO 373 END DO 374 CALL lbc_lnk( zpt, 'T', 1. ) 375 ! 376 ! flux in y-direction 377 DO jj = 1, jpjm1 378 DO ji = 1, fs_jpim1 379 IF( ll_clem ) THEN 380 zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1) 381 ELSE 382 zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * zpt(ji,jj+1) 383 ENDIF 384 END DO 385 END DO 386 387 ! 388 ELSE !== even ice time step: adv_y then adv_x ==! 389 ! flux in y-direction 390 DO jj = 1, jpjm1 391 DO ji = 1, fs_jpim1 392 IF( ll_clem ) THEN 393 zfv_ups(ji,jj) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1) 394 ELSE 395 zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * pt(ji,jj+1) 396 ENDIF 397 END DO 398 END DO 399 400 ! first guess of tracer content from v-flux 401 DO jj = 2, jpjm1 402 DO ji = fs_2, fs_jpim1 ! vector opt. 403 IF( ll_clem ) THEN 404 IF( ll_gurvan ) THEN 405 zpt(ji,jj) = ( pt(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 406 ELSE 407 zpt(ji,jj) = ( pt(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) & 408 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 409 & * tmask(ji,jj,1) 410 ENDIF 411 ELSE 412 zpt(ji,jj) = ( ptc(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) & 413 & * tmask(ji,jj,1) 414 ENDIF 415 !! IF( ji==26 .AND. jj==86) THEN 416 !! WRITE(numout,*) '************************' 417 !! WRITE(numout,*) 'zpt upstream',zpt(ji,jj) 418 !! ENDIF 419 END DO 420 END DO 421 CALL lbc_lnk( zpt, 'T', 1. ) 422 ! 423 ! flux in x-direction 424 DO jj = 1, jpjm1 425 DO ji = 1, fs_jpim1 426 IF( ll_clem ) THEN 427 zfu_ups(ji,jj) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj) 428 ELSE 429 zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( puc(ji,jj), 0._wp ) * zpt(ji+1,jj) 430 ENDIF 431 END DO 432 END DO 433 ! 434 ENDIF 435 436 ENDIF 437 438 IF( ll_clem .AND. kn_limiter /= 1 ) & 439 & CALL ctl_stop( 'STOP', 'icedyn_adv_umx: ll_clem incompatible with limiters other than nonosc' ) 440 441 IF( ll_zeroup2 ) THEN 442 DO jj = 1, jpjm1 443 DO ji = 1, fs_jpim1 ! vector opt. 444 IF( amaxu(ji,jj) == 0._wp ) zfu_ups(ji,jj) = 0._wp 445 IF( amaxv(ji,jj) == 0._wp ) zfv_ups(ji,jj) = 0._wp 446 END DO 447 END DO 448 ENDIF 449 450 ! guess after content field with upstream scheme 451 DO jj = 2, jpjm1 452 DO ji = fs_2, fs_jpim1 453 ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj ) & 454 & + zfv_ups(ji,jj) - zfv_ups(ji ,jj-1) ) * r1_e1e2t(ji,jj) 455 IF( ll_clem ) THEN 456 IF( ll_gurvan ) THEN 457 zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 458 ELSE 459 zt_ups(ji,jj) = ( pt (ji,jj) + pdt * ztra + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 460 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 461 ENDIF 462 ELSE 463 zt_ups(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 464 ENDIF 465 !! IF( ji==26 .AND. jj==86) THEN 466 !! WRITE(numout,*) '**************************' 467 !! WRITE(numout,*) 'zt upstream',zt_ups(ji,jj) 468 !! ENDIF 469 END DO 470 END DO 471 CALL lbc_lnk( zt_ups, 'T', 1. ) 472 186 473 ! High order (_ho) fluxes 187 474 ! ----------------------- 188 SELECT CASE( k_order ) 189 CASE ( 20 ) ! centered second order 475 SELECT CASE( kn_umx ) 476 ! 477 CASE ( 20 ) !== centered second order ==! 478 ! 479 CALL cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho, & 480 & zt_ups, zfu_ups, zfv_ups ) 481 ! 482 CASE ( 1:5 ) !== 1st to 5th order ULTIMATE-MACHO scheme ==! 483 ! 484 CALL macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho, & 485 & zt_ups, zfu_ups, zfv_ups ) 486 ! 487 END SELECT 488 489 IF( ll_thickness ) THEN 490 ! final trend with corrected fluxes 491 ! ------------------------------------ 492 DO jj = 2, jpjm1 493 DO ji = fs_2, fs_jpim1 494 IF( ll_gurvan ) THEN 495 ztra = - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj) 496 ELSE 497 ztra = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) & 498 & + ( pt(ji,jj) * ( pu(ji,jj) - pu(ji-1,jj) ) * (1.-pamsk) ) & 499 & + ( pt(ji,jj) * ( pv(ji,jj) - pv(ji,jj-1) ) * (1.-pamsk) ) ) * r1_e1e2t(ji,jj) 500 ENDIF 501 pt(ji,jj) = ( pt(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 502 503 IF( pt(ji,jj) < -epsi20 ) THEN 504 WRITE(numout,*) 'T<0 ',pt(ji,jj) 505 ENDIF 506 507 IF( pt(ji,jj) < 0._wp .AND. pt(ji,jj) >= -epsi20 ) pt(ji,jj) = 0._wp 508 509 !! IF( ji==26 .AND. jj==86) THEN 510 !! WRITE(numout,*) 'zt high order',pt(ji,jj) 511 !! ENDIF 512 END DO 513 END DO 514 CALL lbc_lnk( pt, 'T', 1. ) 515 ENDIF 516 517 ! Rachid trick 518 ! ------------ 519 IF( ll_clem ) THEN 520 IF( pamsk == 0. ) THEN 521 DO jj = 1, jpjm1 522 DO ji = 1, fs_jpim1 523 IF( ABS( puc(ji,jj) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN 524 zfu_ho (ji,jj) = zfu_ho (ji,jj) * puc(ji,jj) / pu(ji,jj) 525 zfu_ups(ji,jj) = zfu_ups(ji,jj) * puc(ji,jj) / pu(ji,jj) 526 ELSE 527 zfu_ho (ji,jj) = 0._wp 528 zfu_ups(ji,jj) = 0._wp 529 ENDIF 530 ! 531 IF( ABS( pvc(ji,jj) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN 532 zfv_ho (ji,jj) = zfv_ho (ji,jj) * pvc(ji,jj) / pv(ji,jj) 533 zfv_ups(ji,jj) = zfv_ups(ji,jj) * pvc(ji,jj) / pv(ji,jj) 534 ELSE 535 zfv_ho (ji,jj) = 0._wp 536 zfv_ups(ji,jj) = 0._wp 537 ENDIF 538 ENDDO 539 ENDDO 540 ENDIF 541 ENDIF 542 543 IF( ll_zeroup5 ) THEN 544 DO jj = 2, jpjm1 545 DO ji = 2, fs_jpim1 ! vector opt. 546 zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 547 & - ( zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 548 IF( zpt(ji,jj) < 0. ) THEN 549 zfu_ho(ji,jj) = zfu_ups(ji,jj) 550 zfu_ho(ji-1,jj) = zfu_ups(ji-1,jj) 551 zfv_ho(ji,jj) = zfv_ups(ji,jj) 552 zfv_ho(ji,jj-1) = zfv_ups(ji,jj-1) 553 ENDIF 554 END DO 555 END DO 556 CALL lbc_lnk_multi( zfu_ho, 'U', -1., zfv_ho, 'V', -1. ) 557 ENDIF 558 559 ! output upstream trend of concentration and high order fluxes 560 ! ------------------------------------------------------------ 561 IF( ll_dens ) THEN 562 ! high order u*a 190 563 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 ! 564 DO ji = 1, fs_jpim1 565 pua_ho (ji,jj) = zfu_ho (ji,jj) 566 pva_ho (ji,jj) = zfv_ho (ji,jj) 567 END DO 568 END DO 569 !!CALL lbc_lnk( pua_ho, 'U', -1. ) ! clem: not needed I think 570 !!CALL lbc_lnk( pva_ho, 'V', -1. ) 571 ENDIF 572 573 574 IF( .NOT.ll_thickness ) THEN 575 ! final trend with corrected fluxes 576 ! ------------------------------------ 200 577 DO jj = 2, jpjm1 201 DO ji = 1, fs_jpim1 ! vector opt. 202 zfu_ho(ji,jj) = puc(ji,jj) * zt_u(ji,jj) 203 END DO 204 END DO 578 DO ji = fs_2, fs_jpim1 579 ztra = - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) & ! Div(uaH) or Div(ua) 580 & * r1_e1e2t(ji,jj) * pdt 581 582 !!IF( ptc(ji,jj)+ztra < 0._wp ) THEN 583 !! ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) & ! Div(uaH) or Div(ua) 584 !! & * r1_e1e2t(ji,jj) * pdt 585 !!ENDIF 586 !!IF( ptc(ji,jj)+ztra < 0._wp ) THEN 587 !! WRITE(numout,*) 'Tc<0 ',ptc(ji,jj)+ztra 588 !! ztra = 0._wp 589 !!ENDIF 590 591 ptc(ji,jj) = ( ptc(ji,jj) + ztra ) * tmask(ji,jj,1) 592 593 !! IF( ji==26 .AND. jj==86) THEN 594 !! WRITE(numout,*) 'ztc high order',ptc(ji,jj) 595 !! ENDIF 596 597 END DO 598 END DO 599 CALL lbc_lnk( ptc, 'T', 1. ) 600 ENDIF 601 602 ! 603 END SUBROUTINE adv_umx 604 605 SUBROUTINE cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, & 606 & pt_ups, pfu_ups, pfv_ups ) 607 !!--------------------------------------------------------------------- 608 !! *** ROUTINE macho *** 609 !! 610 !! ** Purpose : compute 611 !! 612 !! ** Method : ... ??? 613 !! TIM = transient interpolation Modeling 614 !! 615 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 616 !!---------------------------------------------------------------------- 617 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 618 INTEGER , INTENT(in ) :: kn_limiter ! limiter 619 INTEGER , INTENT(in ) :: jt ! number of sub-iteration 620 INTEGER , INTENT(in ) :: kt ! number of iteration 621 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 622 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields 623 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu, pv ! 2 ice velocity components 624 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc, pvc ! 2 ice velocity * A components 625 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptc ! tracer content at before time step 626 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes 627 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt_ups ! upstream guess of tracer content 628 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pfu_ups, pfv_ups ! upstream fluxes 629 ! 630 INTEGER :: ji, jj ! dummy loop indices 631 LOGICAL :: ll_xy = .TRUE. 632 REAL(wp), DIMENSION(jpi,jpj) :: zzt 633 !!---------------------------------------------------------------------- 634 ! 635 IF( .NOT.ll_xy ) THEN !-- no alternate directions --! 636 ! 205 637 DO jj = 1, jpjm1 206 DO ji = fs_2, fs_jpim1 ! vector opt. 207 zfv_ho(ji,jj) = pvc(ji,jj) * zt_v(ji,jj) 208 END DO 209 END DO 210 ! 211 END SELECT 638 DO ji = 1, fs_jpim1 639 pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) ) 640 pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) ) 641 END DO 642 END DO 643 IF ( kn_limiter == 1 ) THEN 644 IF( ll_clem ) THEN 645 CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 646 ELSE 647 CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 648 ENDIF 649 ELSEIF( kn_limiter == 2 ) THEN 650 CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 651 CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 652 ELSEIF( kn_limiter == 3 ) THEN 653 CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 654 CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 655 ENDIF 656 ! 657 ELSE !-- alternate directions --! 658 ! 659 IF( pamsk == 1. ) z1_ai(:,:) = 1._wp 660 ! 661 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 662 ! 663 ! flux in x-direction 664 DO jj = 1, jpjm1 665 DO ji = 1, fs_jpim1 666 pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) ) 667 END DO 668 END DO 669 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 670 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 671 672 ! first guess of tracer content from u-flux 673 DO jj = 2, jpjm1 674 DO ji = fs_2, fs_jpim1 ! vector opt. 675 IF( ll_clem ) THEN 676 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) * z1_ai(ji,jj) ) * tmask(ji,jj,1) 677 ELSE 678 zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * z1_ai(ji,jj) 679 ENDIF 680 END DO 681 END DO 682 CALL lbc_lnk( zzt, 'T', 1. ) 683 684 ! flux in y-direction 685 DO jj = 1, jpjm1 686 DO ji = 1, fs_jpim1 687 pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( zzt(ji,jj) + zzt(ji,jj+1) ) 688 END DO 689 END DO 690 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 691 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 692 693 ELSE !== even ice time step: adv_y then adv_x ==! 694 ! 695 ! flux in y-direction 696 DO jj = 1, jpjm1 697 DO ji = 1, fs_jpim1 698 pfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) ) 699 END DO 700 END DO 701 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 702 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 703 ! 704 ! first guess of tracer content from v-flux 705 DO jj = 2, jpjm1 706 DO ji = fs_2, fs_jpim1 ! vector opt. 707 IF( ll_clem ) THEN 708 zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) * z1_ai(ji,jj) ) * tmask(ji,jj,1) 709 ELSE 710 zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) * z1_ai(ji,jj) 711 ENDIF 712 END DO 713 END DO 714 CALL lbc_lnk( zzt, 'T', 1. ) 715 ! 716 ! flux in x-direction 717 DO jj = 1, jpjm1 718 DO ji = 1, fs_jpim1 719 pfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( zzt(ji,jj) + zzt(ji+1,jj) ) 720 END DO 721 END DO 722 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 723 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 724 725 ENDIF 726 IF( ll_clem ) THEN 727 IF( kn_limiter == 1 ) CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 728 ELSE 729 IF( kn_limiter == 1 ) CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 730 ENDIF 212 731 213 ! antidiffusive flux : high order minus low order 214 ! -------------------------------------------------- 215 DO jj = 2, jpjm1 216 DO ji = 1, fs_jpim1 ! vector opt. 217 zfu_ho(ji,jj) = zfu_ho(ji,jj) - zfu_ups(ji,jj) 218 END DO 219 END DO 220 DO jj = 1, jpjm1 221 DO ji = fs_2, fs_jpim1 ! vector opt. 222 zfv_ho(ji,jj) = zfv_ho(ji,jj) - zfv_ups(ji,jj) 223 END DO 224 END DO 225 226 ! monotonicity algorithm 227 ! ------------------------- 228 CALL nonosc_2d( ptc, zfu_ho, zfv_ho, zt_ups, pdt ) 229 230 ! final trend with corrected fluxes 231 ! ------------------------------------ 232 DO jj = 2, jpjm1 233 DO ji = fs_2, fs_jpim1 ! vector opt. 234 ztra = ztrd(ji,jj) - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj ) & 235 & + zfv_ho(ji,jj) - zfv_ho(ji ,jj-1) ) * r1_e1e2t(ji,jj) 236 ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 237 END DO 238 END DO 239 CALL lbc_lnk( ptc, 'T', 1. ) 240 ! 241 END SUBROUTINE adv_umx 242 243 244 SUBROUTINE macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v ) 732 ENDIF 733 734 END SUBROUTINE cen2 735 736 737 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, & 738 & pt_ups, pfu_ups, pfv_ups ) 739 !!--------------------------------------------------------------------- 740 !! *** ROUTINE macho *** 741 !! 742 !! ** Purpose : compute 743 !! 744 !! ** Method : ... ??? 745 !! TIM = transient interpolation Modeling 746 !! 747 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 748 !!---------------------------------------------------------------------- 749 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 750 INTEGER , INTENT(in ) :: kn_limiter ! limiter 751 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 752 INTEGER , INTENT(in ) :: jt ! number of sub-iteration 753 INTEGER , INTENT(in ) :: kt ! number of iteration 754 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 755 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields 756 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu, pv ! 2 ice velocity components 757 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc, pvc ! 2 ice velocity * A components 758 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity 759 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptc ! tracer content at before time step 760 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_u, pt_v ! tracer at u- and v-points 761 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes 762 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt_ups ! upstream guess of tracer content 763 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pfu_ups, pfv_ups ! upstream fluxes 764 ! 765 INTEGER :: ji, jj ! dummy loop indices 766 REAL(wp) :: ztra 767 REAL(wp), DIMENSION(jpi,jpj) :: zzt, zzfu_ho, zzfv_ho 768 !!---------------------------------------------------------------------- 769 ! 770 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 771 ! 772 ! !-- ultimate interpolation of pt at u-point --! 773 CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 774 ! !-- limiter in x --! 775 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 776 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 777 ! !-- advective form update in zzt --! 778 779 IF( ll_1stguess_clem ) THEN 780 781 ! first guess of tracer content from u-flux 782 DO jj = 2, jpjm1 783 DO ji = fs_2, fs_jpim1 ! vector opt. 784 IF( ll_clem ) THEN 785 IF( ll_gurvan ) THEN 786 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 787 ELSE 788 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 789 & + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 790 ENDIF 791 ELSE 792 zzt(ji,jj) = ( ptc(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 793 ENDIF 794 END DO 795 END DO 796 CALL lbc_lnk( zzt, 'T', 1. ) 797 798 ELSE 799 800 DO jj = 2, jpjm1 801 DO ji = fs_2, fs_jpim1 ! vector opt. 802 IF( ll_gurvan ) THEN 803 zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj) & 804 & - pt (ji,jj) * pdt * ( pu (ji,jj) - pu (ji-1,jj) ) * r1_e1e2t(ji,jj) 805 ELSE 806 zzt(ji,jj) = pt(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj) & 807 & - pt (ji,jj) * pdt * ( pu (ji,jj) - pu (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk 808 ENDIF 809 zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 810 END DO 811 END DO 812 CALL lbc_lnk( zzt, 'T', 1. ) 813 ENDIF 814 ! 815 ! !-- ultimate interpolation of pt at v-point --! 816 IF( ll_hoxy ) THEN 817 CALL ultimate_y( kn_umx, pdt, zzt, pv, pvc, pt_v, pfv_ho ) 818 ELSE 819 CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 820 ENDIF 821 ! !-- limiter in y --! 822 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 823 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 824 ! 825 ! 826 ELSE !== even ice time step: adv_y then adv_x ==! 827 ! 828 ! !-- ultimate interpolation of pt at v-point --! 829 CALL ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 830 ! !-- limiter in y --! 831 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho ) 832 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 833 ! !-- advective form update in zzt --! 834 IF( ll_1stguess_clem ) THEN 835 836 ! first guess of tracer content from v-flux 837 DO jj = 2, jpjm1 838 DO ji = fs_2, fs_jpim1 ! vector opt. 839 IF( ll_clem ) THEN 840 IF( ll_gurvan ) THEN 841 zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 842 ELSE 843 zzt(ji,jj) = ( pt(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) & 844 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 845 ENDIF 846 ELSE 847 zzt(ji,jj) = ( ptc(ji,jj) - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) & 848 & * tmask(ji,jj,1) 849 ENDIF 850 END DO 851 END DO 852 CALL lbc_lnk( zzt, 'T', 1. ) 853 854 ELSE 855 856 DO jj = 2, jpjm1 857 DO ji = fs_2, fs_jpim1 858 IF( ll_gurvan ) THEN 859 zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj) & 860 & - pt (ji,jj) * pdt * ( pv (ji,jj) - pv (ji,jj-1) ) * r1_e1e2t(ji,jj) 861 ELSE 862 zzt(ji,jj) = pt(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj) & 863 & - pt (ji,jj) * pdt * ( pv (ji,jj) - pv (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk 864 ENDIF 865 zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 866 END DO 867 END DO 868 CALL lbc_lnk( zzt, 'T', 1. ) 869 ENDIF 870 ! 871 ! !-- ultimate interpolation of pt at u-point --! 872 IF( ll_hoxy ) THEN 873 CALL ultimate_x( kn_umx, pdt, zzt, pu, puc, pt_u, pfu_ho ) 874 ELSE 875 CALL ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 876 ENDIF 877 ! !-- limiter in x --! 878 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 879 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 880 ! 881 ! 882 ENDIF 883 884 885 IF( kn_limiter == 1 ) THEN 886 IF( .NOT. ll_limiter_it2 ) THEN 887 IF( ll_clem ) THEN 888 CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 889 ELSE 890 CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 891 ENDIF 892 ELSE 893 zzfu_ho(:,:) = pfu_ho(:,:) 894 zzfv_ho(:,:) = pfv_ho(:,:) 895 ! 1st iteration of nonosc (limit the flux with the upstream solution) 896 IF( ll_clem ) THEN 897 CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho ) 898 ELSE 899 CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho ) 900 ENDIF 901 ! guess after content field with high order 902 DO jj = 2, jpjm1 903 DO ji = fs_2, fs_jpim1 904 ztra = - ( zzfu_ho(ji,jj) - zzfu_ho(ji-1,jj) + zzfv_ho(ji,jj) - zzfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj) 905 zzt(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 906 END DO 907 END DO 908 CALL lbc_lnk( zzt, 'T', 1. ) 909 ! 2nd iteration of nonosc (limit the flux with the limited high order solution) 910 IF( ll_clem ) THEN 911 CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho ) 912 ELSE 913 CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho ) 914 ENDIF 915 ENDIF 916 ENDIF 917 ! 918 END SUBROUTINE macho 919 920 921 SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, puc, pt_u, pfu_ho ) 245 922 !!--------------------------------------------------------------------- 246 923 !! *** ROUTINE ultimate_x *** … … 253 930 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 254 931 !!---------------------------------------------------------------------- 255 INTEGER , INTENT(in ) :: k_order ! order of the ULTIMATE scheme 256 INTEGER , INTENT(in ) :: kt ! number of iteration 257 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 258 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptc ! tracer fields 259 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc, pvc ! 2 ice velocity components 260 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity 261 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_u, pt_v ! tracer at u- and v-points 262 ! 263 INTEGER :: ji, jj ! dummy loop indices 264 REAL(wp) :: zc_box ! - - 265 REAL(wp), DIMENSION(jpi,jpj) :: zzt 266 !!---------------------------------------------------------------------- 267 ! 268 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 269 ! 270 ! !-- ultimate interpolation of pt at u-point --! 271 CALL ultimate_x( k_order, pdt, ptc, puc, pt_u ) 272 ! 273 ! !-- advective form update in zzt --! 274 DO jj = 2, jpjm1 275 DO ji = fs_2, fs_jpim1 ! vector opt. 276 zzt(ji,jj) = ptc(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj) & 277 & - ptc (ji,jj) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e1e2t(ji,jj) 278 zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 279 END DO 280 END DO 281 CALL lbc_lnk( zzt, 'T', 1. ) 282 ! 283 ! !-- ultimate interpolation of pt at v-point --! 284 CALL ultimate_y( k_order, pdt, zzt, pvc, pt_v ) 285 ! 286 ELSE !== even ice time step: adv_y then adv_x ==! 287 ! 288 ! !-- ultimate interpolation of pt at v-point --! 289 CALL ultimate_y( k_order, pdt, ptc, pvc, pt_v ) 290 ! 291 ! !-- advective form update in zzt --! 292 DO jj = 2, jpjm1 293 DO ji = fs_2, fs_jpim1 294 zzt(ji,jj) = ptc(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj) & 295 & - ptc (ji,jj) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e1e2t(ji,jj) 296 zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 297 END DO 298 END DO 299 CALL lbc_lnk( zzt, 'T', 1. ) 300 ! 301 ! !-- ultimate interpolation of pt at u-point --! 302 CALL ultimate_x( k_order, pdt, zzt, puc, pt_u ) 303 ! 304 ENDIF 305 ! 306 END SUBROUTINE macho 307 308 309 SUBROUTINE ultimate_x( k_order, pdt, pt, puc, pt_u ) 310 !!--------------------------------------------------------------------- 311 !! *** ROUTINE ultimate_x *** 312 !! 313 !! ** Purpose : compute 314 !! 315 !! ** Method : ... ??? 316 !! TIM = transient interpolation Modeling 317 !! 318 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 319 !!---------------------------------------------------------------------- 320 INTEGER , INTENT(in ) :: k_order ! ocean time-step index 932 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 321 933 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 322 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc ! ice i-velocity component 934 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu ! ice i-velocity component 935 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc ! ice i-velocity * A component 323 936 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields 324 937 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_u ! tracer at u-point 325 ! 326 INTEGER :: ji, jj ! dummy loop indices 938 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pfu_ho ! high order flux 939 ! 940 INTEGER :: ji, jj ! dummy loop indices 327 941 REAL(wp) :: zcu, zdx2, zdx4 ! - - 328 REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu2, ztu3, ztu4942 REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu2, ztu3, ztu4 329 943 !!---------------------------------------------------------------------- 330 944 ! … … 354 968 ! 355 969 ! 356 SELECT CASE (k _order)970 SELECT CASE (kn_umx ) 357 971 ! 358 972 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 359 973 ! 360 DO jj = 2, jpjm1974 DO jj = 1, jpjm1 361 975 DO ji = 1, fs_jpim1 ! vector opt. 362 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( 363 & - SIGN( 1._wp, pu c(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) )976 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) & 977 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) ) 364 978 END DO 365 979 END DO … … 367 981 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 368 982 ! 369 DO jj = 2, jpjm1983 DO jj = 1, jpjm1 370 984 DO ji = 1, fs_jpim1 ! vector opt. 371 zcu = pu c(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)985 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 372 986 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) & 373 987 & - zcu * ( pt(ji+1,jj) - pt(ji,jj) ) ) … … 377 991 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 378 992 ! 379 DO jj = 2, jpjm1993 DO jj = 1, jpjm1 380 994 DO ji = 1, fs_jpim1 ! vector opt. 381 zcu = pu c(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)995 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 382 996 zdx2 = e1u(ji,jj) * e1u(ji,jj) 383 997 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) … … 391 1005 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 392 1006 ! 393 DO jj = 2, jpjm11007 DO jj = 1, jpjm1 394 1008 DO ji = 1, fs_jpim1 ! vector opt. 395 zcu = pu c(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)1009 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 396 1010 zdx2 = e1u(ji,jj) * e1u(ji,jj) 397 1011 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) … … 405 1019 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 406 1020 ! 407 DO jj = 2, jpjm11021 DO jj = 1, jpjm1 408 1022 DO ji = 1, fs_jpim1 ! vector opt. 409 zcu = pu c(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)1023 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 410 1024 zdx2 = e1u(ji,jj) * e1u(ji,jj) 411 1025 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) … … 421 1035 ! 422 1036 END SELECT 1037 ! !-- High order flux in i-direction --! 1038 IF( ll_neg ) THEN 1039 DO jj = 1, jpjm1 1040 DO ji = 1, fs_jpim1 1041 IF( pt_u(ji,jj) < 0._wp ) THEN 1042 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) & 1043 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) ) 1044 ENDIF 1045 END DO 1046 END DO 1047 ENDIF 1048 1049 DO jj = 1, jpjm1 1050 DO ji = 1, fs_jpim1 ! vector opt. 1051 IF( ll_clem ) THEN 1052 pfu_ho(ji,jj) = pu(ji,jj) * pt_u(ji,jj) 1053 ELSE 1054 pfu_ho(ji,jj) = puc(ji,jj) * pt_u(ji,jj) 1055 ENDIF 1056 END DO 1057 END DO 423 1058 ! 424 1059 END SUBROUTINE ultimate_x 425 1060 426 1061 427 SUBROUTINE ultimate_y( k _order, pdt, pt, pvc, pt_v)1062 SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pvc, pt_v, pfv_ho ) 428 1063 !!--------------------------------------------------------------------- 429 1064 !! *** ROUTINE ultimate_y *** … … 436 1071 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 437 1072 !!---------------------------------------------------------------------- 438 INTEGER , INTENT(in ) :: k _order ! ocean time-step index1073 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 439 1074 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 440 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pvc ! ice j-velocity component 1075 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pv ! ice j-velocity component 1076 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pvc ! ice j-velocity*A component 441 1077 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields 442 1078 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_v ! tracer at v-point 1079 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pfv_ho ! high order flux 443 1080 ! 444 1081 INTEGER :: ji, jj ! dummy loop indices 445 1082 REAL(wp) :: zcv, zdy2, zdy4 ! - - 446 REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv2, ztv3, ztv41083 REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv2, ztv3, ztv4 447 1084 !!---------------------------------------------------------------------- 448 1085 ! … … 474 1111 ! 475 1112 ! 476 SELECT CASE (k _order)1113 SELECT CASE (kn_umx ) 477 1114 ! 478 1115 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 479 1116 DO jj = 1, jpjm1 480 DO ji = fs_2, fs_jpim1481 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( 482 & - SIGN( 1._wp, pv c(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )1117 DO ji = 1, fs_jpim1 1118 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) & 1119 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) 483 1120 END DO 484 1121 END DO … … 486 1123 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 487 1124 DO jj = 1, jpjm1 488 DO ji = fs_2, fs_jpim1489 zcv = pv c(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)1125 DO ji = 1, fs_jpim1 1126 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 490 1127 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) & 491 1128 & - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) ) … … 496 1133 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 497 1134 DO jj = 1, jpjm1 498 DO ji = fs_2, fs_jpim1499 zcv = pv c(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)1135 DO ji = 1, fs_jpim1 1136 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 500 1137 zdy2 = e2v(ji,jj) * e2v(ji,jj) 501 1138 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) … … 509 1146 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 510 1147 DO jj = 1, jpjm1 511 DO ji = fs_2, fs_jpim1512 zcv = pv c(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)1148 DO ji = 1, fs_jpim1 1149 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 513 1150 zdy2 = e2v(ji,jj) * e2v(ji,jj) 514 1151 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) … … 522 1159 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 523 1160 DO jj = 1, jpjm1 524 DO ji = fs_2, fs_jpim1525 zcv = pv c(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)1161 DO ji = 1, fs_jpim1 1162 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 526 1163 zdy2 = e2v(ji,jj) * e2v(ji,jj) 527 1164 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) … … 537 1174 ! 538 1175 END SELECT 1176 ! !-- High order flux in j-direction --! 1177 IF( ll_neg ) THEN 1178 DO jj = 1, jpjm1 1179 DO ji = 1, fs_jpim1 1180 IF( pt_v(ji,jj) < 0._wp ) THEN 1181 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) & 1182 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) 1183 ENDIF 1184 END DO 1185 END DO 1186 ENDIF 1187 1188 DO jj = 1, jpjm1 1189 DO ji = 1, fs_jpim1 ! vector opt. 1190 IF( ll_clem ) THEN 1191 pfv_ho(ji,jj) = pv(ji,jj) * pt_v(ji,jj) 1192 ELSE 1193 pfv_ho(ji,jj) = pvc(ji,jj) * pt_v(ji,jj) 1194 ENDIF 1195 END DO 1196 END DO 539 1197 ! 540 1198 END SUBROUTINE ultimate_y 541 542 543 SUBROUTINE nonosc_2d( p bef, paa, pbb, paft, pdt)1199 1200 1201 SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho ) 544 1202 !!--------------------------------------------------------------------- 545 1203 !! *** ROUTINE nonosc *** 546 1204 !! 547 !! ** Purpose : compute monotonic tracer fluxes from the upstream1205 !! ** Purpose : compute monotonic tracer fluxes from the upstream 548 1206 !! scheme and the before field by a nonoscillatory algorithm 549 1207 !! 550 1208 !! ** Method : ... ??? 551 !! warning : p bef and paftmust be masked, but the boundaries1209 !! warning : pt and pt_low must be masked, but the boundaries 552 1210 !! conditions on the fluxes are not necessary zalezak (1979) 553 1211 !! drange (1995) multi-dimensional forward-in-time and upstream- 554 1212 !! in-space based differencing for fluid 555 1213 !!---------------------------------------------------------------------- 556 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 557 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pbef, paft ! before & after field 558 REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: paa, pbb ! monotonic fluxes in the 2 directions 1214 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 1215 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1216 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pu ! ice i-velocity => u*e2 1217 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: puc ! ice i-velocity *A => u*e2*a 1218 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pv ! ice j-velocity => v*e1 1219 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pvc ! ice j-velocity *A => v*e1*a 1220 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: ptc, pt, pt_low ! before field & upstream guess of after field 1221 REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: pfv_low, pfu_low ! upstream flux 1222 REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: pfv_ho, pfu_ho ! monotonic flux 559 1223 ! 560 1224 INTEGER :: ji, jj ! dummy loop indices 561 INTEGER :: ikm1 ! local integer 562 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zsml, z1_dt ! local scalars 563 REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - 564 REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zmsk, zdiv 565 !!---------------------------------------------------------------------- 566 ! 1225 REAL(wp) :: zpos, zneg, zbig, zsml, z1_dt, zpos2, zneg2 ! local scalars 1226 REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign, zcoef ! - - 1227 REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zti_low, ztj_low, zzt 1228 !!---------------------------------------------------------------------- 567 1229 zbig = 1.e+40_wp 568 zsml = 1.e-15_wp 569 570 ! test on divergence 571 DO jj = 2, jpjm1 572 DO ji = fs_2, fs_jpim1 ! vector opt. 573 zdiv(ji,jj) = - ( paa(ji,jj) - paa(ji-1,jj ) & 574 & + pbb(ji,jj) - pbb(ji ,jj-1) ) 575 END DO 576 END DO 577 CALL lbc_lnk( zdiv, 'T', 1. ) ! Lateral boundary conditions (unchanged sign) 578 579 ! Determine ice masks for before and after tracers 580 WHERE( pbef(:,:) == 0._wp .AND. paft(:,:) == 0._wp .AND. zdiv(:,:) == 0._wp ) ; zmsk(:,:) = 0._wp 581 ELSEWHERE ; zmsk(:,:) = 1._wp * tmask(:,:,1) 582 END WHERE 583 1230 zsml = epsi20 1231 1232 IF( ll_zeroup2 ) THEN 1233 DO jj = 1, jpjm1 1234 DO ji = 1, fs_jpim1 ! vector opt. 1235 IF( amaxu(ji,jj) == 0._wp ) pfu_ho(ji,jj) = 0._wp 1236 IF( amaxv(ji,jj) == 0._wp ) pfv_ho(ji,jj) = 0._wp 1237 END DO 1238 END DO 1239 ENDIF 1240 1241 IF( ll_zeroup4 ) THEN 1242 DO jj = 1, jpjm1 1243 DO ji = 1, fs_jpim1 ! vector opt. 1244 IF( pfu_low(ji,jj) == 0._wp ) pfu_ho(ji,jj) = 0._wp 1245 IF( pfv_low(ji,jj) == 0._wp ) pfv_ho(ji,jj) = 0._wp 1246 END DO 1247 END DO 1248 ENDIF 1249 1250 1251 IF( ll_zeroup1 ) THEN 1252 DO jj = 2, jpjm1 1253 DO ji = fs_2, fs_jpim1 1254 IF( ll_gurvan ) THEN 1255 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 1256 & - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 1257 ELSE 1258 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 1259 & - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) & 1260 & + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1261 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 1262 ENDIF 1263 IF( zzt(ji,jj) < 0._wp ) THEN 1264 pfu_ho(ji,jj) = pfu_low(ji,jj) 1265 pfv_ho(ji,jj) = pfv_low(ji,jj) 1266 WRITE(numout,*) '*** 1 negative high order zzt ***',ji,jj,zzt(ji,jj) 1267 ENDIF 1268 !! IF( ji==26 .AND. jj==86) THEN 1269 !! WRITE(numout,*) 'zzt high order',zzt(ji,jj) 1270 !! WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt 1271 !! WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt 1272 !! WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt 1273 !! WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt 1274 !! ENDIF 1275 IF( ll_gurvan ) THEN 1276 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 1277 & - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 1278 ELSE 1279 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 1280 & - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) & 1281 & + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1282 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 1283 ENDIF 1284 IF( zzt(ji,jj) < 0._wp ) THEN 1285 pfu_ho(ji-1,jj) = pfu_low(ji-1,jj) 1286 pfv_ho(ji,jj-1) = pfv_low(ji,jj-1) 1287 WRITE(numout,*) '*** 2 negative high order zzt ***',ji,jj,zzt(ji,jj) 1288 ENDIF 1289 IF( ll_gurvan ) THEN 1290 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 1291 & - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 1292 ELSE 1293 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 1294 & - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) & 1295 & + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1296 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 1297 ENDIF 1298 IF( zzt(ji,jj) < 0._wp ) THEN 1299 WRITE(numout,*) '*** 3 negative high order zzt ***',ji,jj,zzt(ji,jj) 1300 ENDIF 1301 END DO 1302 END DO 1303 CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) 1304 ENDIF 1305 1306 1307 ! antidiffusive flux : high order minus low order 1308 ! -------------------------------------------------- 1309 DO jj = 1, jpjm1 1310 DO ji = 1, fs_jpim1 ! vector opt. 1311 pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_low(ji,jj) 1312 pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_low(ji,jj) 1313 END DO 1314 END DO 1315 1316 ! extreme case where pfu_ho has to be zero 1317 ! ---------------------------------------- 1318 ! pfu_ho 1319 ! * ---> 1320 ! | | * | | 1321 ! | | | * | 1322 ! | | | | * 1323 ! t_low : i-1 i i+1 i+2 1324 IF( ll_prelimiter_zalesak ) THEN 1325 1326 DO jj = 2, jpjm1 1327 DO ji = fs_2, fs_jpim1 1328 zti_low(ji,jj)= pt_low(ji+1,jj ) 1329 ztj_low(ji,jj)= pt_low(ji ,jj+1) 1330 END DO 1331 END DO 1332 CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. ) 1333 1334 1335 !! this does not work 1336 !! DO jj = 2, jpjm1 1337 !! DO ji = fs_2, fs_jpim1 1338 !! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji+1,jj ) - pt_low (ji ,jj) ) .AND. & 1339 !! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji ,jj+1) - pt_low (ji ,jj) ) & 1340 !! & ) THEN 1341 !! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_low(ji+1,jj ) - zti_low(ji ,jj) ) .AND. & 1342 !! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_low(ji,jj+1 ) - ztj_low(ji ,jj) ) & 1343 !! & ) THEN 1344 !! pfu_ho(ji,jj) = 0. ; pfv_ho(ji,jj) = 0. 1345 !! ENDIF 1346 !! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji ,jj) - pt_low (ji-1,jj ) ) .AND. & 1347 !! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji ,jj) - pt_low (ji ,jj-1) ) & 1348 !! & ) THEN 1349 !! pfu_ho(ji,jj) = 0. ; pfv_ho(ji,jj) = 0. 1350 !! ENDIF 1351 !! ENDIF 1352 !! END DO 1353 !! END DO 1354 1355 DO jj = 2, jpjm1 1356 DO ji = fs_2, fs_jpim1 1357 IF ( pfu_ho(ji,jj) * ( pt_low(ji+1,jj) - pt_low(ji,jj) ) <= 0. .AND. & 1358 & pfv_ho(ji,jj) * ( pt_low(ji,jj+1) - pt_low(ji,jj) ) <= 0. ) THEN 1359 ! 1360 IF( pfu_ho(ji,jj) * ( zti_low(ji+1,jj) - zti_low(ji,jj) ) <= 0 .AND. & 1361 & pfv_ho(ji,jj) * ( ztj_low(ji,jj+1) - ztj_low(ji,jj) ) <= 0) pfu_ho(ji,jj)=0. ; pfv_ho(ji,jj)=0. 1362 ! 1363 IF( pfu_ho(ji,jj) * ( pt_low(ji ,jj) - pt_low(ji-1,jj) ) <= 0 .AND. & 1364 & pfv_ho(ji,jj) * ( pt_low(ji ,jj) - pt_low(ji,jj-1) ) <= 0) pfu_ho(ji,jj)=0. ; pfv_ho(ji,jj)=0. 1365 ! 1366 ENDIF 1367 END DO 1368 END DO 1369 CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond. 1370 1371 ELSEIF( ll_prelimiter_devore ) THEN 1372 DO jj = 2, jpjm1 1373 DO ji = fs_2, fs_jpim1 1374 zti_low(ji,jj)= pt_low(ji+1,jj ) 1375 ztj_low(ji,jj)= pt_low(ji ,jj+1) 1376 END DO 1377 END DO 1378 CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. ) 1379 1380 z1_dt = 1._wp / pdt 1381 DO jj = 2, jpjm1 1382 DO ji = fs_2, fs_jpim1 1383 zsign = SIGN( 1., pt_low(ji+1,jj) - pt_low(ji,jj) ) 1384 pfu_ho(ji,jj) = zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj)) , & 1385 & zsign * ( pt_low (ji ,jj) - pt_low (ji-1,jj) ) * e1e2t(ji ,jj) * z1_dt , & 1386 & zsign * ( zti_low(ji+1,jj) - zti_low(ji ,jj) ) * e1e2t(ji+1,jj) * z1_dt ) ) 1387 1388 zsign = SIGN( 1., pt_low(ji,jj+1) - pt_low(ji,jj) ) 1389 pfv_ho(ji,jj) = zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj)) , & 1390 & zsign * ( pt_low (ji,jj ) - pt_low (ji,jj-1) ) * e1e2t(ji,jj ) * z1_dt , & 1391 & zsign * ( ztj_low(ji,jj+1) - ztj_low(ji,jj ) ) * e1e2t(ji,jj+1) * z1_dt ) ) 1392 END DO 1393 END DO 1394 CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond. 1395 1396 ENDIF 1397 1398 584 1399 ! Search local extrema 585 1400 ! -------------------- 586 ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 587 ! zbup(:,:) = MAX( pbef(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ), & 588 ! & paft(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ) ) 589 ! zbdo(:,:) = MIN( pbef(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ), & 590 ! & paft(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ) ) 591 zbup(:,:) = MAX( pbef(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ), & 592 & paft(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ) ) 593 zbdo(:,:) = MIN( pbef(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ), & 594 & paft(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ) ) 595 1401 ! max/min of pt & pt_low with large negative/positive value (-/+zbig) outside ice cover 1402 DO jj = 1, jpj 1403 DO ji = 1, jpi 1404 IF ( pt(ji,jj) <= 0._wp .AND. pt_low(ji,jj) <= 0._wp ) THEN 1405 zbup(ji,jj) = -zbig 1406 zbdo(ji,jj) = zbig 1407 ELSEIF( pt(ji,jj) <= 0._wp .AND. pt_low(ji,jj) > 0._wp ) THEN 1408 zbup(ji,jj) = pt_low(ji,jj) 1409 zbdo(ji,jj) = pt_low(ji,jj) 1410 ELSEIF( pt(ji,jj) > 0._wp .AND. pt_low(ji,jj) <= 0._wp ) THEN 1411 zbup(ji,jj) = pt(ji,jj) 1412 zbdo(ji,jj) = pt(ji,jj) 1413 ELSE 1414 zbup(ji,jj) = MAX( pt(ji,jj) , pt_low(ji,jj) ) 1415 zbdo(ji,jj) = MIN( pt(ji,jj) , pt_low(ji,jj) ) 1416 ENDIF 1417 END DO 1418 END DO 1419 1420 596 1421 z1_dt = 1._wp / pdt 597 1422 DO jj = 2, jpjm1 598 1423 DO ji = fs_2, fs_jpim1 ! vector opt. 599 1424 ! 600 zup = MAX( zbup(ji,jj), zbup(ji-1,jj ), zbup(ji+1,jj ), & ! search max/min in neighbourhood 601 & zbup(ji ,jj-1), zbup(ji ,jj+1) ) 602 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), & 603 & zbdo(ji ,jj-1), zbdo(ji ,jj+1) ) 604 ! 605 zpos = MAX( 0., paa(ji-1,jj ) ) - MIN( 0., paa(ji ,jj ) ) & ! positive/negative part of the flux 606 & + MAX( 0., pbb(ji ,jj-1) ) - MIN( 0., pbb(ji ,jj ) ) 607 zneg = MAX( 0., paa(ji ,jj ) ) - MIN( 0., paa(ji-1,jj ) ) & 608 & + MAX( 0., pbb(ji ,jj ) ) - MIN( 0., pbb(ji ,jj-1) ) 609 ! 610 zbt = e1e2t(ji,jj) * z1_dt ! up & down beta terms 611 zbetup(ji,jj) = ( zup - paft(ji,jj) ) / ( zpos + zsml ) * zbt 612 zbetdo(ji,jj) = ( paft(ji,jj) - zdo ) / ( zneg + zsml ) * zbt 1425 IF( .NOT. ll_9points ) THEN 1426 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 1427 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), zbdo(ji ,jj-1), zbdo(ji ,jj+1) ) 1428 ! 1429 ELSE 1430 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 1431 & zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1) ) 1432 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), zbdo(ji ,jj-1), zbdo(ji ,jj+1), & 1433 & zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1) ) 1434 ENDIF 1435 ! 1436 zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji ,jj) ) & ! positive/negative part of the flux 1437 & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj ) ) 1438 zneg = MAX( 0., pfu_ho(ji ,jj) ) - MIN( 0., pfu_ho(ji-1,jj) ) & 1439 & + MAX( 0., pfv_ho(ji,jj ) ) - MIN( 0., pfv_ho(ji,jj-1) ) 1440 ! 1441 IF( ll_HgradU .AND. .NOT.ll_gurvan ) THEN 1442 zneg2 = ( pt(ji,jj) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) ) * ( 1. - pamsk ) 1443 zpos2 = ( - pt(ji,jj) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) - pt(ji,jj) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) ) * ( 1. - pamsk ) 1444 ELSE 1445 zneg2 = 0. ; zpos2 = 0. 1446 ENDIF 1447 ! 1448 ! ! up & down beta terms 1449 ! zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt 1450 ! zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt 1451 1452 IF( (zpos+zpos2) > 0. ) THEN ; zbetup(ji,jj) = MAX( 0._wp, zup - pt_low(ji,jj) ) / (zpos+zpos2) * e1e2t(ji,jj) * z1_dt 1453 ELSE ; zbetup(ji,jj) = 0. ! zbig 1454 ENDIF 1455 ! 1456 IF( (zneg+zneg2) > 0. ) THEN ; zbetdo(ji,jj) = MAX( 0._wp, pt_low(ji,jj) - zdo ) / (zneg+zneg2) * e1e2t(ji,jj) * z1_dt 1457 ELSE ; zbetdo(ji,jj) = 0. ! zbig 1458 ENDIF 1459 ! 1460 ! if all the points are outside ice cover 1461 IF( zup == -zbig ) zbetup(ji,jj) = 0. ! zbig 1462 IF( zdo == zbig ) zbetdo(ji,jj) = 0. ! zbig 1463 ! 1464 1465 !! IF( ji==26 .AND. jj==86) THEN 1466 ! WRITE(numout,*) '-----------------' 1467 ! WRITE(numout,*) 'zpos',zpos,zpos2 1468 ! WRITE(numout,*) 'zneg',zneg,zneg2 1469 ! WRITE(numout,*) 'puc/pu',ABS(puc(ji,jj))/MAX(epsi20, ABS(pu(ji,jj))) 1470 ! WRITE(numout,*) 'pvc/pv',ABS(pvc(ji,jj))/MAX(epsi20, ABS(pv(ji,jj))) 1471 ! WRITE(numout,*) 'pucm1/pu',ABS(puc(ji-1,jj))/MAX(epsi20, ABS(pu(ji-1,jj))) 1472 ! WRITE(numout,*) 'pvcm1/pv',ABS(pvc(ji,jj-1))/MAX(epsi20, ABS(pv(ji,jj-1))) 1473 ! WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)+pfu_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt 1474 ! WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)+pfv_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt 1475 ! WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)+pfu_low(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt 1476 ! WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)+pfv_low(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt 1477 ! WRITE(numout,*) 'pfu_low',pfu_low(ji,jj) * r1_e1e2t(ji,jj) * pdt 1478 ! WRITE(numout,*) 'pfv_low',pfv_low(ji,jj) * r1_e1e2t(ji,jj) * pdt 1479 ! WRITE(numout,*) 'pfu_lowm1',pfu_low(ji-1,jj) * r1_e1e2t(ji,jj) * pdt 1480 ! WRITE(numout,*) 'pfv_lowm1',pfv_low(ji,jj-1) * r1_e1e2t(ji,jj) * pdt 1481 ! 1482 ! WRITE(numout,*) 'pt',pt(ji,jj) 1483 ! WRITE(numout,*) 'ptim1',pt(ji-1,jj) 1484 ! WRITE(numout,*) 'ptjm1',pt(ji,jj-1) 1485 ! WRITE(numout,*) 'pt_low',pt_low(ji,jj) 1486 ! WRITE(numout,*) 'zbetup',zbetup(ji,jj) 1487 ! WRITE(numout,*) 'zbetdo',zbetdo(ji,jj) 1488 ! WRITE(numout,*) 'zup',zup 1489 ! WRITE(numout,*) 'zdo',zdo 1490 ! ENDIF 1491 ! 613 1492 END DO 614 1493 END DO 615 1494 CALL lbc_lnk_multi( zbetup, 'T', 1., zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign) 616 1495 617 ! monotonic flux in the i & j direction (paa & pbb) 618 ! ------------------------------------- 619 DO jj = 2, jpjm1 1496 1497 ! monotonic flux in the y direction 1498 ! --------------------------------- 1499 DO jj = 1, jpjm1 620 1500 DO ji = 1, fs_jpim1 ! vector opt. 621 1501 zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) ) 622 1502 zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) ) 623 zcu = 0.5 + SIGN( 0.5 , paa(ji,jj) ) 624 ! 625 paa(ji,jj) = paa(ji,jj) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 626 END DO 627 END DO 628 ! 1503 zcu = 0.5 + SIGN( 0.5 , pfu_ho(ji,jj) ) 1504 ! 1505 zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu ) 1506 1507 pfu_ho(ji,jj) = pfu_ho(ji,jj) * zcoef + pfu_low(ji,jj) 1508 1509 !! IF( ji==26 .AND. jj==86) THEN 1510 !! WRITE(numout,*) 'coefU',zcoef 1511 !! WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt 1512 !! WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt 1513 !! ENDIF 1514 1515 END DO 1516 END DO 1517 629 1518 DO jj = 1, jpjm1 630 DO ji = fs_2, fs_jpim1 ! vector opt.1519 DO ji = 1, fs_jpim1 ! vector opt. 631 1520 zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) ) 632 1521 zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) ) 633 zcv = 0.5 + SIGN( 0.5 , pbb(ji,jj) ) 634 ! 635 pbb(ji,jj) = pbb(ji,jj) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 636 END DO 637 END DO 1522 zcv = 0.5 + SIGN( 0.5 , pfv_ho(ji,jj) ) 1523 ! 1524 zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv ) 1525 1526 pfv_ho(ji,jj) = pfv_ho(ji,jj) * zcoef + pfv_low(ji,jj) 1527 1528 !! IF( ji==26 .AND. jj==86) THEN 1529 !! WRITE(numout,*) 'coefV',zcoef 1530 !! WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)) * r1_e1e2t(ji,jj) * pdt 1531 !! WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt 1532 !! ENDIF 1533 END DO 1534 END DO 1535 1536 ! clem test 1537 DO jj = 2, jpjm1 1538 DO ji = 2, fs_jpim1 ! vector opt. 1539 IF( ll_gurvan ) THEN 1540 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 1541 & - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 1542 ELSE 1543 zzt(ji,jj) = ( pt(ji,jj) - ( pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) & 1544 & - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) & 1545 & + pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1546 & + pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 1547 ENDIF 1548 IF( zzt(ji,jj) < -epsi20 ) THEN 1549 WRITE(numout,*) 'T<0 nonosc',zzt(ji,jj) 1550 ENDIF 1551 END DO 1552 END DO 1553 1554 ! 638 1555 ! 639 1556 END SUBROUTINE nonosc_2d 1557 1558 SUBROUTINE limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 1559 !!--------------------------------------------------------------------- 1560 !! *** ROUTINE limiter_x *** 1561 !! 1562 !! ** Purpose : compute flux limiter 1563 !!---------------------------------------------------------------------- 1564 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1565 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pu ! ice i-velocity => u*e2 1566 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: puc ! ice i-velocity *A => u*e2*a 1567 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pt ! ice tracer 1568 REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: pfu_ho ! high order flux 1569 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ), OPTIONAL :: pfu_ups ! upstream flux 1570 ! 1571 REAL(wp) :: Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr 1572 INTEGER :: ji, jj ! dummy loop indices 1573 REAL(wp), DIMENSION (jpi,jpj) :: zslpx ! tracer slopes 1574 !!---------------------------------------------------------------------- 1575 ! 1576 DO jj = 2, jpjm1 1577 DO ji = fs_2, fs_jpim1 ! vector opt. 1578 zslpx(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * umask(ji,jj,1) 1579 END DO 1580 END DO 1581 CALL lbc_lnk( zslpx, 'U', -1.) ! lateral boundary cond. 1582 1583 DO jj = 2, jpjm1 1584 DO ji = fs_2, fs_jpim1 ! vector opt. 1585 uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj) 1586 1587 Rjm = zslpx(ji-1,jj) 1588 Rj = zslpx(ji ,jj) 1589 Rjp = zslpx(ji+1,jj) 1590 1591 IF( PRESENT(pfu_ups) ) THEN 1592 1593 IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm 1594 ELSE ; Rr = Rjp 1595 ENDIF 1596 1597 zh3 = pfu_ho(ji,jj) - pfu_ups(ji,jj) 1598 IF( Rj > 0. ) THEN 1599 zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(puc(ji,jj)), & 1600 & MIN( 2. * Rr * 0.5 * ABS(puc(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) ) 1601 ELSE 1602 zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(puc(ji,jj)), & 1603 & MIN(-2. * Rr * 0.5 * ABS(puc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(puc(ji,jj)) ) ) ) ) 1604 ENDIF 1605 pfu_ho(ji,jj) = pfu_ups(ji,jj) + zlimiter 1606 1607 ELSE 1608 IF( Rj /= 0. ) THEN 1609 IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj 1610 ELSE ; Cr = Rjp / Rj 1611 ENDIF 1612 ELSE 1613 Cr = 0. 1614 !IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm * 1.e20 1615 !ELSE ; Cr = Rjp * 1.e20 1616 !ENDIF 1617 ENDIF 1618 1619 ! -- superbee -- 1620 zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 1621 ! -- van albada 2 -- 1622 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1623 1624 ! -- sweby (with beta=1) -- 1625 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 1626 ! -- van Leer -- 1627 !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 1628 ! -- ospre -- 1629 !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 1630 ! -- koren -- 1631 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 1632 ! -- charm -- 1633 !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 1634 !ELSE ; zpsi = 0. 1635 !ENDIF 1636 ! -- van albada 1 -- 1637 !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 1638 ! -- smart -- 1639 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 1640 ! -- umist -- 1641 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 1642 1643 ! high order flux corrected by the limiter 1644 pfu_ho(ji,jj) = pfu_ho(ji,jj) - ABS( puc(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5 1645 1646 ENDIF 1647 END DO 1648 END DO 1649 CALL lbc_lnk( pfu_ho, 'U', -1.) ! lateral boundary cond. 1650 ! 1651 END SUBROUTINE limiter_x 1652 1653 SUBROUTINE limiter_y( pdt, pv, pvc, pt, pfv_ho, pfv_ups ) 1654 !!--------------------------------------------------------------------- 1655 !! *** ROUTINE limiter_y *** 1656 !! 1657 !! ** Purpose : compute flux limiter 1658 !!---------------------------------------------------------------------- 1659 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1660 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pv ! ice i-velocity => u*e2 1661 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pvc ! ice i-velocity *A => u*e2*a 1662 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pt ! ice tracer 1663 REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: pfv_ho ! high order flux 1664 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ), OPTIONAL :: pfv_ups ! upstream flux 1665 ! 1666 REAL(wp) :: Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr 1667 INTEGER :: ji, jj ! dummy loop indices 1668 REAL(wp), DIMENSION (jpi,jpj) :: zslpy ! tracer slopes 1669 !!---------------------------------------------------------------------- 1670 ! 1671 DO jj = 2, jpjm1 1672 DO ji = fs_2, fs_jpim1 ! vector opt. 1673 zslpy(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * vmask(ji,jj,1) 1674 END DO 1675 END DO 1676 CALL lbc_lnk( zslpy, 'V', -1.) ! lateral boundary cond. 1677 1678 DO jj = 2, jpjm1 1679 DO ji = fs_2, fs_jpim1 ! vector opt. 1680 vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj) 1681 1682 Rjm = zslpy(ji,jj-1) 1683 Rj = zslpy(ji,jj ) 1684 Rjp = zslpy(ji,jj+1) 1685 1686 IF( PRESENT(pfv_ups) ) THEN 1687 1688 IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm 1689 ELSE ; Rr = Rjp 1690 ENDIF 1691 1692 zh3 = pfv_ho(ji,jj) - pfv_ups(ji,jj) 1693 IF( Rj > 0. ) THEN 1694 zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pvc(ji,jj)), & 1695 & MIN( 2. * Rr * 0.5 * ABS(pvc(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) ) 1696 ELSE 1697 zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pvc(ji,jj)), & 1698 & MIN(-2. * Rr * 0.5 * ABS(pvc(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pvc(ji,jj)) ) ) ) ) 1699 ENDIF 1700 pfv_ho(ji,jj) = pfv_ups(ji,jj) + zlimiter 1701 1702 ELSE 1703 1704 IF( Rj /= 0. ) THEN 1705 IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj 1706 ELSE ; Cr = Rjp / Rj 1707 ENDIF 1708 ELSE 1709 Cr = 0. 1710 !IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm * 1.e20 1711 !ELSE ; Cr = Rjp * 1.e20 1712 !ENDIF 1713 ENDIF 1714 1715 ! -- superbee -- 1716 zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 1717 ! -- van albada 2 -- 1718 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1719 1720 ! -- sweby (with beta=1) -- 1721 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 1722 ! -- van Leer -- 1723 !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 1724 ! -- ospre -- 1725 !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 1726 ! -- koren -- 1727 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 1728 ! -- charm -- 1729 !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 1730 !ELSE ; zpsi = 0. 1731 !ENDIF 1732 ! -- van albada 1 -- 1733 !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 1734 ! -- smart -- 1735 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 1736 ! -- umist -- 1737 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 1738 1739 ! high order flux corrected by the limiter 1740 pfv_ho(ji,jj) = pfv_ho(ji,jj) - ABS( pvc(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5 1741 1742 ENDIF 1743 END DO 1744 END DO 1745 CALL lbc_lnk( pfv_ho, 'V', -1.) ! lateral boundary cond. 1746 ! 1747 END SUBROUTINE limiter_y 640 1748 641 1749 #else
Note: See TracChangeset
for help on using the changeset viewer.
