source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90 @ 8517

Last change on this file since 8517 was 8517, checked in by clem, 3 years ago

changes in style - part6 - one more round

File size: 14.0 KB
Line 
1MODULE icedyn
2   !!======================================================================
3   !!                     ***  MODULE  icedyn  ***
4   !!   Sea-Ice dynamics : 
5   !!======================================================================
6   !! history :  4.0  ! 2017-09  (C. Rousset)  original code
7   !!----------------------------------------------------------------------
8#if defined key_lim3
9   !!----------------------------------------------------------------------
10   !!   'key_lim3'                                       LIM3 sea-ice model
11   !!----------------------------------------------------------------------
12   !!   ice_dyn       : dynamics of sea ice
13   !!   ice_dyn_init  : initialisation of sea-ice dynamics
14   !!----------------------------------------------------------------------
15   USE phycst         ! physical constants
16   USE dom_oce        ! ocean space and time domain
17   USE ice            ! sea-ice: variables
18   USE icerhg         ! sea-ice: rheology
19   USE iceadv         ! sea-ice: advection
20   USE icerdgrft      ! sea-ice: ridging/rafting
21   USE icecor         ! sea-ice: corrections
22   USE icevar         ! sea-ice: operations
23   !
24   USE lbclnk         ! lateral boundary conditions - MPP exchanges
25   USE lib_mpp        ! MPP library
26   USE in_out_manager ! I/O manager
27   USE lib_fortran    ! glob_sum
28   USE timing         ! Timing
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   ice_dyn        ! called by icestp.F90
34   PUBLIC   ice_dyn_init   ! called by icestp.F90
35   
36   INTEGER ::              nice_dyn   ! choice of the type of dynamics
37   !                                        ! associated indices:
38   INTEGER, PARAMETER ::   np_dynFULL    = 1   ! full ice dynamics               (rheology + advection + ridging/rafting + correction)
39   INTEGER, PARAMETER ::   np_dynRHGADV1 = 2   ! no ridging/rafting              (rheology + advection                   + correction)
40   INTEGER, PARAMETER ::   np_dynRHGADV2 = 3   ! pure dynamics                   (rheology + advection)
41   INTEGER, PARAMETER ::   np_dynADV     = 4   ! only advection w prescribed vel.(rn_uvice + advection)
42   !
43   ! ** namelist (namdyn) **
44   REAL(wp) ::   rn_uice          ! prescribed u-vel (case np_dynADV)
45   REAL(wp) ::   rn_vice          ! prescribed v-vel (case np_dynADV)
46   
47   !! * Substitutions
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
51   !! $Id: icedyn.F90 8378 2017-07-26 13:55:59Z clem $
52   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE ice_dyn( kt )
57      !!-------------------------------------------------------------------
58      !!               ***  ROUTINE ice_dyn  ***
59      !!               
60      !! ** Purpose :   this routine manages sea ice dynamics
61      !!
62      !! ** Action : - Initialisation of some variables
63      !!             - call ice_rhg
64      !!--------------------------------------------------------------------
65      INTEGER, INTENT(in) ::   kt     ! ice time step
66      !!
67      INTEGER ::   ji, jj, jl         ! dummy loop indices
68      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhmax
69      !!--------------------------------------------------------------------
70      !
71      IF( nn_timing == 1 )  CALL timing_start('icedyn')
72      !
73      IF( kt == nit000 .AND. lwp ) THEN
74         WRITE(numout,*)
75         WRITE(numout,*)'ice_dyn: sea-ice dynamics'
76         WRITE(numout,*)'~~~~~~~'
77      ENDIF
78
79      !                     
80      IF( ln_landfast ) THEN            !-- Landfast ice parameterization: define max bottom friction
81         tau_icebfr(:,:) = 0._wp
82         DO jl = 1, jpl
83            WHERE( ht_i(:,:,jl) > ht_n(:,:) * rn_gamma )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
84         END DO
85      ENDIF
86
87      zhmax(:,:,:) = ht_i_b(:,:,:)      !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig)
88      DO jl = 1, jpl
89         DO jj = 2, jpjm1
90            DO ji = 2, jpim1
91!!gm use of MAXVAL here is very probably less efficient than expending the 9 values
92               zhmax(ji,jj,jl) = MAX( epsi20, MAXVAL( ht_i_b(ji-1:ji+1,jj-1:jj+1,jl) ) )
93            END DO
94         END DO
95      END DO
96      CALL lbc_lnk( zhmax(:,:,:), 'T', 1. )
97      !
98      !
99      SELECT CASE( nice_dyn )           !-- Set which dynamics is running
100
101      CASE ( np_dynFULL )          !==  all dynamical processes  ==!
102         CALL ice_rhg   ( kt )                            ! -- rheology 
103         CALL ice_adv   ( kt )   ;   CALL Hbig( zhmax )   ! -- advection of ice + correction on ice thickness
104         CALL ice_rdgrft( kt )                            ! -- ridging/rafting
105         CALL ice_cor   ( kt , 1 )                        ! -- Corrections
106
107      CASE ( np_dynRHGADV1 )       !==  no ridge/raft ==!   (mono cat. case 2)
108         CALL ice_rhg   ( kt )                            ! -- rheology 
109         CALL ice_adv   ( kt )                            ! -- advection of ice
110         CALL Hpiling                                     ! -- simple pile-up (replaces ridging/rafting)
111         CALL ice_cor   ( kt , 1 )                        ! -- Corrections
112
113      CASE ( np_dynRHGADV2 )       !==  no ridge/raft & no corrections ==!
114         CALL ice_rhg   ( kt )                            ! -- rheology 
115         CALL ice_adv   ( kt )                            ! -- advection of ice
116         CALL Hpiling                                     ! -- simple pile-up (replaces ridging/rafting)
117
118      CASE ( np_dynADV )           !==  pure advection ==!   (prescribed velocities)
119         u_ice(:,:) = rn_uice * umask(:,:,1)
120         v_ice(:,:) = rn_vice * vmask(:,:,1)
121         !!CALL RANDOM_NUMBER(u_ice(:,:))
122         !!CALL RANDOM_NUMBER(v_ice(:,:))
123         CALL ice_adv   ( kt )                            ! -- advection of ice
124
125      END SELECT
126      !
127      IF( nn_timing == 1 )   CALL timing_stop('icedyn')
128      !
129   END SUBROUTINE ice_dyn
130
131   SUBROUTINE Hbig( phmax )
132      !!-------------------------------------------------------------------
133      !!                  ***  ROUTINE Hbig  ***
134      !!
135      !! ** Purpose : Thickness correction in case advection scheme creates
136      !!              abnormally tick ice
137      !!
138      !! ** Method  : 1- check whether ice thickness resulting from advection is
139      !!                 larger than the surrounding 9-points before advection
140      !!                 and reduce it if a) divergence or b) convergence & at_i>0.8
141      !!              2- bound ice thickness with hi_max (99m)
142      !!
143      !! ** input   : Max thickness of the surrounding 9-points
144      !!-------------------------------------------------------------------
145      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phmax   ! max ice thick from surrounding 9-pts
146      !
147      INTEGER  ::   ji, jj, jl         ! dummy loop indices
148      REAL(wp) ::   zh, zdv
149      !!-------------------------------------------------------------------
150      !
151      CALL ice_var_zapsmall                       !-- zap small areas
152      !
153      DO jl = 1, jpl
154         DO jj = 1, jpj
155            DO ji = 1, jpi
156               IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to hmax
157                  !
158                  zh  = v_i (ji,jj,jl) / a_i(ji,jj,jl)
159                  zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl) 
160                  !
161                  IF ( ( zdv >  0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. &
162                     & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN
163                     a_i (ji,jj,jl) = v_i(ji,jj,jl) / MIN( phmax(ji,jj,jl), hi_max(jpl) )   !-- bound ht_i to hi_max (99 m)
164                  ENDIF
165                  !
166               ENDIF
167            END DO
168         END DO
169      END DO
170           
171      IF ( nn_pnd_scheme > 0 ) THEN               !-- correct pond fraction to avoid a_ip > a_i
172         WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:)
173      ENDIF
174      !
175   END SUBROUTINE Hbig
176
177   SUBROUTINE Hpiling
178      !!-------------------------------------------------------------------
179      !!                  ***  ROUTINE Hpiling  ***
180      !!
181      !! ** Purpose : Simple conservative piling comparable with 1-cat models
182      !!
183      !! ** Method  : pile-up ice when no ridging/rafting
184      !!
185      !! ** input   : a_i
186      !!-------------------------------------------------------------------
187      INTEGER ::   jl         ! dummy loop indices
188      !!-------------------------------------------------------------------
189      !
190      CALL ice_var_zapsmall                       !-- zap small areas
191      !
192      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
193      DO jl = 1, jpl
194         WHERE( at_i(:,:) > epsi20 )
195            a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  )
196         END WHERE
197      END DO
198      !
199   END SUBROUTINE Hpiling
200
201   SUBROUTINE ice_dyn_init
202      !!-------------------------------------------------------------------
203      !!                  ***  ROUTINE ice_dyn_init  ***
204      !!
205      !! ** Purpose : Physical constants and parameters linked to the ice
206      !!      dynamics
207      !!
208      !! ** Method  :  Read the namice_dyn namelist and check the ice-dynamic
209      !!       parameter values called at the first timestep (nit000)
210      !!
211      !! ** input   :   Namelist namice_dyn
212      !!-------------------------------------------------------------------
213      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
214      !!
215      NAMELIST/namice_dyn/ ln_dynFULL, ln_dynRHGADV, ln_dynADV, rn_uice, rn_vice,  &
216         &                 rn_ishlat  , rn_cio   ,                         &
217         &                 ln_landfast, rn_gamma , rn_icebfr, rn_lfrelax
218      !!-------------------------------------------------------------------
219      !
220      REWIND( numnam_ice_ref )         ! Namelist namice_dyn in reference namelist : Ice dynamics
221      READ  ( numnam_ice_ref, namice_dyn, IOSTAT = ios, ERR = 901)
222901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_dyn in reference namelist', lwp )
223      !
224      REWIND( numnam_ice_cfg )         ! Namelist namice_dyn in configuration namelist : Ice dynamics
225      READ  ( numnam_ice_cfg, namice_dyn, IOSTAT = ios, ERR = 902 )
226902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_dyn in configuration namelist', lwp )
227      IF(lwm) WRITE ( numoni, namice_dyn )
228      !
229      IF(lwp) THEN                     ! control print
230         WRITE(numout,*)
231         WRITE(numout,*) 'ice_dyn_init: ice parameters for ice dynamics '
232         WRITE(numout,*) '~~~~~~~~~~~~'
233         WRITE(numout,*) '   Namelist namice_dyn'
234         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynFULL   = ', ln_dynFULL
235         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV = ', ln_dynRHGADV
236         WRITE(numout,*) '      Advection only         (rn_uvice + adv)                 ln_dynADV    = ', ln_dynADV
237         WRITE(numout,*) '           with prescribed velocity given by '
238         WRITE(numout,*) '               a uniform field               (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',', rn_vice,')'
239         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat     = ', rn_ishlat
240         WRITE(numout,*) '      drag coefficient for oceanic stress                    rn_cio        = ', rn_cio
241         WRITE(numout,*) '      Landfast: param (T or F)                               ln_landfast   = ', ln_landfast
242         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_gamma      = ', rn_gamma
243         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr     = ', rn_icebfr
244         WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lfrelax    = ', rn_lfrelax
245      ENDIF
246      !                             !== set the choice of ice dynamics ==!
247      ioptio = 0 
248      !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction)
249      IF( ln_dynFULL .AND. nn_monocat /= 2 ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynFULL      ;   ENDIF
250      !      !--- dynamics without ridging/rafting            (rheology + advection                   + correction)
251      IF( ln_dynFULL .AND. nn_monocat == 2 ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV1   ;   ENDIF
252      !      !--- dynamics without ridging/rafting and corr   (rheology + advection)
253      IF( ln_dynRHGADV                     ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV2   ;   ENDIF
254      !      !--- advection only with prescribed ice velocities (from namelist)
255      IF( ln_dynADV                        ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV       ;   ENDIF
256      !
257      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' )
258      !
259      !                                      !--- Lateral boundary conditions
260      IF     (      rn_ishlat == 0.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  free-slip'
261      ELSEIF (      rn_ishlat == 2.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  no-slip'
262      ELSEIF ( 0. < rn_ishlat .AND. rn_ishlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  partial-slip'
263      ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  strong-slip'
264      ENDIF
265      !                                      !--- NO Landfast ice : set to zero once for all
266      IF( .NOT. ln_landfast )   tau_icebfr(:,:) = 0._wp 
267      !
268   END SUBROUTINE ice_dyn_init
269
270#else
271   !!----------------------------------------------------------------------
272   !!   Default option         Empty module          NO LIM-3 sea-ice model
273   !!----------------------------------------------------------------------
274#endif 
275
276   !!======================================================================
277END MODULE icedyn
Note: See TracBrowser for help on using the repository browser.