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.F90 in NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/ICE – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/ICE/icedyn.F90 @ 11954

Last change on this file since 11954 was 11671, checked in by acc, 5 years ago

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.
  • Property svn:keywords set to Id
File size: 14.3 KB
Line 
1MODULE icedyn
2   !!======================================================================
3   !!                     ***  MODULE  icedyn  ***
4   !!   Sea-Ice dynamics : master routine for sea ice dynamics
5   !!======================================================================
6   !! history :  4.0  ! 2018  (C. Rousset)  original code SI3 [aka Sea Ice cube]
7   !!----------------------------------------------------------------------
8#if defined key_si3
9   !!----------------------------------------------------------------------
10   !!   'key_si3'                                       SI3 sea-ice model
11   !!----------------------------------------------------------------------
12   !!   ice_dyn       : dynamics of sea ice
13   !!   ice_dyn_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     ! sea-ice: rheology
19   USE icedyn_adv     ! sea-ice: advection
20   USE icedyn_rdgrft  ! sea-ice: ridging/rafting
21   USE icecor         ! sea-ice: corrections
22   USE icevar         ! sea-ice: operations
23   USE icectl         ! sea-ice: control prints
24   !
25   USE in_out_manager ! I/O manager
26   USE iom            ! I/O manager library
27   USE lib_mpp        ! MPP library
28   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
29   USE lbclnk         ! lateral boundary conditions (or mpp links)
30   USE timing         ! Timing
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   ice_dyn        ! called by icestp.F90
36   PUBLIC   ice_dyn_init   ! called by icestp.F90
37   
38   INTEGER ::              nice_dyn   ! choice of the type of dynamics
39   !                                        ! associated indices:
40   INTEGER, PARAMETER ::   np_dynALL     = 1   ! full ice dynamics               (rheology + advection + ridging/rafting + correction)
41   INTEGER, PARAMETER ::   np_dynRHGADV  = 2   ! pure dynamics                   (rheology + advection)
42   INTEGER, PARAMETER ::   np_dynADV1D   = 3   ! only advection 1D - test case from Schar & Smolarkiewicz 1996
43   INTEGER, PARAMETER ::   np_dynADV2D   = 4   ! only advection 2D w prescribed vel.(rn_uvice + advection)
44   !
45   ! ** namelist (namdyn) **
46   LOGICAL  ::   ln_dynALL        ! full ice dynamics                      (rheology + advection + ridging/rafting + correction)
47   LOGICAL  ::   ln_dynRHGADV     ! no ridge/raft & no corrections         (rheology + advection)
48   LOGICAL  ::   ln_dynADV1D      ! only advection in 1D w ice convergence (test case from Schar & Smolarkiewicz 1996)
49   LOGICAL  ::   ln_dynADV2D      ! only advection in 2D w prescribed vel. (rn_uvice + advection)
50   REAL(wp) ::   rn_uice          !    prescribed u-vel (case np_dynADV1D & np_dynADV2D)
51   REAL(wp) ::   rn_vice          !    prescribed v-vel (case np_dynADV1D & np_dynADV2D)
52   
53   !! * Substitutions
54#  include "vectopt_loop_substitute.h90"
55   !!----------------------------------------------------------------------
56   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
57   !! $Id$
58   !! Software governed by the CeCILL licence     (./LICENSE)
59   !!----------------------------------------------------------------------
60CONTAINS
61
62   SUBROUTINE ice_dyn( kt )
63      !!-------------------------------------------------------------------
64      !!               ***  ROUTINE ice_dyn  ***
65      !!               
66      !! ** Purpose :   this routine manages sea ice dynamics
67      !!
68      !! ** Action : - calculation of friction in case of landfast ice
69      !!             - call ice_dyn_rhg    = rheology
70      !!             - call ice_dyn_adv    = advection
71      !!             - call ice_dyn_rdgrft = ridging/rafting
72      !!             - call ice_cor        = corrections if fields are out of bounds
73      !!--------------------------------------------------------------------
74      INTEGER, INTENT(in) ::   kt     ! ice time step
75      !!
76      INTEGER  ::   ji, jj        ! dummy loop indices
77      REAL(wp) ::   zcoefu, zcoefv
78      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdivu_i
79      !!--------------------------------------------------------------------
80      !
81      ! controls
82      IF( ln_timing )   CALL timing_start('icedyn')
83      !
84      IF( kt == nit000 .AND. lwp ) THEN
85         WRITE(numout,*)
86         WRITE(numout,*)'ice_dyn: sea-ice dynamics'
87         WRITE(numout,*)'~~~~~~~'
88      ENDIF
89      !                     
90      ! retrieve thickness from volume for landfast param. and UMx advection scheme
91      WHERE( a_i(:,:,:) >= epsi20 )
92         h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:)
93         h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:)
94      ELSEWHERE
95         h_i(:,:,:) = 0._wp
96         h_s(:,:,:) = 0._wp
97      END WHERE
98      !
99      WHERE( a_ip(:,:,:) >= epsi20 )
100         h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)
101      ELSEWHERE
102         h_ip(:,:,:) = 0._wp
103      END WHERE
104      !
105      !
106      SELECT CASE( nice_dyn )          !-- Set which dynamics is running
107
108      CASE ( np_dynALL )           !==  all dynamical processes  ==!
109         !
110         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology 
111         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
112         CALL ice_dyn_rdgrft( kt )                                          ! -- ridging/rafting
113         CALL ice_cor       ( kt , 1 )                                      ! -- Corrections
114         !
115      CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==!
116         !
117         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology 
118         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
119         CALL Hpiling                                                       ! -- simple pile-up (replaces ridging/rafting)
120         CALL ice_var_zapsmall                                              ! -- zap small areas
121         !
122      CASE ( np_dynADV1D )         !==  pure advection ==!   (1D)
123         !
124         ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- !
125         ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length
126         ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s
127         DO jj = 1, jpj
128            DO ji = 1, jpi
129               zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. )
130               zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. )
131               u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1., zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1)
132               v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1., zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1)
133            END DO
134         END DO
135         ! ---
136         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
137         !
138      CASE ( np_dynADV2D )         !==  pure advection ==!   (2D w prescribed velocities)
139         !
140         u_ice(:,:) = rn_uice * umask(:,:,1)
141         v_ice(:,:) = rn_vice * vmask(:,:,1)
142         !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1)
143         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1)
144         ! ---
145         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
146
147      END SELECT
148      !
149      !
150      ! diagnostics: divergence at T points
151      IF( iom_use('icediv') ) THEN
152         !
153         SELECT CASE( nice_dyn )
154
155         CASE ( np_dynADV1D , np_dynADV2D )
156
157            ALLOCATE( zdivu_i(jpi,jpj) )
158            DO jj = 2, jpjm1
159               DO ji = 2, jpim1
160                  zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
161                     &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
162               END DO
163            END DO
164            CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
165            ! output
166            CALL iom_put( 'icediv' , zdivu_i )
167
168            DEALLOCATE( zdivu_i )
169
170         END SELECT
171         !
172      ENDIF
173      !
174      ! controls
175      IF( ln_timing )   CALL timing_stop ('icedyn')
176      !
177   END SUBROUTINE ice_dyn
178
179
180   SUBROUTINE Hpiling
181      !!-------------------------------------------------------------------
182      !!                  ***  ROUTINE Hpiling  ***
183      !!
184      !! ** Purpose : Simple conservative piling comparable with 1-cat models
185      !!
186      !! ** Method  : pile-up ice when no ridging/rafting
187      !!
188      !! ** input   : a_i
189      !!-------------------------------------------------------------------
190      INTEGER ::   jl         ! dummy loop indices
191      !!-------------------------------------------------------------------
192      ! controls
193      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
194      !
195      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
196      DO jl = 1, jpl
197         WHERE( at_i(:,:) > epsi20 )
198            a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  )
199         END WHERE
200      END DO
201      ! controls
202      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
203      !
204   END SUBROUTINE Hpiling
205
206
207   SUBROUTINE ice_dyn_init
208      !!-------------------------------------------------------------------
209      !!                  ***  ROUTINE ice_dyn_init  ***
210      !!
211      !! ** Purpose : Physical constants and parameters linked to the ice
212      !!      dynamics
213      !!
214      !! ** Method  :  Read the namdyn namelist and check the ice-dynamic
215      !!       parameter values called at the first timestep (nit000)
216      !!
217      !! ** input   :   Namelist namdyn
218      !!-------------------------------------------------------------------
219      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
220      !!
221      NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  &
222         &             rn_ishlat ,                                                           &
223         &             ln_landfast_L16, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile
224      !!-------------------------------------------------------------------
225      !
226      READ  ( numnam_ice_ref, namdyn, IOSTAT = ios, ERR = 901)
227901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn in reference namelist' )
228      READ  ( numnam_ice_cfg, namdyn, IOSTAT = ios, ERR = 902 )
229902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn in configuration namelist' )
230      IF(lwm) WRITE( numoni, namdyn )
231      !
232      IF(lwp) THEN                     ! control print
233         WRITE(numout,*)
234         WRITE(numout,*) 'ice_dyn_init: ice parameters for ice dynamics '
235         WRITE(numout,*) '~~~~~~~~~~~~'
236         WRITE(numout,*) '   Namelist namdyn:'
237         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr) ln_dynALL       = ', ln_dynALL
238         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                     ln_dynRHGADV    = ', ln_dynRHGADV
239         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)    ln_dynADV1D     = ', ln_dynADV1D
240         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                ln_dynADV2D     = ', ln_dynADV2D
241         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)   = (', rn_uice,',',rn_vice,')'
242         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat       = ', rn_ishlat
243         WRITE(numout,*) '      Landfast: param from Lemieux 2016                      ln_landfast_L16 = ', ln_landfast_L16
244         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_depfra       = ', rn_depfra
245         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr       = ', rn_icebfr
246         WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lfrelax      = ', rn_lfrelax
247         WRITE(numout,*) '         isotropic tensile strength                          rn_tensile      = ', rn_tensile
248         WRITE(numout,*)
249      ENDIF
250      !                             !== set the choice of ice dynamics ==!
251      ioptio = 0 
252      !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction)
253      IF( ln_dynALL    ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynALL       ;   ENDIF
254      !      !--- dynamics without ridging/rafting and corr   (rheology + advection)
255      IF( ln_dynRHGADV ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV    ;   ENDIF
256      !      !--- advection 1D only - test case from Schar & Smolarkiewicz 1996
257      IF( ln_dynADV1D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV1D     ;   ENDIF
258      !      !--- advection 2D only with prescribed ice velocities (from namelist)
259      IF( ln_dynADV2D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV2D     ;   ENDIF
260      !
261      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' )
262      !
263      !                                      !--- Lateral boundary conditions
264      IF     (      rn_ishlat == 0.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  free-slip'
265      ELSEIF (      rn_ishlat == 2.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  no-slip'
266      ELSEIF ( 0. < rn_ishlat .AND. rn_ishlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  partial-slip'
267      ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  strong-slip'
268      ENDIF
269      !                                      !--- Landfast ice
270      IF( .NOT.ln_landfast_L16 )   tau_icebfr(:,:) = 0._wp
271      !
272      CALL ice_dyn_rdgrft_init          ! set ice ridging/rafting parameters
273      CALL ice_dyn_rhg_init             ! set ice rheology parameters
274      CALL ice_dyn_adv_init             ! set ice advection parameters
275      !
276   END SUBROUTINE ice_dyn_init
277
278#else
279   !!----------------------------------------------------------------------
280   !!   Default option         Empty module           NO SI3 sea-ice model
281   !!----------------------------------------------------------------------
282#endif 
283
284   !!======================================================================
285END MODULE icedyn
Note: See TracBrowser for help on using the repository browser.