New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 15735 – NEMO

Changeset 15735


Ignore:
Timestamp:
2022-03-03T20:02:46+01:00 (2 years ago)
Author:
emmafiedler
Message:

Added max strength as namelist option, ismooth value check, lbc_lnk for zshear and general tidying up

Location:
NEMO/branches/UKMO/NEMO_4.0.4_ice_strength_EAP/src/ICE
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.4_ice_strength_EAP/src/ICE/icedyn_rdgrft.F90

    r15734 r15735  
    7474   LOGICAL  ::   ln_redist_ufm    ! uniform redistribution of ridged ice (Hibler, 1980)  
    7575   LOGICAL  ::   ln_redist_exp    ! exponential redistribution of ridged ice 
     76   REAL(wp) ::   max_strength     !    maximum strength cap    
    7677   REAL(wp) ::   rn_murdg         !    gives e-folding scale of ridged ice (m^.5) for ln_redist_exp 
    7778   REAL(wp) ::   rn_csrdg         ! fraction of shearing energy contributing to ridging             
     
    513514      ! 2) Transfer function 
    514515      !----------------------------------------------------------------- 
     516      ! If assuming ridged ice is uniformly distributed between hrmin and  
     517      ! hrmax (ln_redist_ufm): 
     518      ! 
    515519      ! Compute max and min ridged ice thickness for each ridging category. 
    516       ! Assume ridged ice is uniformly distributed between hrmin and hrmax. 
    517520      !  
    518521      ! This parameterization is a modified version of Hibler (1980). 
     
    528531      ! These modifications have the effect of reducing the ice strength 
    529532      ! (relative to the Hibler formulation) when very thick ice is ridging. 
     533      ! 
     534      !----------------------------------------------------------------- 
     535      ! If assuming ridged ice ITD is a negative exponential  
     536      ! (ln_redist_exp) and following CICE implementation:  
     537      !  
     538      !  g(h) ~ exp[-(h-hrmin)/hrexp], h >= hrmin  
     539      !  
     540      ! where hrmin is the minimum thickness of ridging ice and  
     541      ! hrexp is the e-folding thickness. 
     542      !  
     543      ! Here, assume as above that hrmin = min(2*hi, hi+maxraft). 
     544      ! That is, the minimum ridge thickness results from rafting, 
     545      !  unless the ice is thicker than maxraft. 
     546      ! 
     547      ! Also, assume that hrexp = mu_rdg*sqrt(hi). 
     548      ! The parameter mu_rdg is tuned to give e-folding scales mostly 
     549      !  in the range 2-4 m as observed by upward-looking sonar. 
     550      ! 
     551      ! Values of mu_rdg in the right column give ice strengths 
     552      !  roughly equal to values of Hstar in the left column 
     553      !  (within ~10 kN/m for typical ITDs): 
     554      ! 
     555      !   Hstar     mu_rdg 
     556      ! 
     557      !     25        3.0 
     558      !     50        4.0 
     559      !     75        5.0 
     560      !    100        6.0 
    530561      ! 
    531562      ! zaksum = net area removed/ total area participating 
     
    792823                 ! Compute the fraction of ridged ice area and volume going to thickness category jl2 
    793824                 ! 
    794                  IF ( ln_redist_ufm ) THEN ! Hibler 1980 uniform formulation 
     825                 IF ( ln_redist_ufm ) THEN ! Hibler (1980) uniform formulation 
    795826                 ! 
    796827                  IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN 
     
    806837                  ENDIF 
    807838                  ! 
    808                  ELSE IF ( ln_redist_exp ) THEN ! exponential formulation 
     839                 ELSE IF ( ln_redist_exp ) THEN ! Lipscomb et al. (2007) exponential formulation 
    809840                  ! 
    810841                  IF (jl2 < jpl) THEN 
     
    929960      IF( ln_str_H79 ) THEN          ! Ice strength => Hibler (1979) method             ! 
    930961      !                              !--------------------------------------------------! 
     962      ! Recommend using ismooth = 1 (spatial smoothing) 
    931963        strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 
    932       ! Recommend ismooth = 1 
    933964      !                              !--------------------------------------------------! 
    934965      ELSE IF( ln_str_R75 ) THEN     ! Ice strength => Rothrock (1975) method           ! 
    935966      !                              !--------------------------------------------------! 
    936       ! 
    937         ! Initialise the strength field at each timestep 
     967      ! Recommend using ismooth = 0 (no smoothing needed) 
     968      ! Initialise the strength field at each timestep 
    938969        strength(:,:) = 0._wp 
    939       ! Recommend ismooth = 0 
    940  
     970      ! 
    941971        !-------------------------------- 
    942972        ! 0) Identify grid cells with ice 
     
    9791009              IF ( apartf(ji,jl) > 0._wp ) THEN 
    9801010                ! 
    981                 IF ( ln_redist_ufm ) THEN   ! Uniform distribution of ridged ice                
     1011                IF ( ln_redist_ufm ) THEN   ! Uniform redistribution of ridged ice                
    9821012                  h2rdg = (1._wp/3._wp) * (hrmax(ji,jl)**3._wp - hrmin(ji,jl)**3._wp)  & 
    9831013                           / (hrmax(ji,jl) - hrmin(ji,jl)) 
    9841014                ! 
    985                 ELSE IF ( ln_redist_exp ) THEN   ! Exponential distribution of ridged ice 
     1015                ELSE IF ( ln_redist_exp ) THEN   ! Exponential redistribution of ridged ice 
    9861016        h2rdg = hrmin(ji,jl) * hrmin(ji,jl) &  
    9871017                + 2._wp * hrmin(ji,jl) * hrexp(ji,jl) &  
     
    10041034 
    10051035          DO ji = 1, npti 
    1006             IF ( zaksum(ji) > epsi10 ) THEN  ! Unclear why is this check should be needed? Occasional x10E-24... 
     1036            IF ( zaksum(ji) > epsi10 ) THEN 
    10071037              strength_1d(ji) = rn_Cf * Cp * strength_pe_sum(ji) / zaksum(ji) 
    10081038            ELSE 
     
    10131043 
    10141044          ! Enforce a maximum for strength 
    1015           WHERE( strength_1d(:) > 2.0e5) ; strength_1d(:) = 2.0e5 
     1045          ! Richter-Menge and Elder (1998) estimate maximum in Beaufort Sea in wintertime of the order 150 kN/m 
     1046          WHERE( strength_1d(:) > max_strength ) ; strength_1d(:) = max_strength 
    10161047          END WHERE 
    10171048 
     
    10351066      CASE( 0 )               !--- No smoothing 
    10361067 
    1037       CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. ) 
     1068         CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. ) 
    10381069 
    10391070      CASE( 1 )               !--- Spatial smoothing 
     
    11741205      INTEGER :: ios                 ! Local integer output status for namelist read 
    11751206      !! 
    1176       NAMELIST/namdyn_rdgrft/ ismooth, ln_str_H79, rn_pstar,          & 
    1177          &                    ln_str_R75, rn_Cf,                      & 
    1178          &                    ln_redist_ufm, ln_redist_exp, rn_murdg, & 
     1207      NAMELIST/namdyn_rdgrft/ ismooth, ln_str_H79, rn_pstar, & 
     1208         &                    ln_str_R75, rn_Cf,             & 
     1209         &                    ln_redist_ufm, ln_redist_exp,  & 
     1210         &                    max_strength, rn_murdg,        & 
    11791211         &                    rn_crhg, rn_csrdg,             & 
    11801212         &                    ln_partf_lin, rn_gstar,        & 
     
    12031235         WRITE(numout,*) '            ratio of ridging work to PE change in ridging      rn_Cf        = ', rn_Cf 
    12041236         WRITE(numout,*) '      smoothing the resistance to deformation (strength)       ismooth      = ', ismooth 
     1237         WRITE(numout,*) '      maximum strength cap                                     max_strength = ', max_strength 
    12051238         WRITE(numout,*) '      uniform redistribution of ridged ice (Hibler, 1980)     ln_redist_ufm = ', ln_redist_ufm 
    12061239         WRITE(numout,*) '      exponential redistribution of ridged ice                ln_redist_exp = ', ln_redist_exp 
     
    12351268      ENDIF 
    12361269      ! 
    1237       !! Add check for ismooth too, =0, 1 or 2 
    1238  
     1270      IF ( ( ismooth .NE. 0 ) .AND. ( ismooth .NE. 1 ) .AND. (ismooth .NE. 2 ) ) THEN 
     1271         CALL ctl_stop( 'ice_dyn_rdgrft_init: ismooth must be one of 0 (no smoothing), 1 (spatial), or 2 (temporal)' ) 
     1272      ENDIF 
     1273      ! 
     1274      ! 
    12391275      IF( .NOT. ln_icethd ) THEN 
    12401276         rn_porordg = 0._wp 
  • NEMO/branches/UKMO/NEMO_4.0.4_ice_strength_EAP/src/ICE/icedyn_rhg_evp.F90

    r15734 r15735  
    792792      END DO 
    793793      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1., zten_i, 'T', 1., & 
    794          &                                  zs1     , 'T', 1., zs2    , 'T', 1., zs12    , 'F', 1. ) 
     794         &                                  zs1     , 'T', 1., zs2    , 'T', 1., zs12    , 'F', 1., zshear, 'T', 1. ) 
    795795       
    796796      ! --- Store the stress tensor for the next time step --- ! 
     
    885885         END DO                
    886886         ! 
    887          CALL iom_put( 'sig1_pnorm' , zsig1_p(:,:) * zmsk00(:,:) )  
    888          CALL iom_put( 'sig2_pnorm' , zsig2_p(:,:) * zmsk00(:,:) )  
     887         IF( iom_use('sig1_pnorm') )  CALL iom_put( 'sig1_pnorm' , zsig1_p(:,:) * zmsk00(:,:) )  
     888         IF( iom_use('sig2_pnorm') )  CALL iom_put( 'sig2_pnorm' , zsig2_p(:,:) * zmsk00(:,:) )  
    889889       
    890890         DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 
     
    10091009 
    10101010      ! time 
    1011       it = ( kt - 1 ) * kitermax + kiter 
     1011      it = ( kt - nit000 ) * kitermax + kiter 
    10121012       
    10131013      ! convergence 
     
    10291029         istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 
    10301030         ! close file 
    1031          IF( kt == nitend - nn_fsbc + 1 )   istatus = NF90_CLOSE(ncvgid) 
     1031         IF( kt == nitend - nn_fsbc + 1 .AND. kiter == kitermax )   istatus = NF90_CLOSE(ncvgid) 
    10321032      ENDIF 
    10331033       
Note: See TracChangeset for help on using the changeset viewer.