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.
limdyn.F90 in branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90 @ 6515

Last change on this file since 6515 was 6515, checked in by clem, 8 years ago

implement several developments for LIM3: new advection scheme (ultimate-macho, not yet perfect) ; lateral ice melt ; enabling/disabling thermo and dyn with namelist options ; simplifications (including a clarified namelist)
  • Property svn:keywords set to Id
File size: 12.1 KB
Line 
1MODULE limdyn
2   !!======================================================================
3   !!                     ***  MODULE  limdyn  ***
4   !!   Sea-Ice dynamics : 
5   !!======================================================================
6   !! history :  1.0  ! 2002-08  (C. Ethe, G. Madec)  original VP code
7   !!            3.0  ! 2007-03  (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle)  LIM3: EVP-Cgrid
8   !!            3.5  ! 2011-02  (G. Madec) dynamical allocation
9   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3' :                                 LIM3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!    lim_dyn      : computes ice velocities
15   !!    lim_dyn_init : initialization and namelist read
16   !!----------------------------------------------------------------------
17   USE phycst           ! physical constants
18   USE dom_oce          ! ocean space and time domain
19   USE sbc_ice          ! Surface boundary condition: ice   fields
20   USE ice              ! LIM-3 variables
21   USE dom_ice          ! LIM-3 domain
22   USE limrhg           ! LIM-3 rheology
23   USE lbclnk           ! lateral boundary conditions - MPP exchanges
24   USE lib_mpp          ! MPP library
25   USE wrk_nemo         ! work arrays
26   USE in_out_manager   ! I/O manager
27   USE prtctl           ! Print control
28   USE lib_fortran      ! glob_sum
29   USE timing          ! Timing
30   USE limcons        ! conservation tests
31   USE limvar
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   lim_dyn   ! routine called by ice_step
37
38   !! * Substitutions
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
42   !! $Id$
43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE lim_dyn( kt )
48      !!-------------------------------------------------------------------
49      !!               ***  ROUTINE lim_dyn  ***
50      !!               
51      !! ** Purpose :   compute ice velocity
52      !!               
53      !! ** Method  :
54      !!
55      !! ** Action  : - Initialisation
56      !!              - Call of the dynamic routine for each hemisphere
57      !!------------------------------------------------------------------------------------
58      INTEGER, INTENT(in) ::   kt     ! number of iteration
59      !!
60      INTEGER  ::   ji, jj, jl, ja    ! dummy loop indices
61      INTEGER  ::   i_j1, i_jpj       ! Starting/ending j-indices for rheology
62      REAL(wp), POINTER, DIMENSION(:)   ::   zswitch        ! i-averaged indicator of sea-ice
63      REAL(wp), POINTER, DIMENSION(:)   ::   zmsk           ! i-averaged of tmask
64      !
65      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
66     !!---------------------------------------------------------------------
67
68      IF( nn_timing == 1 )  CALL timing_start('limdyn')
69
70      CALL wrk_alloc( jpj, zswitch, zmsk )
71
72      CALL lim_var_agg(1)             ! aggregate ice categories
73
74      IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only)
75      !
76      ! conservation test
77      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
78     
79      u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1)
80      v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1)
81     
82      ! Rheology (ice dynamics)
83      ! ========
84     
85      !  Define the j-limits where ice rheology is computed
86      ! ---------------------------------------------------
87     
88      IF( lk_mpp .OR. lk_mpp_rep ) THEN                    ! mpp: compute over the whole domain
89         i_j1 = 1
90         i_jpj = jpj
91         IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
92         CALL lim_rhg( i_j1, i_jpj )
93      ELSE                                 ! optimization of the computational area
94         !
95         DO jj = 1, jpj
96            zswitch(jj) = SUM( 1.0 - at_i(:,jj) )   ! = REAL(jpj) if ocean everywhere on a j-line
97            zmsk   (jj) = SUM( tmask(:,jj,1)    )   ! = 0         if land  everywhere on a j-line
98         END DO
99         
100         IF( l_jeq ) THEN                     ! local domain include both hemisphere
101            !                                 ! Rheology is computed in each hemisphere
102            !                                 ! only over the ice cover latitude strip
103            ! Northern hemisphere
104            i_j1  = njeq
105            i_jpj = jpj
106            DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
107               i_j1 = i_j1 + 1
108            END DO
109            i_j1 = MAX( 1, i_j1-2 )
110            IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : NH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
111            CALL lim_rhg( i_j1, i_jpj )
112            !
113            ! Southern hemisphere
114            i_j1  =  1
115            i_jpj = njeq
116            DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
117               i_jpj = i_jpj - 1
118            END DO
119            i_jpj = MIN( jpj, i_jpj+1 )
120            IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : SH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
121            !
122            CALL lim_rhg( i_j1, i_jpj )
123            !
124         ELSE                                 ! local domain extends over one hemisphere only
125            !                                 ! Rheology is computed only over the ice cover
126            !                                 ! latitude strip
127            i_j1  = 1
128            DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
129               i_j1 = i_j1 + 1
130            END DO
131            i_j1 = MAX( 1, i_j1-2 )
132           
133            i_jpj  = jpj
134            DO WHILE ( i_jpj >= 1  .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
135               i_jpj = i_jpj - 1
136            END DO
137            i_jpj = MIN( jpj, i_jpj+1)
138            !
139            IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : one hemisphere:  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
140            !
141            CALL lim_rhg( i_j1, i_jpj )
142            !
143         ENDIF
144         !
145      ENDIF
146      !
147      IF(ln_ctl) THEN   ! Control print
148         CALL prt_ctl_info(' ')
149         CALL prt_ctl_info(' - Cell values : ')
150         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
151         CALL prt_ctl(tab2d_1=divu_i    , clinfo1=' lim_dyn  : divu_i    :')
152         CALL prt_ctl(tab2d_1=delta_i   , clinfo1=' lim_dyn  : delta_i   :')
153         CALL prt_ctl(tab2d_1=strength  , clinfo1=' lim_dyn  : strength  :')
154         CALL prt_ctl(tab2d_1=e12t      , clinfo1=' lim_dyn  : cell area :')
155         CALL prt_ctl(tab2d_1=at_i      , clinfo1=' lim_dyn  : at_i      :')
156         CALL prt_ctl(tab2d_1=vt_i      , clinfo1=' lim_dyn  : vt_i      :')
157         CALL prt_ctl(tab2d_1=vt_s      , clinfo1=' lim_dyn  : vt_s      :')
158         CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' lim_dyn  : stress1_i :')
159         CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' lim_dyn  : stress2_i :')
160         CALL prt_ctl(tab2d_1=stress12_i, clinfo1=' lim_dyn  : stress12_i:')
161         DO jl = 1, jpl
162            CALL prt_ctl_info(' ')
163            CALL prt_ctl_info(' - Category : ', ivar1=jl)
164            CALL prt_ctl_info('   ~~~~~~~~~~')
165            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_dyn  : a_i      : ')
166            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_i     : ')
167            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_s     : ')
168            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_dyn  : v_i      : ')
169            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_dyn  : v_s      : ')
170            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : e_s      : ')
171            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_dyn  : t_su     : ')
172            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : t_snow   : ')
173            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_dyn  : sm_i     : ')
174            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_dyn  : smv_i    : ')
175            DO ja = 1, nlay_i
176               CALL prt_ctl_info(' ')
177               CALL prt_ctl_info(' - Layer : ', ivar1=ja)
178               CALL prt_ctl_info('   ~~~~~~~')
179               CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : t_i      : ')
180               CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : e_i      : ')
181            END DO
182         END DO
183      ENDIF
184      !
185      ! conservation test
186      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
187      !
188      CALL wrk_dealloc( jpj, zswitch, zmsk )
189      !
190      IF( nn_timing == 1 )  CALL timing_stop('limdyn')
191
192   END SUBROUTINE lim_dyn
193
194
195   SUBROUTINE lim_dyn_init
196      !!-------------------------------------------------------------------
197      !!                  ***  ROUTINE lim_dyn_init  ***
198      !!
199      !! ** Purpose : Physical constants and parameters linked to the ice
200      !!      dynamics
201      !!
202      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
203      !!       parameter values called at the first timestep (nit000)
204      !!
205      !! ** input   :   Namelist namicedyn
206      !!-------------------------------------------------------------------
207      INTEGER  ::   ios                 ! Local integer output status for namelist read
208      NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, &
209         &                nn_nevp, rn_relast
210      !!-------------------------------------------------------------------
211
212      REWIND( numnam_ice_ref )              ! Namelist namicedyn in reference namelist : Ice dynamics
213      READ  ( numnam_ice_ref, namicedyn, IOSTAT = ios, ERR = 901)
214901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in reference namelist', lwp )
215
216      REWIND( numnam_ice_cfg )              ! Namelist namicedyn in configuration namelist : Ice dynamics
217      READ  ( numnam_ice_cfg, namicedyn, IOSTAT = ios, ERR = 902 )
218902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in configuration namelist', lwp )
219      IF(lwm) WRITE ( numoni, namicedyn )
220     
221      IF(lwp) THEN                        ! control print
222         WRITE(numout,*)
223         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
224         WRITE(numout,*) '~~~~~~~~~~~~'
225         WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)  nn_icestr     = ', nn_icestr 
226         WRITE(numout,*)'    Including brine volume in ice strength comp.         ln_icestr_bvf = ', ln_icestr_bvf
227         WRITE(numout,*)'    Ratio of ridging work to PotEner change in ridging   rn_pe_rdg     = ', rn_pe_rdg 
228         WRITE(numout,*) '   drag coefficient for oceanic stress                  rn_cio        = ', rn_cio
229         WRITE(numout,*) '   first bulk-rheology parameter                        rn_pstar      = ', rn_pstar
230         WRITE(numout,*) '   second bulk-rhelogy parameter                        rn_crhg       = ', rn_crhg
231         WRITE(numout,*) '   creep limit                                          rn_creepl     = ', rn_creepl
232         WRITE(numout,*) '   eccentricity of the elliptical yield curve           rn_ecc        = ', rn_ecc
233         WRITE(numout,*) '   number of iterations for subcycling                  nn_nevp       = ', nn_nevp
234         WRITE(numout,*) '   ratio of elastic timescale over ice time step        rn_relast     = ', rn_relast
235      ENDIF
236      !
237      usecc2 = 1._wp / ( rn_ecc * rn_ecc )
238      rhoco  = rau0  * rn_cio
239      !
240   END SUBROUTINE lim_dyn_init
241
242#else
243   !!----------------------------------------------------------------------
244   !!   Default option          Empty module           NO LIM sea-ice model
245   !!----------------------------------------------------------------------
246CONTAINS
247   SUBROUTINE lim_dyn         ! Empty routine
248   END SUBROUTINE lim_dyn
249#endif 
250
251   !!======================================================================
252END MODULE limdyn
Note: See TracBrowser for help on using the repository browser.