Changeset 10439
- Timestamp:
- 2018-12-21T12:18:00+01:00 (6 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/ICE/icedyn_adv_umx.F90 (modified) (54 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icedyn_adv_umx.F90
r10425 r10439 34 34 REAL(wp) :: z1_6 = 1._wp / 6._wp ! =1/6 35 35 REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120 36 37 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: amaxu, amaxv38 36 39 ! advect H all the way (and get V=H*A at the end)40 LOGICAL :: ll_thickness = .FALSE.41 42 ! look for 9 points around in nonosc limiter43 LOGICAL :: ll_9points = .FALSE. ! false=better h?44 45 ! use HgradU in nonosc limiter46 LOGICAL :: ll_HgradU = .TRUE. ! no effect?47 48 37 ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme 49 LOGICAL :: ll_neg = .TRUE. ! keep TRUE 50 51 ! limit the fluxes 52 LOGICAL :: ll_zeroup1 = .FALSE. ! false ok if Hbig otherwise needed for 2D sinon on a des valeurs de H trop fortes !! 53 LOGICAL :: ll_zeroup2 = .FALSE. ! false ok for 1D, 2D, 3D 54 LOGICAL :: ll_zeroup4 = .FALSE. ! false ok for 1D, 2D, 3D 55 LOGICAL :: ll_zeroup5 = .FALSE. ! false ok for 1D, 2D 56 57 ! fluxes that are limited are u*H, then (u*H)*(ua/u) is used for V (only for nonosc) 58 LOGICAL :: ll_clem = .TRUE. ! simpler than rachid and works 59 60 ! First advect H as H*=Hdiv(u), then use H* for H(n+1)=H(n)-div(uH*) 61 LOGICAL :: ll_gurvan = .FALSE. ! must be false for 1D case !! 62 63 ! First guess as div(uH) (-true-) or Hdiv(u)+ugradH (-false-) 64 LOGICAL :: ll_1stguess_clem = .FALSE. ! better negative values but less good h 65 66 ! advect (or not) open water. If not, retrieve it from advection of A 67 LOGICAL :: ll_ADVopw = .FALSE. ! 38 LOGICAL :: ll_neg = .TRUE. 68 39 69 40 ! alternate directions for upstream … … 74 45 75 46 ! prelimiter: use it to avoid overshoot in H 76 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? 77 LOGICAL :: ll_prelimiter_devore = .FALSE. ! from: Devore eq. 11 (worth than zalesak) 78 79 ! iterate on the limiter (only nonosc) 80 LOGICAL :: ll_limiter_it2 = .FALSE. 81 47 LOGICAL :: ll_prelimiter_zalesak = .TRUE. ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D 48 82 49 83 50 !! * Substitutions … … 120 87 REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers 121 88 REAL(wp) :: zdt 122 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! send zcflnow and receive zcflprv 123 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx 124 REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2 125 126 127 128 REAL(wp), DIMENSION(jpi,jpj) :: zcu_box, zcv_box 89 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! send zcflnow and receive zcflprv 90 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box 91 REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2 129 92 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zua_ho, zva_ho 130 93 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_ai, z1_aip 131 94 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhvar 132 133 95 !!---------------------------------------------------------------------- 134 96 ! 135 97 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 136 !137 98 ! 138 99 ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! … … 170 131 END DO 171 132 172 IF( ll_zeroup2 ) THEN173 IF(.NOT. ALLOCATED(amaxu)) ALLOCATE(amaxu (jpi,jpj,jpl))174 IF(.NOT. ALLOCATED(amaxv)) ALLOCATE(amaxv (jpi,jpj,jpl))175 ENDIF176 133 !---------------! 177 134 !== advection ==! … … 179 136 DO jt = 1, icycle 180 137 181 !!$ IF( ll_ADVopw ) THEN 182 !!$ zamsk = 1._wp 183 !!$ CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:), pato_i(:,:) ) ! Open water area 184 !!$ zamsk = 0._wp 185 !!$ ELSE 186 zati1(:,:) = SUM( pa_i(:,:,:), dim=3 ) 187 !!$ ENDIF 138 ! record at_i before advection (for open water) 139 zati1(:,:) = SUM( pa_i(:,:,:), dim=3 ) 188 140 141 ! inverse of A and Ap 189 142 WHERE( pa_i(:,:,:) >= epsi20 ) ; z1_ai(:,:,:) = 1._wp / pa_i(:,:,:) 190 143 ELSEWHERE ; z1_ai(:,:,:) = 0. 191 144 END WHERE 192 !193 145 WHERE( pa_ip(:,:,:) >= epsi20 ) ; z1_aip(:,:,:) = 1._wp / pa_ip(:,:,:) 194 146 ELSEWHERE ; z1_aip(:,:,:) = 0. 195 147 END WHERE 196 148 ! 197 IF( ll_zeroup2 ) THEN 198 DO jl = 1, jpl 199 DO jj = 2, jpjm1 200 DO ji = fs_2, fs_jpim1 201 amaxu(ji,jj,jl)=MAX( pa_i(ji,jj,jl), pa_i(ji,jj-1,jl), pa_i(ji,jj+1,jl), & 202 & pa_i(ji+1,jj,jl), pa_i(ji+1,jj-1,jl), pa_i(ji+1,jj+1,jl) ) 203 amaxv(ji,jj,jl)=MAX( pa_i(ji,jj,jl), pa_i(ji-1,jj,jl), pa_i(ji+1,jj,jl), & 204 & pa_i(ji,jj+1,jl), pa_i(ji-1,jj+1,jl), pa_i(ji+1,jj+1,jl) ) 205 END DO 206 END DO 207 END DO 208 CALL lbc_lnk_multi('icedyn_adv_umx', amaxu, 'T', 1., amaxv, 'T', 1.) 209 ENDIF 210 ! 149 ! set u*a=u for advection of A only 211 150 DO jl = 1, jpl 212 151 zua_ho(:,:,jl) = zudy(:,:) … … 218 157 zamsk = 0._wp 219 158 ! 220 IF( ll_thickness ) THEN221 DO jl = 1, jpl222 zua_ho(:,:,jl) = zudy(:,:)223 zva_ho(:,:,jl) = zvdx(:,:)224 END DO225 ENDIF226 !227 159 zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 228 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_i ) ! Ice volume 229 IF( ll_thickness ) pv_i(:,:,:) = zhvar(:,:,:) * pa_i(:,:,:) 160 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_i ) ! Ice volume 230 161 ! 231 162 zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 232 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_s ) ! Snw volume 233 IF( ll_thickness ) pv_s(:,:,:) = zhvar(:,:,:) * pa_i(:,:,:) 163 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_s ) ! Snw volume 234 164 ! 235 165 zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:) 236 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, psv_i )! Salt content166 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, psv_i ) ! Salt content 237 167 ! 238 168 zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:) 239 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, poa_i )! Age content169 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i ) ! Age content 240 170 ! 241 171 DO jk = 1, nlay_i 242 172 zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:) 243 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) ) ! Ice heat content173 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) ) ! Ice heat content 244 174 END DO 245 175 ! 246 176 DO jk = 1, nlay_s 247 177 zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:) 248 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) ) ! Snw heat content249 END DO 250 178 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) ) ! Snw heat content 179 END DO 180 ! 251 181 IF ( ln_pnd_H12 ) THEN 252 ! 182 ! set u*a=u for advection of Ap only 253 183 DO jl = 1, jpl 254 184 zua_ho(:,:,jl) = zudy(:,:) … … 261 191 ! 262 192 zhvar(:,:,:) = pv_ip(:,:,:) * z1_ai(:,:,:) 263 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, zhvar, pv_ip )! mp volume193 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_ip ) ! mp volume 264 194 ENDIF 265 195 ! 266 ! 267 !!$ IF( .NOT. ll_ADVopw ) THEN 268 zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 269 DO jj = 2, jpjm1 270 DO ji = fs_2, fs_jpim1 271 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & ! Open water area 272 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) )*r1_e1e2t(ji,jj)*zdt 273 END DO 274 END DO 275 CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T', 1. ) 276 !!$ ENDIF 196 zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 197 DO jj = 2, jpjm1 198 DO ji = fs_2, fs_jpim1 199 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & ! Open water area 200 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 201 END DO 202 END DO 203 CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T', 1. ) 277 204 ! 278 205 END DO … … 294 221 !! ** Action : - pt the after advective tracer 295 222 !!---------------------------------------------------------------------- 296 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)297 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)298 INTEGER , INTENT(in ) :: jt ! number of sub-iteration299 INTEGER , INTENT(in ) :: kt ! number of iteration300 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step301 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2302 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 or u*a*e2u303 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity304 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pt ! tracer field305 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: ptc ! tracer content field223 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 224 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 225 INTEGER , INTENT(in ) :: jt ! number of sub-iteration 226 INTEGER , INTENT(in ) :: kt ! number of iteration 227 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 228 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2 229 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 or u*a*e2u 230 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity 231 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pt ! tracer field 232 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: ptc ! tracer content field 306 233 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes 307 234 ! … … 310 237 INTEGER :: kn_limiter = 1 ! 1=nonosc ; 2=superbee ; 3=h3(rachid) 311 238 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zfu_ho , zfv_ho , zt_u, zt_v, zpt 312 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zfu_ups, zfv_ups, zt_ups ! only for nonosc239 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zfu_ups, zfv_ups, zt_ups 313 240 !!---------------------------------------------------------------------- 314 241 ! 315 242 ! upstream (_ups) advection with initial mass fluxes 316 243 ! --------------------------------------------------- 317 318 IF( ll_gurvan .AND. pamsk==0. ) THEN 319 DO jl = 1, jpl 320 DO jj = 2, jpjm1 321 DO ji = fs_2, fs_jpim1 322 pt(ji,jj,jl) = ( pt (ji,jj,jl) + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) & 323 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) ) & 324 & * tmask(ji,jj,1) 325 END DO 326 END DO 327 END DO 328 CALL lbc_lnk( 'icedyn_adv_umx', pt, 'T', 1. ) 329 ENDIF 330 331 332 IF( .NOT. ll_upsxy ) THEN 333 334 ! fluxes in both x-y directions 244 ! 245 IF( .NOT. ll_upsxy ) THEN !** no alternate directions **! 335 246 DO jl = 1, jpl 336 247 DO jj = 1, jpjm1 337 248 DO ji = 1, fs_jpim1 338 IF( ll_clem ) THEN 339 zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 340 zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 341 ELSE 342 zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * pt(ji+1,jj,jl) 343 zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj+1,jl) 344 ENDIF 345 END DO 346 END DO 347 END DO 348 349 ELSE 249 zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 250 zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 251 END DO 252 END DO 253 END DO 254 255 ELSE !** alternate directions **! 350 256 ! 351 257 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 352 ! flux in x-direction353 DO jl = 1, jpl 258 ! 259 DO jl = 1, jpl !-- flux in x-direction 354 260 DO jj = 1, jpjm1 355 261 DO ji = 1, fs_jpim1 356 IF( ll_clem ) THEN 357 zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 358 ELSE 359 zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * pt(ji+1,jj,jl) 360 ENDIF 262 zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 361 263 END DO 362 264 END DO 363 265 END DO 364 365 ! first guess of tracer content from u-flux 366 DO jl = 1, jpl 266 ! 267 DO jl = 1, jpl !-- first guess of tracer from u-flux 367 268 DO jj = 2, jpjm1 368 DO ji = fs_2, fs_jpim1 ! vector opt. 369 IF( ll_clem ) THEN 370 IF( ll_gurvan ) THEN 371 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 372 ELSE 373 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 374 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 375 & ) * tmask(ji,jj,1) 376 ENDIF 377 ELSE 378 zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) & 379 & * tmask(ji,jj,1) 380 ENDIF 381 !! IF( ji==26 .AND. jj==86) THEN 382 !! WRITE(numout,*) '************************' 383 !! WRITE(numout,*) 'zpt upstream',zpt(ji,jj) 384 !! ENDIF 269 DO ji = fs_2, fs_jpim1 270 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 271 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 272 & ) * tmask(ji,jj,1) 385 273 END DO 386 274 END DO … … 388 276 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 389 277 ! 390 ! flux in y-direction 391 DO jl = 1, jpl 278 DO jl = 1, jpl !-- flux in y-direction 392 279 DO jj = 1, jpjm1 393 280 DO ji = 1, fs_jpim1 394 IF( ll_clem ) THEN 395 zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl) 396 ELSE 397 zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * zpt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * zpt(ji,jj+1,jl) 398 ENDIF 281 zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl) 399 282 END DO 400 283 END DO 401 284 END DO 402 !285 ! 403 286 ELSE !== even ice time step: adv_y then adv_x ==! 404 ! flux in y-direction405 DO jl = 1, jpl 287 ! 288 DO jl = 1, jpl !-- flux in y-direction 406 289 DO jj = 1, jpjm1 407 290 DO ji = 1, fs_jpim1 408 IF( ll_clem ) THEN 409 zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 410 ELSE 411 zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj,jl) + MIN( pvc(ji,jj,jl), 0._wp ) * pt(ji,jj+1,jl) 412 ENDIF 291 zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 413 292 END DO 414 293 END DO 415 294 END DO 416 417 ! first guess of tracer content from v-flux 418 DO jl = 1, jpl 295 ! 296 DO jl = 1, jpl !-- first guess of tracer from v-flux 419 297 DO jj = 2, jpjm1 420 DO ji = fs_2, fs_jpim1 ! vector opt. 421 IF( ll_clem ) THEN 422 IF( ll_gurvan ) THEN 423 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 424 ELSE 425 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 426 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 427 & * tmask(ji,jj,1) 428 ENDIF 429 ELSE 430 zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) & 431 & * tmask(ji,jj,1) 432 ENDIF 433 !! IF( ji==26 .AND. jj==86) THEN 434 !! WRITE(numout,*) '************************' 435 !! WRITE(numout,*) 'zpt upstream',zpt(ji,jj) 436 !! ENDIF 298 DO ji = fs_2, fs_jpim1 299 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 300 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 301 & * tmask(ji,jj,1) 437 302 END DO 438 303 END DO … … 440 305 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 441 306 ! 442 ! flux in x-direction 443 DO jl = 1, jpl 307 DO jl = 1, jpl !-- flux in x-direction 444 308 DO jj = 1, jpjm1 445 309 DO ji = 1, fs_jpim1 446 IF( ll_clem ) THEN 447 zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl) 448 ELSE 449 zfu_ups(ji,jj,jl) = MAX( puc(ji,jj,jl), 0._wp ) * zpt(ji,jj,jl) + MIN( puc(ji,jj,jl), 0._wp ) * zpt(ji+1,jj,jl) 450 ENDIF 310 zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl) 451 311 END DO 452 312 END DO … … 456 316 457 317 ENDIF 458 459 IF( ll_clem .AND. kn_limiter /= 1 ) & 460 & CALL ctl_stop( 'STOP', 'icedyn_adv_umx: ll_clem incompatible with limiters other than nonosc' ) 461 462 IF( ll_zeroup2 ) THEN 463 DO jl = 1, jpl 464 DO jj = 1, jpjm1 465 DO ji = 1, fs_jpim1 ! vector opt. 466 IF( amaxu(ji,jj,jl) == 0._wp ) zfu_ups(ji,jj,jl) = 0._wp 467 IF( amaxv(ji,jj,jl) == 0._wp ) zfv_ups(ji,jj,jl) = 0._wp 468 END DO 469 END DO 470 END DO 471 ENDIF 472 473 ! guess after content field with upstream scheme 474 DO jl = 1, jpl 318 ! 319 DO jl = 1, jpl !-- after tracer with upstream scheme 475 320 DO jj = 2, jpjm1 476 321 DO ji = fs_2, fs_jpim1 477 322 ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj ,jl) & 478 323 & + zfv_ups(ji,jj,jl) - zfv_ups(ji ,jj-1,jl) ) * r1_e1e2t(ji,jj) 479 IF( ll_clem ) THEN 480 IF( ll_gurvan ) THEN 481 zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) 482 ELSE 483 zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra + ( pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) & 484 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) ) & 485 & * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 486 ENDIF 487 ELSE 488 zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) 489 ENDIF 490 !! IF( ji==26 .AND. jj==86) THEN 491 !! WRITE(numout,*) '**************************' 492 !! WRITE(numout,*) 'zt upstream',zt_ups(ji,jj) 493 !! ENDIF 324 zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra + ( pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) & 325 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) ) & 326 & * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 494 327 END DO 495 328 END DO … … 512 345 ! 513 346 END SELECT 514 515 IF( ll_thickness ) THEN 516 ! final trend with corrected fluxes 517 ! ------------------------------------ 518 DO jl = 1, jpl 519 DO jj = 2, jpjm1 520 DO ji = fs_2, fs_jpim1 521 IF( ll_gurvan ) THEN 522 ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) 347 ! 348 ! --ho --ho 349 ! new fluxes = u*H * u*a / u 350 ! ---------------------------- 351 IF( pamsk == 0. ) THEN 352 DO jl = 1, jpl 353 DO jj = 1, jpjm1 354 DO ji = 1, fs_jpim1 355 IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN 356 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 357 zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 523 358 ELSE 524 ztra = ( - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) & 525 & + ( pt(ji,jj,jl) * ( pu(ji,jj) - pu(ji-1,jj) ) * (1.-pamsk) ) & 526 & + ( pt(ji,jj,jl) * ( pv(ji,jj) - pv(ji,jj-1) ) * (1.-pamsk) ) ) * r1_e1e2t(ji,jj) 359 zfu_ho (ji,jj,jl) = 0._wp 360 zfu_ups(ji,jj,jl) = 0._wp 527 361 ENDIF 528 pt(ji,jj,jl) = ( pt(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) 529 530 IF( pt(ji,jj,jl) < -epsi20 ) THEN 531 WRITE(numout,*) 'T<0 ',pt(ji,jj,jl) 362 ! 363 IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN 364 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 365 zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 366 ELSE 367 zfv_ho (ji,jj,jl) = 0._wp 368 zfv_ups(ji,jj,jl) = 0._wp 532 369 ENDIF 533 534 IF( pt(ji,jj,jl) < 0._wp .AND. pt(ji,jj,jl) >= -epsi20 ) pt(ji,jj,jl) = 0._wp 535 536 !! IF( ji==26 .AND. jj==86) THEN 537 !! WRITE(numout,*) 'zt high order',pt(ji,jj) 538 !! ENDIF 539 END DO 540 END DO 541 END DO 542 CALL lbc_lnk( 'icedyn_adv_umx', pt, 'T', 1. ) 370 END DO 371 END DO 372 END DO 543 373 ENDIF 544 545 ! Rachid trick 546 ! ------------ 547 IF( ll_clem ) THEN 548 IF( pamsk == 0. ) THEN 549 DO jl = 1, jpl 550 DO jj = 1, jpjm1 551 DO ji = 1, fs_jpim1 552 IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN 553 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 554 zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 555 ELSE 556 zfu_ho (ji,jj,jl) = 0._wp 557 zfu_ups(ji,jj,jl) = 0._wp 558 ENDIF 559 ! 560 IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN 561 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 562 zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 563 ELSE 564 zfv_ho (ji,jj,jl) = 0._wp 565 zfv_ups(ji,jj,jl) = 0._wp 566 ENDIF 567 END DO 568 END DO 569 END DO 570 ENDIF 571 ENDIF 572 573 IF( ll_zeroup5 ) THEN 574 DO jl = 1, jpl 575 DO jj = 2, jpjm1 576 DO ji = 2, fs_jpim1 ! vector opt. 577 zpt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 578 & - ( zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 579 IF( zpt(ji,jj,jl) < 0. ) THEN 580 zfu_ho(ji ,jj,jl) = zfu_ups(ji ,jj,jl) 581 zfu_ho(ji-1,jj,jl) = zfu_ups(ji-1,jj,jl) 582 zfv_ho(ji ,jj,jl) = zfv_ups(ji ,jj,jl) 583 zfv_ho(ji,jj-1,jl) = zfv_ups(ji,jj-1,jl) 584 ENDIF 585 END DO 586 END DO 587 END DO 588 CALL lbc_lnk_multi( 'icedyn_adv_umx', zfu_ho, 'U', -1., zfv_ho, 'V', -1. ) 589 ENDIF 590 591 ! output high order fluxes u*a 592 ! ---------------------------- 374 ! 375 ! in case of advection of A: output u*a(ho) 376 ! ----------------------------------------- 593 377 IF( PRESENT( pua_ho ) ) THEN 594 378 DO jl = 1, jpl … … 601 385 END DO 602 386 ENDIF 603 604 605 IF( .NOT.ll_thickness ) THEN 606 ! final trend with corrected fluxes 607 ! ------------------------------------ 608 DO jl = 1, jpl 609 DO jj = 2, jpjm1 610 DO ji = fs_2, fs_jpim1 611 ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) * pdt 612 613 ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra ) * tmask(ji,jj,1) 614 615 !! IF( ji==26 .AND. jj==86) THEN 616 !! WRITE(numout,*) 'ztc high order',ptc(ji,jj) 617 !! ENDIF 618 619 END DO 620 END DO 621 END DO 622 CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1. ) 623 ENDIF 624 387 ! 388 ! final trend with corrected fluxes 389 ! --------------------------------- 390 DO jl = 1, jpl 391 DO jj = 2, jpjm1 392 DO ji = fs_2, fs_jpim1 393 ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) * pdt 394 395 ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra ) * tmask(ji,jj,1) 396 END DO 397 END DO 398 END DO 399 CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1. ) 625 400 ! 626 401 END SUBROUTINE adv_umx 402 627 403 628 404 SUBROUTINE cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, & … … 638 414 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 639 415 !!---------------------------------------------------------------------- 640 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)641 INTEGER , INTENT(in ) :: kn_limiter ! limiter642 INTEGER , INTENT(in ) :: jt ! number of sub-iteration643 INTEGER , INTENT(in ) :: kt ! number of iteration644 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step645 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields646 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu, pv ! 2 ice velocity components647 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc, pvc ! 2 ice velocity * A components648 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ptc ! tracer content at before time step416 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 417 INTEGER , INTENT(in ) :: kn_limiter ! limiter 418 INTEGER , INTENT(in ) :: jt ! number of sub-iteration 419 INTEGER , INTENT(in ) :: kt ! number of iteration 420 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 421 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields 422 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu, pv ! 2 ice velocity components 423 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc, pvc ! 2 ice velocity * A components 424 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ptc ! tracer content at before time step 649 425 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes 650 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt_ups ! upstream guess of tracer content651 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pfu_ups, pfv_ups ! upstream fluxes426 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt_ups ! upstream guess of tracer content 427 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pfu_ups, pfv_ups ! upstream fluxes 652 428 ! 653 429 INTEGER :: ji, jj, jl ! dummy loop indices 654 LOGICAL :: ll_xy = .TRUE.655 430 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zzt 656 431 !!---------------------------------------------------------------------- 657 432 ! 658 IF( .NOT.ll_ xy ) THEN !-- no alternate directions --!433 IF( .NOT.ll_hoxy ) THEN !** no alternate directions **! 659 434 ! 660 435 DO jl = 1, jpl 661 436 DO jj = 1, jpjm1 662 437 DO ji = 1, fs_jpim1 663 IF( ll_clem ) THEN 664 pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 665 pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 666 ELSE 667 pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 668 pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 669 ENDIF 438 pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 439 pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 670 440 END DO 671 441 END DO 672 442 END DO 673 443 IF ( kn_limiter == 1 ) THEN 674 IF( ll_clem ) THEN 675 CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 676 ELSE 677 CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 678 ENDIF 444 CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 679 445 ELSEIF( kn_limiter == 2 ) THEN 680 CALL limiter_x( pdt, pu, p uc, pt, pfu_ho )681 CALL limiter_y( pdt, pv, p vc, pt, pfv_ho )446 CALL limiter_x( pdt, pu, pt, pfu_ho ) 447 CALL limiter_y( pdt, pv, pt, pfv_ho ) 682 448 ELSEIF( kn_limiter == 3 ) THEN 683 CALL limiter_x( pdt, pu, p uc, pt, pfu_ho, pfu_ups )684 CALL limiter_y( pdt, pv, p vc, pt, pfv_ho, pfv_ups )449 CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 450 CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 685 451 ENDIF 686 452 ! 687 ELSE !-- alternate directions --!453 ELSE !** alternate directions **! 688 454 ! 689 455 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 690 456 ! 691 ! flux in x-direction 692 DO jl = 1, jpl 457 DO jl = 1, jpl !-- flux in x-direction 693 458 DO jj = 1, jpjm1 694 459 DO ji = 1, fs_jpim1 695 IF( ll_clem ) THEN 696 pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 697 ELSE 698 pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 699 ENDIF 460 pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 700 461 END DO 701 462 END DO 702 463 END DO 703 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho ) 704 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, puc, pt, pfu_ho, pfu_ups ) 705 706 ! first guess of tracer content from u-flux 707 DO jl = 1, jpl 464 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, pt, pfu_ho ) 465 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 466 467 DO jl = 1, jpl !-- first guess of tracer from u-flux 708 468 DO jj = 2, jpjm1 709 DO ji = fs_2, fs_jpim1 ! vector opt. 710 IF( ll_clem ) THEN 711 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 712 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 713 & * tmask(ji,jj,1) 714 ELSE 715 zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 716 ENDIF 469 DO ji = fs_2, fs_jpim1 470 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 471 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 472 & * tmask(ji,jj,1) 717 473 END DO 718 474 END DO … … 720 476 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 721 477 722 ! flux in y-direction 723 DO jl = 1, jpl 478 DO jl = 1, jpl !-- flux in y-direction 724 479 DO jj = 1, jpjm1 725 480 DO ji = 1, fs_jpim1 726 IF( ll_clem ) THEN 727 pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) ) 728 ELSE 729 pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) ) 730 ENDIF 481 pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji,jj+1,jl) ) 731 482 END DO 732 483 END DO 733 484 END DO 734 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, p vc, pt, pfv_ho )735 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, p vc, pt, pfv_ho, pfv_ups )485 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pt, pfv_ho ) 486 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 736 487 737 488 ELSE !== even ice time step: adv_y then adv_x ==! 738 489 ! 739 ! flux in y-direction 740 DO jl = 1, jpl 490 DO jl = 1, jpl !-- flux in y-direction 741 491 DO jj = 1, jpjm1 742 492 DO ji = 1, fs_jpim1 743 IF( ll_clem ) THEN 744 pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 745 ELSE 746 pfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj,jl) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 747 ENDIF 493 pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 748 494 END DO 749 495 END DO 750 496 END DO 751 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, p vc, pt, pfv_ho )752 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, p vc, pt, pfv_ho, pfv_ups )497 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pt, pfv_ho ) 498 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 753 499 ! 754 ! first guess of tracer content from v-flux 755 DO jl = 1, jpl 500 DO jl = 1, jpl !-- first guess of tracer from v-flux 756 501 DO jj = 2, jpjm1 757 DO ji = fs_2, fs_jpim1 ! vector opt. 758 IF( ll_clem ) THEN 759 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 760 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 761 & * tmask(ji,jj,1) 762 ELSE 763 zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 764 ENDIF 502 DO ji = fs_2, fs_jpim1 503 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 504 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 505 & * tmask(ji,jj,1) 765 506 END DO 766 507 END DO … … 768 509 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 769 510 ! 770 ! flux in x-direction 771 DO jl = 1, jpl 511 DO jl = 1, jpl !-- flux in x-direction 772 512 DO jj = 1, jpjm1 773 513 DO ji = 1, fs_jpim1 774 IF( ll_clem ) THEN 775 pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) ) 776 ELSE 777 pfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj,jl) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) ) 778 ENDIF 514 pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( zzt(ji,jj,jl) + zzt(ji+1,jj,jl) ) 779 515 END DO 780 516 END DO 781 517 END DO 782 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, p uc, pt, pfu_ho )783 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, p uc, pt, pfu_ho, pfu_ups )518 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, pt, pfu_ho ) 519 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 784 520 785 521 ENDIF 786 IF( ll_clem ) THEN 787 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 ) 788 ELSE 789 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 ) 790 ENDIF 522 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 ) 791 523 792 524 ENDIF … … 807 539 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 808 540 !!---------------------------------------------------------------------- 809 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0)810 INTEGER , INTENT(in ) :: kn_limiter ! limiter811 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2)812 INTEGER , INTENT(in ) :: jt ! number of sub-iteration813 INTEGER , INTENT(in ) :: kt ! number of iteration814 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step815 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields816 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu, pv ! 2 ice velocity components817 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc, pvc ! 2 ice velocity * A components818 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity819 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ptc ! tracer content at before time step541 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 542 INTEGER , INTENT(in ) :: kn_limiter ! limiter 543 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 544 INTEGER , INTENT(in ) :: jt ! number of sub-iteration 545 INTEGER , INTENT(in ) :: kt ! number of iteration 546 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 547 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields 548 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu, pv ! 2 ice velocity components 549 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc, pvc ! 2 ice velocity * A components 550 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity 551 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ptc ! tracer content at before time step 820 552 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_u, pt_v ! tracer at u- and v-points 821 553 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes 822 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt_ups ! upstream guess of tracer content823 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pfu_ups, pfv_ups ! upstream fluxes554 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt_ups ! upstream guess of tracer content 555 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pfu_ups, pfv_ups ! upstream fluxes 824 556 ! 825 557 INTEGER :: ji, jj, jl ! dummy loop indices … … 831 563 ! 832 564 ! !-- ultimate interpolation of pt at u-point --! 833 CALL ultimate_x( kn_umx, pdt, pt, pu, p uc, pt_u, pfu_ho )565 CALL ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 834 566 ! !-- limiter in x --! 835 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, p uc, pt, pfu_ho )836 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, p uc, pt, pfu_ho, pfu_ups )567 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, pt, pfu_ho ) 568 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 837 569 ! !-- advective form update in zzt --! 838 839 IF( ll_1stguess_clem ) THEN 840 841 ! first guess of tracer content from u-flux 842 DO jl = 1, jpl 843 DO jj = 2, jpjm1 844 DO ji = fs_2, fs_jpim1 ! vector opt. 845 IF( ll_clem ) THEN 846 IF( ll_gurvan ) THEN 847 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 848 ELSE 849 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 850 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 851 & ) * tmask(ji,jj,1) 852 ENDIF 853 ELSE 854 zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 855 ENDIF 856 END DO 857 END DO 858 END DO 859 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 860 861 ELSE 862 863 DO jl = 1, jpl 864 DO jj = 2, jpjm1 865 DO ji = fs_2, fs_jpim1 ! vector opt. 866 IF( ll_gurvan ) THEN 867 zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj) & 868 & - pt (ji,jj,jl) * pdt * ( pu (ji,jj) - pu (ji-1,jj) ) * r1_e1e2t(ji,jj) 869 ELSE 870 zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj) & 871 & - pt (ji,jj,jl) * pdt * ( pu (ji,jj) - pu (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk 872 ENDIF 873 zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 874 END DO 875 END DO 876 END DO 877 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 878 ENDIF 570 DO jl = 1, jpl 571 DO jj = 2, jpjm1 572 DO ji = fs_2, fs_jpim1 573 zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj) & 574 & - pt (ji,jj,jl) * pdt * ( pu (ji,jj) - pu (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk 575 zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 576 END DO 577 END DO 578 END DO 579 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 879 580 ! 880 581 ! !-- ultimate interpolation of pt at v-point --! 881 582 IF( ll_hoxy ) THEN 882 CALL ultimate_y( kn_umx, pdt, zzt, pv, p vc, pt_v, pfv_ho )583 CALL ultimate_y( kn_umx, pdt, zzt, pv, pt_v, pfv_ho ) 883 584 ELSE 884 CALL ultimate_y( kn_umx, pdt, pt, pv, p vc, pt_v, pfv_ho )585 CALL ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 885 586 ENDIF 886 587 ! !-- limiter in y --! 887 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, p vc, pt, pfv_ho )888 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, p vc, pt, pfv_ho, pfv_ups )588 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pt, pfv_ho ) 589 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 889 590 ! 890 591 ! … … 892 593 ! 893 594 ! !-- ultimate interpolation of pt at v-point --! 894 CALL ultimate_y( kn_umx, pdt, pt, pv, p vc, pt_v, pfv_ho )595 CALL ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 895 596 ! !-- limiter in y --! 896 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, p vc, pt, pfv_ho )897 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, p vc, pt, pfv_ho, pfv_ups )597 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pt, pfv_ho ) 598 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 898 599 ! !-- advective form update in zzt --! 899 IF( ll_1stguess_clem ) THEN 900 901 ! first guess of tracer content from v-flux 902 DO jl = 1, jpl 903 DO jj = 2, jpjm1 904 DO ji = fs_2, fs_jpim1 ! vector opt. 905 IF( ll_clem ) THEN 906 IF( ll_gurvan ) THEN 907 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 908 ELSE 909 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 910 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 911 & ) * tmask(ji,jj,1) 912 ENDIF 913 ELSE 914 zzt(ji,jj,jl) = ( ptc(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) & 915 & * tmask(ji,jj,1) 916 ENDIF 917 END DO 918 END DO 919 END DO 920 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 921 922 ELSE 923 924 DO jl = 1, jpl 925 DO jj = 2, jpjm1 926 DO ji = fs_2, fs_jpim1 927 IF( ll_gurvan ) THEN 928 zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj) & 929 & - pt (ji,jj,jl) * pdt * ( pv (ji,jj) - pv (ji,jj-1) ) * r1_e1e2t(ji,jj) 930 ELSE 931 zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj) & 932 & - pt (ji,jj,jl) * pdt * ( pv (ji,jj) - pv (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk 933 ENDIF 934 zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 935 END DO 936 END DO 937 END DO 938 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 939 ENDIF 600 DO jl = 1, jpl 601 DO jj = 2, jpjm1 602 DO ji = fs_2, fs_jpim1 603 zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj) & 604 & - pt (ji,jj,jl) * pdt * ( pv (ji,jj) - pv (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk 605 zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 606 END DO 607 END DO 608 END DO 609 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 940 610 ! 941 611 ! !-- ultimate interpolation of pt at u-point --! 942 612 IF( ll_hoxy ) THEN 943 CALL ultimate_x( kn_umx, pdt, zzt, pu, p uc, pt_u, pfu_ho )613 CALL ultimate_x( kn_umx, pdt, zzt, pu, pt_u, pfu_ho ) 944 614 ELSE 945 CALL ultimate_x( kn_umx, pdt, pt, pu, p uc, pt_u, pfu_ho )615 CALL ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 946 616 ENDIF 947 617 ! !-- limiter in x --! 948 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, p uc, pt, pfu_ho )949 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, p uc, pt, pfu_ho, pfu_ups )618 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, pt, pfu_ho ) 619 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 950 620 ! 951 621 ! 952 622 ENDIF 953 623 954 955 624 IF( kn_limiter == 1 ) THEN 956 IF( .NOT. ll_limiter_it2 ) THEN 957 IF( ll_clem ) THEN 958 CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 959 ELSE 960 CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 961 ENDIF 962 ELSE 963 zzfu_ho(:,:,:) = pfu_ho(:,:,:) 964 zzfv_ho(:,:,:) = pfv_ho(:,:,:) 965 ! 1st iteration of nonosc (limit the flux with the upstream solution) 966 IF( ll_clem ) THEN 967 CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho ) 968 ELSE 969 CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, pt_ups, pfu_ups, pfv_ups, zzfu_ho, zzfv_ho ) 970 ENDIF 971 ! guess after content field with high order 972 DO jl = 1, jpl 973 DO jj = 2, jpjm1 974 DO ji = fs_2, fs_jpim1 975 ztra = - ( zzfu_ho(ji,jj,jl) - zzfu_ho(ji-1,jj,jl) + zzfv_ho(ji,jj,jl) - zzfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) 976 zzt(ji,jj,jl) = ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) 977 END DO 978 END DO 979 END DO 980 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 981 ! 2nd iteration of nonosc (limit the flux with the limited high order solution) 982 IF( ll_clem ) THEN 983 CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho ) 984 ELSE 985 CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, ptc, zzt, zzfu_ho, zzfv_ho, pfu_ho, pfv_ho ) 986 ENDIF 987 ENDIF 625 CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 988 626 ENDIF 989 627 ! … … 991 629 992 630 993 SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, p uc, pt_u, pfu_ho )631 SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 994 632 !!--------------------------------------------------------------------- 995 633 !! *** ROUTINE ultimate_x *** … … 1002 640 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 1003 641 !!---------------------------------------------------------------------- 1004 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 1005 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1006 REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pu ! ice i-velocity component 1007 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: puc ! ice i-velocity * A component 1008 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt ! tracer fields 642 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 643 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 644 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu ! ice i-velocity component 645 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields 1009 646 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_u ! tracer at u-point 1010 647 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho ! high order flux 1011 648 ! 1012 649 INTEGER :: ji, jj, jl ! dummy loop indices 1013 REAL(wp) :: zcu, zdx2, zdx4 ! - -650 REAL(wp) :: zcu, zdx2, zdx4 ! - - 1014 651 REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztu1, ztu2, ztu3, ztu4 1015 652 !!---------------------------------------------------------------------- … … 1076 713 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 1077 714 zdx2 = e1u(ji,jj) * e1u(ji,jj) 1078 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)715 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 1079 716 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 1080 717 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & … … 1092 729 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 1093 730 zdx2 = e1u(ji,jj) * e1u(ji,jj) 1094 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)731 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 1095 732 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 1096 733 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & … … 1108 745 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 1109 746 zdx2 = e1u(ji,jj) * e1u(ji,jj) 1110 !!rachidzdx2 = e1u(ji,jj) * e1t(ji,jj)747 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 1111 748 zdx4 = zdx2 * zdx2 1112 749 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & … … 1121 758 ! 1122 759 END SELECT 1123 ! !-- High order flux in i-direction --! 760 ! 761 ! if pt at u-point is negative then use the upstream value 762 ! this should not be necessary if a proper sea-ice mask is set in Ultimate 763 ! to degrade the order of the scheme when necessary (for ex. at the ice edge) 1124 764 IF( ll_neg ) THEN 1125 765 DO jl = 1, jpl … … 1134 774 END DO 1135 775 ENDIF 1136 776 ! !-- High order flux in i-direction --! 1137 777 DO jl = 1, jpl 1138 778 DO jj = 1, jpjm1 1139 779 DO ji = 1, fs_jpim1 ! vector opt. 1140 IF( ll_clem ) THEN 1141 pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 1142 ELSE 1143 pfu_ho(ji,jj,jl) = puc(ji,jj,jl) * pt_u(ji,jj,jl) 1144 ENDIF 780 pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 1145 781 END DO 1146 782 END DO … … 1150 786 1151 787 1152 SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, p vc, pt_v, pfv_ho )788 SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 1153 789 !!--------------------------------------------------------------------- 1154 790 !! *** ROUTINE ultimate_y *** … … 1161 797 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 1162 798 !!---------------------------------------------------------------------- 1163 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 1164 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1165 REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pv ! ice j-velocity component 1166 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pvc ! ice j-velocity*A component 1167 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt ! tracer fields 799 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 800 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 801 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pv ! ice j-velocity component 802 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields 1168 803 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_v ! tracer at v-point 1169 804 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfv_ho ! high order flux 1170 805 ! 1171 INTEGER :: ji, jj, jl ! dummy loop indices806 INTEGER :: ji, jj, jl ! dummy loop indices 1172 807 REAL(wp) :: zcv, zdy2, zdy4 ! - - 1173 808 REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztv1, ztv2, ztv3, ztv4 … … 1235 870 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1236 871 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1237 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)872 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 1238 873 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1239 874 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & … … 1250 885 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1251 886 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1252 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)887 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 1253 888 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1254 889 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & … … 1265 900 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1266 901 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1267 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)902 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 1268 903 zdy4 = zdy2 * zdy2 1269 904 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & … … 1278 913 ! 1279 914 END SELECT 1280 ! !-- High order flux in j-direction --! 915 ! 916 ! if pt at v-point is negative then use the upstream value 917 ! this should not be necessary if a proper sea-ice mask is set in Ultimate 918 ! to degrade the order of the scheme when necessary (for ex. at the ice edge) 1281 919 IF( ll_neg ) THEN 1282 920 DO jl = 1, jpl … … 1291 929 END DO 1292 930 ENDIF 1293 931 ! !-- High order flux in j-direction --! 1294 932 DO jl = 1, jpl 1295 933 DO jj = 1, jpjm1 1296 934 DO ji = 1, fs_jpim1 ! vector opt. 1297 IF( ll_clem ) THEN 1298 pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 1299 ELSE 1300 pfv_ho(ji,jj,jl) = pvc(ji,jj,jl) * pt_v(ji,jj,jl) 1301 ENDIF 935 pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 1302 936 END DO 1303 937 END DO … … 1307 941 1308 942 1309 SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ low, pfu_low, pfv_low, pfu_ho, pfv_ho )943 SUBROUTINE nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 1310 944 !!--------------------------------------------------------------------- 1311 945 !! *** ROUTINE nonosc *** … … 1315 949 !! 1316 950 !! ** Method : ... ??? 1317 !! warning : pt and pt_ lowmust be masked, but the boundaries951 !! warning : pt and pt_ups must be masked, but the boundaries 1318 952 !! conditions on the fluxes are not necessary zalezak (1979) 1319 953 !! drange (1995) multi-dimensional forward-in-time and upstream- 1320 954 !! in-space based differencing for fluid 1321 955 !!---------------------------------------------------------------------- 1322 REAL(wp) 1323 REAL(wp) 956 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 957 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1324 958 REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pu ! ice i-velocity => u*e2 1325 959 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: puc ! ice i-velocity *A => u*e2*a 1326 960 REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pv ! ice j-velocity => v*e1 1327 961 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pvc ! ice j-velocity *A => v*e1*a 1328 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: ptc, pt, pt_ low! before field & upstream guess of after field1329 REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ low, pfu_low! upstream flux962 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: ptc, pt, pt_ups ! before field & upstream guess of after field 963 REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ups, pfu_ups ! upstream flux 1330 964 REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ho, pfu_ho ! monotonic flux 1331 965 ! 1332 966 INTEGER :: ji, jj, jl ! dummy loop indices 1333 REAL(wp) :: zpos, zneg, zbig, zsml, z 1_dt, zpos2, zneg2! local scalars1334 REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, z up, zdo, zsign, zcoef! - -967 REAL(wp) :: zpos, zneg, zbig, zsml, zup, zdo, z1_dt ! local scalars 968 REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zsign, zcoef, zzt ! - - 1335 969 REAL(wp), DIMENSION(jpi,jpj ) :: zbup, zbdo 1336 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_ low, ztj_low, zzt970 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zbetup, zbetdo, zti_ups, ztj_ups 1337 971 !!---------------------------------------------------------------------- 1338 972 zbig = 1.e+40_wp 1339 973 zsml = epsi20 1340 1341 IF( ll_zeroup2 ) THEN1342 DO jl = 1, jpl1343 DO jj = 1, jpjm11344 DO ji = 1, fs_jpim1 ! vector opt.1345 IF( amaxu(ji,jj,jl) == 0._wp ) pfu_ho(ji,jj,jl) = 0._wp1346 IF( amaxv(ji,jj,jl) == 0._wp ) pfv_ho(ji,jj,jl) = 0._wp1347 END DO1348 END DO1349 END DO1350 ENDIF1351 1352 IF( ll_zeroup4 ) THEN1353 DO jl = 1, jpl1354 DO jj = 1, jpjm11355 DO ji = 1, fs_jpim1 ! vector opt.1356 IF( pfu_low(ji,jj,jl) == 0._wp ) pfu_ho(ji,jj,jl) = 0._wp1357 IF( pfv_low(ji,jj,jl) == 0._wp ) pfv_ho(ji,jj,jl) = 0._wp1358 END DO1359 END DO1360 END DO1361 ENDIF1362 1363 1364 IF( ll_zeroup1 ) THEN1365 DO jl = 1, jpl1366 DO jj = 2, jpjm11367 DO ji = fs_2, fs_jpim11368 IF( ll_gurvan ) THEN1369 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &1370 & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)1371 ELSE1372 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &1373 & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &1374 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &1375 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &1376 & ) * tmask(ji,jj,1)1377 ENDIF1378 IF( zzt(ji,jj,jl) < 0._wp ) THEN1379 pfu_ho(ji,jj,jl) = pfu_low(ji,jj,jl)1380 pfv_ho(ji,jj,jl) = pfv_low(ji,jj,jl)1381 WRITE(numout,*) '*** 1 negative high order zzt ***',ji,jj,zzt(ji,jj,jl)1382 ENDIF1383 !! IF( ji==26 .AND. jj==86) THEN1384 !! WRITE(numout,*) 'zzt high order',zzt(ji,jj)1385 !! WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt1386 !! WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt1387 !! WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj,jl)) * r1_e1e2t(ji,jj) * pdt1388 !! WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1,jl)) * r1_e1e2t(ji,jj) * pdt1389 !! ENDIF1390 IF( ll_gurvan ) THEN1391 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &1392 & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)1393 ELSE1394 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &1395 & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &1396 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &1397 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &1398 & ) * tmask(ji,jj,1)1399 ENDIF1400 IF( zzt(ji,jj,jl) < 0._wp ) THEN1401 pfu_ho(ji-1,jj,jl) = pfu_low(ji-1,jj,jl)1402 pfv_ho(ji,jj-1,jl) = pfv_low(ji,jj-1,jl)1403 WRITE(numout,*) '*** 2 negative high order zzt ***',ji,jj,zzt(ji,jj,jl)1404 ENDIF1405 IF( ll_gurvan ) THEN1406 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &1407 & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)1408 ELSE1409 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &1410 & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &1411 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &1412 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &1413 & ) * tmask(ji,jj,1)1414 ENDIF1415 IF( zzt(ji,jj,jl) < 0._wp ) THEN1416 WRITE(numout,*) '*** 3 negative high order zzt ***',ji,jj,zzt(ji,jj,jl)1417 ENDIF1418 END DO1419 END DO1420 END DO1421 CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )1422 ENDIF1423 1424 974 1425 975 ! antidiffusive flux : high order minus low order … … 1428 978 DO jj = 1, jpjm1 1429 979 DO ji = 1, fs_jpim1 ! vector opt. 1430 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ low(ji,jj,jl)1431 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ low(ji,jj,jl)980 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 981 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1432 982 END DO 1433 983 END DO … … 1441 991 ! | | | * | 1442 992 ! | | | | * 1443 ! t_ low: i-1 i i+1 i+2993 ! t_ups : i-1 i i+1 i+2 1444 994 IF( ll_prelimiter_zalesak ) THEN 1445 995 … … 1447 997 DO jj = 2, jpjm1 1448 998 DO ji = fs_2, fs_jpim1 1449 zti_ low(ji,jj,jl)= pt_low(ji+1,jj ,jl)1450 ztj_ low(ji,jj,jl)= pt_low(ji ,jj+1,jl)1451 END DO 1452 END DO 1453 END DO 1454 CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ low, 'T', 1., ztj_low, 'T', 1. )999 zti_ups(ji,jj,jl)= pt_ups(ji+1,jj ,jl) 1000 ztj_ups(ji,jj,jl)= pt_ups(ji ,jj+1,jl) 1001 END DO 1002 END DO 1003 END DO 1004 CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1., ztj_ups, 'T', 1. ) 1455 1005 1456 1006 !! this does not work ?? 1457 1007 !! DO jj = 2, jpjm1 1458 1008 !! DO ji = fs_2, fs_jpim1 1459 !! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_ low (ji+1,jj ) - pt_low(ji ,jj) ) .AND. &1460 !! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_ low (ji ,jj+1) - pt_low(ji ,jj) ) &1009 !! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_ups (ji+1,jj ) - pt_ups (ji ,jj) ) .AND. & 1010 !! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_ups (ji ,jj+1) - pt_ups (ji ,jj) ) & 1461 1011 !! & ) THEN 1462 !! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_ low(ji+1,jj ) - zti_low(ji ,jj) ) .AND. &1463 !! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_ low(ji,jj+1 ) - ztj_low(ji ,jj) ) &1012 !! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_ups(ji+1,jj ) - zti_ups(ji ,jj) ) .AND. & 1013 !! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_ups(ji,jj+1 ) - ztj_ups(ji ,jj) ) & 1464 1014 !! & ) THEN 1465 1015 !! pfu_ho(ji,jj) = 0. ; pfv_ho(ji,jj) = 0. 1466 1016 !! ENDIF 1467 !! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_ low (ji ,jj) - pt_low(ji-1,jj ) ) .AND. &1468 !! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_ low (ji ,jj) - pt_low(ji ,jj-1) ) &1017 !! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_ups (ji ,jj) - pt_ups (ji-1,jj ) ) .AND. & 1018 !! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_ups (ji ,jj) - pt_ups (ji ,jj-1) ) & 1469 1019 !! & ) THEN 1470 1020 !! pfu_ho(ji,jj) = 0. ; pfv_ho(ji,jj) = 0. … … 1477 1027 DO jj = 2, jpjm1 1478 1028 DO ji = fs_2, fs_jpim1 1479 IF ( pfu_ho(ji,jj,jl) * ( pt_ low(ji+1,jj,jl) - pt_low(ji,jj,jl) ) <= 0. .AND. &1480 & pfv_ho(ji,jj,jl) * ( pt_ low(ji,jj+1,jl) - pt_low(ji,jj,jl) ) <= 0. ) THEN1029 IF ( pfu_ho(ji,jj,jl) * ( pt_ups(ji+1,jj,jl) - pt_ups(ji,jj,jl) ) <= 0. .AND. & 1030 & pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0. ) THEN 1481 1031 ! 1482 IF( pfu_ho(ji,jj,jl) * ( zti_low(ji+1,jj,jl) - zti_low(ji,jj,jl) ) <= 0 .AND. & 1483 & pfv_ho(ji,jj,jl) * ( ztj_low(ji,jj+1,jl) - ztj_low(ji,jj,jl) ) <= 0) pfu_ho(ji,jj,jl)=0. ; pfv_ho(ji,jj,jl)=0. 1032 IF( pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj,jl) - zti_ups(ji,jj,jl) ) <= 0 .AND. & 1033 & pfv_ho(ji,jj,jl) * ( ztj_ups(ji,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0) THEN 1034 pfu_ho(ji,jj,jl)=0. 1035 pfv_ho(ji,jj,jl)=0. 1036 ENDIF 1484 1037 ! 1485 IF( pfu_ho(ji,jj,jl) * ( pt_low(ji ,jj,jl) - pt_low(ji-1,jj,jl) ) <= 0 .AND. & 1486 & pfv_ho(ji,jj,jl) * ( pt_low(ji ,jj,jl) - pt_low(ji,jj-1,jl) ) <= 0) pfu_ho(ji,jj,jl)=0. ; pfv_ho(ji,jj,jl)=0. 1038 IF( pfu_ho(ji,jj,jl) * ( pt_ups(ji ,jj,jl) - pt_ups(ji-1,jj,jl) ) <= 0 .AND. & 1039 & pfv_ho(ji,jj,jl) * ( pt_ups(ji ,jj,jl) - pt_ups(ji,jj-1,jl) ) <= 0) THEN 1040 pfu_ho(ji,jj,jl)=0. 1041 pfv_ho(ji,jj,jl)=0. 1042 ENDIF 1487 1043 ! 1488 1044 ENDIF … … 1492 1048 CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond. 1493 1049 1494 ELSEIF( ll_prelimiter_devore ) THEN1495 DO jl = 1, jpl1496 DO jj = 2, jpjm11497 DO ji = fs_2, fs_jpim11498 zti_low(ji,jj,jl)= pt_low(ji+1,jj ,jl)1499 ztj_low(ji,jj,jl)= pt_low(ji ,jj+1,jl)1500 END DO1501 END DO1502 END DO1503 CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_low, 'T', 1., ztj_low, 'T', 1. )1504 1505 z1_dt = 1._wp / pdt1506 DO jl = 1, jpl1507 DO jj = 2, jpjm11508 DO ji = fs_2, fs_jpim11509 zsign = SIGN( 1., pt_low(ji+1,jj,jl) - pt_low(ji,jj,jl) )1510 pfu_ho(ji,jj,jl) = zsign * MAX( 0. , MIN( ABS(pfu_ho(ji,jj,jl)) , &1511 & zsign * ( pt_low (ji ,jj,jl) - pt_low (ji-1,jj,jl) ) * e1e2t(ji ,jj) * z1_dt , &1512 & zsign * ( zti_low(ji+1,jj,jl) - zti_low(ji ,jj,jl) ) * e1e2t(ji+1,jj) * z1_dt ) )1513 1514 zsign = SIGN( 1., pt_low(ji,jj+1,jl) - pt_low(ji,jj,jl) )1515 pfv_ho(ji,jj,jl) = zsign * MAX( 0. , MIN( ABS(pfv_ho(ji,jj,jl)) , &1516 & zsign * ( pt_low (ji,jj ,jl) - pt_low (ji,jj-1,jl) ) * e1e2t(ji,jj ) * z1_dt , &1517 & zsign * ( ztj_low(ji,jj+1,jl) - ztj_low(ji,jj ,jl) ) * e1e2t(ji,jj+1) * z1_dt ) )1518 END DO1519 END DO1520 END DO1521 CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond.1522 1523 1050 ENDIF 1524 1051 … … 1526 1053 ! Search local extrema 1527 1054 ! -------------------- 1528 ! max/min of pt & pt_ lowwith large negative/positive value (-/+zbig) outside ice cover1055 ! max/min of pt & pt_ups with large negative/positive value (-/+zbig) outside ice cover 1529 1056 z1_dt = 1._wp / pdt 1530 1057 DO jl = 1, jpl … … 1532 1059 DO jj = 1, jpj 1533 1060 DO ji = 1, jpi 1534 IF ( pt(ji,jj,jl) <= 0._wp .AND. pt_ low(ji,jj,jl) <= 0._wp ) THEN1061 IF ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN 1535 1062 zbup(ji,jj) = -zbig 1536 1063 zbdo(ji,jj) = zbig 1537 ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ low(ji,jj,jl) > 0._wp ) THEN1538 zbup(ji,jj) = pt_ low(ji,jj,jl)1539 zbdo(ji,jj) = pt_ low(ji,jj,jl)1540 ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ low(ji,jj,jl) <= 0._wp ) THEN1064 ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) > 0._wp ) THEN 1065 zbup(ji,jj) = pt_ups(ji,jj,jl) 1066 zbdo(ji,jj) = pt_ups(ji,jj,jl) 1067 ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN 1541 1068 zbup(ji,jj) = pt(ji,jj,jl) 1542 1069 zbdo(ji,jj) = pt(ji,jj,jl) 1543 1070 ELSE 1544 zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ low(ji,jj,jl) )1545 zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ low(ji,jj,jl) )1071 zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ups(ji,jj,jl) ) 1072 zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ups(ji,jj,jl) ) 1546 1073 ENDIF 1547 1074 END DO … … 1551 1078 DO ji = fs_2, fs_jpim1 ! vector opt. 1552 1079 ! 1553 IF( .NOT. ll_9points ) THEN 1554 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 1555 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), zbdo(ji ,jj-1), zbdo(ji ,jj+1) ) 1556 ! 1557 ELSE 1558 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 1559 & zbup(ji-1,jj-1), zbup(ji+1,jj+1), zbup(ji+1,jj-1), zbup(ji-1,jj+1) ) 1560 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), zbdo(ji ,jj-1), zbdo(ji ,jj+1), & 1561 & zbdo(ji-1,jj-1), zbdo(ji+1,jj+1), zbdo(ji+1,jj-1), zbdo(ji-1,jj+1) ) 1562 ENDIF 1080 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 1081 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj ), zbdo(ji+1,jj ), zbdo(ji ,jj-1), zbdo(ji ,jj+1) ) 1563 1082 ! 1564 1083 zpos = MAX( 0., pfu_ho(ji-1,jj,jl) ) - MIN( 0., pfu_ho(ji ,jj,jl) ) & ! positive/negative part of the flux … … 1567 1086 & + MAX( 0., pfv_ho(ji,jj ,jl) ) - MIN( 0., pfv_ho(ji,jj-1,jl) ) 1568 1087 ! 1569 IF( ll_HgradU .AND. .NOT.ll_gurvan ) THEN 1570 zneg2 = ( pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) & 1571 & ) * ( 1. - pamsk ) 1572 zpos2 = ( - pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) - pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) & 1573 & ) * ( 1. - pamsk ) 1574 ELSE 1575 zneg2 = 0. ; zpos2 = 0. 1088 zpos = zpos - ( pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) & 1089 & ) * ( 1. - pamsk ) 1090 zneg = zneg + ( pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) & 1091 & ) * ( 1. - pamsk ) 1092 ! 1093 ! ! up & down beta terms 1094 IF( zpos > 0. ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 1095 ELSE ; zbetup(ji,jj,jl) = 0. ! zbig 1576 1096 ENDIF 1577 1097 ! 1578 ! ! up & down beta terms 1579 IF( (zpos+zpos2) > 0. ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_low(ji,jj,jl) ) / (zpos+zpos2) * e1e2t(ji,jj) * z1_dt 1580 ELSE ; zbetup(ji,jj,jl) = 0. ! zbig 1581 ENDIF 1582 ! 1583 IF( (zneg+zneg2) > 0. ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_low(ji,jj,jl) - zdo ) / (zneg+zneg2) * e1e2t(ji,jj) * z1_dt 1584 ELSE ; zbetdo(ji,jj,jl) = 0. ! zbig 1098 IF( zneg > 0. ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 1099 ELSE ; zbetdo(ji,jj,jl) = 0. ! zbig 1585 1100 ENDIF 1586 1101 ! … … 1589 1104 IF( zdo == zbig ) zbetdo(ji,jj,jl) = 0. ! zbig 1590 1105 ! 1591 1592 !! IF( ji==26 .AND. jj==86) THEN1593 ! WRITE(numout,*) '-----------------'1594 ! WRITE(numout,*) 'zpos',zpos,zpos21595 ! WRITE(numout,*) 'zneg',zneg,zneg21596 ! WRITE(numout,*) 'puc/pu',ABS(puc(ji,jj))/MAX(epsi20, ABS(pu(ji,jj)))1597 ! WRITE(numout,*) 'pvc/pv',ABS(pvc(ji,jj))/MAX(epsi20, ABS(pv(ji,jj)))1598 ! WRITE(numout,*) 'pucm1/pu',ABS(puc(ji-1,jj))/MAX(epsi20, ABS(pu(ji-1,jj)))1599 ! WRITE(numout,*) 'pvcm1/pv',ABS(pvc(ji,jj-1))/MAX(epsi20, ABS(pv(ji,jj-1)))1600 ! WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj)+pfu_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt1601 ! WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj)+pfv_low(ji,jj)) * r1_e1e2t(ji,jj) * pdt1602 ! WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj)+pfu_low(ji-1,jj)) * r1_e1e2t(ji,jj) * pdt1603 ! WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1)+pfv_low(ji,jj-1)) * r1_e1e2t(ji,jj) * pdt1604 ! WRITE(numout,*) 'pfu_low',pfu_low(ji,jj) * r1_e1e2t(ji,jj) * pdt1605 ! WRITE(numout,*) 'pfv_low',pfv_low(ji,jj) * r1_e1e2t(ji,jj) * pdt1606 ! WRITE(numout,*) 'pfu_lowm1',pfu_low(ji-1,jj) * r1_e1e2t(ji,jj) * pdt1607 ! WRITE(numout,*) 'pfv_lowm1',pfv_low(ji,jj-1) * r1_e1e2t(ji,jj) * pdt1608 !1609 ! WRITE(numout,*) 'pt',pt(ji,jj)1610 ! WRITE(numout,*) 'ptim1',pt(ji-1,jj)1611 ! WRITE(numout,*) 'ptjm1',pt(ji,jj-1)1612 ! WRITE(numout,*) 'pt_low',pt_low(ji,jj)1613 ! WRITE(numout,*) 'zbetup',zbetup(ji,jj)1614 ! WRITE(numout,*) 'zbetdo',zbetdo(ji,jj)1615 ! WRITE(numout,*) 'zup',zup1616 ! WRITE(numout,*) 'zdo',zdo1617 ! ENDIF1618 !1619 1106 END DO 1620 1107 END DO … … 1633 1120 ! 1634 1121 zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu ) 1635 1636 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_low(ji,jj,jl) 1637 1638 !! IF( ji==26 .AND. jj==86) THEN 1639 !! WRITE(numout,*) 'coefU',zcoef 1640 !! WRITE(numout,*) 'pfu_ho',(pfu_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt 1641 !! WRITE(numout,*) 'pfu_hom1',(pfu_ho(ji-1,jj,jl)) * r1_e1e2t(ji,jj) * pdt 1642 !! ENDIF 1643 1122 ! 1123 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_ups(ji,jj,jl) 1124 ! 1644 1125 END DO 1645 1126 END DO … … 1652 1133 ! 1653 1134 zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv ) 1654 1655 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_low(ji,jj,jl) 1656 1657 !! IF( ji==26 .AND. jj==86) THEN 1658 !! WRITE(numout,*) 'coefV',zcoef 1659 !! WRITE(numout,*) 'pfv_ho',(pfv_ho(ji,jj,jl)) * r1_e1e2t(ji,jj) * pdt 1660 !! WRITE(numout,*) 'pfv_hom1',(pfv_ho(ji,jj-1,jl)) * r1_e1e2t(ji,jj) * pdt 1661 !! ENDIF 1135 ! 1136 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_ups(ji,jj,jl) 1137 ! 1662 1138 END DO 1663 1139 END DO 1664 1140 1665 1141 ! clem test 1666 DO jj = 2, jpjm1 1667 DO ji = 2, fs_jpim1 ! vector opt. 1668 IF( ll_gurvan ) THEN 1669 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 1670 & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 1671 ELSE 1672 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 1673 & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 1674 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1675 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1676 & ) * tmask(ji,jj,1) 1677 ENDIF 1678 IF( zzt(ji,jj,jl) < -epsi20 ) THEN 1679 WRITE(numout,*) 'T<0 nonosc',zzt(ji,jj,jl) 1680 ENDIF 1681 END DO 1682 END DO 1142 !! DO jj = 2, jpjm1 1143 !! DO ji = 2, fs_jpim1 ! vector opt. 1144 !! zzt = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 1145 !! & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 1146 !! & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1147 !! & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 1148 !! & ) * tmask(ji,jj,1) 1149 !! IF( zzt < -epsi20 ) THEN 1150 !! WRITE(numout,*) 'T<0 nonosc',zzt 1151 !! ENDIF 1152 !! END DO 1153 !! END DO 1683 1154 1684 1155 END DO … … 1688 1159 END SUBROUTINE nonosc_2d 1689 1160 1690 SUBROUTINE limiter_x( pdt, pu, p uc, pt, pfu_ho, pfu_ups )1161 SUBROUTINE limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 1691 1162 !!--------------------------------------------------------------------- 1692 1163 !! *** ROUTINE limiter_x *** … … 1696 1167 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1697 1168 REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pu ! ice i-velocity => u*e2 1698 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: puc ! ice i-velocity *A => u*e2*a1699 1169 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pt ! ice tracer 1700 1170 REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfu_ho ! high order flux … … 1747 1217 ELSE 1748 1218 Cr = 0. 1749 !IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm * 1.e201750 !ELSE ; Cr = Rjp * 1.e201751 !ENDIF1752 1219 ENDIF 1753 1220 … … 1787 1254 END SUBROUTINE limiter_x 1788 1255 1789 SUBROUTINE limiter_y( pdt, pv, p vc, pt, pfv_ho, pfv_ups )1256 SUBROUTINE limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 1790 1257 !!--------------------------------------------------------------------- 1791 1258 !! *** ROUTINE limiter_y *** … … 1795 1262 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1796 1263 REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pv ! ice i-velocity => u*e2 1797 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pvc ! ice i-velocity *A => u*e2*a1798 1264 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pt ! ice tracer 1799 1265 REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ho ! high order flux … … 1847 1313 ELSE 1848 1314 Cr = 0. 1849 !IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm * 1.e201850 !ELSE ; Cr = Rjp * 1.e201851 !ENDIF1852 1315 ENDIF 1853 1316
Note: See TracChangeset
for help on using the changeset viewer.
