source: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn.F90 @ 12475

Last change on this file since 12475 was 12475, checked in by dancopsey, 14 months ago
  • Add more print statements.
  • Move away from using snow to ice diagnostics and use a new snow to pond one instead.
File size: 16.2 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         lh_ip(:,:,:) = 0._wp
104      END WHERE
105      !
106
107      write(numout,*)'ice_dyn 0:sv_i, e_i, t_i = ',sv_i(3,4,1), ' ', e_i(3,4,1,1), ' ', t_i(3,4,1,1)
108      write(numout,*)'ice_dyn 0:a_i, h_i, v_i = ',a_i(3,4,1), ' ', h_i(3,4,1), ' ', v_i(3,4,1)
109      write(numout,*)'ice_dyn 0:s_i = ',s_i(3,4,1)
110      write(numout,*)'ice_dyn 0: MEAN(sv_i) = ',SUM(sv_i)/SIZE(sv_i)
111      write(numout,*)'ice_dyn 0: MEAN(e_i) = ',SUM(e_i)/SIZE(e_i)
112      write(numout,*)'ice_dyn 0: MEAN(t_i) = ',SUM(t_i)/SIZE(t_i)
113      write(numout,*)'ice_dyn 0: MEAN(a_i) = ',SUM(a_i)/SIZE(a_i)
114      write(numout,*)'ice_dyn 0: MEAN(h_i) = ',SUM(h_i)/SIZE(h_i)
115      write(numout,*)'ice_dyn 0: MEAN(v_i) = ',SUM(v_i)/SIZE(v_i)
116      write(numout,*)'ice_dyn 0: MEAN(s_i) = ',SUM(s_i)/SIZE(s_i)
117      write(numout,*)'ice_dyn 0: MEAN(u_ice) = ',SUM(u_ice)/SIZE(u_ice)
118      write(numout,*)'ice_dyn 0: MEAN(e2u) = ',SUM(e2u)/SIZE(e2u)
119      write(numout,*)'ice_dyn 0: u_ice, e2u = ',u_ice(3,4), e2u(3,4)
120     
121      !
122      SELECT CASE( nice_dyn )          !-- Set which dynamics is running
123
124      CASE ( np_dynALL )           !==  all dynamical processes  ==!
125         !
126         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology
127          write(numout,*)'ice_dyn 1:sv_i, e_s = ',sv_i(3,4,1), ' ', e_i (3,4,1,1)
128          write(numout,*)'ice_dyn 1: u_ice, e2u = ',u_ice(3,4), e2u(3,4)
129         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
130          write(numout,*)'ice_dyn 2:sv_i, e_s = ',sv_i(3,4,1), ' ', e_i (3,4,1,1)
131         CALL ice_dyn_rdgrft( kt )                                          ! -- ridging/rafting
132          write(numout,*)'ice_dyn 3:sv_i, e_s = ',sv_i(3,4,1), ' ', e_i (3,4,1,1)
133         CALL ice_cor       ( kt , 1 )                                      ! -- Corrections
134         !
135      CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==!
136         !
137         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology 
138         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
139         CALL Hpiling                                                       ! -- simple pile-up (replaces ridging/rafting)
140         CALL ice_var_zapsmall                                              ! -- zap small areas
141         !
142      CASE ( np_dynADV1D )         !==  pure advection ==!   (1D)
143         !
144         ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- !
145         ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length
146         ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s
147         DO jj = 1, jpj
148            DO ji = 1, jpi
149               zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. )
150               zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. )
151               u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1., zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1)
152               v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1., zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1)
153            END DO
154         END DO
155         ! ---
156         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
157         !
158      CASE ( np_dynADV2D )         !==  pure advection ==!   (2D w prescribed velocities)
159         !
160         u_ice(:,:) = rn_uice * umask(:,:,1)
161         v_ice(:,:) = rn_vice * vmask(:,:,1)
162         !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1)
163         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1)
164         ! ---
165         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
166
167      END SELECT
168
169      write(numout,*)'ice_dyn 4: e_i = ',e_i (42,26,1,1)
170      !
171      !
172      ! diagnostics: divergence at T points
173      IF( iom_use('icediv') ) THEN
174         !
175         SELECT CASE( nice_dyn )
176
177         CASE ( np_dynADV1D , np_dynADV2D )
178
179            ALLOCATE( zdivu_i(jpi,jpj) )
180            DO jj = 2, jpjm1
181               DO ji = 2, jpim1
182                  zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
183                     &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
184               END DO
185            END DO
186            CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
187            CALL iom_put( "icediv" , zdivu_i(:,:) )
188            DEALLOCATE( zdivu_i )
189
190         END SELECT
191         !
192      ENDIF
193      write(numout,*)'ice_dyn 5: e_i = ',e_i (42,26,1,1)
194      !
195      ! controls
196      IF( ln_timing )   CALL timing_stop ('icedyn')
197      !
198   END SUBROUTINE ice_dyn
199
200
201   SUBROUTINE Hpiling
202      !!-------------------------------------------------------------------
203      !!                  ***  ROUTINE Hpiling  ***
204      !!
205      !! ** Purpose : Simple conservative piling comparable with 1-cat models
206      !!
207      !! ** Method  : pile-up ice when no ridging/rafting
208      !!
209      !! ** input   : a_i
210      !!-------------------------------------------------------------------
211      INTEGER ::   jl         ! dummy loop indices
212      !!-------------------------------------------------------------------
213      ! controls
214      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
215      !
216      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
217      DO jl = 1, jpl
218         WHERE( at_i(:,:) > epsi20 )
219            a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  )
220         END WHERE
221      END DO
222      ! controls
223      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
224      !
225   END SUBROUTINE Hpiling
226
227
228   SUBROUTINE ice_dyn_init
229      !!-------------------------------------------------------------------
230      !!                  ***  ROUTINE ice_dyn_init  ***
231      !!
232      !! ** Purpose : Physical constants and parameters linked to the ice
233      !!      dynamics
234      !!
235      !! ** Method  :  Read the namdyn namelist and check the ice-dynamic
236      !!       parameter values called at the first timestep (nit000)
237      !!
238      !! ** input   :   Namelist namdyn
239      !!-------------------------------------------------------------------
240      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
241      !!
242      NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  &
243         &             rn_ishlat ,                                                           &
244         &             ln_landfast_L16, ln_landfast_home, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile
245      !!-------------------------------------------------------------------
246      !
247      REWIND( numnam_ice_ref )         ! Namelist namdyn in reference namelist : Ice dynamics
248      READ  ( numnam_ice_ref, namdyn, IOSTAT = ios, ERR = 901)
249901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn in reference namelist', lwp )
250      REWIND( numnam_ice_cfg )         ! Namelist namdyn in configuration namelist : Ice dynamics
251      READ  ( numnam_ice_cfg, namdyn, IOSTAT = ios, ERR = 902 )
252902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn in configuration namelist', lwp )
253      IF(lwm) WRITE( numoni, namdyn )
254      !
255      IF(lwp) THEN                     ! control print
256         WRITE(numout,*)
257         WRITE(numout,*) 'ice_dyn_init: ice parameters for ice dynamics '
258         WRITE(numout,*) '~~~~~~~~~~~~'
259         WRITE(numout,*) '   Namelist namdyn:'
260         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr) ln_dynALL       = ', ln_dynALL
261         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                     ln_dynRHGADV    = ', ln_dynRHGADV
262         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)    ln_dynADV1D     = ', ln_dynADV1D
263         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                ln_dynADV2D     = ', ln_dynADV2D
264         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)   = (', rn_uice,',',rn_vice,')'
265         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat       = ', rn_ishlat
266         WRITE(numout,*) '      Landfast: param from Lemieux 2016                      ln_landfast_L16 = ', ln_landfast_L16
267         WRITE(numout,*) '      Landfast: param from home made                         ln_landfast_home= ', ln_landfast_home
268         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_depfra       = ', rn_depfra
269         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr       = ', rn_icebfr
270         WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lfrelax      = ', rn_lfrelax
271         WRITE(numout,*) '         isotropic tensile strength                          rn_tensile      = ', rn_tensile
272         WRITE(numout,*)
273      ENDIF
274      !                             !== set the choice of ice dynamics ==!
275      ioptio = 0 
276      !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction)
277      IF( ln_dynALL    ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynALL       ;   ENDIF
278      !      !--- dynamics without ridging/rafting and corr   (rheology + advection)
279      IF( ln_dynRHGADV ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV    ;   ENDIF
280      !      !--- advection 1D only - test case from Schar & Smolarkiewicz 1996
281      IF( ln_dynADV1D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV1D     ;   ENDIF
282      !      !--- advection 2D only with prescribed ice velocities (from namelist)
283      IF( ln_dynADV2D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV2D     ;   ENDIF
284      !
285      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' )
286      !
287      !                                      !--- Lateral boundary conditions
288      IF     (      rn_ishlat == 0.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  free-slip'
289      ELSEIF (      rn_ishlat == 2.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  no-slip'
290      ELSEIF ( 0. < rn_ishlat .AND. rn_ishlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  partial-slip'
291      ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  strong-slip'
292      ENDIF
293      !                                      !--- Landfast ice
294      IF( .NOT.ln_landfast_L16 .AND. .NOT.ln_landfast_home )   tau_icebfr(:,:) = 0._wp
295      !
296      IF ( ln_landfast_L16 .AND. ln_landfast_home ) THEN
297         CALL ctl_stop( 'ice_dyn_init: choose one and only one landfast parameterization (ln_landfast_L16 or ln_landfast_home)' )
298      ENDIF
299      !
300      CALL ice_dyn_rdgrft_init          ! set ice ridging/rafting parameters
301      CALL ice_dyn_rhg_init             ! set ice rheology parameters
302      CALL ice_dyn_adv_init             ! set ice advection parameters
303      !
304   END SUBROUTINE ice_dyn_init
305
306#else
307   !!----------------------------------------------------------------------
308   !!   Default option         Empty module           NO SI3 sea-ice model
309   !!----------------------------------------------------------------------
310#endif 
311
312   !!======================================================================
313END MODULE icedyn
Note: See TracBrowser for help on using the repository browser.