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.
icedyn_rhg.F90 in NEMO/branches/2020/SI3-03_VP_rheology/src/ICE – NEMO

source: NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg.F90 @ 13658

Last change on this file since 13658 was 13544, checked in by vancop, 4 years ago

now VP rheology compiles
  • Property svn:keywords set to Id
File size: 9.8 KB
Line 
1MODULE icedyn_rhg
2   !!======================================================================
3   !!                     ***  MODULE  icedyn_rhg  ***
4   !!   Sea-Ice dynamics : master routine for rheology
5   !!======================================================================
6   !! history :  4.0  !  2018     (C. Rousset)      Original code
7   !!----------------------------------------------------------------------
8#if defined key_si3
9   !!----------------------------------------------------------------------
10   !!   'key_si3'                                       SI3 sea-ice model
11   !!----------------------------------------------------------------------
12   !!    ice_dyn_rhg      : computes ice velocities
13   !!    ice_dyn_rhg_init : initialization and namelist read
14   !!----------------------------------------------------------------------
15   USE phycst         ! physical constants
16   USE dom_oce        ! ocean space and time domain
17   USE ice            ! sea-ice: variables
18   USE icedyn_rhg_evp ! sea-ice: EVP rheology
19   USE icedyn_rhg_vp  ! sea-ice: VP  rheology
20   USE icectl         ! sea-ice: control prints
21   !
22   USE in_out_manager ! I/O manager
23   USE lib_mpp        ! MPP library
24   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
25   USE timing         ! Timing
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   ice_dyn_rhg        ! called by icestp.F90
31   PUBLIC   ice_dyn_rhg_init   ! called by icestp.F90
32
33   INTEGER ::              nice_rhg   ! choice of the type of rheology
34   !                                        ! associated indices:
35   INTEGER, PARAMETER ::   np_rhgEVP = 1   ! EVP rheology
36!! INTEGER, PARAMETER ::   np_rhgEAP = 2   ! EAP rheology
37   INTEGER, PARAMETER ::   np_rhgVP  = 3   ! VP rheology
38
39   ! ** namelist (namrhg) **
40   LOGICAL ::   ln_rhg_EVP       ! EVP rheology ! MV: why is it declared here while all other parameters are declared in ice.F90
41   LOGICAL ::   ln_rhg_VP        ! EVP rheology
42   !
43   !! * Substitutions
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (./LICENSE)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE ice_dyn_rhg( kt )
53      !!-------------------------------------------------------------------
54      !!               ***  ROUTINE ice_dyn_rhg  ***
55      !!               
56      !! ** Purpose :   compute ice velocity
57      !!
58      !! ** Action  : comupte - ice velocity (u_ice, v_ice)
59      !!                      - 3 components of the stress tensor (stress1_i, stress2_i, stress12_i) - EVP only
60      !!                      - shear, divergence and delta (shear_i, divu_i, delta_i)
61      !!--------------------------------------------------------------------
62      INTEGER, INTENT(in) ::   kt     ! ice time step
63      !
64      INTEGER  ::   jl   ! dummy loop indices
65      !!--------------------------------------------------------------------
66      ! controls
67      IF( ln_timing    )   CALL timing_start('icedyn_rhg')                                                             ! timing
68      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_rhg', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
69      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icedyn_rhg',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
70      !
71      IF( kt == nit000 .AND. lwp ) THEN
72         WRITE(numout,*)
73         WRITE(numout,*)'ice_dyn_rhg: sea-ice rheology'
74         WRITE(numout,*)'~~~~~~~~~~~'
75      ENDIF
76      !
77      !--------------!
78      !== Rheology ==!
79      !--------------!   
80      SELECT CASE( nice_rhg )
81      !                                !---------------------------------!
82      CASE( np_rhgEVP )                ! Elasto-Viscous-Plastic rheology !
83         !                             !---------------------------------!
84         CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i )
85      !                                !---------------------------------------------!
86      CASE( np_rhgVP  )                ! Viscous-Plastic rheology                    !
87         !                             !---------------------------------------------!
88         CALL ice_dyn_rhg_vp ( kt, shear_i, divu_i, delta_i )
89         !         
90      END SELECT
91      !
92      IF( lrst_ice ) THEN                       !* write EVP fields in the restart file
93         IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'WRITE', kt )
94         ! MV note: no restart needed for VP as there is no time equation for stress tensor
95      ENDIF
96      !
97      ! controls
98      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rhg')                                                             ! prints
99      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rhg', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
100      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icedyn_rhg',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
101      IF( ln_timing    )   CALL timing_stop ('icedyn_rhg')                                                             ! timing
102      !
103   END SUBROUTINE ice_dyn_rhg
104
105
106   SUBROUTINE ice_dyn_rhg_init
107      !!-------------------------------------------------------------------
108      !!                  ***  ROUTINE ice_dyn_rhg_init  ***
109      !!
110      !! ** Purpose : Physical constants and parameters linked to the ice
111      !!      dynamics
112      !!
113      !! ** Method  :  Read the namdyn_rhg namelist and check the ice-dynamic
114      !!       parameter values called at the first timestep (nit000)
115      !!
116      !! ** input   :   Namelist namdyn_rhg
117      !!-------------------------------------------------------------------
118      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
119      !!
120      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, ln_rhg_chkcvg, &           !- evp
121    &                       ln_rhg_VP, nn_nout_vp, nn_ninn_vp, ln_zebra_vp, rn_relaxu_vp, rn_relaxv_vp, rn_uerr_max_vp, rn_uerr_min_vp, nn_cvgchk_vp  !-vp
122      !!-------------------------------------------------------------------
123     
124      !
125      REWIND( numnam_ice_ref )         ! Namelist namdyn_rhg in reference namelist : Ice dynamics
126      READ  ( numnam_ice_ref, namdyn_rhg, IOSTAT = ios, ERR = 901)
127901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_rhg in reference namelist' )
128      REWIND( numnam_ice_cfg )         ! Namelist namdyn_rhg in configuration namelist : Ice dynamics
129      READ  ( numnam_ice_cfg, namdyn_rhg, IOSTAT = ios, ERR = 902 )
130902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namdyn_rhg in configuration namelist' )
131      IF(lwm) WRITE ( numoni, namdyn_rhg )
132      !
133      IF(lwp) THEN                     ! control print
134         WRITE(numout,*)
135         WRITE(numout,*) 'ice_dyn_rhg_init: ice parameters for ice dynamics '
136         WRITE(numout,*) '~~~~~~~~~~~~~~~'
137         WRITE(numout,*) '   Namelist : namdyn_rhg:'
138         WRITE(numout,*) '      rheology EVP (icedyn_rhg_evp)                        ln_rhg_EVP    = ', ln_rhg_EVP
139         WRITE(numout,*) '         use adaptive EVP (aEVP)                           ln_aEVP       = ', ln_aEVP
140         WRITE(numout,*) '         creep limit                                       rn_creepl     = ', rn_creepl     ! common with VP
141         WRITE(numout,*) '         eccentricity of the elliptical yield curve        rn_ecc        = ', rn_ecc        ! common with VP
142         WRITE(numout,*) '         number of iterations for subcycling               nn_nevp       = ', nn_nevp
143         WRITE(numout,*) '         ratio of elastic timescale over ice time step     rn_relast     = ', rn_relast
144         WRITE(numout,*) '      check convergence of rheology                        ln_rhg_chkcvg = ', ln_rhg_chkcvg ! common with VP
145         WRITE(numout,*) ''   
146         WRITE(numout,*) '      rheology VP     (icedyn_rhg_VP)                      ln_rhg_VP     = ', ln_rhg_VP
147         WRITE(numout,*) '         number of outer iterations                        nn_nout_vp    = ', nn_nout_vp
148         WRITE(numout,*) '         number of inner iterations                        nn_ninn_vp    = ', nn_ninn_vp
149         WRITE(numout,*) '         activate zebra solver                             ln_zebra_vp   = ', ln_zebra_vp
150         WRITE(numout,*) '         relaxation factor for u                           rn_relaxu_vp  = ', rn_relaxu_vp
151         WRITE(numout,*) '         relaxation factor for v                           rn_relaxv_vp  = ', rn_relaxv_vp
152         WRITE(numout,*) '         maximum error on velocity                         rn_uerr_max_vp = ', rn_uerr_max_vp
153         WRITE(numout,*) '         velocity to decide convergence                    rn_uerr_min_vp = ', rn_uerr_min_vp
154         WRITE(numout,*) '         iteration step for convergence check              nn_cvgchk_vp   = ', nn_cvgchk_vp
155                           
156      ENDIF
157      !
158      !                             !== set the choice of ice advection ==!
159      ioptio = 0 
160      IF( ln_rhg_EVP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEVP    ;   ENDIF
161!!    IF( ln_rhg_EAP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEAP    ;   ENDIF
162      IF( ln_rhg_VP  ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgVP     ;   ENDIF
163      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_dyn_rhg_init: choose one and only one ice rheology' )
164      !
165      IF( ln_rhg_EVP  )   CALL rhg_evp_rst( 'READ' )  !* read or initialize all required files
166      !
167   END SUBROUTINE ice_dyn_rhg_init
168
169#else
170   !!----------------------------------------------------------------------
171   !!   Default option         Empty module           NO SI3 sea-ice model
172   !!----------------------------------------------------------------------
173#endif 
174
175   !!======================================================================
176END MODULE icedyn_rhg
Note: See TracBrowser for help on using the repository browser.