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/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90 @ 7963

Last change on this file since 7963 was 5870, checked in by acc, 9 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
  • Property svn:keywords set to Id
File size: 16.2 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_oce          ! Surface boundary condition: ocean fields
20   USE sbc_ice          ! Surface boundary condition: ice   fields
21   USE ice              ! LIM-3 variables
22   USE dom_ice          ! LIM-3 domain
23   USE limrhg           ! LIM-3 rheology
24   USE lbclnk           ! lateral boundary conditions - MPP exchanges
25   USE lib_mpp          ! MPP library
26   USE wrk_nemo         ! work arrays
27   USE in_out_manager   ! I/O manager
28   USE prtctl           ! Print control
29   USE lib_fortran      ! glob_sum
30   USE timing          ! Timing
31   USE limcons        ! conservation tests
32   USE limvar
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   lim_dyn   ! routine called by ice_step
38
39   !! * Substitutions
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
43   !! $Id$
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE lim_dyn( kt )
49      !!-------------------------------------------------------------------
50      !!               ***  ROUTINE lim_dyn  ***
51      !!               
52      !! ** Purpose :   compute ice velocity and ocean-ice stress
53      !!               
54      !! ** Method  :
55      !!
56      !! ** Action  : - Initialisation
57      !!              - Call of the dynamic routine for each hemisphere
58      !!              - computation of the stress at the ocean surface         
59      !!              - treatment of the case if no ice dynamic
60      !!------------------------------------------------------------------------------------
61      INTEGER, INTENT(in) ::   kt     ! number of iteration
62      !!
63      INTEGER  ::   ji, jj, jl, ja    ! dummy loop indices
64      INTEGER  ::   i_j1, i_jpj       ! Starting/ending j-indices for rheology
65      REAL(wp) ::   zcoef             ! local scalar
66      REAL(wp), POINTER, DIMENSION(:)   ::   zswitch        ! i-averaged indicator of sea-ice
67      REAL(wp), POINTER, DIMENSION(:)   ::   zmsk           ! i-averaged of tmask
68      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_io, zv_io   ! ice-ocean velocity
69      !
70      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
71     !!---------------------------------------------------------------------
72
73      IF( nn_timing == 1 )  CALL timing_start('limdyn')
74
75      CALL wrk_alloc( jpi, jpj, zu_io, zv_io )
76      CALL wrk_alloc( jpj, zswitch, zmsk )
77
78      CALL lim_var_agg(1)             ! aggregate ice categories
79
80      IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only)
81
82      IF( ln_limdyn ) THEN
83         !
84         ! conservation test
85         IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
86
87         u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1)
88         v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1)
89
90         ! Rheology (ice dynamics)
91         ! ========
92
93         !  Define the j-limits where ice rheology is computed
94         ! ---------------------------------------------------
95
96         IF( lk_mpp .OR. lk_mpp_rep ) THEN                    ! mpp: compute over the whole domain
97            i_j1 = 1
98            i_jpj = jpj
99            IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
100            CALL lim_rhg( i_j1, i_jpj )
101         ELSE                                 ! optimization of the computational area
102            !
103            DO jj = 1, jpj
104               zswitch(jj) = SUM( 1.0 - at_i(:,jj) )   ! = REAL(jpj) if ocean everywhere on a j-line
105               zmsk   (jj) = SUM( tmask(:,jj,1)    )   ! = 0         if land  everywhere on a j-line
106            END DO
107
108            IF( l_jeq ) THEN                     ! local domain include both hemisphere
109               !                                 ! Rheology is computed in each hemisphere
110               !                                 ! only over the ice cover latitude strip
111               ! Northern hemisphere
112               i_j1  = njeq
113               i_jpj = jpj
114               DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
115                  i_j1 = i_j1 + 1
116               END DO
117               i_j1 = MAX( 1, i_j1-2 )
118               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : NH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
119               CALL lim_rhg( i_j1, i_jpj )
120               !
121               ! Southern hemisphere
122               i_j1  =  1
123               i_jpj = njeq
124               DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
125                  i_jpj = i_jpj - 1
126               END DO
127               i_jpj = MIN( jpj, i_jpj+1 )
128               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : SH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
129               !
130               CALL lim_rhg( i_j1, i_jpj )
131               !
132            ELSE                                 ! local domain extends over one hemisphere only
133               !                                 ! Rheology is computed only over the ice cover
134               !                                 ! latitude strip
135               i_j1  = 1
136               DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
137                  i_j1 = i_j1 + 1
138               END DO
139               i_j1 = MAX( 1, i_j1-2 )
140
141               i_jpj  = jpj
142               DO WHILE ( i_jpj >= 1  .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
143                  i_jpj = i_jpj - 1
144               END DO
145               i_jpj = MIN( jpj, i_jpj+1)
146               !
147               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : one hemisphere:  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
148               !
149               CALL lim_rhg( i_j1, i_jpj )
150               !
151            ENDIF
152            !
153         ENDIF
154
155         ! computation of friction velocity
156         ! --------------------------------
157         ! ice-ocean velocity at U & V-points (u_ice v_ice at U- & V-points ; ssu_m, ssv_m at U- & V-points)
158         zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
159         zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
160         ! frictional velocity at T-point
161         zcoef = 0.5_wp * rn_cio
162         DO jj = 2, jpjm1 
163            DO ji = fs_2, fs_jpim1   ! vector opt.
164               ust2s(ji,jj) = zcoef * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
165                  &                    + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) * tmask(ji,jj,1)
166            END DO
167         END DO
168         !
169         ! conservation test
170         IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
171         !
172      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
173         !
174         zcoef = SQRT( 0.5_wp ) * r1_rau0
175         DO jj = 2, jpjm1
176            DO ji = fs_2, fs_jpim1   ! vector opt.
177               ust2s(ji,jj) = zcoef * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
178                  &                        + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) * tmask(ji,jj,1)
179            END DO
180         END DO
181         !
182      ENDIF
183      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
184
185      IF(ln_ctl) THEN   ! Control print
186         CALL prt_ctl_info(' ')
187         CALL prt_ctl_info(' - Cell values : ')
188         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
189         CALL prt_ctl(tab2d_1=ust2s     , clinfo1=' lim_dyn  : ust2s     :')
190         CALL prt_ctl(tab2d_1=divu_i    , clinfo1=' lim_dyn  : divu_i    :')
191         CALL prt_ctl(tab2d_1=delta_i   , clinfo1=' lim_dyn  : delta_i   :')
192         CALL prt_ctl(tab2d_1=strength  , clinfo1=' lim_dyn  : strength  :')
193         CALL prt_ctl(tab2d_1=e1e2t     , clinfo1=' lim_dyn  : cell area :')
194         CALL prt_ctl(tab2d_1=at_i      , clinfo1=' lim_dyn  : at_i      :')
195         CALL prt_ctl(tab2d_1=vt_i      , clinfo1=' lim_dyn  : vt_i      :')
196         CALL prt_ctl(tab2d_1=vt_s      , clinfo1=' lim_dyn  : vt_s      :')
197         CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' lim_dyn  : stress1_i :')
198         CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' lim_dyn  : stress2_i :')
199         CALL prt_ctl(tab2d_1=stress12_i, clinfo1=' lim_dyn  : stress12_i:')
200         DO jl = 1, jpl
201            CALL prt_ctl_info(' ')
202            CALL prt_ctl_info(' - Category : ', ivar1=jl)
203            CALL prt_ctl_info('   ~~~~~~~~~~')
204            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_dyn  : a_i      : ')
205            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_i     : ')
206            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_s     : ')
207            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_dyn  : v_i      : ')
208            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_dyn  : v_s      : ')
209            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : e_s      : ')
210            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_dyn  : t_su     : ')
211            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : t_snow   : ')
212            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_dyn  : sm_i     : ')
213            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_dyn  : smv_i    : ')
214            DO ja = 1, nlay_i
215               CALL prt_ctl_info(' ')
216               CALL prt_ctl_info(' - Layer : ', ivar1=ja)
217               CALL prt_ctl_info('   ~~~~~~~')
218               CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : t_i      : ')
219               CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : e_i      : ')
220            END DO
221         END DO
222      ENDIF
223      !
224      CALL wrk_dealloc( jpi, jpj, zu_io, zv_io )
225      CALL wrk_dealloc( jpj, zswitch, zmsk )
226      !
227      IF( nn_timing == 1 )  CALL timing_stop('limdyn')
228
229   END SUBROUTINE lim_dyn
230
231
232   SUBROUTINE lim_dyn_init
233      !!-------------------------------------------------------------------
234      !!                  ***  ROUTINE lim_dyn_init  ***
235      !!
236      !! ** Purpose : Physical constants and parameters linked to the ice
237      !!      dynamics
238      !!
239      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
240      !!       parameter values called at the first timestep (nit000)
241      !!
242      !! ** input   :   Namelist namicedyn
243      !!-------------------------------------------------------------------
244      INTEGER  ::   ios                 ! Local integer output status for namelist read
245      NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, &
246         &                nn_nevp, rn_relast, nn_ahi0, rn_ahi0_ref
247      INTEGER  ::   ji, jj
248      REAL(wp) ::   za00, zd_max
249      !!-------------------------------------------------------------------
250
251      REWIND( numnam_ice_ref )              ! Namelist namicedyn in reference namelist : Ice dynamics
252      READ  ( numnam_ice_ref, namicedyn, IOSTAT = ios, ERR = 901)
253901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in reference namelist', lwp )
254
255      REWIND( numnam_ice_cfg )              ! Namelist namicedyn in configuration namelist : Ice dynamics
256      READ  ( numnam_ice_cfg, namicedyn, IOSTAT = ios, ERR = 902 )
257902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in configuration namelist', lwp )
258      IF(lwm) WRITE ( numoni, namicedyn )
259     
260      IF(lwp) THEN                        ! control print
261         WRITE(numout,*)
262         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
263         WRITE(numout,*) '~~~~~~~~~~~~'
264         WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)  nn_icestr     = ', nn_icestr 
265         WRITE(numout,*)'    Including brine volume in ice strength comp.         ln_icestr_bvf = ', ln_icestr_bvf
266         WRITE(numout,*)'   Ratio of ridging work to PotEner change in ridging    rn_pe_rdg     = ', rn_pe_rdg 
267         WRITE(numout,*) '   drag coefficient for oceanic stress                  rn_cio        = ', rn_cio
268         WRITE(numout,*) '   first bulk-rheology parameter                        rn_pstar      = ', rn_pstar
269         WRITE(numout,*) '   second bulk-rhelogy parameter                        rn_crhg       = ', rn_crhg
270         WRITE(numout,*) '   creep limit                                          rn_creepl     = ', rn_creepl
271         WRITE(numout,*) '   eccentricity of the elliptical yield curve           rn_ecc        = ', rn_ecc
272         WRITE(numout,*) '   number of iterations for subcycling                  nn_nevp       = ', nn_nevp
273         WRITE(numout,*) '   ratio of elastic timescale over ice time step        rn_relast     = ', rn_relast
274         WRITE(numout,*) '   horizontal diffusivity calculation                   nn_ahi0       = ', nn_ahi0
275         WRITE(numout,*) '   horizontal diffusivity coeff. (orca2 grid)           rn_ahi0_ref   = ', rn_ahi0_ref
276      ENDIF
277      !
278      usecc2 = 1._wp / ( rn_ecc * rn_ecc )
279      rhoco  = rau0  * rn_cio
280      !
281      !  Diffusion coefficients
282      SELECT CASE( nn_ahi0 )
283
284      CASE( 0 )
285         ahiu(:,:) = rn_ahi0_ref
286         ahiv(:,:) = rn_ahi0_ref
287
288         IF(lwp) WRITE(numout,*) ''
289         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim constant = rn_ahi0_ref'
290
291      CASE( 1 ) 
292
293         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
294         IF( lk_mpp )   CALL mpp_max( zd_max )          ! max over the global domain
295         
296         ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp   ! 1.e05 = 100km = max grid space at 60° latitude in orca2
297                                                        !                    (60° = min latitude for ice cover) 
298         ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp
299
300         IF(lwp) WRITE(numout,*) ''
301         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')'
302         IF(lwp) WRITE(numout,*) '   value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp 
303         
304      CASE( 2 ) 
305
306         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
307         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
308         
309         za00 = rn_ahi0_ref * 1.e-05_wp          ! 1.e05 = 100km = max grid space at 60° latitude in orca2
310                                                 !                    (60° = min latitude for ice cover) 
311         DO jj = 1, jpj
312            DO ji = 1, jpi
313               ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1)
314               ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1)
315            END DO
316         END DO
317         !
318         IF(lwp) WRITE(numout,*) ''
319         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to e1'
320         IF(lwp) WRITE(numout,*) '   maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max
321         
322      END SELECT
323
324   END SUBROUTINE lim_dyn_init
325
326#else
327   !!----------------------------------------------------------------------
328   !!   Default option          Empty module           NO LIM sea-ice model
329   !!----------------------------------------------------------------------
330CONTAINS
331   SUBROUTINE lim_dyn         ! Empty routine
332   END SUBROUTINE lim_dyn
333#endif 
334
335   !!======================================================================
336END MODULE limdyn
Note: See TracBrowser for help on using the repository browser.