Changeset 10930
- Timestamp:
- 2019-05-05T21:37:21+02:00 (4 years ago)
- Location:
- NEMO/trunk/src/ICE
- Files:
-
- 2 edited
-
icedyn_adv_umx.F90 (modified) (19 diffs)
-
icevar.F90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icedyn_adv_umx.F90
r10911 r10930 36 36 REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120 37 37 38 ! limiter: 1=nonosc_ice, 2=superbee, 3=h3(rachid) 39 INTEGER :: kn_limiter = 1 40 41 ! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge), 42 ! interpolated T at u/v points can be non-zero while it should 43 ! (because of the high order of the advection scheme). Thus set it to 0 in this case 44 LOGICAL :: ll_icedge = .TRUE. 45 38 ! advection for S and T: dVS/dt = -div( uA * uHS / u ) => ll_advS = F 39 ! or dVS/dt = -div( uV * uS / u ) => ll_advS = T 40 LOGICAL :: ll_advS = .FALSE. 41 ! 42 ! alternate directions for upstream 43 LOGICAL :: ll_upsxy = .TRUE. 44 ! 45 ! alternate directions for high order 46 LOGICAL :: ll_hoxy = .TRUE. 47 ! 46 48 ! if T interpolated at u/v points is negative or v_i < 1.e-6 47 49 ! then interpolate T at u/v points using the upstream scheme 48 50 LOGICAL :: ll_neg = .TRUE. 49 50 ! alternate directions for upstream 51 LOGICAL :: ll_upsxy = .TRUE. 52 53 ! alternate directions for high order 54 LOGICAL :: ll_hoxy = .TRUE. 55 51 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zmski_small, zmskj_small 52 ! 53 ! limiter: 1=nonosc_ice, 2=superbee, 3=h3(rachid) 54 INTEGER :: kn_limiter = 1 55 ! 56 56 ! prelimiter: use it to avoid overshoot in H 57 57 LOGICAL :: ll_prelimiter_zalesak = .FALSE. ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D 58 59 ! advection for S and T: dVS/dt = -div( uA * uHS / u ) => ll_advS = F 60 ! or dVS/dt = -div( uV * uS / u ) => ll_advS = T 61 LOGICAL :: ll_advS = .FALSE. 62 58 ! 63 59 !! * Substitutions 64 60 # include "vectopt_loop_substitute.h90" … … 102 98 INTEGER :: icycle ! number of sub-timestep for the advection 103 99 REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers 104 REAL(wp) :: zdt 100 REAL(wp) :: zdt, zvi_cen 105 101 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 106 102 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box … … 134 130 END DO 135 131 END DO 136 CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max (:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. )132 CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. ) 137 133 ! 138 134 ! … … 185 181 ELSEWHERE ; z1_aip(:,:,:) = 0. 186 182 END WHERE 183 ! 184 ! setup a mask where advection will be upstream 185 IF( ll_neg ) THEN 186 IF( .NOT. ALLOCATED(zmski_small) ) ALLOCATE( zmski_small(jpi,jpj,jpl) ) 187 IF( .NOT. ALLOCATED(zmskj_small) ) ALLOCATE( zmskj_small(jpi,jpj,jpl) ) 188 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 ; zmski_small(ji,jj,jl) = 0. 193 ELSE ; zmski_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 ; zmskj_small(ji,jj,jl) = 0. 196 ELSE ; zmskj_small(ji,jj,jl) = 1. ; ENDIF 197 END DO 198 END DO 199 END DO 200 ENDIF 201 ! 202 ! ----------------------- ! 203 ! ==> start advection <== ! 204 ! ----------------------- ! 187 205 ! 188 206 ! set u*A=u for advection of A only … … 319 337 END DO 320 338 END DO 321 CALL lbc_lnk( 'icedyn_adv_umx', pato_i (:,:), 'T', 1. )339 CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1. ) 322 340 ! 323 341 ! … … 326 344 ! Remove negative values (conservation is ensured) 327 345 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 328 CALL ice_var_zapneg( pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )346 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 ) 329 347 ! 330 348 ! Make sure ice thickness is not too big 331 349 ! (because ice thickness can be too large where ice concentration is very small) 332 CALL Hbig( z hi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )350 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 ) 333 351 334 352 END DO … … 762 780 & * pamsk & 763 781 & ) * pdt ) * tmask(ji,jj,1) 764 !!clem test765 !!zpt(ji,jj,jl) = MAX( 0._wp, zpt(ji,jj,jl) )766 !!clem test767 782 END DO 768 783 END DO … … 794 809 & * pamsk & 795 810 & ) * pdt ) * tmask(ji,jj,1) 796 !!clem test797 !!zpt(ji,jj,jl) = MAX( 0._wp, zpt(ji,jj,jl) )798 !!clem test799 811 END DO 800 812 END DO … … 837 849 ! 838 850 INTEGER :: ji, jj, jl ! dummy loop indices 839 REAL(wp) :: zcu, zdx2, zdx4 , zvi_cen2! - -851 REAL(wp) :: zcu, zdx2, zdx4 ! - - 840 852 REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztu1, ztu2, ztu3, ztu4 841 853 !!---------------------------------------------------------------------- … … 948 960 END SELECT 949 961 ! 950 ! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge),951 ! interpolated T at u/v points can be non-zero while it should952 ! (because of the high order of the advection scheme). Thus set it to 0 in this case953 IF( ll_icedge ) THEN954 DO jl = 1, jpl955 DO jj = 1, jpjm1956 DO ji = 1, fs_jpim1957 IF( pt(ji,jj,jl) <= 0._wp .AND. pu(ji,jj) >= 0._wp ) THEN958 pt_u(ji,jj,jl) = 0._wp959 ENDIF960 END DO961 END DO962 END DO963 ENDIF964 !965 962 ! if pt at u-point is negative then use the upstream value 966 963 ! this should not be necessary if a proper sea-ice mask is set in Ultimate … … 970 967 DO jj = 1, jpjm1 971 968 DO ji = 1, fs_jpim1 972 zvi_cen2 = 0.5_wp * ( v_i(ji+1,jj,jl) + v_i(ji,jj,jl) ) 973 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( zvi_cen2 < epsi06 .AND. pamsk == 0._wp ) ) THEN 969 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( zmski_small(ji,jj,jl) == 0. .AND. pamsk == 0. ) ) THEN 974 970 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 975 971 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) … … 1010 1006 ! 1011 1007 INTEGER :: ji, jj, jl ! dummy loop indices 1012 REAL(wp) :: zcv, zdy2, zdy4 , zvi_cen2! - -1008 REAL(wp) :: zcv, zdy2, zdy4 ! - - 1013 1009 REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztv1, ztv2, ztv3, ztv4 1014 1010 !!---------------------------------------------------------------------- … … 1118 1114 END SELECT 1119 1115 ! 1120 ! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge),1121 ! interpolated T at u/v points can be non-zero while it should1122 ! (because of the high order of the advection scheme). Thus set it to 0 in this case1123 IF( ll_icedge ) THEN1124 DO jl = 1, jpl1125 DO jj = 1, jpjm11126 DO ji = 1, fs_jpim11127 IF( pt(ji,jj,jl) <= 0._wp .AND. pv(ji,jj) >= 0._wp ) THEN1128 pt_v(ji,jj,jl) = 0._wp1129 ENDIF1130 END DO1131 END DO1132 END DO1133 ENDIF1134 !1135 1116 ! if pt at v-point is negative then use the upstream value 1136 1117 ! this should not be necessary if a proper sea-ice mask is set in Ultimate … … 1140 1121 DO jj = 1, jpjm1 1141 1122 DO ji = 1, fs_jpim1 1142 zvi_cen2 = 0.5_wp * ( v_i(ji,jj+1,jl) + v_i(ji,jj,jl) ) 1143 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( zvi_cen2 < epsi06 .AND. pamsk == 0._wp ) ) THEN 1123 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( zmskj_small(ji,jj,jl) == 0. .AND. pamsk == 0. ) ) THEN 1144 1124 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & 1145 1125 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) … … 1528 1508 1529 1509 1530 SUBROUTINE Hbig( p hi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )1510 SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 1531 1511 !!------------------------------------------------------------------- 1532 1512 !! *** ROUTINE Hbig *** … … 1545 1525 !! ** input : Max thickness of the surrounding 9-points 1546 1526 !!------------------------------------------------------------------- 1527 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1547 1528 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts 1548 1529 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip … … 1551 1532 ! 1552 1533 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1553 REAL(wp) :: z hip, zhi, zhs, zvs_excess, zfra1534 REAL(wp) :: z1_dt, zhip, zhi, zhs, zvs_excess, zfra 1554 1535 REAL(wp), DIMENSION(jpi,jpj) :: zswitch 1555 1536 !!------------------------------------------------------------------- 1556 1537 ! 1538 z1_dt = 1._wp / pdt 1557 1539 ! 1558 1540 DO jl = 1, jpl … … 1584 1566 zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 1585 1567 ! 1586 wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice1587 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice! W.m-2 <01568 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 1569 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 1588 1570 ! 1589 1571 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra … … 1598 1580 IF( zvs_excess > 0._wp ) THEN 1599 1581 zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 1600 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice1601 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice! W.m-2 <01582 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 1583 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 1602 1584 ! 1603 1585 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra -
NEMO/trunk/src/ICE/icevar.F90
r10589 r10930 229 229 IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area 230 230 ! 231 ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i 231 ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3] 232 232 ztmelts = - sz_i(ji,jj,jk,jl) * rTmlt ! Ice layer melt temperature [C] 233 233 ! Conversion q(S,T) -> T (second order equation) … … 236 236 t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts 237 237 ! 238 ELSE !--- no ice238 ELSE !--- no ice 239 239 t_i(ji,jj,jk,jl) = rt0 240 240 ENDIF … … 537 537 538 538 539 SUBROUTINE ice_var_zapneg( p ato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )539 SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 540 540 !!------------------------------------------------------------------- 541 541 !! *** ROUTINE ice_var_zapneg *** … … 543 543 !! ** Purpose : Remove negative sea ice fields and correct fluxes 544 544 !!------------------------------------------------------------------- 545 INTEGER :: ji, jj, jl, jk ! dummy loop indices 546 ! 545 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 547 546 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area 548 547 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume … … 555 554 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content 556 555 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content 557 !!------------------------------------------------------------------- 558 ! 556 ! 557 INTEGER :: ji, jj, jl, jk ! dummy loop indices 558 REAL(wp) :: z1_dt 559 !!------------------------------------------------------------------- 560 ! 561 z1_dt = 1._wp / pdt 559 562 ! 560 563 DO jl = 1, jpl !== loop over the categories ==! … … 567 570 DO ji = 1 , jpi 568 571 IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 569 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * r1_rdtice! W.m-2 >0572 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0 570 573 pe_i(ji,jj,jk,jl) = 0._wp 571 574 ENDIF … … 578 581 DO ji = 1 , jpi 579 582 IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 580 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * r1_rdtice! W.m-2 <0583 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0 581 584 pe_s(ji,jj,jk,jl) = 0._wp 582 585 ENDIF … … 591 594 DO ji = 1 , jpi 592 595 IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 593 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * r1_rdtice596 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt 594 597 pv_i (ji,jj,jl) = 0._wp 595 598 ENDIF 596 599 IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 597 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * r1_rdtice600 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt 598 601 pv_s (ji,jj,jl) = 0._wp 599 602 ENDIF 600 603 IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 601 sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * r1_rdtice604 sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt 602 605 psv_i (ji,jj,jl) = 0._wp 603 606 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
