Index: NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/ICE/icedyn_rhg.F90
===================================================================
--- NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/ICE/icedyn_rhg.F90	(revision 13824)
+++ NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/ICE/icedyn_rhg.F90	(revision 14009)
@@ -17,4 +17,6 @@
    USE ice            ! sea-ice: variables
    USE icedyn_rhg_evp ! sea-ice: EVP rheology
+   USE icedyn_rhg_eap ! sea-ice: EAP rheology
+   USE icedyn_rhg_vp  ! sea-ice: VP  rheology
    USE icectl         ! sea-ice: control prints
    !
@@ -33,8 +35,7 @@
    !                                        ! associated indices:
    INTEGER, PARAMETER ::   np_rhgEVP = 1   ! EVP rheology
-!! INTEGER, PARAMETER ::   np_rhgEAP = 2   ! EAP rheology
+   INTEGER, PARAMETER ::   np_rhgEAP = 2   ! EAP rheology
+   INTEGER, PARAMETER ::   np_rhgVP  = 3   ! VP rheology
 
-   ! ** namelist (namrhg) **
-   LOGICAL ::   ln_rhg_EVP       ! EVP rheology
    !
    !!----------------------------------------------------------------------
@@ -77,9 +78,20 @@
          !                             !------------------------!
          CALL ice_dyn_rhg_evp( kt, Kmm, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i )
-         !         
+         !        
+         !                             !------------------------!
+      CASE( np_rhgVP  )                ! Viscous-Plastic        !
+         !                             !------------------------!
+         CALL ice_dyn_rhg_vp ( kt, 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( lrst_ice ) THEN
+         IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'WRITE', kt ) !* write EVP fields in the restart file
+         IF( ln_rhg_EAP )   CALL rhg_eap_rst( 'WRITE', kt ) !* write EAP fields in the restart file
+         ! MV note: no restart needed for VP as there is no time equation for stress tensor
       ENDIF
       !
@@ -108,5 +120,6 @@
       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, &  !-- evp
+         &                  ln_rhg_VP, nn_vp_nout, nn_vp_ninn, nn_vp_chkcvg                                            !-- vp 
       !!-------------------------------------------------------------------
       !
@@ -124,13 +137,20 @@
          WRITE(numout,*) '      rheology EVP (icedyn_rhg_evp)                        ln_rhg_EVP    = ', ln_rhg_EVP
          WRITE(numout,*) '         use adaptive EVP (aEVP)                           ln_aEVP       = ', ln_aEVP
-         WRITE(numout,*) '         creep limit                                       rn_creepl     = ', rn_creepl
-         WRITE(numout,*) '         eccentricity of the elliptical yield curve        rn_ecc        = ', rn_ecc
+         WRITE(numout,*) '         creep limit                                       rn_creepl     = ', rn_creepl ! also used by vp
+         WRITE(numout,*) '         eccentricity of the elliptical yield curve        rn_ecc        = ', rn_ecc    ! also used by vp
          WRITE(numout,*) '         number of iterations for subcycling               nn_nevp       = ', nn_nevp
          WRITE(numout,*) '         ratio of elastic timescale over ice time step     rn_relast     = ', rn_relast
-         WRITE(numout,*) '      check convergence of rheology                        nn_rhg_chkcvg = ', nn_rhg_chkcvg
-         IF    ( nn_rhg_chkcvg == 0 ) THEN   ;   WRITE(numout,*) '         no check'
-         ELSEIF( nn_rhg_chkcvg == 1 ) THEN   ;   WRITE(numout,*) '         check cvg at the main time step'
-         ELSEIF( nn_rhg_chkcvg == 2 ) THEN   ;   WRITE(numout,*) '         check cvg at both main and rheology time steps'
+         WRITE(numout,*) '         check convergence of rheology                     nn_rhg_chkcvg = ', nn_rhg_chkcvg
+         WRITE(numout,*) '      rheology VP   (icedyn_rhg_VP)                        ln_rhg_VP     = ', ln_rhg_VP
+         WRITE(numout,*) '         number of outer iterations                        nn_vp_nout    = ', nn_vp_nout
+         WRITE(numout,*) '         number of inner iterations                        nn_vp_ninn    = ', nn_vp_ninn
+         WRITE(numout,*) '         iteration step for convergence check              nn_vp_chkcvg  = ', nn_vp_chkcvg
+         IF( ln_rhg_EVP ) THEN
+            IF    ( nn_rhg_chkcvg == 0 ) THEN   ;   WRITE(numout,*) '         no check cvg'
+            ELSEIF( nn_rhg_chkcvg == 1 ) THEN   ;   WRITE(numout,*) '         check cvg at the main time step'
+            ELSEIF( nn_rhg_chkcvg == 2 ) THEN   ;   WRITE(numout,*) '         check cvg at both main and rheology time steps'
+            ENDIF
          ENDIF
+         WRITE(numout,*) '      rheology EAP (icedyn_rhg_eap)                        ln_rhg_EAP = ', ln_rhg_EAP
       ENDIF
       !
@@ -138,8 +158,11 @@
       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( ln_rhg_VP  ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgVP     ;   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
+      ! no restart for VP as there is no explicit time dependency in the equation
       !
    END SUBROUTINE ice_dyn_rhg_init
