Changeset 12377 for NEMO/trunk/src/ICE/icedyn_adv_umx.F90
- Timestamp:
- 2020-02-12T15:39:06+01:00 (4 years ago)
- Location:
- NEMO/trunk
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_adv_umx.F90 (modified) (36 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEAD ext/AGRIF5 ^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL
-
- Property svn:externals
-
NEMO/trunk/src/ICE/icedyn_adv_umx.F90
r12197 r12377 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. ) … … 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 ! … … 449 441 IF( pamsk == 0._wp ) THEN 450 442 DO jl = 1, jpl 451 DO jj = 1, jpjm1 452 DO ji = 1, fs_jpim1 453 IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 454 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj) 455 zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj) 456 ELSE 457 zfu_ho (ji,jj,jl) = 0._wp 458 zfu_ups(ji,jj,jl) = 0._wp 459 ENDIF 460 ! 461 IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 462 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj) 463 zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj) 464 ELSE 465 zfv_ho (ji,jj,jl) = 0._wp 466 zfv_ups(ji,jj,jl) = 0._wp 467 ENDIF 468 END DO 469 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 470 460 END DO 471 461 … … 473 463 ! thus we calculate the upstream solution and apply a limiter again 474 464 DO jl = 1, jpl 475 DO jj = 2, jpjm1 476 DO ji = fs_2, fs_jpim1 477 ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) 478 ! 479 zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 480 END DO 481 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 482 470 END DO 483 471 CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. ) … … 496 484 IF( PRESENT( pua_ho ) ) THEN 497 485 DO jl = 1, jpl 498 DO jj = 1, jpjm1 499 DO ji = 1, fs_jpim1 500 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 501 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 502 END DO 503 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 504 490 END DO 505 491 ENDIF … … 508 494 ! --------------------------------- 509 495 DO jl = 1, jpl 510 DO jj = 2, jpjm1 511 DO ji = fs_2, fs_jpim1 512 ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) 513 ! 514 ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 515 END DO 516 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 517 501 END DO 518 502 CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1. ) … … 544 528 ! 545 529 DO jl = 1, jpl 546 DO jj = 1, jpjm1 547 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 548 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 549 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) 550 END DO 551 END DO 552 END DO 553 ! 554 ELSE !** alternate directions **! 555 ! 556 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. ) 557 579 ! 558 580 DO jl = 1, jpl !-- flux in x-direction 559 DO jj = 1, jpjm1 560 DO ji = 1, fs_jpim1 561 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) 562 END DO 563 END DO 564 END DO 565 ! 566 DO jl = 1, jpl !-- first guess of tracer from u-flux 567 DO jj = 2, jpjm1 568 DO ji = fs_2, fs_jpim1 569 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) ) & 570 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) 571 ! 572 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 573 END DO 574 END DO 575 END DO 576 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 577 ! 578 DO jl = 1, jpl !-- flux in y-direction 579 DO jj = 1, jpjm1 580 DO ji = 1, fs_jpim1 581 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) 582 END DO 583 END DO 584 END DO 585 ! 586 ELSE !== even ice time step: adv_y then adv_x ==! 587 ! 588 DO jl = 1, jpl !-- flux in y-direction 589 DO jj = 1, jpjm1 590 DO ji = 1, fs_jpim1 591 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) 592 END DO 593 END DO 594 END DO 595 ! 596 DO jl = 1, jpl !-- first guess of tracer from v-flux 597 DO jj = 2, jpjm1 598 DO ji = fs_2, fs_jpim1 599 ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) ) & 600 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) 601 ! 602 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 603 END DO 604 END DO 605 END DO 606 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 607 ! 608 DO jl = 1, jpl !-- flux in x-direction 609 DO jj = 1, jpjm1 610 DO ji = 1, fs_jpim1 611 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) 612 END DO 613 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 614 584 END DO 615 585 ! … … 619 589 ! 620 590 DO jl = 1, jpl !-- after tracer with upstream scheme 621 DO jj = 2, jpjm1 622 DO ji = fs_2, fs_jpim1 623 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj ,jl) & 624 & + pfv_ups(ji,jj,jl) - pfv_ups(ji ,jj-1,jl) ) & 625 & + ( pu (ji,jj ) - pu (ji-1,jj ) & 626 & + pv (ji,jj ) - pv (ji ,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) 627 ! 628 pt_ups(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 629 END DO 630 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 631 599 END DO 632 600 CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1. ) … … 660 628 ! 661 629 DO jl = 1, jpl 662 DO jj = 1, jpjm1 663 DO ji = 1, fs_jpim1 664 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj ,jl) ) 665 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji ,jj+1,jl) ) 666 END DO 667 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 668 634 END DO 669 635 ! … … 680 646 ! 681 647 DO jl = 1, jpl !-- flux in x-direction 682 DO jj = 1, jpjm1 683 DO ji = 1, fs_jpim1 684 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 685 END DO 686 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 687 651 END DO 688 652 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 689 653 690 654 DO jl = 1, jpl !-- first guess of tracer from u-flux 691 DO jj = 2, jpjm1 692 DO ji = fs_2, fs_jpim1 693 ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) & 694 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) 695 ! 696 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 697 END DO 698 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 699 661 END DO 700 662 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 701 663 702 664 DO jl = 1, jpl !-- flux in y-direction 703 DO jj = 1, jpjm1 704 DO ji = 1, fs_jpim1 705 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 706 END DO 707 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 708 668 END DO 709 669 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) … … 712 672 ! 713 673 DO jl = 1, jpl !-- flux in y-direction 714 DO jj = 1, jpjm1 715 DO ji = 1, fs_jpim1 716 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 717 END DO 718 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 719 677 END DO 720 678 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 721 679 ! 722 680 DO jl = 1, jpl !-- first guess of tracer from v-flux 723 DO jj = 2, jpjm1 724 DO ji = fs_2, fs_jpim1 725 ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) & 726 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) 727 ! 728 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 729 END DO 730 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 731 687 END DO 732 688 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 733 689 ! 734 690 DO jl = 1, jpl !-- flux in x-direction 735 DO jj = 1, jpjm1 736 DO ji = 1, fs_jpim1 737 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 738 END DO 739 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 740 694 END DO 741 695 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) … … 783 737 ! !-- advective form update in zpt --! 784 738 DO jl = 1, jpl 785 DO jj = 2, jpjm1 786 DO ji = fs_2, fs_jpim1 787 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) & 788 & + pt (ji,jj,jl) * ( pu (ji,jj ) - pu (ji-1,jj ) ) * r1_e1e2t(ji,jj) & 789 & * pamsk & 790 & ) * pdt ) * tmask(ji,jj,1) 791 END DO 792 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 793 745 END DO 794 746 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) … … 812 764 ! !-- advective form update in zpt --! 813 765 DO jl = 1, jpl 814 DO jj = 2, jpjm1 815 DO ji = fs_2, fs_jpim1 816 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) & 817 & + pt (ji,jj,jl) * ( pv (ji,jj ) - pv (ji,jj-1 ) ) * r1_e1e2t(ji,jj) & 818 & * pamsk & 819 & ) * pdt ) * tmask(ji,jj,1) 820 END DO 821 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 822 772 END DO 823 773 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) … … 865 815 DO jl = 1, jpl 866 816 DO jj = 2, jpjm1 ! First derivative (gradient) 867 DO ji = 1, fs_jpim1817 DO ji = 1, jpim1 868 818 ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 869 819 END DO 870 820 ! ! Second derivative (Laplacian) 871 DO ji = fs_2, fs_jpim1821 DO ji = 2, jpim1 872 822 ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj) 873 823 END DO … … 879 829 DO jl = 1, jpl 880 830 DO jj = 2, jpjm1 ! Third derivative 881 DO ji = 1, fs_jpim1831 DO ji = 1, jpim1 882 832 ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 883 833 END DO 884 834 ! ! Fourth derivative 885 DO ji = fs_2, fs_jpim1835 DO ji = 2, jpim1 886 836 ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj) 887 837 END DO … … 896 846 ! 897 847 DO jl = 1, jpl 898 DO jj = 1, jpjm1 899 DO ji = 1, fs_jpim1 ! vector opt. 900 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 901 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 902 END DO 903 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 904 852 END DO 905 853 ! … … 907 855 ! 908 856 DO jl = 1, jpl 909 DO jj = 1, jpjm1 910 DO ji = 1, fs_jpim1 ! vector opt. 911 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 912 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 913 & - zcu * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 914 END DO 915 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 916 862 END DO 917 863 ! … … 919 865 ! 920 866 DO jl = 1, jpl 921 DO jj = 1, jpjm1 922 DO ji = 1, fs_jpim1 ! vector opt. 923 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 924 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) 925 870 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 926 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 927 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 928 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 929 & - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 930 END DO 931 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 932 876 END DO 933 877 ! … … 935 879 ! 936 880 DO jl = 1, jpl 937 DO jj = 1, jpjm1 938 DO ji = 1, fs_jpim1 ! vector opt. 939 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 940 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) 941 884 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 942 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 943 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 944 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 945 & - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 946 END DO 947 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 948 890 END DO 949 891 ! … … 951 893 ! 952 894 DO jl = 1, jpl 953 DO jj = 1, jpjm1 954 DO ji = 1, fs_jpim1 ! vector opt. 955 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 956 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) 957 898 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 958 zdx4 = zdx2 * zdx2 959 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 960 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 961 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 962 & - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) & 963 & + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) & 964 & - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) ) 965 END DO 966 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 967 907 END DO 968 908 ! … … 974 914 IF( ll_neg ) THEN 975 915 DO jl = 1, jpl 976 DO jj = 1, jpjm1 977 DO ji = 1, fs_jpim1 978 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 979 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 980 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 981 ENDIF 982 END DO 983 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 984 922 END DO 985 923 ENDIF 986 924 ! !-- High order flux in i-direction --! 987 925 DO jl = 1, jpl 988 DO jj = 1, jpjm1 989 DO ji = 1, fs_jpim1 ! vector opt. 990 pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 991 END DO 992 END DO 926 DO_2D_10_10 927 pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 928 END_2D 993 929 END DO 994 930 ! … … 1021 957 ! !-- Laplacian in j-direction --! 1022 958 DO jl = 1, jpl 1023 DO jj = 1, jpjm1 ! First derivative (gradient) 1024 DO ji = fs_2, fs_jpim1 1025 ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 1026 END DO 1027 END DO 1028 DO jj = 2, jpjm1 ! Second derivative (Laplacian) 1029 DO ji = fs_2, fs_jpim1 1030 ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj) 1031 END DO 1032 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 1033 965 END DO 1034 966 CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1. ) … … 1036 968 ! !-- BiLaplacian in j-direction --! 1037 969 DO jl = 1, jpl 1038 DO jj = 1, jpjm1 ! First derivative 1039 DO ji = fs_2, fs_jpim1 1040 ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 1041 END DO 1042 END DO 1043 DO jj = 2, jpjm1 ! Second derivative 1044 DO ji = fs_2, fs_jpim1 1045 ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj) 1046 END DO 1047 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 1048 976 END DO 1049 977 CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1. ) … … 1054 982 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 1055 983 DO jl = 1, jpl 1056 DO jj = 1, jpjm1 1057 DO ji = 1, fs_jpim1 1058 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & 1059 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 1060 END DO 1061 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 1062 988 END DO 1063 989 ! 1064 990 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 1065 991 DO jl = 1, jpl 1066 DO jj = 1, jpjm1 1067 DO ji = 1, fs_jpim1 1068 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1069 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & 1070 & - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 1071 END DO 1072 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 1073 997 END DO 1074 998 ! 1075 999 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 1076 1000 DO jl = 1, jpl 1077 DO jj = 1, jpjm1 1078 DO ji = 1, fs_jpim1 1079 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1080 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) 1081 1004 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 1082 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1083 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 1084 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 1085 & - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 1086 END DO 1087 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 1088 1010 END DO 1089 1011 ! 1090 1012 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 1091 1013 DO jl = 1, jpl 1092 DO jj = 1, jpjm1 1093 DO ji = 1, fs_jpim1 1094 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1095 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) 1096 1017 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 1097 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1098 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 1099 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 1100 & - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 1101 END DO 1102 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 1103 1023 END DO 1104 1024 ! 1105 1025 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 1106 1026 DO jl = 1, jpl 1107 DO jj = 1, jpjm1 1108 DO ji = 1, fs_jpim1 1109 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1110 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) 1111 1030 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 1112 zdy4 = zdy2 * zdy2 1113 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 1114 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 1115 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 1116 & - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) & 1117 & + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) & 1118 & - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) ) 1119 END DO 1120 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 1121 1039 END DO 1122 1040 ! … … 1128 1046 IF( ll_neg ) THEN 1129 1047 DO jl = 1, jpl 1130 DO jj = 1, jpjm1 1131 DO ji = 1, fs_jpim1 1132 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 1133 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & 1134 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 1135 ENDIF 1136 END DO 1137 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 1138 1054 END DO 1139 1055 ENDIF 1140 1056 ! !-- High order flux in j-direction --! 1141 1057 DO jl = 1, jpl 1142 DO jj = 1, jpjm1 1143 DO ji = 1, fs_jpim1 ! vector opt. 1144 pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 1145 END DO 1146 END DO 1058 DO_2D_10_10 1059 pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 1060 END_2D 1147 1061 END DO 1148 1062 ! … … 1178 1092 ! -------------------------------------------------- 1179 1093 DO jl = 1, jpl 1180 DO jj = 1, jpjm1 1181 DO ji = 1, fs_jpim1 ! vector opt. 1182 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 1183 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1184 END DO 1185 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 1186 1098 END DO 1187 1099 … … 1197 1109 1198 1110 DO jl = 1, jpl 1199 DO jj = 2, jpjm1 1200 DO ji = fs_2, fs_jpim1 1201 zti_ups(ji,jj,jl)= pt_ups(ji+1,jj ,jl) 1202 ztj_ups(ji,jj,jl)= pt_ups(ji ,jj+1,jl) 1203 END DO 1204 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 1205 1115 END DO 1206 1116 CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1., ztj_ups, 'T', 1. ) 1207 1117 1208 1118 DO jl = 1, jpl 1209 DO jj = 2, jpjm1 1210 DO ji = fs_2, fs_jpim1 1211 IF ( pfu_ho(ji,jj,jl) * ( pt_ups(ji+1,jj ,jl) - pt_ups(ji,jj,jl) ) <= 0._wp .AND. & 1212 & pfv_ho(ji,jj,jl) * ( pt_ups(ji ,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0._wp ) THEN 1213 ! 1214 IF( pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj ,jl) - zti_ups(ji,jj,jl) ) <= 0._wp .AND. & 1215 & pfv_ho(ji,jj,jl) * ( ztj_ups(ji ,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0._wp ) THEN 1216 pfu_ho(ji,jj,jl)=0._wp 1217 pfv_ho(ji,jj,jl)=0._wp 1218 ENDIF 1219 ! 1220 IF( pfu_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji-1,jj ,jl) ) <= 0._wp .AND. & 1221 & pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji ,jj-1,jl) ) <= 0._wp ) THEN 1222 pfu_ho(ji,jj,jl)=0._wp 1223 pfv_ho(ji,jj,jl)=0._wp 1224 ENDIF 1225 ! 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 1226 1127 ENDIF 1227 END DO 1228 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 1229 1137 END DO 1230 1138 CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. ) ! lateral boundary cond. … … 1238 1146 DO jl = 1, jpl 1239 1147 1240 DO jj = 1, jpj 1241 DO ji = 1, jpi 1242 IF ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN 1243 zbup(ji,jj) = -zbig 1244 zbdo(ji,jj) = zbig 1245 ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) > 0._wp ) THEN 1246 zbup(ji,jj) = pt_ups(ji,jj,jl) 1247 zbdo(ji,jj) = pt_ups(ji,jj,jl) 1248 ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN 1249 zbup(ji,jj) = pt(ji,jj,jl) 1250 zbdo(ji,jj) = pt(ji,jj,jl) 1251 ELSE 1252 zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ups(ji,jj,jl) ) 1253 zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ups(ji,jj,jl) ) 1254 ENDIF 1255 END DO 1256 END DO 1257 1258 DO jj = 2, jpjm1 1259 DO ji = fs_2, fs_jpim1 ! vector opt. 1260 ! 1261 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 1262 zdo = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) ) 1263 ! 1264 zpos = MAX( 0._wp, pfu_ho(ji-1,jj ,jl) ) - MIN( 0._wp, pfu_ho(ji ,jj ,jl) ) & ! positive/negative part of the flux 1265 & + MAX( 0._wp, pfv_ho(ji ,jj-1,jl) ) - MIN( 0._wp, pfv_ho(ji ,jj ,jl) ) 1266 zneg = MAX( 0._wp, pfu_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfu_ho(ji-1,jj ,jl) ) & 1267 & + MAX( 0._wp, pfv_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfv_ho(ji ,jj-1,jl) ) 1268 ! 1269 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) ) & 1270 & ) * ( 1. - pamsk ) 1271 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) ) & 1272 & ) * ( 1. - pamsk ) 1273 ! 1274 ! ! up & down beta terms 1275 ! 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 !!!) 1276 IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 1277 ELSE ; zbetup(ji,jj,jl) = 0._wp ! zbig 1278 ENDIF 1279 ! 1280 IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 1281 ELSE ; zbetdo(ji,jj,jl) = 0._wp ! zbig 1282 ENDIF 1283 ! 1284 ! if all the points are outside ice cover 1285 IF( zup == -zbig ) zbetup(ji,jj,jl) = 0._wp ! zbig 1286 IF( zdo == zbig ) zbetdo(ji,jj,jl) = 0._wp ! zbig 1287 ! 1288 END DO 1289 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 1290 1194 END DO 1291 1195 CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign) … … 1295 1199 ! --------------------------------- 1296 1200 DO jl = 1, jpl 1297 DO jj = 1, jpjm1 1298 DO ji = 1, fs_jpim1 ! vector opt. 1299 zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 1300 zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) 1301 zcu = 0.5_wp + SIGN( 0.5_wp , pfu_ho(ji,jj,jl) ) 1302 ! 1303 zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu ) 1304 ! 1305 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_ups(ji,jj,jl) 1306 ! 1307 END DO 1308 END DO 1309 1310 DO jj = 1, jpjm1 1311 DO ji = 1, fs_jpim1 ! vector opt. 1312 zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 1313 zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) 1314 zcv = 0.5_wp + SIGN( 0.5_wp , pfv_ho(ji,jj,jl) ) 1315 ! 1316 zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv ) 1317 ! 1318 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_ups(ji,jj,jl) 1319 ! 1320 END DO 1321 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 1322 1222 1323 1223 END DO … … 1344 1244 ! 1345 1245 DO jl = 1, jpl 1346 DO jj = 2, jpjm1 1347 DO ji = fs_2, fs_jpim1 ! vector opt. 1348 zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1) 1349 END DO 1350 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 1351 1249 END DO 1352 1250 CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.) ! lateral boundary cond. 1353 1251 1354 1252 DO jl = 1, jpl 1355 DO jj = 2, jpjm1 1356 DO ji = fs_2, fs_jpim1 ! vector opt. 1357 uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj) 1358 1359 Rjm = zslpx(ji-1,jj,jl) 1360 Rj = zslpx(ji ,jj,jl) 1361 Rjp = zslpx(ji+1,jj,jl) 1362 1363 IF( np_limiter == 3 ) THEN 1364 1365 IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm 1366 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 1367 1280 ENDIF 1368 1369 zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 1370 IF( Rj > 0. ) THEN 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 ELSE 1374 zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)), & 1375 & MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 1376 ENDIF 1377 pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter 1378 1379 ELSEIF( np_limiter == 2 ) THEN 1380 IF( Rj /= 0. ) THEN 1381 IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj 1382 ELSE ; Cr = Rjp / Rj 1383 ENDIF 1384 ELSE 1385 Cr = 0. 1386 ENDIF 1387 1388 ! -- superbee -- 1389 zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 1390 ! -- van albada 2 -- 1391 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1392 ! -- sweby (with beta=1) -- 1393 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 1394 ! -- van Leer -- 1395 !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 1396 ! -- ospre -- 1397 !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 1398 ! -- koren -- 1399 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 1400 ! -- charm -- 1401 !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 1402 !ELSE ; zpsi = 0. 1403 !ENDIF 1404 ! -- van albada 1 -- 1405 !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 1406 ! -- smart -- 1407 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 1408 ! -- umist -- 1409 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 1410 1411 ! high order flux corrected by the limiter 1412 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5 1413 1281 ELSE 1282 Cr = 0. 1414 1283 ENDIF 1415 END DO 1416 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 1417 1313 END DO 1418 1314 CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.) ! lateral boundary cond. … … 1439 1335 ! 1440 1336 DO jl = 1, jpl 1441 DO jj = 2, jpjm1 1442 DO ji = fs_2, fs_jpim1 ! vector opt. 1443 zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1) 1444 END DO 1445 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 1446 1340 END DO 1447 1341 CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.) ! lateral boundary cond. 1448 1342 1449 1343 DO jl = 1, jpl 1450 DO jj = 2, jpjm1 1451 DO ji = fs_2, fs_jpim1 ! vector opt. 1452 vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj) 1453 1454 Rjm = zslpy(ji,jj-1,jl) 1455 Rj = zslpy(ji,jj ,jl) 1456 Rjp = zslpy(ji,jj+1,jl) 1457 1458 IF( np_limiter == 3 ) THEN 1459 1460 IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm 1461 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 1462 1372 ENDIF 1463 1464 zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1465 IF( Rj > 0. ) THEN 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 ELSE 1469 zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)), & 1470 & MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 1471 ENDIF 1472 pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter 1473 1474 ELSEIF( np_limiter == 2 ) THEN 1475 1476 IF( Rj /= 0. ) THEN 1477 IF( pv(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj 1478 ELSE ; Cr = Rjp / Rj 1479 ENDIF 1480 ELSE 1481 Cr = 0. 1482 ENDIF 1483 1484 ! -- superbee -- 1485 zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 1486 ! -- van albada 2 -- 1487 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1488 ! -- sweby (with beta=1) -- 1489 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 1490 ! -- van Leer -- 1491 !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 1492 ! -- ospre -- 1493 !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 1494 ! -- koren -- 1495 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 1496 ! -- charm -- 1497 !IF( Cr > 0. ) THEN ; zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 1498 !ELSE ; zpsi = 0. 1499 !ENDIF 1500 ! -- van albada 1 -- 1501 !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 1502 ! -- smart -- 1503 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 1504 ! -- umist -- 1505 !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 1506 1507 ! high order flux corrected by the limiter 1508 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5 1509 1373 ELSE 1374 Cr = 0. 1510 1375 ENDIF 1511 END DO 1512 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 1513 1405 END DO 1514 1406 CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.) ! lateral boundary cond. … … 1544 1436 DO jl = 1, jpl 1545 1437 1546 DO jj = 1, jpj 1547 DO ji = 1, jpi 1548 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 ) 1549 1462 ! 1550 ! ! -- check h_ip -- ! 1551 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 1552 IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 1553 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 1554 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 1555 pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 1556 ENDIF 1557 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 1558 1465 ! 1559 ! ! -- check h_i -- ! 1560 ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 1561 zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 1562 IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1563 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) 1564 ENDIF 1565 ! 1566 ! ! -- check h_s -- ! 1567 ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 1568 zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 1569 IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1570 zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 1571 ! 1572 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 1573 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 1574 ! 1575 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 1576 pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 1577 ENDIF 1578 ! 1579 ENDIF 1580 END DO 1581 END DO 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 1582 1472 END DO 1583 1473 ! … … 1612 1502 ! -- check snow load -- ! 1613 1503 DO jl = 1, jpl 1614 DO jj = 1, jpj 1615 DO ji = 1, jpi 1616 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 1617 ! 1618 zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 1619 ! 1620 IF( zvs_excess > 0._wp ) THEN ! snow-ice interface deplets below the ocean surface 1621 ! put snow excess in the ocean 1622 zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 1623 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 1624 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 1625 ! correct snow volume and heat content 1626 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 1627 pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess 1628 ENDIF 1629 ! 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) * (rau0-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 1630 1517 ENDIF 1631 END DO 1632 END DO 1518 ! 1519 ENDIF 1520 END_2D 1633 1521 END DO 1634 1522 !
Note: See TracChangeset
for help on using the changeset viewer.
