Changeset 10930


Ignore:
Timestamp:
2019-05-05T21:37:21+02:00 (18 months ago)
Author:
clem
Message:

solve a problem of conservation when ice advection is called twice (because of a cfl that is too large)

Location:
NEMO/trunk/src/ICE
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/icedyn_adv_umx.F90

    r10911 r10930  
    3636   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120 
    3737    
    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   ! 
    4648   ! if T interpolated at u/v points is negative or v_i < 1.e-6 
    4749   !    then interpolate T at u/v points using the upstream scheme 
    4850   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   ! 
    5656   ! prelimiter: use it to avoid overshoot in H 
    5757   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   ! 
    6359   !! * Substitutions 
    6460#  include "vectopt_loop_substitute.h90" 
     
    10298      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
    10399      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers 
    104       REAL(wp) ::   zdt 
     100      REAL(wp) ::   zdt, zvi_cen 
    105101      REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! for global communication 
    106102      REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box  
     
    134130         END DO 
    135131      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. ) 
    137133      ! 
    138134      ! 
     
    185181         ELSEWHERE                        ;   z1_aip(:,:,:) = 0. 
    186182         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         ! ----------------------- ! 
    187205         ! 
    188206         ! set u*A=u for advection of A only  
     
    319337            END DO 
    320338         END DO 
    321          CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T',  1. ) 
     339         CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1. ) 
    322340         ! 
    323341         ! 
     
    326344         ! Remove negative values (conservation is ensured) 
    327345         !    (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 ) 
    329347         ! 
    330348         ! Make sure ice thickness is not too big 
    331349         !    (because ice thickness can be too large where ice concentration is very small) 
    332          CALL Hbig( zhi_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 ) 
    333351 
    334352      END DO 
     
    762780                     &                                                                                        * pamsk           & 
    763781                     &                             ) * pdt ) * tmask(ji,jj,1) 
    764                   !!clem test 
    765                   !!zpt(ji,jj,jl) = MAX( 0._wp, zpt(ji,jj,jl) ) 
    766                   !!clem test 
    767782               END DO 
    768783            END DO 
     
    794809                     &                                                                                        * pamsk           & 
    795810                     &                             ) * pdt ) * tmask(ji,jj,1)  
    796                   !!clem test 
    797                   !!zpt(ji,jj,jl) = MAX( 0._wp, zpt(ji,jj,jl) ) 
    798                   !!clem test 
    799811               END DO 
    800812            END DO 
     
    837849      ! 
    838850      INTEGER  ::   ji, jj, jl             ! dummy loop indices 
    839       REAL(wp) ::   zcu, zdx2, zdx4, zvi_cen2        !   -      - 
     851      REAL(wp) ::   zcu, zdx2, zdx4        !   -      - 
    840852      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztu1, ztu2, ztu3, ztu4 
    841853      !!---------------------------------------------------------------------- 
     
    948960      END SELECT 
    949961      ! 
    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 should 
    952       !    (because of the high order of the advection scheme). Thus set it to 0 in this case 
    953       IF( ll_icedge ) THEN 
    954          DO jl = 1, jpl 
    955             DO jj = 1, jpjm1 
    956                DO ji = 1, fs_jpim1 
    957                   IF( pt(ji,jj,jl) <= 0._wp .AND. pu(ji,jj) >= 0._wp ) THEN 
    958                      pt_u(ji,jj,jl) = 0._wp 
    959                   ENDIF 
    960                END DO 
    961             END DO 
    962          END DO 
    963       ENDIF 
    964       ! 
    965962      ! if pt at u-point is negative then use the upstream value 
    966963      !    this should not be necessary if a proper sea-ice mask is set in Ultimate 
     
    970967            DO jj = 1, jpjm1 
    971968               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 
    974970                     pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    975971                        &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
     
    10101006      ! 
    10111007      INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    1012       REAL(wp) ::   zcv, zdy2, zdy4, zvi_cen2    !   -      - 
     1008      REAL(wp) ::   zcv, zdy2, zdy4    !   -      - 
    10131009      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztv1, ztv2, ztv3, ztv4 
    10141010      !!---------------------------------------------------------------------- 
     
    11181114      END SELECT 
    11191115      ! 
    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 should 
    1122       !    (because of the high order of the advection scheme). Thus set it to 0 in this case 
    1123       IF( ll_icedge ) THEN 
    1124          DO jl = 1, jpl 
    1125             DO jj = 1, jpjm1 
    1126                DO ji = 1, fs_jpim1 
    1127                   IF( pt(ji,jj,jl) <= 0._wp .AND. pv(ji,jj) >= 0._wp ) THEN 
    1128                      pt_v(ji,jj,jl) = 0._wp 
    1129                   ENDIF 
    1130                END DO 
    1131             END DO 
    1132          END DO 
    1133       ENDIF 
    1134       ! 
    11351116      ! if pt at v-point is negative then use the upstream value 
    11361117      !    this should not be necessary if a proper sea-ice mask is set in Ultimate 
     
    11401121            DO jj = 1, jpjm1 
    11411122               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 
    11441124                     pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
    11451125                        &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
     
    15281508 
    15291509 
    1530    SUBROUTINE Hbig( phi_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 ) 
    15311511      !!------------------------------------------------------------------- 
    15321512      !!                  ***  ROUTINE Hbig  *** 
     
    15451525      !! ** input   : Max thickness of the surrounding 9-points 
    15461526      !!------------------------------------------------------------------- 
     1527      REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step 
    15471528      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
    15481529      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip 
     
    15511532      ! 
    15521533      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices 
    1553       REAL(wp) ::   zhip, zhi, zhs, zvs_excess, zfra 
     1534      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zvs_excess, zfra 
    15541535      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch 
    15551536      !!------------------------------------------------------------------- 
    15561537      ! 
     1538      z1_dt = 1._wp / pdt 
    15571539      ! 
    15581540      DO jl = 1, jpl 
     
    15841566                     zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 
    15851567                     ! 
    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_rdtice 
    1587                      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 <0 
     1568                     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 
    15881570                     ! 
    15891571                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
     
    15981580                  IF( zvs_excess > 0._wp ) THEN 
    15991581                     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_rdtice 
    1601                      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 <0 
     1582                     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 
    16021584                     ! 
    16031585                     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  
    229229                  IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area  
    230230                     ! 
    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] 
     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] 
    232232                     ztmelts          = - sz_i(ji,jj,jk,jl) * rTmlt                                 ! Ice layer melt temperature [C] 
    233233                     ! Conversion q(S,T) -> T (second order equation) 
     
    236236                     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 
    237237                     ! 
    238                   ELSE                                !--- no ice 
     238                  ELSE                                   !--- no ice 
    239239                     t_i(ji,jj,jk,jl) = rt0 
    240240                  ENDIF 
     
    537537 
    538538 
    539    SUBROUTINE ice_var_zapneg( pato_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 ) 
    540540      !!------------------------------------------------------------------- 
    541541      !!                   ***  ROUTINE ice_var_zapneg *** 
     
    543543      !! ** Purpose :   Remove negative sea ice fields and correct fluxes 
    544544      !!------------------------------------------------------------------- 
    545       INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices 
    546       ! 
     545      REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step 
    547546      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area 
    548547      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume 
     
    555554      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    556555      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 
    559562      ! 
    560563      DO jl = 1, jpl       !==  loop over the categories  ==! 
     
    567570               DO ji = 1 , jpi 
    568571                  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 >0 
     572                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0 
    570573                     pe_i(ji,jj,jk,jl) = 0._wp 
    571574                  ENDIF 
     
    578581               DO ji = 1 , jpi 
    579582                  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 <0 
     583                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0 
    581584                     pe_s(ji,jj,jk,jl) = 0._wp 
    582585                  ENDIF 
     
    591594            DO ji = 1 , jpi 
    592595               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_rdtice 
     596                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt 
    594597                  pv_i   (ji,jj,jl) = 0._wp 
    595598               ENDIF 
    596599               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_rdtice 
     600                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt 
    598601                  pv_s   (ji,jj,jl) = 0._wp 
    599602               ENDIF 
    600603               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_rdtice 
     604                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt 
    602605                  psv_i  (ji,jj,jl) = 0._wp 
    603606               ENDIF 
Note: See TracChangeset for help on using the changeset viewer.