- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_adv_umx.F90 (modified) (41 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/ICE/icedyn_adv_umx.F90
r12178 r12928 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) … … 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 109 DO_2D_00_00 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 125 123 END DO 126 124 CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. ) … … 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_00_00 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_10_10 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 … … 338 332 !== Open water area ==! 339 333 zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 340 DO jj = 2, jpjm1 341 DO ji = fs_2, fs_jpim1 342 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & 343 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 344 END DO 345 END DO 334 DO_2D_00_00 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 346 338 CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1. ) 347 339 ! … … 352 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 ) 353 345 ! 354 ! Make sure ice thickness is not too big 355 ! (because ice thickness can be too large where ice concentration is very small) 356 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 ) 357 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 ! 358 353 END DO 359 354 ! … … 446 441 IF( pamsk == 0._wp ) THEN 447 442 DO jl = 1, jpl 448 DO jj = 1, jpjm1 449 DO ji = 1, fs_jpim1 450 IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 451 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj) 452 zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj) 453 ELSE 454 zfu_ho (ji,jj,jl) = 0._wp 455 zfu_ups(ji,jj,jl) = 0._wp 456 ENDIF 457 ! 458 IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 459 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj) 460 zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj) 461 ELSE 462 zfv_ho (ji,jj,jl) = 0._wp 463 zfv_ups(ji,jj,jl) = 0._wp 464 ENDIF 465 END DO 466 END DO 443 DO_2D_10_10 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 467 460 END DO 468 461 … … 470 463 ! thus we calculate the upstream solution and apply a limiter again 471 464 DO jl = 1, jpl 472 DO jj = 2, jpjm1 473 DO ji = fs_2, fs_jpim1 474 ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) 475 ! 476 zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 477 END DO 478 END DO 465 DO_2D_00_00 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 479 470 END DO 480 471 CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. ) … … 493 484 IF( PRESENT( pua_ho ) ) THEN 494 485 DO jl = 1, jpl 495 DO jj = 1, jpjm1 496 DO ji = 1, fs_jpim1 497 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 498 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 499 END DO 500 END DO 486 DO_2D_10_10 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 501 490 END DO 502 491 ENDIF … … 505 494 ! --------------------------------- 506 495 DO jl = 1, jpl 507 DO jj = 2, jpjm1 508 DO ji = fs_2, fs_jpim1 509 ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) 510 ! 511 ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 512 END DO 513 END DO 496 DO_2D_00_00 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 514 501 END DO 515 502 CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1. ) … … 541 528 ! 542 529 DO jl = 1, jpl 543 DO jj = 1, jpjm1 544 DO ji = 1, fs_jpim1 530 DO_2D_10_10 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_10_10 545 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_00_00 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. ) 555 ! 556 DO jl = 1, jpl !-- flux in y-direction 557 DO_2D_10_10 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_10_10 546 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) 547 END DO 548 END DO 549 END DO 550 ! 551 ELSE !** alternate directions **! 552 ! 553 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_00_00 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. ) 554 579 ! 555 580 DO jl = 1, jpl !-- flux in x-direction 556 DO jj = 1, jpjm1 557 DO ji = 1, fs_jpim1 558 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) 559 END DO 560 END DO 561 END DO 562 ! 563 DO jl = 1, jpl !-- first guess of tracer from u-flux 564 DO jj = 2, jpjm1 565 DO ji = fs_2, fs_jpim1 566 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) ) & 567 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) 568 ! 569 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 570 END DO 571 END DO 572 END DO 573 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 574 ! 575 DO jl = 1, jpl !-- flux in y-direction 576 DO jj = 1, jpjm1 577 DO ji = 1, fs_jpim1 578 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) 579 END DO 580 END DO 581 END DO 582 ! 583 ELSE !== even ice time step: adv_y then adv_x ==! 584 ! 585 DO jl = 1, jpl !-- flux in y-direction 586 DO jj = 1, jpjm1 587 DO ji = 1, fs_jpim1 588 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) 589 END DO 590 END DO 591 END DO 592 ! 593 DO jl = 1, jpl !-- first guess of tracer from v-flux 594 DO jj = 2, jpjm1 595 DO ji = fs_2, fs_jpim1 596 ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) ) & 597 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) 598 ! 599 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 600 END DO 601 END DO 602 END DO 603 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 604 ! 605 DO jl = 1, jpl !-- flux in x-direction 606 DO jj = 1, jpjm1 607 DO ji = 1, fs_jpim1 608 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) 609 END DO 610 END DO 581 DO_2D_10_10 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 611 584 END DO 612 585 ! … … 616 589 ! 617 590 DO jl = 1, jpl !-- after tracer with upstream scheme 618 DO jj = 2, jpjm1 619 DO ji = fs_2, fs_jpim1 620 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj ,jl) & 621 & + pfv_ups(ji,jj,jl) - pfv_ups(ji ,jj-1,jl) ) & 622 & + ( pu (ji,jj ) - pu (ji-1,jj ) & 623 & + pv (ji,jj ) - pv (ji ,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) 624 ! 625 pt_ups(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 626 END DO 627 END DO 591 DO_2D_00_00 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 628 599 END DO 629 600 CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1. ) … … 657 628 ! 658 629 DO jl = 1, jpl 659 DO jj = 1, jpjm1 660 DO ji = 1, fs_jpim1 661 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj ,jl) ) 662 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji ,jj+1,jl) ) 663 END DO 664 END DO 630 DO_2D_10_10 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 665 634 END DO 666 635 ! … … 677 646 ! 678 647 DO jl = 1, jpl !-- flux in x-direction 679 DO jj = 1, jpjm1 680 DO ji = 1, fs_jpim1 681 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 682 END DO 683 END DO 648 DO_2D_10_10 649 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 650 END_2D 684 651 END DO 685 652 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 686 653 687 654 DO jl = 1, jpl !-- first guess of tracer from u-flux 688 DO jj = 2, jpjm1 689 DO ji = fs_2, fs_jpim1 690 ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) & 691 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) 692 ! 693 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 694 END DO 695 END DO 655 DO_2D_00_00 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 696 661 END DO 697 662 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 698 663 699 664 DO jl = 1, jpl !-- flux in y-direction 700 DO jj = 1, jpjm1 701 DO ji = 1, fs_jpim1 702 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 703 END DO 704 END DO 665 DO_2D_10_10 666 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 667 END_2D 705 668 END DO 706 669 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) … … 709 672 ! 710 673 DO jl = 1, jpl !-- flux in y-direction 711 DO jj = 1, jpjm1 712 DO ji = 1, fs_jpim1 713 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 714 END DO 715 END DO 674 DO_2D_10_10 675 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 676 END_2D 716 677 END DO 717 678 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 718 679 ! 719 680 DO jl = 1, jpl !-- first guess of tracer from v-flux 720 DO jj = 2, jpjm1 721 DO ji = fs_2, fs_jpim1 722 ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) & 723 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) 724 ! 725 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 726 END DO 727 END DO 681 DO_2D_00_00 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 728 687 END DO 729 688 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 730 689 ! 731 690 DO jl = 1, jpl !-- flux in x-direction 732 DO jj = 1, jpjm1 733 DO ji = 1, fs_jpim1 734 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 735 END DO 736 END DO 691 DO_2D_10_10 692 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 693 END_2D 737 694 END DO 738 695 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) … … 780 737 ! !-- advective form update in zpt --! 781 738 DO jl = 1, jpl 782 DO jj = 2, jpjm1 783 DO ji = fs_2, fs_jpim1 784 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) & 785 & + pt (ji,jj,jl) * ( pu (ji,jj ) - pu (ji-1,jj ) ) * r1_e1e2t(ji,jj) & 786 & * pamsk & 787 & ) * pdt ) * tmask(ji,jj,1) 788 END DO 789 END DO 739 DO_2D_00_00 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 790 745 END DO 791 746 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) … … 809 764 ! !-- advective form update in zpt --! 810 765 DO jl = 1, jpl 811 DO jj = 2, jpjm1 812 DO ji = fs_2, fs_jpim1 813 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) & 814 & + pt (ji,jj,jl) * ( pv (ji,jj ) - pv (ji,jj-1 ) ) * r1_e1e2t(ji,jj) & 815 & * pamsk & 816 & ) * pdt ) * tmask(ji,jj,1) 817 END DO 818 END DO 766 DO_2D_00_00 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 819 772 END DO 820 773 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) … … 862 815 DO jl = 1, jpl 863 816 DO jj = 2, jpjm1 ! First derivative (gradient) 864 DO ji = 1, fs_jpim1817 DO ji = 1, jpim1 865 818 ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 866 819 END DO 867 820 ! ! Second derivative (Laplacian) 868 DO ji = fs_2, fs_jpim1821 DO ji = 2, jpim1 869 822 ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj) 870 823 END DO … … 876 829 DO jl = 1, jpl 877 830 DO jj = 2, jpjm1 ! Third derivative 878 DO ji = 1, fs_jpim1831 DO ji = 1, jpim1 879 832 ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 880 833 END DO 881 834 ! ! Fourth derivative 882 DO ji = fs_2, fs_jpim1835 DO ji = 2, jpim1 883 836 ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj) 884 837 END DO … … 893 846 ! 894 847 DO jl = 1, jpl 895 DO jj = 1, jpjm1 896 DO ji = 1, fs_jpim1 ! vector opt. 897 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 898 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 899 END DO 900 END DO 848 DO_2D_10_10 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 901 852 END DO 902 853 ! … … 904 855 ! 905 856 DO jl = 1, jpl 906 DO jj = 1, jpjm1 907 DO ji = 1, fs_jpim1 ! vector opt. 908 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 909 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 910 & - zcu * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 911 END DO 912 END DO 857 DO_2D_10_10 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 913 862 END DO 914 863 ! … … 916 865 ! 917 866 DO jl = 1, jpl 918 DO jj = 1, jpjm1 919 DO ji = 1, fs_jpim1 ! vector opt. 920 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 921 zdx2 = e1u(ji,jj) * e1u(ji,jj) 867 DO_2D_10_10 868 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 869 zdx2 = e1u(ji,jj) * e1u(ji,jj) 922 870 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 923 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 924 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 925 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 926 & - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 927 END DO 928 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 929 876 END DO 930 877 ! … … 932 879 ! 933 880 DO jl = 1, jpl 934 DO jj = 1, jpjm1 935 DO ji = 1, fs_jpim1 ! vector opt. 936 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 937 zdx2 = e1u(ji,jj) * e1u(ji,jj) 881 DO_2D_10_10 882 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 883 zdx2 = e1u(ji,jj) * e1u(ji,jj) 938 884 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 939 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 940 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 941 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 942 & - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 943 END DO 944 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 945 890 END DO 946 891 ! … … 948 893 ! 949 894 DO jl = 1, jpl 950 DO jj = 1, jpjm1 951 DO ji = 1, fs_jpim1 ! vector opt. 952 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 953 zdx2 = e1u(ji,jj) * e1u(ji,jj) 895 DO_2D_10_10 896 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 897 zdx2 = e1u(ji,jj) * e1u(ji,jj) 954 898 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 955 zdx4 = zdx2 * zdx2 956 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 957 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 958 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 959 & - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) & 960 & + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) & 961 & - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) ) 962 END DO 963 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 964 907 END DO 965 908 ! … … 971 914 IF( ll_neg ) THEN 972 915 DO jl = 1, jpl 973 DO jj = 1, jpjm1 974 DO ji = 1, fs_jpim1 975 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 976 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 977 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 978 ENDIF 979 END DO 980 END DO 916 DO_2D_10_10 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 981 922 END DO 982 923 ENDIF 983 924 ! !-- High order flux in i-direction --! 984 925 DO jl = 1, jpl 985 DO jj = 1, jpjm1 986 DO ji = 1, fs_jpim1 ! vector opt. 987 pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 988 END DO 989 END DO 926 DO_2D_10_10 927 pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 928 END_2D 990 929 END DO 991 930 ! … … 1018 957 ! !-- Laplacian in j-direction --! 1019 958 DO jl = 1, jpl 1020 DO jj = 1, jpjm1 ! First derivative (gradient) 1021 DO ji = fs_2, fs_jpim1 1022 ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 1023 END DO 1024 END DO 1025 DO jj = 2, jpjm1 ! Second derivative (Laplacian) 1026 DO ji = fs_2, fs_jpim1 1027 ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj) 1028 END DO 1029 END DO 959 DO_2D_10_00 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_00_00 963 ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj) 964 END_2D 1030 965 END DO 1031 966 CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1. ) … … 1033 968 ! !-- BiLaplacian in j-direction --! 1034 969 DO jl = 1, jpl 1035 DO jj = 1, jpjm1 ! First derivative 1036 DO ji = fs_2, fs_jpim1 1037 ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 1038 END DO 1039 END DO 1040 DO jj = 2, jpjm1 ! Second derivative 1041 DO ji = fs_2, fs_jpim1 1042 ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj) 1043 END DO 1044 END DO 970 DO_2D_10_00 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_00_00 974 ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj) 975 END_2D 1045 976 END DO 1046 977 CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1. ) … … 1051 982 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 1052 983 DO jl = 1, jpl 1053 DO jj = 1, jpjm1 1054 DO ji = 1, fs_jpim1 1055 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & 1056 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 1057 END DO 1058 END DO 984 DO_2D_10_10 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 1059 988 END DO 1060 989 ! 1061 990 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 1062 991 DO jl = 1, jpl 1063 DO jj = 1, jpjm1 1064 DO ji = 1, fs_jpim1 1065 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1066 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & 1067 & - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 1068 END DO 1069 END DO 992 DO_2D_10_10 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 1070 997 END DO 1071 998 ! 1072 999 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 1073 1000 DO jl = 1, jpl 1074 DO jj = 1, jpjm1 1075 DO ji = 1, fs_jpim1 1076 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1077 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1001 DO_2D_10_10 1002 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1003 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1078 1004 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 1079 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1080 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 1081 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 1082 & - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 1083 END DO 1084 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 1085 1010 END DO 1086 1011 ! 1087 1012 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 1088 1013 DO jl = 1, jpl 1089 DO jj = 1, jpjm1 1090 DO ji = 1, fs_jpim1 1091 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1092 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1014 DO_2D_10_10 1015 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1016 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1093 1017 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 1094 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1095 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 1096 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 1097 & - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 1098 END DO 1099 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 1100 1023 END DO 1101 1024 ! 1102 1025 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 1103 1026 DO jl = 1, jpl 1104 DO jj = 1, jpjm1 1105 DO ji = 1, fs_jpim1 1106 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1107 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1027 DO_2D_10_10 1028 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1029 zdy2 = e2v(ji,jj) * e2v(ji,jj) 1108 1030 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 1109 zdy4 = zdy2 * zdy2 1110 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1111 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 1112 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 1113 & - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) & 1114 & + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) & 1115 & - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) ) 1116 END DO 1117 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 1118 1039 END DO 1119 1040 ! … … 1125 1046 IF( ll_neg ) THEN 1126 1047 DO jl = 1, jpl 1127 DO jj = 1, jpjm1 1128 DO ji = 1, fs_jpim1 1129 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 1130 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & 1131 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 1132 ENDIF 1133 END DO 1134 END DO 1048 DO_2D_10_10 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 1135 1054 END DO 1136 1055 ENDIF 1137 1056 ! !-- High order flux in j-direction --! 1138 1057 DO jl = 1, jpl 1139 DO jj = 1, jpjm1 1140 DO ji = 1, fs_jpim1 ! vector opt. 1141 pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 1142 END DO 1143 END DO 1058 DO_2D_10_10 1059 pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 1060 END_2D 1144 1061 END DO 1145 1062 ! … … 1175 1092 ! -------------------------------------------------- 1176 1093 DO jl = 1, jpl 1177 DO jj = 1, jpjm1 1178 DO ji = 1, fs_jpim1 ! vector opt. 1179 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 1180 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1181 END DO 1182 END DO 1094 DO_2D_10_10 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 1183 1098 END DO 1184 1099 … … 1194 1109 1195 1110 DO jl = 1, jpl 1196 DO jj = 2, jpjm1 1197 DO ji = fs_2, fs_jpim1 1198 zti_ups(ji,jj,jl)= pt_ups(ji+1,jj ,jl) 1199 ztj_ups(ji,jj,jl)= pt_ups(ji ,jj+1,jl) 1200 END DO 1201 END DO 1111 DO_2D_00_00 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 1202 1115 END DO 1203 1116 CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1., ztj_ups, 'T', 1. ) 1204 1117 1205 1118 DO jl = 1, jpl 1206 DO jj = 2, jpjm1 1207 DO ji = fs_2, fs_jpim1 1208 IF ( pfu_ho(ji,jj,jl) * ( pt_ups(ji+1,jj ,jl) - pt_ups(ji,jj,jl) ) <= 0._wp .AND. & 1209 & pfv_ho(ji,jj,jl) * ( pt_ups(ji ,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0._wp ) THEN 1210 ! 1211 IF( pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj ,jl) - zti_ups(ji,jj,jl) ) <= 0._wp .AND. & 1212 & pfv_ho(ji,jj,jl) * ( ztj_ups(ji ,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0._wp ) THEN 1213 pfu_ho(ji,jj,jl)=0._wp 1214 pfv_ho(ji,jj,jl)=0._wp 1215 ENDIF 1216 ! 1217 IF( pfu_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji-1,jj ,jl) ) <= 0._wp .AND. & 1218 & pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji ,jj-1,jl) ) <= 0._wp ) THEN 1219 pfu_ho(ji,jj,jl)=0._wp 1220 pfv_ho(ji,jj,jl)=0._wp 1221 ENDIF 1222 ! 1119 DO_2D_00_00 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 1223 1127 ENDIF 1224 END DO 1225 END DO 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 1226 1137 END DO 1227 1138 CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond. … … 1235 1146 DO jl = 1, jpl 1236 1147 1237 DO jj = 1, jpj 1238 DO ji = 1, jpi 1239 IF ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN 1240 zbup(ji,jj) = -zbig 1241 zbdo(ji,jj) = zbig 1242 ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) > 0._wp ) THEN 1243 zbup(ji,jj) = pt_ups(ji,jj,jl) 1244 zbdo(ji,jj) = pt_ups(ji,jj,jl) 1245 ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN 1246 zbup(ji,jj) = pt(ji,jj,jl) 1247 zbdo(ji,jj) = pt(ji,jj,jl) 1248 ELSE 1249 zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ups(ji,jj,jl) ) 1250 zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ups(ji,jj,jl) ) 1251 ENDIF 1252 END DO 1253 END DO 1254 1255 DO jj = 2, jpjm1 1256 DO ji = fs_2, fs_jpim1 ! vector opt. 1257 ! 1258 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 1259 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) ) 1260 ! 1261 zpos = MAX( 0._wp, pfu_ho(ji-1,jj ,jl) ) - MIN( 0._wp, pfu_ho(ji ,jj ,jl) ) & ! positive/negative part of the flux 1262 & + MAX( 0._wp, pfv_ho(ji ,jj-1,jl) ) - MIN( 0._wp, pfv_ho(ji ,jj ,jl) ) 1263 zneg = MAX( 0._wp, pfu_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfu_ho(ji-1,jj ,jl) ) & 1264 & + MAX( 0._wp, pfv_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfv_ho(ji ,jj-1,jl) ) 1265 ! 1266 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) ) & 1267 & ) * ( 1. - pamsk ) 1268 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) ) & 1269 & ) * ( 1. - pamsk ) 1270 ! 1271 ! ! up & down beta terms 1272 ! 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 !!!) 1273 IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 1274 ELSE ; zbetup(ji,jj,jl) = 0._wp ! zbig 1275 ENDIF 1276 ! 1277 IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 1278 ELSE ; zbetdo(ji,jj,jl) = 0._wp ! zbig 1279 ENDIF 1280 ! 1281 ! if all the points are outside ice cover 1282 IF( zup == -zbig ) zbetup(ji,jj,jl) = 0._wp ! zbig 1283 IF( zdo == zbig ) zbetdo(ji,jj,jl) = 0._wp ! zbig 1284 ! 1285 END DO 1286 END DO 1148 DO_2D_11_11 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_00_00 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 1287 1194 END DO 1288 1195 CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign) … … 1292 1199 ! --------------------------------- 1293 1200 DO jl = 1, jpl 1294 DO jj = 1, jpjm1 1295 DO ji = 1, fs_jpim1 ! vector opt. 1296 zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 1297 zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) 1298 zcu = 0.5_wp + SIGN( 0.5_wp , pfu_ho(ji,jj,jl) ) 1299 ! 1300 zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu ) 1301 ! 1302 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_ups(ji,jj,jl) 1303 ! 1304 END DO 1305 END DO 1306 1307 DO jj = 1, jpjm1 1308 DO ji = 1, fs_jpim1 ! vector opt. 1309 zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 1310 zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) 1311 zcv = 0.5_wp + SIGN( 0.5_wp , pfv_ho(ji,jj,jl) ) 1312 ! 1313 zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv ) 1314 ! 1315 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_ups(ji,jj,jl) 1316 ! 1317 END DO 1318 END DO 1201 DO_2D_10_10 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_10_10 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 1319 1222 1320 1223 END DO … … 1341 1244 ! 1342 1245 DO jl = 1, jpl 1343 DO jj = 2, jpjm1 1344 DO ji = fs_2, fs_jpim1 ! vector opt. 1345 zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1) 1346 END DO 1347 END DO 1246 DO_2D_00_00 1247 zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1) 1248 END_2D 1348 1249 END DO 1349 1250 CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.) ! lateral boundary cond. 1350 1251 1351 1252 DO jl = 1, jpl 1352 DO jj = 2, jpjm1 1353 DO ji = fs_2, fs_jpim1 ! vector opt. 1354 uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj) 1355 1356 Rjm = zslpx(ji-1,jj,jl) 1357 Rj = zslpx(ji ,jj,jl) 1358 Rjp = zslpx(ji+1,jj,jl) 1359 1360 IF( np_limiter == 3 ) THEN 1361 1362 IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm 1363 ELSE ; Rr = Rjp 1253 DO_2D_00_00 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 1364 1280 ENDIF 1365 1366 zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 1367 IF( Rj > 0. ) THEN 1368 zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)), & 1369 & MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 1370 ELSE 1371 zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)), & 1372 & MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 1373 ENDIF 1374 pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter 1375 1376 ELSEIF( np_limiter == 2 ) THEN 1377 IF( Rj /= 0. ) THEN 1378 IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj 1379 ELSE ; Cr = Rjp / Rj 1380 ENDIF 1381 ELSE 1382 Cr = 0. 1383 ENDIF 1384 1385 ! -- superbee -- 1386 zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 1387 ! -- van albada 2 -- 1388 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1389 ! -- sweby (with beta=1) -- 1390 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 1391 ! -- van Leer -- 1392 !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 1393 ! -- ospre -- 1394 !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 1395 ! -- koren -- 1396 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 1397 ! -- charm -- 1398 !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 1399 !ELSE ; zpsi = 0. 1400 !ENDIF 1401 ! -- van albada 1 -- 1402 !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 1403 ! -- smart -- 1404 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 1405 ! -- umist -- 1406 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 1407 1408 ! high order flux corrected by the limiter 1409 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5 1410 1281 ELSE 1282 Cr = 0. 1411 1283 ENDIF 1412 END DO 1413 END DO 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 1414 1313 END DO 1415 1314 CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.) ! lateral boundary cond. … … 1436 1335 ! 1437 1336 DO jl = 1, jpl 1438 DO jj = 2, jpjm1 1439 DO ji = fs_2, fs_jpim1 ! vector opt. 1440 zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1) 1441 END DO 1442 END DO 1337 DO_2D_00_00 1338 zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1) 1339 END_2D 1443 1340 END DO 1444 1341 CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.) ! lateral boundary cond. 1445 1342 1446 1343 DO jl = 1, jpl 1447 DO jj = 2, jpjm1 1448 DO ji = fs_2, fs_jpim1 ! vector opt. 1449 vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj) 1450 1451 Rjm = zslpy(ji,jj-1,jl) 1452 Rj = zslpy(ji,jj ,jl) 1453 Rjp = zslpy(ji,jj+1,jl) 1454 1455 IF( np_limiter == 3 ) THEN 1456 1457 IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm 1458 ELSE ; Rr = Rjp 1344 DO_2D_00_00 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 1459 1372 ENDIF 1460 1461 zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1462 IF( Rj > 0. ) THEN 1463 zlimiter = MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)), & 1464 & MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)), zh3, 1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 1465 ELSE 1466 zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)), & 1467 & MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 1468 ENDIF 1469 pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter 1470 1471 ELSEIF( np_limiter == 2 ) THEN 1472 1473 IF( Rj /= 0. ) THEN 1474 IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj 1475 ELSE ; Cr = Rjp / Rj 1476 ENDIF 1477 ELSE 1478 Cr = 0. 1479 ENDIF 1480 1481 ! -- superbee -- 1482 zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 1483 ! -- van albada 2 -- 1484 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1485 ! -- sweby (with beta=1) -- 1486 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 1487 ! -- van Leer -- 1488 !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 1489 ! -- ospre -- 1490 !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 1491 ! -- koren -- 1492 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 1493 ! -- charm -- 1494 !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 1495 !ELSE ; zpsi = 0. 1496 !ENDIF 1497 ! -- van albada 1 -- 1498 !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 1499 ! -- smart -- 1500 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 1501 ! -- umist -- 1502 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 1503 1504 ! high order flux corrected by the limiter 1505 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5 1506 1373 ELSE 1374 Cr = 0. 1507 1375 ENDIF 1508 END DO 1509 END DO 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 1510 1405 END DO 1511 1406 CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.) ! lateral boundary cond. … … 1514 1409 1515 1410 1516 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 ) 1517 1412 !!------------------------------------------------------------------- 1518 1413 !! *** ROUTINE Hbig *** … … 1525 1420 !! 2- check whether snow thickness is larger than the surrounding 9-points 1526 1421 !! (before advection) and reduce it by sending the excess in the ocean 1527 !! 3- check whether snow load deplets the snow-ice interface below sea level$1528 !! and reduce it by sending the excess in the ocean1529 !! 4- correct pond concentration to avoid a_ip > a_i1530 1422 !! 1531 1423 !! ** input : Max thickness of the surrounding 9-points … … 1533 1425 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1534 1426 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts 1535 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 1536 1428 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s 1537 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i 1538 ! 1539 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1540 REAL(wp) :: z1_dt, zhip, zhi, zhs, zvs_excess, zfra 1541 REAL(wp), DIMENSION(jpi,jpj) :: zswitch 1429 ! 1430 INTEGER :: ji, jj, jl ! dummy loop indices 1431 REAL(wp) :: z1_dt, zhip, zhi, zhs, zfra 1542 1432 !!------------------------------------------------------------------- 1543 1433 ! … … 1546 1436 DO jl = 1, jpl 1547 1437 1548 DO jj = 1, jpj 1549 DO ji = 1, jpi 1550 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 1438 DO_2D_11_11 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 ) 1551 1462 ! 1552 ! ! -- check h_ip -- ! 1553 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 1554 IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 1555 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 1556 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 1557 pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 1558 ENDIF 1559 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 1560 1465 ! 1561 ! ! -- check h_i -- ! 1562 ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 1563 zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 1564 IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1565 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) 1566 ENDIF 1567 ! 1568 ! ! -- check h_s -- ! 1569 ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 1570 zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 1571 IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1572 zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 1573 ! 1574 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 1575 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 1576 ! 1577 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 1578 pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 1579 ENDIF 1580 ! 1581 ! ! -- check snow load -- ! 1582 ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 1583 ! this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 1584 ! this imposed mini can artificially make the snow very thick (if concentration decreases drastically) 1585 zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 1586 IF( zvs_excess > 0._wp ) THEN 1587 zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 1588 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 1589 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 1590 ! 1591 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 1592 pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess 1593 ENDIF 1594 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_11_11 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 1595 1517 ENDIF 1596 END DO 1597 END DO 1598 END DO 1599 ! !-- correct pond concentration 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 -- ! 1600 1524 WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:) 1601 1525 ! 1602 !1603 END SUBROUTINE Hbig 1604 1526 END SUBROUTINE Hsnow 1527 1528 1605 1529 #else 1606 1530 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
