Changeset 10315 for NEMO/branches
- Timestamp:
- 2018-11-15T17:31:50+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_adv_umx.F90
r10267 r10315 35 35 REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120 36 36 37 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z1_a, z1_a_ups, zua_ups, zva_ups 38 39 ! alternate directions for upstream 40 ! clem: it gives results for Lipscomb test that are the same as "ll_upsxy=false" 41 ! clem: needs to be set to true in 2D when using prelimiter (otherwise "wavy solutions" are created) 42 LOGICAL :: ll_upsxy = .TRUE. 43 44 ! prelimiter 45 ! clem: use it to avoid overshoot in H 46 LOGICAL :: ll_prelimiter = .TRUE. 47 LOGICAL :: ll_prelimiter_zalesak = .TRUE. ! from: Zalesak(1979) eq. 14 => better than Devore especially in 2D 48 LOGICAL :: ll_prelimiter_devore = .FALSE. ! from: Devore eq. 11 49 50 ! iterate on the limiter (only nonosc) 51 ! clem: useless 52 LOGICAL :: ll_limiter_it2 = .FALSE. 53 54 37 55 !! * Substitutions 38 56 # include "vectopt_loop_substitute.h90" … … 74 92 REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers 75 93 REAL(wp) :: zcfl , zdt 76 REAL(wp) :: zeps = 0. 1_wp ! shift in concentration to avoid division by 094 REAL(wp) :: zeps = 0.0_wp ! shift in concentration to avoid division by 0 77 95 ! ! must be >= 0.01 and the best seems to be 0.1 78 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box, z fu, zfv79 REAL(wp), DIMENSION(jpi,jpj) :: z 1_a, zh_i, zh_s, zs_i, zo_i, ze_i, ze_s, z1_ap, zh_ip96 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box, zua_ho, zva_ho, za_ups 97 REAL(wp), DIMENSION(jpi,jpj) :: zh_i, zh_s, zs_i, zo_i, ze_i, ze_s, zh_ip 80 98 !!---------------------------------------------------------------------- 81 99 ! … … 112 130 END DO 113 131 132 IF(.NOT. ALLOCATED(z1_a) ) ALLOCATE(z1_a (jpi,jpj)) 133 IF(.NOT. ALLOCATED(z1_a_ups)) ALLOCATE(z1_a_ups(jpi,jpj)) 134 IF(.NOT. ALLOCATED(zua_ups) ) ALLOCATE(zua_ups (jpi,jpj)) 135 IF(.NOT. ALLOCATED(zva_ups) ) ALLOCATE(zva_ups (jpi,jpj)) 114 136 !---------------! 115 137 !== advection ==! … … 117 139 DO jt = 1, icycle 118 140 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 area141 zamsk = 1._wp ; zua_ups(:,:) = zudy(:,:) ; zva_ups(:,:) = zvdx(:,:) 142 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:), zua_ups, zva_ups, za_ups, zua_ho, zva_ho ) ! Open water area 121 143 122 144 DO jl = 1, jpl … … 125 147 ! 126 148 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.149 ELSEWHERE ; z1_a(:,:) = 0. 128 150 END WHERE 129 151 ! 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) ) 152 zamsk = 1._wp ; zua_ups(:,:) = zudy(:,:) ; zva_ups(:,:) = zvdx(:,:) 153 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl), pa_i(:,:,jl), zua_ups, zva_ups, za_ups, zua_ho, zva_ho ) ! Ice area 154 155 ! 1/A_ups 156 WHERE( za_ups(:,:) > epsi20 ) ; z1_a_ups(:,:) = 1._wp / za_ups(:,:) 157 ELSEWHERE ; z1_a_ups(:,:) = 0. 158 END WHERE 159 160 !!zua_ho = zudy ; zva_ho = zvdx 161 !!CALL adv_umx( kn_umx, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, pv_i(:,:,jl), pv_i(:,:,jl) ) 134 162 zamsk = 0._wp ; zh_i(:,:) = pv_i (:,:,jl) * z1_a(:,:) 135 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, z fu , zfv , zcu_box, zcv_box, zh_i(:,:), pv_i (:,:,jl)) ! Ice volume163 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zh_i(:,:), pv_i (:,:,jl), zua_ups, zva_ups ) ! Ice volume 136 164 zamsk = 0._wp ; zh_s(:,:) = pv_s (:,:,jl) * z1_a(:,:) 137 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, z fu , zfv , zcu_box, zcv_box, zh_s(:,:), pv_s (:,:,jl)) ! Snw volume165 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zh_s(:,:), pv_s (:,:,jl), zua_ups, zva_ups ) ! Snw volume 138 166 zamsk = 0._wp ; zs_i(:,:) = psv_i(:,:,jl) * z1_a(:,:) 139 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, z fu , zfv , zcu_box, zcv_box, zs_i(:,:), psv_i(:,:,jl)) ! Salt content167 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zs_i(:,:), psv_i(:,:,jl), zua_ups, zva_ups ) ! Salt content 140 168 zamsk = 0._wp ; zo_i(:,:) = poa_i(:,:,jl) * z1_a(:,:) 141 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, z fu , zfv , zcu_box, zcv_box, zo_i(:,:), poa_i(:,:,jl)) ! Age content169 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zo_i(:,:), poa_i(:,:,jl), zua_ups, zva_ups ) ! Age content 142 170 DO jk = 1, nlay_i 143 171 zamsk = 0._wp ; ze_i(:,:) = pe_i(:,:,jk,jl) * z1_a(:,:) 144 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, z fu, zfv, zcu_box, zcv_box, ze_i(:,:), pe_i(:,:,jk,jl)) ! Ice heat content172 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, ze_i(:,:), pe_i(:,:,jk,jl), zua_ups, zva_ups ) ! Ice heat content 145 173 END DO 146 174 DO jk = 1, nlay_s 147 175 zamsk = 0._wp ; ze_s(:,:) = pe_s(:,:,jk,jl) * z1_a(:,:) 148 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, z fu, zfv, zcu_box, zcv_box, ze_s(:,:), pe_s(:,:,jk,jl)) ! Snw heat content149 END DO 150 ! 151 IF ( ln_pnd_H12 ) THEN ! melt ponds 176 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, ze_s(:,:), pe_s(:,:,jk,jl), zua_ups, zva_ups ) ! Snw heat content 177 END DO 178 ! 179 IF ( ln_pnd_H12 ) THEN ! melt ponds (must be the last ones to be advected because of z1_a...) 152 180 ! to avoid a problem with the limiter nonosc when A gets close to 0 153 181 pa_ip(:,:,jl) = pa_ip(:,:,jl) + zeps * tmask(:,:,1) 154 182 ! 155 WHERE( pa_ip(:,:,jl) > epsi20 ) ; z1_a p(:,:) = 1._wp / pa_ip(:,:,jl)156 ELSEWHERE ; z1_a p(:,:) = 0.183 WHERE( pa_ip(:,:,jl) > epsi20 ) ; z1_a(:,:) = 1._wp / pa_ip(:,:,jl) 184 ELSEWHERE ; z1_a(:,:) = 0. 157 185 END WHERE 158 186 ! 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 187 zamsk = 1._wp ; zua_ups(:,:) = zudy(:,:) ; zva_ups(:,:) = zvdx(:,:) 188 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl), pa_ip(:,:,jl), zua_ups, zva_ups, za_ups, zua_ho, zva_ho ) ! mp fraction 189 190 ! 1/A_ups 191 WHERE( za_ups(:,:) > epsi20 ) ; z1_a_ups(:,:) = 1._wp / za_ups(:,:) 192 ELSEWHERE ; z1_a_ups(:,:) = 0. 193 END WHERE 194 161 195 zamsk = 0._wp ; zh_ip(:,:) = pv_ip(:,:,jl) * z1_a(:,:) 162 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, z fu , zfv , zcu_box, zcv_box, zh_ip(:,:), pv_ip(:,:,jl)) ! mp volume196 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zh_ip(:,:), pv_ip(:,:,jl), zua_ups, zva_ups ) ! mp volume 163 197 ENDIF 164 198 ! … … 177 211 END DO 178 212 END DO 179 CALL lbc_lnk ( pa_i(:,:,jl), 'T', 1. )213 CALL lbc_lnk_multi( pa_i(:,:,jl), 'T', 1., pa_ip(:,:,jl), 'T', 1. ) 180 214 ! 181 215 END DO … … 186 220 187 221 188 SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, p fu, pfv)222 SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ups, pva_ups, pa_ups, pua_ho, pva_ho ) 189 223 !!---------------------------------------------------------------------- 190 224 !! *** ROUTINE adv_umx *** … … 209 243 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer field 210 244 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 245 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua_ups, pva_ups ! upstream u*a fluxes or u 246 REAL(wp), DIMENSION(jpi,jpj), INTENT( out), OPTIONAL :: pa_ups ! concentration advected upstream 247 REAL(wp), DIMENSION(jpi,jpj), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes 212 248 ! 213 249 INTEGER :: ji, jj ! dummy loop indices 214 250 REAL(wp) :: ztra ! local scalar 215 !!clem REAL(wp) :: zeps = 1.e-02 ! local scalar216 251 INTEGER :: kn_limiter = 1 ! 1=nonosc ; 2=superbee ; 3=h3(rachid) 217 252 REAL(wp), DIMENSION(jpi,jpj) :: zfu_ho , zfv_ho , zt_u, zt_v, zpt … … 219 254 !!---------------------------------------------------------------------- 220 255 ! 221 ! add a constant value to avoid problems with zeros222 DO jj = 1, jpj223 DO ji = 1, jpi224 zpt(ji,jj) = pt(ji,jj) !!clem + zeps * tmask(ji,jj,1)225 END DO226 END DO227 228 256 ! upstream (_ups) advection with initial mass fluxes 229 257 ! --------------------------------------------------- 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 258 IF( .NOT. ll_upsxy ) THEN 259 260 ! fluxes 261 DO jj = 1, jpjm1 262 DO ji = 1, fs_jpim1 263 zfu_ups(ji,jj) = MAX( pua_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pua_ups(ji,jj), 0._wp ) * pt(ji+1,jj) 264 zfv_ups(ji,jj) = MAX( pva_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pva_ups(ji,jj), 0._wp ) * pt(ji,jj+1) 265 END DO 266 END DO 267 268 ELSE 269 ! 1 if advection of A 270 ! z1_a_ups already defined IF advection of other tracers 271 IF( pamsk == 1. ) z1_a_ups(:,:) = 1._wp 272 ! 273 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 274 ! flux in x-direction 275 DO jj = 1, jpjm1 276 DO ji = 1, fs_jpim1 277 zfu_ups(ji,jj) = MAX( pua_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pua_ups(ji,jj), 0._wp ) * pt(ji+1,jj) 278 END DO 279 END DO 280 ! first guess of tracer content from u-flux 281 DO jj = 2, jpjm1 282 DO ji = fs_2, fs_jpim1 ! vector opt. 283 zpt(ji,jj) = ( ptc(ji,jj) - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) * pdt * r1_e1e2t(ji,jj) ) & 284 & * tmask(ji,jj,1) * z1_a_ups(ji,jj) 285 END DO 286 END DO 287 CALL lbc_lnk( zpt, 'T', 1. ) 288 ! 289 ! flux in y-direction 290 DO jj = 1, jpjm1 291 DO ji = 1, fs_jpim1 292 zfv_ups(ji,jj) = MAX( pva_ups(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pva_ups(ji,jj), 0._wp ) * zpt(ji,jj+1) 293 END DO 294 END DO 295 ! 296 ELSE !== even ice time step: adv_y then adv_x ==! 297 ! flux in y-direction 298 DO jj = 1, jpjm1 299 DO ji = 1, fs_jpim1 300 zfv_ups(ji,jj) = MAX( pva_ups(ji,jj), 0._wp ) * pt(ji,jj) + MIN( pva_ups(ji,jj), 0._wp ) * pt(ji,jj+1) 301 END DO 302 END DO 303 ! first guess of tracer content from v-flux 304 DO jj = 2, jpjm1 305 DO ji = fs_2, fs_jpim1 ! vector opt. 306 zpt(ji,jj) = ( ptc(ji,jj) - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) * pdt * r1_e1e2t(ji,jj) ) & 307 & * tmask(ji,jj,1) * z1_a_ups(ji,jj) 308 END DO 309 END DO 310 CALL lbc_lnk( zpt, 'T', 1. ) 311 ! 312 ! flux in x-direction 313 DO jj = 1, jpjm1 314 DO ji = 1, fs_jpim1 315 zfu_ups(ji,jj) = MAX( pua_ups(ji,jj), 0._wp ) * zpt(ji,jj) + MIN( pua_ups(ji,jj), 0._wp ) * zpt(ji+1,jj) 316 END DO 317 END DO 318 ! 319 ENDIF 320 321 ENDIF 237 322 ! guess after content field with upstream scheme 238 323 DO jj = 2, jpjm1 … … 244 329 END DO 245 330 CALL lbc_lnk( zt_ups, 'T', 1. ) 246 331 247 332 ! High order (_ho) fluxes 248 333 ! ----------------------- … … 251 336 CASE ( 20 ) !== centered second order ==! 252 337 ! 253 CALL cen2( kn_limiter, jt, kt, pdt, zpt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho, &338 CALL cen2( kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho, & 254 339 & zt_ups, zfu_ups, zfv_ups ) 255 340 ! 256 341 CASE ( 1:5 ) !== 1st to 5th order ULTIMATE-MACHO scheme ==! 257 342 ! 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, &343 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, & 259 344 & zt_ups, zfu_ups, zfv_ups ) 260 345 ! 261 346 END SELECT 262 347 263 ! output high order fluxes 264 ! ------------------------ 265 IF( PRESENT(pfu) ) THEN 348 ! output upstream trend of concentration and high order fluxes 349 ! ------------------------------------------------------------ 350 IF( pamsk == 1. ) THEN 351 ! upstream trend of concentration 352 pa_ups(:,:) = zt_ups(:,:) 353 ! upstream and high order u*a 266 354 DO jj = 1, jpjm1 267 355 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. ) 356 pua_ups(ji,jj) = zfu_ups(ji,jj) 357 pva_ups(ji,jj) = zfv_ups(ji,jj) 358 pua_ho (ji,jj) = zfu_ho (ji,jj) 359 pva_ho (ji,jj) = zfv_ho (ji,jj) 360 END DO 361 END DO 362 !!CALL lbc_lnk( pua_ho, 'U', -1. ) ! clem: not needed I think 363 !!CALL lbc_lnk( pva_ho, 'V', -1. ) 274 364 ENDIF 275 365 … … 278 368 DO jj = 2, jpjm1 279 369 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) 370 ztra = ( - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj) + zfv_ho(ji,jj) - zfv_ho(ji,jj-1) ) & ! Div(uaH) or Div(ua) 282 371 & ) * r1_e1e2t(ji,jj) 283 372 ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) … … 430 519 ! 431 520 INTEGER :: ji, jj ! dummy loop indices 432 REAL(wp), DIMENSION(jpi,jpj) :: zzt 521 REAL(wp) :: ztra 522 REAL(wp), DIMENSION(jpi,jpj) :: zzt, zzfu_ho, zzfv_ho 433 523 !!---------------------------------------------------------------------- 434 524 ! … … 478 568 ! 479 569 ENDIF 480 IF( kn_limiter == 1 ) CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 570 IF( kn_limiter == 1 ) THEN 571 IF( .NOT. ll_limiter_it2 ) THEN 572 CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 573 ELSE 574 zzfu_ho(:,:) = pfu_ho(:,:) 575 zzfv_ho(:,:) = pfv_ho(:,:) 576 ! 1st iteration of nonosc (limit the flux with the upstream solution) 577 CALL nonosc_2d ( pdt, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho ) 578 ! guess after content field with high order 579 DO jj = 2, jpjm1 580 DO ji = fs_2, fs_jpim1 581 ztra = - ( zzfu_ho(ji,jj) - zzfu_ho(ji-1,jj) + zzfv_ho(ji,jj) - zzfv_ho(ji,jj-1) ) * r1_e1e2t(ji,jj) 582 zzt(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 583 END DO 584 END DO 585 CALL lbc_lnk( zzt, 'T', 1. ) 586 ! 2nd iteration of nonosc (limit the flux with the limited high order solution) 587 CALL nonosc_2d ( pdt, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho ) 588 ENDIF 589 ENDIF 481 590 ! 482 591 END SUBROUTINE macho … … 733 842 734 843 735 SUBROUTINE nonosc_2d( pdt, ptc, pt_ ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )844 SUBROUTINE nonosc_2d( pdt, ptc, pt_low, pfu_low, pfv_low, pfu_ho, pfv_ho ) 736 845 !!--------------------------------------------------------------------- 737 846 !! *** ROUTINE nonosc *** 738 847 !! 739 !! ** Purpose : compute monotonic tracer fluxes from the upstream848 !! ** Purpose : compute monotonic tracer fluxes from the upstream 740 849 !! scheme and the before field by a nonoscillatory algorithm 741 850 !! 742 851 !! ** Method : ... ??? 743 !! warning : ptc and pt_ upsmust be masked, but the boundaries852 !! warning : ptc and pt_low must be masked, but the boundaries 744 853 !! conditions on the fluxes are not necessary zalezak (1979) 745 854 !! drange (1995) multi-dimensional forward-in-time and upstream- … … 747 856 !!---------------------------------------------------------------------- 748 857 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 field750 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pfv_ ups, pfu_ups! upstream flux858 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: ptc, pt_low ! before field & upstream guess of after field 859 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pfv_low, pfu_low ! upstream flux 751 860 REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: pfv_ho, pfu_ho ! monotonic flux 752 861 ! 753 862 INTEGER :: ji, jj ! dummy loop indices 754 863 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, z div864 REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo, zsign ! - - 865 REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zti_low, ztj_low 757 866 !!---------------------------------------------------------------------- 758 867 zbig = 1.e+40_wp 759 zsml = 1.e-15_wp 760 761 ! test on divergence 762 DO jj = 2, jpjm1 763 DO ji = fs_2, fs_jpim1 ! vector opt. 764 zdiv(ji,jj) = - ( pfv_ho(ji,jj) - pfv_ho(ji,jj-1) + pfu_ho(ji,jj) - pfu_ho(ji-1,jj) ) 765 END DO 766 END DO 767 CALL lbc_lnk( zdiv, 'T', 1. ) ! Lateral boundary conditions (unchanged sign) 868 zsml = epsi20 768 869 769 870 ! antidiffusive flux : high order minus low order … … 771 872 DO jj = 1, jpjm1 772 873 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 777 874 pfu_ho(ji,jj) = pfu_ho(ji,jj) - pfu_low(ji,jj) 875 pfv_ho(ji,jj) = pfv_ho(ji,jj) - pfv_low(ji,jj) 876 END DO 877 END DO 878 879 ! extreme case where pfu_ho has to be zero 880 ! ---------------------------------------- 881 ! pfu_ho 882 ! * ---> 883 ! | | * | | 884 ! | | | * | 885 ! | | | | * 886 ! t_low : i-1 i i+1 i+2 887 IF( ll_prelimiter ) THEN 888 889 DO jj = 2, jpjm1 890 DO ji = fs_2, fs_jpim1 891 zti_low(ji,jj)= pt_low(ji+1,jj ) 892 ztj_low(ji,jj)= pt_low(ji ,jj+1) 893 END DO 894 END DO 895 CALL lbc_lnk_multi( zti_low, 'T', 1., ztj_low, 'T', 1. ) 896 897 IF( ll_prelimiter_zalesak ) THEN 898 899 !! this does not work 900 !! DO jj = 2, jpjm1 901 !! DO ji = fs_2, fs_jpim1 902 !! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji+1,jj ) - pt_low (ji ,jj) ) .AND. & 903 !! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji ,jj+1) - pt_low (ji ,jj) ) & 904 !! & ) THEN 905 !! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_low(ji+1,jj ) - zti_low(ji ,jj) ) .AND. & 906 !! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_low(ji,jj+1 ) - ztj_low(ji ,jj) ) & 907 !! & ) THEN 908 !! pfu_ho(ji,jj) = 0. ; pfv_ho(ji,jj) = 0. 909 !! ENDIF 910 !! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_low (ji ,jj) - pt_low (ji-1,jj ) ) .AND. & 911 !! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_low (ji ,jj) - pt_low (ji ,jj-1) ) & 912 !! & ) THEN 913 !! pfu_ho(ji,jj) = 0. ; pfv_ho(ji,jj) = 0. 914 !! ENDIF 915 !! ENDIF 916 !! END DO 917 !! END DO 918 919 DO jj = 2, jpjm1 920 DO ji = fs_2, fs_jpim1 921 IF ( pfu_ho(ji,jj) * ( pt_low(ji+1,jj) - pt_low(ji,jj) ) <= 0. .AND. & 922 & pfv_ho(ji,jj) * ( pt_low(ji,jj+1) - pt_low(ji,jj) ) <= 0. ) THEN 923 ! 924 IF( pfu_ho(ji,jj) * ( zti_low(ji+1,jj) - zti_low(ji,jj) ) <= 0 .AND. & 925 & pfv_ho(ji,jj) * ( ztj_low(ji,jj+1) - ztj_low(ji,jj) ) <= 0) pfu_ho(ji,jj)=0. ; pfv_ho(ji,jj)=0. 926 ! 927 IF( pfu_ho(ji,jj) * ( pt_low(ji ,jj) - pt_low(ji-1,jj) ) <= 0 .AND. & 928 & pfv_ho(ji,jj) * ( pt_low(ji ,jj) - pt_low(ji,jj-1) ) <= 0) pfu_ho(ji,jj)=0. ; pfv_ho(ji,jj)=0. 929 ! 930 ENDIF 931 END DO 932 END DO 933 CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond. 934 935 ELSEIF( ll_prelimiter_devore ) THEN 936 z1_dt = 1._wp / pdt 937 DO jj = 2, jpjm1 938 DO ji = fs_2, fs_jpim1 939 zsign = SIGN( 1., pt_low(ji+1,jj) - pt_low(ji,jj) ) 940 pfu_ho(ji,jj) = zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj)) , & 941 & zsign * ( pt_low (ji ,jj) - pt_low (ji-1,jj) ) * e1e2t(ji ,jj) * z1_dt , & 942 & zsign * ( zti_low(ji+1,jj) - zti_low(ji ,jj) ) * e1e2t(ji+1,jj) * z1_dt ) ) 943 944 zsign = SIGN( 1., pt_low(ji,jj+1) - pt_low(ji,jj) ) 945 pfv_ho(ji,jj) = zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj)) , & 946 & zsign * ( pt_low (ji,jj ) - pt_low (ji,jj-1) ) * e1e2t(ji,jj ) * z1_dt , & 947 & zsign * ( ztj_low(ji,jj+1) - ztj_low(ji,jj ) ) * e1e2t(ji,jj+1) * z1_dt ) ) 948 END DO 949 END DO 950 CALL lbc_lnk_multi( pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond. 951 952 ENDIF 953 954 ENDIF 955 778 956 ! Search local extrema 779 957 ! -------------------- 780 ! max/min of ptc & pt_ upswith large negative/positive value (-/+zbig) outside ice cover958 ! max/min of ptc & pt_low with large negative/positive value (-/+zbig) outside ice cover 781 959 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 960 DO ji = 1, jpi 961 !!clem IF ( ptc(ji,jj) == 0._wp .AND. pt_low(ji,jj) == 0._wp ) THEN 962 IF ( ptc(ji,jj) <= epsi20 .AND. pt_low(ji,jj) <= epsi20 ) THEN 784 963 zbup(ji,jj) = -zbig 785 964 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) ) 965 !!clem ELSEIF( ptc(ji,jj) == 0._wp .AND. pt_low(ji,jj) /= 0._wp ) THEN 966 ELSEIF( ptc(ji,jj) <= epsi20 .AND. pt_low(ji,jj) > epsi20 ) THEN 967 zbup(ji,jj) = pt_low(ji,jj) 968 zbdo(ji,jj) = pt_low(ji,jj) 969 !!clem ELSEIF( ptc(ji,jj) /= 0._wp .AND. pt_low(ji,jj) == 0._wp ) THEN 970 ELSEIF( ptc(ji,jj) > epsi20 .AND. pt_low(ji,jj) <= epsi20 ) THEN 971 zbup(ji,jj) = ptc(ji,jj) 972 zbdo(ji,jj) = ptc(ji,jj) 973 ELSE 974 zbup(ji,jj) = MAX( ptc(ji,jj) , pt_low(ji,jj) ) 975 zbdo(ji,jj) = MIN( ptc(ji,jj) , pt_low(ji,jj) ) 789 976 ENDIF 790 977 END DO … … 795 982 DO ji = fs_2, fs_jpim1 ! vector opt. 796 983 ! 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 984 !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 985 !zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), zbdo(ji ,jj-1), zbdo(ji ,jj+1) ) 986 ! 987 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 988 & zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1) ) 989 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), zbdo(ji ,jj-1), zbdo(ji ,jj+1), & 990 & zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1) ) 991 ! 992 zpos = MAX( 0., pfu_ho(ji-1,jj) ) - MIN( 0., pfu_ho(ji ,jj) ) & ! positive/negative part of the flux 993 & + MAX( 0., pfv_ho(ji,jj-1) ) - MIN( 0., pfv_ho(ji,jj ) ) 802 994 zneg = MAX( 0., pfu_ho(ji ,jj) ) - MIN( 0., pfu_ho(ji-1,jj) ) & 803 995 & + MAX( 0., pfv_ho(ji,jj ) ) - MIN( 0., pfv_ho(ji,jj-1) ) 804 996 ! 805 997 ! ! 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 998 ! zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / ( zpos + zsml ) * e1e2t(ji,jj) * z1_dt 999 ! zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / ( zneg + zsml ) * e1e2t(ji,jj) * z1_dt 1000 IF( zpos >= epsi10 ) THEN ; zbetup(ji,jj) = ( zup - pt_low(ji,jj) ) / zpos * e1e2t(ji,jj) * z1_dt 1001 ELSE ; zbetup(ji,jj) = 0. 812 1002 ENDIF 813 1003 ! 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 1004 IF( zneg >= epsi10 ) THEN ; zbetdo(ji,jj) = ( pt_low(ji,jj) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 1005 ELSE ; zbetdo(ji,jj) = 0. 818 1006 ENDIF 1007 ! 1008 IF( zbetdo(ji,jj) < 0._wp ) zbetdo(ji,jj)=0. 1009 IF( zbetup(ji,jj) < 0._wp ) zbetup(ji,jj)=0. 1010 ! 819 1011 ! 820 1012 END DO … … 830 1022 zcu = 0.5 + SIGN( 0.5 , pfu_ho(ji,jj) ) 831 1023 ! 832 pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_ ups(ji,jj)1024 pfu_ho(ji,jj) = pfu_ho(ji,jj) * ( zcu * zau + ( 1._wp - zcu ) * zbu ) + pfu_low(ji,jj) 833 1025 END DO 834 1026 END DO … … 840 1032 zcv = 0.5 + SIGN( 0.5 , pfv_ho(ji,jj) ) 841 1033 ! 842 pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_ ups(ji,jj)1034 pfv_ho(ji,jj) = pfv_ho(ji,jj) * ( zcv * zav + ( 1._wp - zcv ) * zbv ) + pfv_low(ji,jj) 843 1035 END DO 844 1036 END DO
Note: See TracChangeset
for help on using the changeset viewer.
