Changeset 14004 for NEMO/branches/2020/SI3_martin_ponds/src

Timestamp:

2020-12-02T15:12:26+01:00 (4 years ago)

Author:

clem

Message:

update to trunk

Location:

NEMO/branches/2020/SI3_martin_ponds/src/ICE

Files:

• ice.F90 (modified) (3 diffs)
• icedyn_rdgrft.F90 (modified) (4 diffs)
• icedyn_rhg.F90 (modified) (6 diffs)
• icedyn_rhg_eap.F90 (copied) (copied from NEMO/trunk/src/ICE/icedyn_rhg_eap.F90)

NEMO/branches/2020/SI3_martin_ponds/src/ICE/ice.F90

@@ -150 +150 @@
    !
    !                                     !!** ice-rheology namelist (namdyn_rhg) **
+   LOGICAL , PUBLIC ::   ln_rhg_EVP       ! EVP rheology switch, used for rdgrft and rheology
+   LOGICAL , PUBLIC ::   ln_rhg_EAP       ! EAP rheology switch, used for rdgrft and rheology
    LOGICAL , PUBLIC ::   ln_aEVP          !: using adaptive EVP (T or F) 
    REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9
@@ -248 +250 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   divu_i          !: Divergence of the velocity field             [s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   shear_i         !: Shear of the velocity field                  [s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   aniso_11, aniso_12   !: structure tensor elements
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdg_conv
    !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   t_bo            !: Sea-Ice bottom temperature [Kelvin]     
@@ -445 +449 @@
       ALLOCATE( u_oce    (jpi,jpj) , v_oce    (jpi,jpj) , ht_i_new  (jpi,jpj) , strength(jpi,jpj) ,  &
          &      stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,                      &
-         &      delta_i  (jpi,jpj) , divu_i   (jpi,jpj) , shear_i   (jpi,jpj) , STAT=ierr(ii) )
+         &      delta_i  (jpi,jpj) , divu_i   (jpi,jpj) , shear_i   (jpi,jpj) ,                      &
+         &      aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , rdg_conv  (jpi,jpj) , STAT=ierr(ii) )
 
       ii = ii + 1

NEMO/branches/2020/SI3_martin_ponds/src/ICE/icedyn_rdgrft.F90

@@ -140 +140 @@
       INTEGER , DIMENSION(jpij) ::   iptidx        ! compute ridge/raft or not
       REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i
+      REAL(wp), DIMENSION(jpij) ::   zconv         ! 1D rdg_conv (if EAP rheology)
       !
       INTEGER, PARAMETER ::   jp_itermax = 20    
@@ -175 +176 @@
          ! just needed here
          CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i )
+         CALL tab_2d_1d( npti, nptidx(1:npti), zconv   (1:npti)      , rdg_conv )
          ! needed here and in the iteration loop
          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)
@@ -184 +186 @@
             ! closing_net = rate at which open water area is removed + ice area removed by ridging 
             !                                                        - ice area added in new ridges
-            closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp )
+            IF( ln_rhg_EVP )  closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp )
+            IF( ln_rhg_EAP )  closing_net(ji) = zconv(ji)
             !
             IF( zdivu(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) )   ! make sure the closing rate is large enough
@@ -776 +779 @@
       !                              !--------------------------------------------------!
          strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
-         ismooth = 1
+         ismooth = 1    ! original code
+!        ismooth = 0    ! try for EAP stability
          !                           !--------------------------------------------------!
       ELSE                           ! Zero strength                                    !

NEMO/branches/2020/SI3_martin_ponds/src/ICE/icedyn_rhg.F90

@@ -17 +17 @@
    USE ice            ! sea-ice: variables
    USE icedyn_rhg_evp ! sea-ice: EVP rheology
+   USE icedyn_rhg_eap ! sea-ice: EAP rheology
    USE icectl         ! sea-ice: control prints
    !
@@ -33 +34 @@
    !                                        ! associated indices:
    INTEGER, PARAMETER ::   np_rhgEVP = 1   ! EVP rheology
-!! INTEGER, PARAMETER ::   np_rhgEAP = 2   ! EAP rheology
+   INTEGER, PARAMETER ::   np_rhgEAP = 2   ! EAP rheology
 
    ! ** namelist (namrhg) **
-   LOGICAL ::   ln_rhg_EVP       ! EVP rheology
    !
    !!----------------------------------------------------------------------
@@ -78 +78 @@
          CALL ice_dyn_rhg_evp( kt, Kmm, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i )
          !         
+         !                             !----------------------------!
+      CASE( np_rhgEAP )                ! Elasto-Anisotropic-Plastic !
+         !                             !----------------------------!
+         CALL ice_dyn_rhg_eap( kt, Kmm, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i, aniso_11, aniso_12, rdg_conv )
       END SELECT
       !
       IF( lrst_ice ) THEN                       !* write EVP fields in the restart file
          IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'WRITE', kt )
+         IF( ln_rhg_EAP )   CALL rhg_eap_rst( 'WRITE', kt ) !* write EAP fields in the restart file
       ENDIF
       !
@@ -108 +113 @@
       INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, nn_rhg_chkcvg
+      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, ln_rhg_EAP, rn_creepl, rn_ecc , nn_nevp, rn_relast, nn_rhg_chkcvg
       !!-------------------------------------------------------------------
       !
@@ -133 +138 @@
          ELSEIF( nn_rhg_chkcvg == 2 ) THEN   ;   WRITE(numout,*) '         check cvg at both main and rheology time steps'
          ENDIF
+         WRITE(numout,*) '      rheology EAP (icedyn_rhg_eap)                        ln_rhg_EAP = ', ln_rhg_EAP
       ENDIF
       !
@@ -138 +144 @@
       ioptio = 0 
       IF( ln_rhg_EVP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEVP    ;   ENDIF
-!!    IF( ln_rhg_EAP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEAP    ;   ENDIF
+      IF( ln_rhg_EAP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEAP    ;   ENDIF
       IF( ioptio /= 1 )   CALL ctl_stop( 'ice_dyn_rhg_init: choose one and only one ice rheology' )
       !
       IF( ln_rhg_EVP  )   CALL rhg_evp_rst( 'READ' )  !* read or initialize all required files
+      IF( ln_rhg_EAP  )   CALL rhg_eap_rst( 'READ' )  !* read or initialize all required files
       !
    END SUBROUTINE ice_dyn_rhg_init