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 14998 – NEMO

Changeset 14998


Ignore:
Timestamp:
2021-06-16T10:06:19+02:00 (3 years ago)
Author:
annkeen
Message:

Add EAP sea ice rheology code (from NOCS)

Location:
NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology
Files:
20 added
6 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology/cfgs/SHARED/field_def_nemo-ice.xml

    r13648 r14998  
    8181          <field id="sig1_pnorm"   long_name="P-normalized 1st principal stress component"                                                                       unit=""     /> 
    8282          <field id="sig2_pnorm"   long_name="P-normalized 2nd principal stress component"                                                                       unit=""     /> 
     83          <field id="icedlt"       long_name="delta"                                                   standard_name="delta"                                     unit=""     /> 
    8384          <field id="normstr"      long_name="Average normal stress in sea ice"                        standard_name="average_normal_stress"                     unit="N/m"  /> 
    8485          <field id="sheastr"      long_name="Maximum shear stress in sea ice"                         standard_name="maximum_shear_stress"                      unit="N/m"  /> 
     
    8687          <field id="icediv"       long_name="Divergence of the sea-ice velocity field"                standard_name="divergence_of_sea_ice_velocity"            unit="s-1"  /> 
    8788          <field id="iceshe"       long_name="Maximum shear of sea-ice velocity field"                 standard_name="maximum_shear_of_sea_ice_velocity"         unit="s-1"  /> 
     89          <field id="aniso"        long_name="anisotropy of sea ice floe orientation (0.5 - 1)"        standard_name="anisotropy"                                unit=""     /> 
     90          <field id="yield11"      long_name="yield surface tensor component 11"                       standard_name="yield11"                                   unit="N/m"  /> 
     91          <field id="yield22"      long_name="yield surface tensor component 22"                       standard_name="yield22"                                   unit="N/m"  /> 
     92          <field id="yield12"      long_name="yield surface tensor component 12"                       standard_name="yield12"                                   unit="N/m"  /> 
    8893          <field id="beta_evp"     long_name="Relaxation parameter of ice rheology (beta)"             standard_name="relaxation_parameter_of_ice_rheology"      unit=""  />    
     94  
    8995  
    9096     <!-- surface heat fluxes --> 
     
    167173          <!-- diags --> 
    168174          <field id="hfxdhc"       long_name="Heat content variation in snow and ice (neg = ice cooling)"   unit="W/m2" /> 
     175          <field id="icedrift_mass" long_name="Ice mass drift"   unit="kg/m2/s" /> 
     176          <field id="icedrift_salt" long_name="Ice salt drift"   unit="g/m2/s"  /> 
     177          <field id="icedrift_heat" long_name="Ice heat drift"   unit="W/m2"    /> 
    169178 
    170179     <!-- sbcssm variables --> 
     
    401410     <field field_ref="normstr"          name="normstr" /> 
    402411     <field field_ref="sheastr"          name="sheastr" /> 
     412          <field field_ref="sig1_pnorm"       name="sig1_pnorm"/> 
     413          <field field_ref="sig2_pnorm"       name="sig2_pnorm"/> 
     414     <field field_ref="icedlt"           name="sidelta" /> 
    403415      
    404416     <!-- heat fluxes --> 
  • NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology/src/ICE/ice.F90

    r14075 r14998  
    150150   ! 
    151151   !                                     !!** ice-rheology namelist (namdyn_rhg) ** 
     152   LOGICAL , PUBLIC ::   ln_rhg_EVP       ! EVP rheology switch, used for rdgrft and rheology 
     153   LOGICAL , PUBLIC ::   ln_rhg_EAP       ! EAP rheology switch, used for rdgrft and rheology 
    152154   LOGICAL , PUBLIC ::   ln_aEVP          !: using adaptive EVP (T or F)  
    153155   REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9 
     
    246248   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   divu_i          !: Divergence of the velocity field             [s-1] 
    247249   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   shear_i         !: Shear of the velocity field                  [s-1] 
     250   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   aniso_11, aniso_12   !: structure tensor elements 
     251   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdg_conv 
    248252   ! 
    249253   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   t_bo            !: Sea-Ice bottom temperature [Kelvin]      
     
    436440      ALLOCATE( u_oce    (jpi,jpj) , v_oce    (jpi,jpj) , ht_i_new  (jpi,jpj) , strength(jpi,jpj) ,  & 
    437441         &      stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,                      & 
    438          &      delta_i  (jpi,jpj) , divu_i   (jpi,jpj) , shear_i   (jpi,jpj) , STAT=ierr(ii) ) 
     442         &      delta_i  (jpi,jpj) , divu_i   (jpi,jpj) , shear_i   (jpi,jpj) ,                      & 
     443         &      aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , rdg_conv  (jpi,jpj) , STAT=ierr(ii) ) 
    439444 
    440445      ii = ii + 1 
  • NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology/src/ICE/icedyn_rdgrft.F90

    r14075 r14998  
    2323   USE icevar         ! sea-ice: operations 
    2424   USE icectl         ! sea-ice: control prints 
     25   !USE icedyn_rhg     ! sea-ice: rheology choice 
    2526   ! 
    2627   USE in_out_manager ! I/O manager 
     
    138139      INTEGER , DIMENSION(jpij) ::   iptidx        ! compute ridge/raft or not 
    139140      REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i 
     141      REAL(wp), DIMENSION(jpij) ::   zconv         ! 1D rdg_conv (if EAP rheology) 
    140142      ! 
    141143      INTEGER, PARAMETER ::   jp_itermax = 20     
     
    174176         
    175177         ! just needed here 
    176          CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i ) 
     178         CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i  ) 
     179         CALL tab_2d_1d( npti, nptidx(1:npti), zconv   (1:npti)      , rdg_conv ) 
    177180         ! needed here and in the iteration loop 
    178181         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i) ! zdivu is used as a work array here (no change in divu_i) 
     
    185188            !                                                        - ice area added in new ridges 
    186189            closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 
     190            IF( ln_rhg_EVP )  closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 
     191            IF( ln_rhg_EAP )  closing_net(ji) = zconv(ji) 
    187192            ! 
    188193            IF( zdivu(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) )   ! make sure the closing rate is large enough 
     
    768773      !                              !--------------------------------------------------! 
    769774         strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 
    770          ismooth = 1 
     775!         ismooth = 1    ! original code 
     776         ismooth = 0    ! try for EAP stability 
    771777         !                           !--------------------------------------------------! 
    772778      ELSE                           ! Zero strength                                    ! 
  • NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology/src/ICE/icedyn_rhg.F90

    r14075 r14998  
    1717   USE ice            ! sea-ice: variables 
    1818   USE icedyn_rhg_evp ! sea-ice: EVP rheology 
     19   USE icedyn_rhg_eap ! sea-ice: EAP rheology 
    1920   USE icectl         ! sea-ice: control prints 
    2021   ! 
     
    3334   !                                        ! associated indices: 
    3435   INTEGER, PARAMETER ::   np_rhgEVP = 1   ! EVP rheology 
    35 !! INTEGER, PARAMETER ::   np_rhgEAP = 2   ! EAP rheology 
     36  INTEGER, PARAMETER ::   np_rhgEAP = 2   ! EAP rheology 
    3637 
    3738   ! ** namelist (namrhg) ** 
    38    LOGICAL ::   ln_rhg_EVP       ! EVP rheology 
     39!   LOGICAL ::   ln_rhg_EVP       ! EVP rheology, used by icedyn_rdgrft 
     40!   LOGICAL ::   ln_rhg_EAP       ! EAP rheology, used by icedyn_rdgrft 
     41   !LOGICAL, PUBLIC ::   ln_rhg_EVP       ! EVP rheology, used by icedyn_rdgrft 
     42   !LOGICAL, PUBLIC ::   ln_rhg_EAP       ! EAP rheology, used by icedyn_rdgrft 
    3943   ! 
    4044   !! * Substitutions 
     
    8185         CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i ) 
    8286         !          
     87         !                             !----------------------------! 
     88      CASE( np_rhgEAP )                ! Elasto-Anisotropic-Plastic ! 
     89         !                             !----------------------------! 
     90         CALL ice_dyn_rhg_eap( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i, aniso_11, aniso_12, rdg_conv ) 
    8391      END SELECT 
    8492      ! 
    85       IF( lrst_ice ) THEN                       !* write EVP fields in the restart file 
    86          IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'WRITE', kt ) 
     93      IF( lrst_ice ) THEN                       
     94         IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'WRITE', kt ) !* write EVP fields in the restart file 
     95         IF( ln_rhg_EAP )   CALL rhg_eap_rst( 'WRITE', kt ) !* write EAP fields in the restart file 
    8796      ENDIF 
    8897      ! 
     
    110119      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read 
    111120      !! 
    112       NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, nn_rhg_chkcvg 
     121      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, ln_rhg_EAP, rn_creepl, rn_ecc , nn_nevp, rn_relast, nn_rhg_chkcvg 
    113122      !!------------------------------------------------------------------- 
    114123      ! 
     
    137146         ELSEIF( nn_rhg_chkcvg == 2 ) THEN   ;   WRITE(numout,*) '         check cvg at both main and rheology time steps' 
    138147         ENDIF 
     148         WRITE(numout,*) '      rheology EAP (icedyn_rhg_eap)                        ln_rhg_EAP = ', ln_rhg_EAP 
    139149      ENDIF 
    140150      ! 
     
    142152      ioptio = 0  
    143153      IF( ln_rhg_EVP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEVP    ;   ENDIF 
    144 !!    IF( ln_rhg_EAP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEAP    ;   ENDIF 
     154      IF( ln_rhg_EAP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEAP    ;   ENDIF 
    145155      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_dyn_rhg_init: choose one and only one ice rheology' ) 
    146156      ! 
    147157      IF( ln_rhg_EVP  )   CALL rhg_evp_rst( 'READ' )  !* read or initialize all required files 
     158      IF( ln_rhg_EAP  )   CALL rhg_eap_rst( 'READ' )  !* read or initialize all required files 
    148159      ! 
    149160   END SUBROUTINE ice_dyn_rhg_init 
  • NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology/src/ICE/icedyn_rhg_evp.F90

    r14075 r14998  
    818818      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence 
    819819      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear 
     820      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta 
    820821      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
    821822 
     
    856857      ! Recommendation 1 : we use ice strength, not replacement pressure 
    857858      ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 
    858 !!$      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 
    859 !!$         ! 
    860 !!$         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )          
    861 !!$         !          
    862 !!$         DO jj = 1, jpj 
    863 !!$            DO ji = 1, jpi 
    864 !!$             
    865 !!$               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates  
    866 !!$               !                        and **deformations** at current iterates 
    867 !!$               !                        following Lemieux & Dupont (2020) 
    868 !!$               zfac             =   zp_delt(ji,jj) 
    869 !!$               zsig1            =   zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) 
    870 !!$               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
    871 !!$               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
    872 !!$                
    873 !!$               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 
    874 !!$               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure 
    875 !!$               zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress 
    876 !!$       
    877 !!$               ! Normalized  principal stresses (used to display the ellipse) 
    878 !!$               z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) ) 
    879 !!$               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 
    880 !!$               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 
    881 !!$            END DO 
    882 !!$         END DO                
    883 !!$         ! 
    884 !!$         CALL iom_put( 'sig1_pnorm' , zsig1_p )  
    885 !!$         CALL iom_put( 'sig2_pnorm' , zsig2_p )  
    886 !!$       
    887 !!$         DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 
    888 !!$          
    889 !!$      ENDIF 
     859      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 
     860         ! 
     861         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )          
     862         !          
     863         DO jj = 1, jpj 
     864            DO ji = 1, jpi 
     865             
     866               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates  
     867               !                        and **deformations** at current iterates 
     868               !                        following Lemieux & Dupont (2020) 
     869               zfac             =   zp_delt(ji,jj) 
     870               zsig1            =   zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) 
     871               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
     872               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     873                
     874               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 
     875               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure 
     876               zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress 
     877       
     878               ! Normalized  principal stresses (used to display the ellipse) 
     879               z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) ) 
     880               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 
     881               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 
     882            END DO 
     883         END DO                
     884         ! 
     885         CALL iom_put( 'sig1_pnorm' , zsig1_p )  
     886         CALL iom_put( 'sig2_pnorm' , zsig2_p )  
     887       
     888         DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 
     889          
     890      ENDIF 
    890891 
    891892      ! --- SIMIP --- ! 
  • NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology/src/ICE/icewri.F90

    r14075 r14998  
    254254      CALL iom_rstput( 0, 0, kid, 'sivolu', vt_i         )   ! Ice volume 
    255255      CALL iom_rstput( 0, 0, kid, 'sidive', divu_i*1.0e8 )   ! Ice divergence 
     256      CALL iom_rstput( 0, 0, kid, 'sishea', shear_i*1.0e8 )   ! Ice shear 
    256257      CALL iom_rstput( 0, 0, kid, 'si_amp', at_ip        )   ! Melt pond fraction 
    257258      CALL iom_rstput( 0, 0, kid, 'si_vmp', vt_ip        )   ! Melt pond volume 
Note: See TracChangeset for help on using the changeset viewer.