Changeset 10929 for NEMO/releases/release-4.0/src/ICE/icedyn_adv_umx.F90
- Timestamp:
- 2019-05-05T21:35:24+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/releases/release-4.0/src/ICE/icedyn_adv_umx.F90
r10910 r10929 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
Note: See TracChangeset
for help on using the changeset viewer.