- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_adv_umx.F90 (modified) (40 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/ICE/icedyn_adv_umx.F90
r10945 r13463 51 51 ! 52 52 !! * Substitutions 53 # include " vectopt_loop_substitute.h90"53 # include "do_loop_substitute.h90" 54 54 !!---------------------------------------------------------------------- 55 55 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 83 83 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content 84 84 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration 85 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction85 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond concentration 86 86 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume 87 87 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content … … 107 107 ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 108 108 DO jl = 1, jpl 109 DO jj = 2, jpjm1 110 DO ji = fs_2, fs_jpim1 111 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 112 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 113 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 114 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 115 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 116 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 117 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 118 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 119 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 120 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 121 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 122 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 123 END DO 124 END DO 125 END DO 126 CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. ) 109 DO_2D( 0, 0, 0, 0 ) 110 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 111 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 112 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 113 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 114 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 115 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 116 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 117 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 118 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 119 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 120 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 121 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 122 END_2D 123 END DO 124 CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1.0_wp, zhs_max, 'T', 1.0_wp, zhip_max, 'T', 1.0_wp ) 127 125 ! 128 126 ! … … 130 128 ! Note: the advection split is applied at the next time-step in order to avoid blocking global comm. 131 129 ! this should not affect too much the stability 132 zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * r dt_ice * r1_e1u(:,:) )133 zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * r dt_ice * r1_e2v(:,:) ) )130 zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rDt_ice * r1_e1u(:,:) ) 131 zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rDt_ice * r1_e2v(:,:) ) ) 134 132 135 133 ! non-blocking global communication send zcflnow and receive zcflprv … … 139 137 ELSE ; icycle = 1 140 138 ENDIF 141 zdt = r dt_ice / REAL(icycle)139 zdt = rDt_ice / REAL(icycle) 142 140 143 141 ! --- transport --- ! … … 152 150 ! 153 151 ! --- define velocity for advection: u*grad(H) --- ! 154 DO jj = 2, jpjm1 155 DO ji = fs_2, fs_jpim1 156 IF ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN ; zcu_box(ji,jj) = 0._wp 157 ELSEIF( pu_ice(ji,jj) > 0._wp ) THEN ; zcu_box(ji,jj) = pu_ice(ji-1,jj) 158 ELSE ; zcu_box(ji,jj) = pu_ice(ji ,jj) 159 ENDIF 160 161 IF ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN ; zcv_box(ji,jj) = 0._wp 162 ELSEIF( pv_ice(ji,jj) > 0._wp ) THEN ; zcv_box(ji,jj) = pv_ice(ji,jj-1) 163 ELSE ; zcv_box(ji,jj) = pv_ice(ji,jj ) 164 ENDIF 165 END DO 166 END DO 152 DO_2D( 0, 0, 0, 0 ) 153 IF ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN ; zcu_box(ji,jj) = 0._wp 154 ELSEIF( pu_ice(ji,jj) > 0._wp ) THEN ; zcu_box(ji,jj) = pu_ice(ji-1,jj) 155 ELSE ; zcu_box(ji,jj) = pu_ice(ji ,jj) 156 ENDIF 157 158 IF ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN ; zcv_box(ji,jj) = 0._wp 159 ELSEIF( pv_ice(ji,jj) > 0._wp ) THEN ; zcv_box(ji,jj) = pv_ice(ji,jj-1) 160 ELSE ; zcv_box(ji,jj) = pv_ice(ji,jj ) 161 ENDIF 162 END_2D 167 163 168 164 !---------------! … … 187 183 IF( .NOT. ALLOCATED(jmsk_small) ) ALLOCATE( jmsk_small(jpi,jpj,jpl) ) 188 184 DO jl = 1, jpl 189 DO jj = 1, jpjm1 190 DO ji = 1, jpim1 191 zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) ) 192 IF( zvi_cen < epsi06) THEN ; imsk_small(ji,jj,jl) = 0 193 ELSE ; imsk_small(ji,jj,jl) = 1 ; ENDIF 194 zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) ) 195 IF( zvi_cen < epsi06) THEN ; jmsk_small(ji,jj,jl) = 0 196 ELSE ; jmsk_small(ji,jj,jl) = 1 ; ENDIF 197 END DO 198 END DO 185 DO_2D( 1, 0, 1, 0 ) 186 zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) ) 187 IF( zvi_cen < epsi06) THEN ; imsk_small(ji,jj,jl) = 0 188 ELSE ; imsk_small(ji,jj,jl) = 1 ; ENDIF 189 zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) ) 190 IF( zvi_cen < epsi06) THEN ; jmsk_small(ji,jj,jl) = 0 191 ELSE ; jmsk_small(ji,jj,jl) = 1 ; ENDIF 192 END_2D 199 193 END DO 200 194 ENDIF … … 319 313 ! 320 314 !== Ice age ==! 321 IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN 322 zamsk = 1._wp 323 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 324 & poa_i, poa_i ) 325 ENDIF 315 zamsk = 1._wp 316 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 317 & poa_i, poa_i ) 326 318 ! 327 319 !== melt ponds ==! 328 320 IF ( ln_pnd_H12 ) THEN 329 ! fraction321 ! concentration 330 322 zamsk = 1._wp 331 323 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat , zv_cat , zcu_box, zcv_box, & … … 340 332 !== Open water area ==! 341 333 zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 342 DO jj = 2, jpjm1 343 DO ji = fs_2, fs_jpim1 344 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & 345 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 346 END DO 347 END DO 348 CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1. ) 334 DO_2D( 0, 0, 0, 0 ) 335 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & 336 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 337 END_2D 338 CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1.0_wp ) 349 339 ! 350 340 ! … … 354 344 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 355 345 ! 356 ! Make sure ice thickness is not too big 357 ! (because ice thickness can be too large where ice concentration is very small) 358 CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 359 346 ! --- Make sure ice thickness is not too big --- ! 347 ! (because ice thickness can be too large where ice concentration is very small) 348 CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 349 ! 350 ! --- Ensure snow load is not too big --- ! 351 CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 352 ! 360 353 END DO 361 354 ! … … 448 441 IF( pamsk == 0._wp ) THEN 449 442 DO jl = 1, jpl 450 DO jj = 1, jpjm1 451 DO ji = 1, fs_jpim1 452 IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 453 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj) 454 zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj) 455 ELSE 456 zfu_ho (ji,jj,jl) = 0._wp 457 zfu_ups(ji,jj,jl) = 0._wp 458 ENDIF 459 ! 460 IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 461 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj) 462 zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj) 463 ELSE 464 zfv_ho (ji,jj,jl) = 0._wp 465 zfv_ups(ji,jj,jl) = 0._wp 466 ENDIF 467 END DO 468 END DO 443 DO_2D( 1, 0, 1, 0 ) 444 IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 445 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj) 446 zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj) 447 ELSE 448 zfu_ho (ji,jj,jl) = 0._wp 449 zfu_ups(ji,jj,jl) = 0._wp 450 ENDIF 451 ! 452 IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 453 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj) 454 zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj) 455 ELSE 456 zfv_ho (ji,jj,jl) = 0._wp 457 zfv_ups(ji,jj,jl) = 0._wp 458 ENDIF 459 END_2D 469 460 END DO 470 461 … … 472 463 ! thus we calculate the upstream solution and apply a limiter again 473 464 DO jl = 1, jpl 474 DO jj = 2, jpjm1 475 DO ji = fs_2, fs_jpim1 476 ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) 477 ! 478 zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 479 END DO 480 END DO 481 END DO 482 CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. ) 465 DO_2D( 0, 0, 0, 0 ) 466 ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) 467 ! 468 zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 469 END_2D 470 END DO 471 CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1.0_wp ) 483 472 ! 484 473 IF ( np_limiter == 1 ) THEN … … 495 484 IF( PRESENT( pua_ho ) ) THEN 496 485 DO jl = 1, jpl 497 DO jj = 1, jpjm1 498 DO ji = 1, fs_jpim1 499 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 500 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 501 END DO 502 END DO 486 DO_2D( 1, 0, 1, 0 ) 487 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 488 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 489 END_2D 503 490 END DO 504 491 ENDIF … … 507 494 ! --------------------------------- 508 495 DO jl = 1, jpl 509 DO jj = 2, jpjm1 510 DO ji = fs_2, fs_jpim1 511 ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) 512 ! 513 ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 514 END DO 515 END DO 516 END DO 517 CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1. ) 496 DO_2D( 0, 0, 0, 0 ) 497 ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) 498 ! 499 ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 500 END_2D 501 END DO 502 CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1.0_wp ) 518 503 ! 519 504 END SUBROUTINE adv_umx … … 543 528 ! 544 529 DO jl = 1, jpl 545 DO jj = 1, jpjm1 546 DO ji = 1, fs_jpim1 530 DO_2D( 1, 0, 1, 0 ) 531 pfu_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) 532 pfv_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) 533 END_2D 534 END DO 535 ! 536 ELSE !** alternate directions **! 537 ! 538 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 539 ! 540 DO jl = 1, jpl !-- flux in x-direction 541 DO_2D( 1, 0, 1, 0 ) 547 542 pfu_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) 543 END_2D 544 END DO 545 ! 546 DO jl = 1, jpl !-- first guess of tracer from u-flux 547 DO_2D( 0, 0, 0, 0 ) 548 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) ) & 549 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) 550 ! 551 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 552 END_2D 553 END DO 554 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 555 ! 556 DO jl = 1, jpl !-- flux in y-direction 557 DO_2D( 1, 0, 1, 0 ) 558 pfv_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) 559 END_2D 560 END DO 561 ! 562 ELSE !== even ice time step: adv_y then adv_x ==! 563 ! 564 DO jl = 1, jpl !-- flux in y-direction 565 DO_2D( 1, 0, 1, 0 ) 548 566 pfv_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) 549 END DO 550 END DO 551 END DO 552 ! 553 ELSE !** alternate directions **! 554 ! 555 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 567 END_2D 568 END DO 569 ! 570 DO jl = 1, jpl !-- first guess of tracer from v-flux 571 DO_2D( 0, 0, 0, 0 ) 572 ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) ) & 573 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) 574 ! 575 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 576 END_2D 577 END DO 578 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 556 579 ! 557 580 DO jl = 1, jpl !-- flux in x-direction 558 DO jj = 1, jpjm1 559 DO ji = 1, fs_jpim1 560 pfu_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) 561 END DO 562 END DO 563 END DO 564 ! 565 DO jl = 1, jpl !-- first guess of tracer from u-flux 566 DO jj = 2, jpjm1 567 DO ji = fs_2, fs_jpim1 568 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) ) & 569 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) 570 ! 571 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 572 END DO 573 END DO 574 END DO 575 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 576 ! 577 DO jl = 1, jpl !-- flux in y-direction 578 DO jj = 1, jpjm1 579 DO ji = 1, fs_jpim1 580 pfv_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) 581 END DO 582 END DO 583 END DO 584 ! 585 ELSE !== even ice time step: adv_y then adv_x ==! 586 ! 587 DO jl = 1, jpl !-- flux in y-direction 588 DO jj = 1, jpjm1 589 DO ji = 1, fs_jpim1 590 pfv_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) 591 END DO 592 END DO 593 END DO 594 ! 595 DO jl = 1, jpl !-- first guess of tracer from v-flux 596 DO jj = 2, jpjm1 597 DO ji = fs_2, fs_jpim1 598 ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) ) & 599 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) 600 ! 601 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 602 END DO 603 END DO 604 END DO 605 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 606 ! 607 DO jl = 1, jpl !-- flux in x-direction 608 DO jj = 1, jpjm1 609 DO ji = 1, fs_jpim1 610 pfu_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) 611 END DO 612 END DO 581 DO_2D( 1, 0, 1, 0 ) 582 pfu_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) 583 END_2D 613 584 END DO 614 585 ! … … 618 589 ! 619 590 DO jl = 1, jpl !-- after tracer with upstream scheme 620 DO jj = 2, jpjm1 621 DO ji = fs_2, fs_jpim1 622 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj ,jl) & 623 & + pfv_ups(ji,jj,jl) - pfv_ups(ji ,jj-1,jl) ) & 624 & + ( pu (ji,jj ) - pu (ji-1,jj ) & 625 & + pv (ji,jj ) - pv (ji ,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) 626 ! 627 pt_ups(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 628 END DO 629 END DO 630 END DO 631 CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1. ) 591 DO_2D( 0, 0, 0, 0 ) 592 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj ,jl) & 593 & + pfv_ups(ji,jj,jl) - pfv_ups(ji ,jj-1,jl) ) & 594 & + ( pu (ji,jj ) - pu (ji-1,jj ) & 595 & + pv (ji,jj ) - pv (ji ,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) 596 ! 597 pt_ups(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 598 END_2D 599 END DO 600 CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1.0_wp ) 632 601 633 602 END SUBROUTINE upstream … … 659 628 ! 660 629 DO jl = 1, jpl 661 DO jj = 1, jpjm1 662 DO ji = 1, fs_jpim1 663 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj ,jl) ) 664 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji ,jj+1,jl) ) 665 END DO 666 END DO 630 DO_2D( 1, 0, 1, 0 ) 631 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj ,jl) ) 632 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji ,jj+1,jl) ) 633 END_2D 667 634 END DO 668 635 ! … … 679 646 ! 680 647 DO jl = 1, jpl !-- flux in x-direction 681 DO jj = 1, jpjm1 682 DO ji = 1, fs_jpim1 683 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 684 END DO 685 END DO 648 DO_2D( 1, 0, 1, 0 ) 649 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 650 END_2D 686 651 END DO 687 652 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 688 653 689 654 DO jl = 1, jpl !-- first guess of tracer from u-flux 690 DO jj = 2, jpjm1 691 DO ji = fs_2, fs_jpim1 692 ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) & 693 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) 694 ! 695 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 696 END DO 697 END DO 698 END DO 699 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 655 DO_2D( 0, 0, 0, 0 ) 656 ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) & 657 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) 658 ! 659 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 660 END_2D 661 END DO 662 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 700 663 701 664 DO jl = 1, jpl !-- flux in y-direction 702 DO jj = 1, jpjm1 703 DO ji = 1, fs_jpim1 704 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 705 END DO 706 END DO 665 DO_2D( 1, 0, 1, 0 ) 666 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 667 END_2D 707 668 END DO 708 669 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) … … 711 672 ! 712 673 DO jl = 1, jpl !-- flux in y-direction 713 DO jj = 1, jpjm1 714 DO ji = 1, fs_jpim1 715 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 716 END DO 717 END DO 674 DO_2D( 1, 0, 1, 0 ) 675 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 676 END_2D 718 677 END DO 719 678 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 720 679 ! 721 680 DO jl = 1, jpl !-- first guess of tracer from v-flux 722 DO jj = 2, jpjm1 723 DO ji = fs_2, fs_jpim1 724 ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) & 725 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) 726 ! 727 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 728 END DO 729 END DO 730 END DO 731 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 681 DO_2D( 0, 0, 0, 0 ) 682 ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) & 683 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) 684 ! 685 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 686 END_2D 687 END DO 688 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 732 689 ! 733 690 DO jl = 1, jpl !-- flux in x-direction 734 DO jj = 1, jpjm1 735 DO ji = 1, fs_jpim1 736 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 737 END DO 738 END DO 691 DO_2D( 1, 0, 1, 0 ) 692 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 693 END_2D 739 694 END DO 740 695 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) … … 782 737 ! !-- advective form update in zpt --! 783 738 DO jl = 1, jpl 784 DO jj = 2, jpjm1 785 DO ji = fs_2, fs_jpim1 786 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pubox(ji,jj ) * ( zt_u(ji,jj,jl) - zt_u(ji-1,jj,jl) ) * r1_e1t (ji,jj) & 787 & + pt (ji,jj,jl) * ( pu (ji,jj ) - pu (ji-1,jj ) ) * r1_e1e2t(ji,jj) & 788 & * pamsk & 789 & ) * pdt ) * tmask(ji,jj,1) 790 END DO 791 END DO 792 END DO 793 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 739 DO_2D( 0, 0, 0, 0 ) 740 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pubox(ji,jj ) * ( zt_u(ji,jj,jl) - zt_u(ji-1,jj,jl) ) * r1_e1t (ji,jj) & 741 & + pt (ji,jj,jl) * ( pu (ji,jj ) - pu (ji-1,jj ) ) * r1_e1e2t(ji,jj) & 742 & * pamsk & 743 & ) * pdt ) * tmask(ji,jj,1) 744 END_2D 745 END DO 746 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 794 747 ! 795 748 ! !-- ultimate interpolation of pt at v-point --! … … 811 764 ! !-- advective form update in zpt --! 812 765 DO jl = 1, jpl 813 DO jj = 2, jpjm1 814 DO ji = fs_2, fs_jpim1 815 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pvbox(ji,jj ) * ( zt_v(ji,jj,jl) - zt_v(ji,jj-1,jl) ) * r1_e2t (ji,jj) & 816 & + pt (ji,jj,jl) * ( pv (ji,jj ) - pv (ji,jj-1 ) ) * r1_e1e2t(ji,jj) & 817 & * pamsk & 818 & ) * pdt ) * tmask(ji,jj,1) 819 END DO 820 END DO 821 END DO 822 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 766 DO_2D( 0, 0, 0, 0 ) 767 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pvbox(ji,jj ) * ( zt_v(ji,jj,jl) - zt_v(ji,jj-1,jl) ) * r1_e2t (ji,jj) & 768 & + pt (ji,jj,jl) * ( pv (ji,jj ) - pv (ji,jj-1 ) ) * r1_e1e2t(ji,jj) & 769 & * pamsk & 770 & ) * pdt ) * tmask(ji,jj,1) 771 END_2D 772 END DO 773 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 823 774 ! 824 775 ! !-- ultimate interpolation of pt at u-point --! … … 864 815 DO jl = 1, jpl 865 816 DO jj = 2, jpjm1 ! First derivative (gradient) 866 DO ji = 1, fs_jpim1817 DO ji = 1, jpim1 867 818 ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 868 819 END DO 869 820 ! ! Second derivative (Laplacian) 870 DO ji = fs_2, fs_jpim1821 DO ji = 2, jpim1 871 822 ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj) 872 823 END DO 873 824 END DO 874 825 END DO 875 CALL lbc_lnk( 'icedyn_adv_umx', ztu2, 'T', 1. )826 CALL lbc_lnk( 'icedyn_adv_umx', ztu2, 'T', 1.0_wp ) 876 827 ! 877 828 ! !-- BiLaplacian in i-direction --! 878 829 DO jl = 1, jpl 879 830 DO jj = 2, jpjm1 ! Third derivative 880 DO ji = 1, fs_jpim1831 DO ji = 1, jpim1 881 832 ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 882 833 END DO 883 834 ! ! Fourth derivative 884 DO ji = fs_2, fs_jpim1835 DO ji = 2, jpim1 885 836 ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj) 886 837 END DO 887 838 END DO 888 839 END DO 889 CALL lbc_lnk( 'icedyn_adv_umx', ztu4, 'T', 1. )840 CALL lbc_lnk( 'icedyn_adv_umx', ztu4, 'T', 1.0_wp ) 890 841 ! 891 842 ! … … 895 846 ! 896 847 DO jl = 1, jpl 897 DO jj = 1, jpjm1 898 DO ji = 1, fs_jpim1 ! vector opt. 899 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 900 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 901 END DO 902 END DO 848 DO_2D( 1, 0, 1, 0 ) 849 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 850 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 851 END_2D 903 852 END DO 904 853 ! … … 906 855 ! 907 856 DO jl = 1, jpl 908 DO jj = 1, jpjm1 909 DO ji = 1, fs_jpim1 ! vector opt. 910 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 911 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 912 & - zcu * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 913 END DO 914 END DO 857 DO_2D( 1, 0, 1, 0 ) 858 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 859 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 860 & - zcu * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 861 END_2D 915 862 END DO 916 863 ! … … 918 865 ! 919 866 DO jl = 1, jpl 920 DO jj = 1, jpjm1 921 DO ji = 1, fs_jpim1 ! vector opt. 922 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 923 zdx2 = e1u(ji,jj) * e1u(ji,jj) 867 DO_2D( 1, 0, 1, 0 ) 868 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 869 zdx2 = e1u(ji,jj) * e1u(ji,jj) 924 870 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 925 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 926 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 927 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 928 & - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 929 END DO 930 END DO 871 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 872 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 873 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 874 & - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 875 END_2D 931 876 END DO 932 877 ! … … 934 879 ! 935 880 DO jl = 1, jpl 936 DO jj = 1, jpjm1 937 DO ji = 1, fs_jpim1 ! vector opt. 938 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 939 zdx2 = e1u(ji,jj) * e1u(ji,jj) 881 DO_2D( 1, 0, 1, 0 ) 882 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 883 zdx2 = e1u(ji,jj) * e1u(ji,jj) 940 884 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 941 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 942 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 943 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 944 & - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 945 END DO 946 END DO 885 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 886 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 887 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 888 & - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 889 END_2D 947 890 END DO 948 891 ! … … 950 893 ! 951 894 DO jl = 1, jpl 952 DO jj = 1, jpjm1 953 DO ji = 1, fs_jpim1 ! vector opt. 954 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 955 zdx2 = e1u(ji,jj) * e1u(ji,jj) 895 DO_2D( 1, 0, 1, 0 ) 896 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 897 zdx2 = e1u(ji,jj) * e1u(ji,jj) 956 898 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 957 zdx4 = zdx2 * zdx2 958 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 959 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 960 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 961 & - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) & 962 & + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) & 963 & - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) ) 964 END DO 965 END DO 899 zdx4 = zdx2 * zdx2 900 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 901 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 902 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 903 & - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) & 904 & + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) & 905 & - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) ) 906 END_2D 966 907 END DO 967 908 ! … … 973 914 IF( ll_neg ) THEN 974 915 DO jl = 1, jpl 975 DO jj = 1, jpjm1 976 DO ji = 1, fs_jpim1 977 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 978 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 979 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 980 ENDIF 981 END DO 982 END DO 916 DO_2D( 1, 0, 1, 0 ) 917 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 918 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 919 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 920 ENDIF 921 END_2D 983 922 END DO 984 923 ENDIF 985 924 ! !-- High order flux in i-direction --! 986 925 DO jl = 1, jpl 987 DO jj = 1, jpjm1 988 DO ji = 1, fs_jpim1 ! vector opt. 989 pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 990 END DO 991 END DO 926 DO_2D( 1, 0, 1, 0 ) 927 pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 928 END_2D 992 929 END DO 993 930 ! … … 1020 957 ! !-- Laplacian in j-direction --! 1021 958 DO jl = 1, jpl 1022 DO jj = 1, jpjm1 ! First derivative (gradient) 1023 DO ji = fs_2, fs_jpim1 1024 ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 1025 END DO 1026 END DO 1027 DO jj = 2, jpjm1 ! Second derivative (Laplacian) 1028 DO ji = fs_2, fs_jpim1 1029 ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj) 1030 END DO 1031 END DO 1032 END DO 1033 CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1. ) 959 DO_2D( 1, 0, 0, 0 ) 960 ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 961 END_2D 962 DO_2D( 0, 0, 0, 0 ) 963 ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj) 964 END_2D 965 END DO 966 CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1.0_wp ) 1034 967 ! 1035 968 ! !-- BiLaplacian in j-direction --! 1036 969 DO jl = 1, jpl 1037 DO jj = 1, jpjm1 ! First derivative 1038 DO ji = fs_2, fs_jpim1 1039 ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 1040 END DO 1041 END DO 1042 DO jj = 2, jpjm1 ! Second derivative 1043 DO ji = fs_2, fs_jpim1 1044 ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj) 1045 END DO 1046 END DO 1047 END DO 1048 CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1. ) 970 DO_2D( 1, 0, 0, 0 ) 971 ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 972 END_2D 973 DO_2D( 0, 0, 0, 0 ) 974 ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj) 975 END_2D 976 END DO 977 CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1.0_wp ) 1049 978 ! 1050 979 ! … … 1053 982 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 1054 983 DO jl = 1, jpl 1055 DO jj = 1, jpjm1 1056 DO ji = 1, fs_jpim1 1057 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & 1058 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 1059 END DO 1060 END DO 984 DO_2D( 1, 0, 1, 0 ) 985 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & 986 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 987 END_2D 1061 988 END DO 1062 989 ! 1063 990 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 1064 991 DO jl = 1, jpl 1065 DO jj = 1, jpjm1 1066 DO ji = 1, fs_jpim1 1067 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1068 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & 1069 & - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 1070 END DO 1071 END DO 992 DO_2D( 1, 0, 1, 0 ) 993 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 994 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & 995 & - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 996 END_2D 1072 997 END DO 1073 998 ! 1074 999 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 1075 1000 DO jl = 1, jpl 1076 DO jj = 1, jpjm1 1077 DO ji = 1, fs_jpim1 1078 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1079 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1001 DO_2D( 1, 0, 1, 0 ) 1002 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1003 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1080 1004 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 1081 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1082 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 1083 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 1084 & - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 1085 END DO 1086 END DO 1005 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1006 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 1007 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 1008 & - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 1009 END_2D 1087 1010 END DO 1088 1011 ! 1089 1012 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 1090 1013 DO jl = 1, jpl 1091 DO jj = 1, jpjm1 1092 DO ji = 1, fs_jpim1 1093 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1094 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1014 DO_2D( 1, 0, 1, 0 ) 1015 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1016 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1095 1017 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 1096 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1097 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 1098 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 1099 & - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 1100 END DO 1101 END DO 1018 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1019 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 1020 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 1021 & - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 1022 END_2D 1102 1023 END DO 1103 1024 ! 1104 1025 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 1105 1026 DO jl = 1, jpl 1106 DO jj = 1, jpjm1 1107 DO ji = 1, fs_jpim1 1108 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1109 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1027 DO_2D( 1, 0, 1, 0 ) 1028 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1029 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1110 1030 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 1111 zdy4 = zdy2 * zdy2 1112 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1113 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 1114 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 1115 & - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) & 1116 & + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) & 1117 & - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) ) 1118 END DO 1119 END DO 1031 zdy4 = zdy2 * zdy2 1032 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1033 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 1034 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 1035 & - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) & 1036 & + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) & 1037 & - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) ) 1038 END_2D 1120 1039 END DO 1121 1040 ! … … 1127 1046 IF( ll_neg ) THEN 1128 1047 DO jl = 1, jpl 1129 DO jj = 1, jpjm1 1130 DO ji = 1, fs_jpim1 1131 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 1132 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & 1133 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 1134 ENDIF 1135 END DO 1136 END DO 1048 DO_2D( 1, 0, 1, 0 ) 1049 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 1050 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & 1051 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 1052 ENDIF 1053 END_2D 1137 1054 END DO 1138 1055 ENDIF 1139 1056 ! !-- High order flux in j-direction --! 1140 1057 DO jl = 1, jpl 1141 DO jj = 1, jpjm1 1142 DO ji = 1, fs_jpim1 ! vector opt. 1143 pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 1144 END DO 1145 END DO 1058 DO_2D( 1, 0, 1, 0 ) 1059 pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 1060 END_2D 1146 1061 END DO 1147 1062 ! … … 1177 1092 ! -------------------------------------------------- 1178 1093 DO jl = 1, jpl 1179 DO jj = 1, jpjm1 1180 DO ji = 1, fs_jpim1 ! vector opt. 1181 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 1182 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1183 END DO 1184 END DO 1094 DO_2D( 1, 0, 1, 0 ) 1095 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 1096 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1097 END_2D 1185 1098 END DO 1186 1099 … … 1196 1109 1197 1110 DO jl = 1, jpl 1198 DO jj = 2, jpjm1 1199 DO ji = fs_2, fs_jpim1 1200 zti_ups(ji,jj,jl)= pt_ups(ji+1,jj ,jl) 1201 ztj_ups(ji,jj,jl)= pt_ups(ji ,jj+1,jl) 1202 END DO 1203 END DO 1204 END DO 1205 CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1., ztj_ups, 'T', 1. ) 1206 1207 DO jl = 1, jpl 1208 DO jj = 2, jpjm1 1209 DO ji = fs_2, fs_jpim1 1210 IF ( pfu_ho(ji,jj,jl) * ( pt_ups(ji+1,jj ,jl) - pt_ups(ji,jj,jl) ) <= 0._wp .AND. & 1211 & pfv_ho(ji,jj,jl) * ( pt_ups(ji ,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0._wp ) THEN 1212 ! 1213 IF( pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj ,jl) - zti_ups(ji,jj,jl) ) <= 0._wp .AND. & 1214 & pfv_ho(ji,jj,jl) * ( ztj_ups(ji ,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0._wp ) THEN 1215 pfu_ho(ji,jj,jl)=0._wp 1216 pfv_ho(ji,jj,jl)=0._wp 1217 ENDIF 1218 ! 1219 IF( pfu_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji-1,jj ,jl) ) <= 0._wp .AND. & 1220 & pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji ,jj-1,jl) ) <= 0._wp ) THEN 1221 pfu_ho(ji,jj,jl)=0._wp 1222 pfv_ho(ji,jj,jl)=0._wp 1223 ENDIF 1224 ! 1111 DO_2D( 0, 0, 0, 0 ) 1112 zti_ups(ji,jj,jl)= pt_ups(ji+1,jj ,jl) 1113 ztj_ups(ji,jj,jl)= pt_ups(ji ,jj+1,jl) 1114 END_2D 1115 END DO 1116 CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1.0_wp, ztj_ups, 'T', 1.0_wp ) 1117 1118 DO jl = 1, jpl 1119 DO_2D( 0, 0, 0, 0 ) 1120 IF ( pfu_ho(ji,jj,jl) * ( pt_ups(ji+1,jj ,jl) - pt_ups(ji,jj,jl) ) <= 0._wp .AND. & 1121 & pfv_ho(ji,jj,jl) * ( pt_ups(ji ,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0._wp ) THEN 1122 ! 1123 IF( pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj ,jl) - zti_ups(ji,jj,jl) ) <= 0._wp .AND. & 1124 & pfv_ho(ji,jj,jl) * ( ztj_ups(ji ,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0._wp ) THEN 1125 pfu_ho(ji,jj,jl)=0._wp 1126 pfv_ho(ji,jj,jl)=0._wp 1225 1127 ENDIF 1226 END DO 1227 END DO 1228 END DO 1229 CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond. 1128 ! 1129 IF( pfu_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji-1,jj ,jl) ) <= 0._wp .AND. & 1130 & pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji ,jj-1,jl) ) <= 0._wp ) THEN 1131 pfu_ho(ji,jj,jl)=0._wp 1132 pfv_ho(ji,jj,jl)=0._wp 1133 ENDIF 1134 ! 1135 ENDIF 1136 END_2D 1137 END DO 1138 CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1.0_wp, pfv_ho, 'V', -1.0_wp ) ! lateral boundary cond. 1230 1139 1231 1140 ENDIF … … 1237 1146 DO jl = 1, jpl 1238 1147 1239 DO jj = 1, jpj 1240 DO ji = 1, jpi 1241 IF ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN 1242 zbup(ji,jj) = -zbig 1243 zbdo(ji,jj) = zbig 1244 ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) > 0._wp ) THEN 1245 zbup(ji,jj) = pt_ups(ji,jj,jl) 1246 zbdo(ji,jj) = pt_ups(ji,jj,jl) 1247 ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN 1248 zbup(ji,jj) = pt(ji,jj,jl) 1249 zbdo(ji,jj) = pt(ji,jj,jl) 1250 ELSE 1251 zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ups(ji,jj,jl) ) 1252 zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ups(ji,jj,jl) ) 1253 ENDIF 1254 END DO 1255 END DO 1256 1257 DO jj = 2, jpjm1 1258 DO ji = fs_2, fs_jpim1 ! vector opt. 1259 ! 1260 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 1261 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) ) 1262 ! 1263 zpos = MAX( 0._wp, pfu_ho(ji-1,jj ,jl) ) - MIN( 0._wp, pfu_ho(ji ,jj ,jl) ) & ! positive/negative part of the flux 1264 & + MAX( 0._wp, pfv_ho(ji ,jj-1,jl) ) - MIN( 0._wp, pfv_ho(ji ,jj ,jl) ) 1265 zneg = MAX( 0._wp, pfu_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfu_ho(ji-1,jj ,jl) ) & 1266 & + MAX( 0._wp, pfv_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfv_ho(ji ,jj-1,jl) ) 1267 ! 1268 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) ) & 1269 & ) * ( 1. - pamsk ) 1270 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) ) & 1271 & ) * ( 1. - pamsk ) 1272 ! 1273 ! ! up & down beta terms 1274 ! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!) 1275 IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 1276 ELSE ; zbetup(ji,jj,jl) = 0._wp ! zbig 1277 ENDIF 1278 ! 1279 IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 1280 ELSE ; zbetdo(ji,jj,jl) = 0._wp ! zbig 1281 ENDIF 1282 ! 1283 ! if all the points are outside ice cover 1284 IF( zup == -zbig ) zbetup(ji,jj,jl) = 0._wp ! zbig 1285 IF( zdo == zbig ) zbetdo(ji,jj,jl) = 0._wp ! zbig 1286 ! 1287 END DO 1288 END DO 1289 END DO 1290 CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign) 1148 DO_2D( 1, 1, 1, 1 ) 1149 IF ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN 1150 zbup(ji,jj) = -zbig 1151 zbdo(ji,jj) = zbig 1152 ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) > 0._wp ) THEN 1153 zbup(ji,jj) = pt_ups(ji,jj,jl) 1154 zbdo(ji,jj) = pt_ups(ji,jj,jl) 1155 ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN 1156 zbup(ji,jj) = pt(ji,jj,jl) 1157 zbdo(ji,jj) = pt(ji,jj,jl) 1158 ELSE 1159 zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ups(ji,jj,jl) ) 1160 zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ups(ji,jj,jl) ) 1161 ENDIF 1162 END_2D 1163 1164 DO_2D( 0, 0, 0, 0 ) 1165 ! 1166 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 1167 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) ) 1168 ! 1169 zpos = MAX( 0._wp, pfu_ho(ji-1,jj ,jl) ) - MIN( 0._wp, pfu_ho(ji ,jj ,jl) ) & ! positive/negative part of the flux 1170 & + MAX( 0._wp, pfv_ho(ji ,jj-1,jl) ) - MIN( 0._wp, pfv_ho(ji ,jj ,jl) ) 1171 zneg = MAX( 0._wp, pfu_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfu_ho(ji-1,jj ,jl) ) & 1172 & + MAX( 0._wp, pfv_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfv_ho(ji ,jj-1,jl) ) 1173 ! 1174 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) ) & 1175 & ) * ( 1. - pamsk ) 1176 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) ) & 1177 & ) * ( 1. - pamsk ) 1178 ! 1179 ! ! up & down beta terms 1180 ! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!) 1181 IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 1182 ELSE ; zbetup(ji,jj,jl) = 0._wp ! zbig 1183 ENDIF 1184 ! 1185 IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 1186 ELSE ; zbetdo(ji,jj,jl) = 0._wp ! zbig 1187 ENDIF 1188 ! 1189 ! if all the points are outside ice cover 1190 IF( zup == -zbig ) zbetup(ji,jj,jl) = 0._wp ! zbig 1191 IF( zdo == zbig ) zbetdo(ji,jj,jl) = 0._wp ! zbig 1192 ! 1193 END_2D 1194 END DO 1195 CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1.0_wp, zbetdo, 'T', 1.0_wp ) ! lateral boundary cond. (unchanged sign) 1291 1196 1292 1197 … … 1294 1199 ! --------------------------------- 1295 1200 DO jl = 1, jpl 1296 DO jj = 1, jpjm1 1297 DO ji = 1, fs_jpim1 ! vector opt. 1298 zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 1299 zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) 1300 zcu = 0.5_wp + SIGN( 0.5_wp , pfu_ho(ji,jj,jl) ) 1301 ! 1302 zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu ) 1303 ! 1304 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_ups(ji,jj,jl) 1305 ! 1306 END DO 1307 END DO 1308 1309 DO jj = 1, jpjm1 1310 DO ji = 1, fs_jpim1 ! vector opt. 1311 zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 1312 zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) 1313 zcv = 0.5_wp + SIGN( 0.5_wp , pfv_ho(ji,jj,jl) ) 1314 ! 1315 zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv ) 1316 ! 1317 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_ups(ji,jj,jl) 1318 ! 1319 END DO 1320 END DO 1201 DO_2D( 1, 0, 1, 0 ) 1202 zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 1203 zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) 1204 zcu = 0.5_wp + SIGN( 0.5_wp , pfu_ho(ji,jj,jl) ) 1205 ! 1206 zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu ) 1207 ! 1208 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_ups(ji,jj,jl) 1209 ! 1210 END_2D 1211 1212 DO_2D( 1, 0, 1, 0 ) 1213 zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 1214 zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) 1215 zcv = 0.5_wp + SIGN( 0.5_wp , pfv_ho(ji,jj,jl) ) 1216 ! 1217 zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv ) 1218 ! 1219 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_ups(ji,jj,jl) 1220 ! 1221 END_2D 1321 1222 1322 1223 END DO … … 1343 1244 ! 1344 1245 DO jl = 1, jpl 1345 DO jj = 2, jpjm1 1346 DO ji = fs_2, fs_jpim1 ! vector opt. 1347 zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1) 1348 END DO 1349 END DO 1350 END DO 1351 CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.) ! lateral boundary cond. 1246 DO_2D( 0, 0, 0, 0 ) 1247 zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1) 1248 END_2D 1249 END DO 1250 CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.0_wp) ! lateral boundary cond. 1352 1251 1353 1252 DO jl = 1, jpl 1354 DO jj = 2, jpjm1 1355 DO ji = fs_2, fs_jpim1 ! vector opt. 1356 uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj) 1357 1358 Rjm = zslpx(ji-1,jj,jl) 1359 Rj = zslpx(ji ,jj,jl) 1360 Rjp = zslpx(ji+1,jj,jl) 1361 1362 IF( np_limiter == 3 ) THEN 1363 1364 IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm 1365 ELSE ; Rr = Rjp 1253 DO_2D( 0, 0, 0, 0 ) 1254 uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj) 1255 1256 Rjm = zslpx(ji-1,jj,jl) 1257 Rj = zslpx(ji ,jj,jl) 1258 Rjp = zslpx(ji+1,jj,jl) 1259 1260 IF( np_limiter == 3 ) THEN 1261 1262 IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm 1263 ELSE ; Rr = Rjp 1264 ENDIF 1265 1266 zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 1267 IF( Rj > 0. ) THEN 1268 zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)), & 1269 & MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 1270 ELSE 1271 zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)), & 1272 & MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 1273 ENDIF 1274 pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter 1275 1276 ELSEIF( np_limiter == 2 ) THEN 1277 IF( Rj /= 0. ) THEN 1278 IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj 1279 ELSE ; Cr = Rjp / Rj 1366 1280 ENDIF 1367 1368 zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 1369 IF( Rj > 0. ) THEN 1370 zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)), & 1371 & MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 1372 ELSE 1373 zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)), & 1374 & MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 1375 ENDIF 1376 pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter 1377 1378 ELSEIF( np_limiter == 2 ) THEN 1379 IF( Rj /= 0. ) THEN 1380 IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj 1381 ELSE ; Cr = Rjp / Rj 1382 ENDIF 1383 ELSE 1384 Cr = 0. 1385 ENDIF 1386 1387 ! -- superbee -- 1388 zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 1389 ! -- van albada 2 -- 1390 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1391 ! -- sweby (with beta=1) -- 1392 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 1393 ! -- van Leer -- 1394 !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 1395 ! -- ospre -- 1396 !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 1397 ! -- koren -- 1398 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 1399 ! -- charm -- 1400 !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 1401 !ELSE ; zpsi = 0. 1402 !ENDIF 1403 ! -- van albada 1 -- 1404 !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 1405 ! -- smart -- 1406 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 1407 ! -- umist -- 1408 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 1409 1410 ! high order flux corrected by the limiter 1411 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5 1412 1281 ELSE 1282 Cr = 0. 1413 1283 ENDIF 1414 END DO 1415 END DO 1416 END DO 1417 CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.) ! lateral boundary cond. 1284 1285 ! -- superbee -- 1286 zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 1287 ! -- van albada 2 -- 1288 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1289 ! -- sweby (with beta=1) -- 1290 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 1291 ! -- van Leer -- 1292 !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 1293 ! -- ospre -- 1294 !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 1295 ! -- koren -- 1296 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 1297 ! -- charm -- 1298 !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 1299 !ELSE ; zpsi = 0. 1300 !ENDIF 1301 ! -- van albada 1 -- 1302 !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 1303 ! -- smart -- 1304 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 1305 ! -- umist -- 1306 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 1307 1308 ! high order flux corrected by the limiter 1309 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5 1310 1311 ENDIF 1312 END_2D 1313 END DO 1314 CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.0_wp) ! lateral boundary cond. 1418 1315 ! 1419 1316 END SUBROUTINE limiter_x … … 1438 1335 ! 1439 1336 DO jl = 1, jpl 1440 DO jj = 2, jpjm1 1441 DO ji = fs_2, fs_jpim1 ! vector opt. 1442 zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1) 1443 END DO 1444 END DO 1445 END DO 1446 CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.) ! lateral boundary cond. 1447 1448 DO jl = 1, jpl 1449 DO jj = 2, jpjm1 1450 DO ji = fs_2, fs_jpim1 ! vector opt. 1451 vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj) 1452 1453 Rjm = zslpy(ji,jj-1,jl) 1454 Rj = zslpy(ji,jj ,jl) 1455 Rjp = zslpy(ji,jj+1,jl) 1456 1457 IF( np_limiter == 3 ) THEN 1458 1459 IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm 1460 ELSE ; Rr = Rjp 1337 DO_2D( 0, 0, 0, 0 ) 1338 zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1) 1339 END_2D 1340 END DO 1341 CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.0_wp) ! lateral boundary cond. 1342 1343 DO jl = 1, jpl 1344 DO_2D( 0, 0, 0, 0 ) 1345 vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj) 1346 1347 Rjm = zslpy(ji,jj-1,jl) 1348 Rj = zslpy(ji,jj ,jl) 1349 Rjp = zslpy(ji,jj+1,jl) 1350 1351 IF( np_limiter == 3 ) THEN 1352 1353 IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm 1354 ELSE ; Rr = Rjp 1355 ENDIF 1356 1357 zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1358 IF( Rj > 0. ) THEN 1359 zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)), & 1360 & MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 1361 ELSE 1362 zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)), & 1363 & MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 1364 ENDIF 1365 pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter 1366 1367 ELSEIF( np_limiter == 2 ) THEN 1368 1369 IF( Rj /= 0. ) THEN 1370 IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj 1371 ELSE ; Cr = Rjp / Rj 1461 1372 ENDIF 1462 1463 zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1464 IF( Rj > 0. ) THEN 1465 zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)), & 1466 & MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 1467 ELSE 1468 zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)), & 1469 & MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 1470 ENDIF 1471 pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter 1472 1473 ELSEIF( np_limiter == 2 ) THEN 1474 1475 IF( Rj /= 0. ) THEN 1476 IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj 1477 ELSE ; Cr = Rjp / Rj 1478 ENDIF 1479 ELSE 1480 Cr = 0. 1481 ENDIF 1482 1483 ! -- superbee -- 1484 zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 1485 ! -- van albada 2 -- 1486 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1487 ! -- sweby (with beta=1) -- 1488 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 1489 ! -- van Leer -- 1490 !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 1491 ! -- ospre -- 1492 !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 1493 ! -- koren -- 1494 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 1495 ! -- charm -- 1496 !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 1497 !ELSE ; zpsi = 0. 1498 !ENDIF 1499 ! -- van albada 1 -- 1500 !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 1501 ! -- smart -- 1502 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 1503 ! -- umist -- 1504 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 1505 1506 ! high order flux corrected by the limiter 1507 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5 1508 1373 ELSE 1374 Cr = 0. 1509 1375 ENDIF 1510 END DO 1511 END DO 1512 END DO 1513 CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.) ! lateral boundary cond. 1376 1377 ! -- superbee -- 1378 zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 1379 ! -- van albada 2 -- 1380 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1381 ! -- sweby (with beta=1) -- 1382 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 1383 ! -- van Leer -- 1384 !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 1385 ! -- ospre -- 1386 !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 1387 ! -- koren -- 1388 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 1389 ! -- charm -- 1390 !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 1391 !ELSE ; zpsi = 0. 1392 !ENDIF 1393 ! -- van albada 1 -- 1394 !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 1395 ! -- smart -- 1396 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 1397 ! -- umist -- 1398 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 1399 1400 ! high order flux corrected by the limiter 1401 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5 1402 1403 ENDIF 1404 END_2D 1405 END DO 1406 CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.0_wp) ! lateral boundary cond. 1514 1407 ! 1515 1408 END SUBROUTINE limiter_y 1516 1409 1517 1410 1518 SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, p sv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i)1411 SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 1519 1412 !!------------------------------------------------------------------- 1520 1413 !! *** ROUTINE Hbig *** … … 1527 1420 !! 2- check whether snow thickness is larger than the surrounding 9-points 1528 1421 !! (before advection) and reduce it by sending the excess in the ocean 1529 !! 3- check whether snow load deplets the snow-ice interface below sea level$1530 !! and reduce it by sending the excess in the ocean1531 !! 4- correct pond fraction to avoid a_ip > a_i1532 1422 !! 1533 1423 !! ** input : Max thickness of the surrounding 9-points … … 1535 1425 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1536 1426 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts 1537 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, p sv_i, poa_i, pa_i, pa_ip, pv_ip1427 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip, pv_ip 1538 1428 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s 1539 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i 1540 ! 1541 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1542 REAL(wp) :: z1_dt, zhip, zhi, zhs, zvs_excess, zfra 1543 REAL(wp), DIMENSION(jpi,jpj) :: zswitch 1429 ! 1430 INTEGER :: ji, jj, jl ! dummy loop indices 1431 REAL(wp) :: z1_dt, zhip, zhi, zhs, zfra 1544 1432 !!------------------------------------------------------------------- 1545 1433 ! … … 1548 1436 DO jl = 1, jpl 1549 1437 1550 DO jj = 1, jpj 1551 DO ji = 1, jpi 1552 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 1438 DO_2D( 1, 1, 1, 1 ) 1439 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 1440 ! 1441 ! ! -- check h_ip -- ! 1442 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 1443 IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 1444 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 1445 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 1446 pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 1447 ENDIF 1448 ENDIF 1449 ! 1450 ! ! -- check h_i -- ! 1451 ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 1452 zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 1453 IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1454 pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m) 1455 ENDIF 1456 ! 1457 ! ! -- check h_s -- ! 1458 ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 1459 zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 1460 IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1461 zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 1553 1462 ! 1554 ! ! -- check h_ip -- ! 1555 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 1556 IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 1557 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 1558 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 1559 pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 1560 ENDIF 1561 ENDIF 1463 wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt 1464 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 1562 1465 ! 1563 ! ! -- check h_i -- ! 1564 ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 1565 zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 1566 IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1567 pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m) 1568 ENDIF 1569 ! 1570 ! ! -- check h_s -- ! 1571 ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 1572 zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 1573 IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1574 zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 1575 ! 1576 wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt 1577 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 1578 ! 1579 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 1580 pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 1581 ENDIF 1582 ! 1583 ! ! -- check snow load -- ! 1584 ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 1585 ! this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 1586 ! this imposed mini can artificially make the snow very thick (if concentration decreases drastically) 1587 zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 1588 IF( zvs_excess > 0._wp ) THEN 1589 zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 1590 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 1591 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 1592 ! 1593 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 1594 pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess 1595 ENDIF 1596 1466 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 1467 pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 1468 ENDIF 1469 ! 1470 ENDIF 1471 END_2D 1472 END DO 1473 ! 1474 END SUBROUTINE Hbig 1475 1476 1477 SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 1478 !!------------------------------------------------------------------- 1479 !! *** ROUTINE Hsnow *** 1480 !! 1481 !! ** Purpose : 1- Check snow load after advection 1482 !! 2- Correct pond concentration to avoid a_ip > a_i 1483 !! 1484 !! ** Method : If snow load makes snow-ice interface to deplet below the ocean surface 1485 !! then put the snow excess in the ocean 1486 !! 1487 !! ** Notes : This correction is crucial because of the call to routine icecor afterwards 1488 !! which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially 1489 !! make the snow very thick (if concentration decreases drastically) 1490 !! This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather 1491 !!------------------------------------------------------------------- 1492 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1493 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip 1494 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s 1495 ! 1496 INTEGER :: ji, jj, jl ! dummy loop indices 1497 REAL(wp) :: z1_dt, zvs_excess, zfra 1498 !!------------------------------------------------------------------- 1499 ! 1500 z1_dt = 1._wp / pdt 1501 ! 1502 ! -- check snow load -- ! 1503 DO jl = 1, jpl 1504 DO_2D( 1, 1, 1, 1 ) 1505 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 1506 ! 1507 zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rho0-rhoi) * r1_rhos ) 1508 ! 1509 IF( zvs_excess > 0._wp ) THEN ! snow-ice interface deplets below the ocean surface 1510 ! put snow excess in the ocean 1511 zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 1512 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 1513 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 1514 ! correct snow volume and heat content 1515 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 1516 pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess 1597 1517 ENDIF 1598 END DO 1599 END DO 1600 END DO 1601 ! !-- correct pond fraction to avoid a_ip > a_i 1518 ! 1519 ENDIF 1520 END_2D 1521 END DO 1522 ! 1523 !-- correct pond concentration to avoid a_ip > a_i -- ! 1602 1524 WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:) 1603 1525 ! 1604 !1605 END SUBROUTINE Hbig 1606 1526 END SUBROUTINE Hsnow 1527 1528 1607 1529 #else 1608 1530 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
