source: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90 @ 5086

Last change on this file since 5086 was 5086, checked in by timgraham, 7 years ago

Merged head of trunk into branch in preparation for putting code back onto the trunk
In working copy ran the command:
svn merge svn+sshtimgraham@…/ipsl/forge/projets/nemo/svn/trunk

Also recompiled NEMO_book.pdf with merged input files

  • Property svn:keywords set to Id
File size: 14.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   !!            4.0  ! 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 par_ice          ! LIM-3 parameters
23   USE dom_ice          ! LIM-3 domain
24   USE limrhg           ! LIM-3 rheology
25   USE lbclnk           ! lateral boundary conditions - MPP exchanges
26   USE lib_mpp          ! MPP library
27   USE wrk_nemo         ! work arrays
28   USE in_out_manager   ! I/O manager
29   USE prtctl           ! Print control
30   USE lib_fortran      ! glob_sum
31   USE timing          ! Timing
32   USE limcons        ! conservation tests
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      IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only)
79
80      IF( ln_limdyn ) THEN
81         !
82         ! conservation test
83         IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
84
85         u_ice_b(:,:) = u_ice(:,:) * tmu(:,:)
86         v_ice_b(:,:) = v_ice(:,:) * tmv(:,:)
87
88         ! Rheology (ice dynamics)
89         ! ========
90
91         !  Define the j-limits where ice rheology is computed
92         ! ---------------------------------------------------
93
94         IF( lk_mpp .OR. lk_mpp_rep ) THEN                    ! mpp: compute over the whole domain
95            i_j1 = 1
96            i_jpj = jpj
97            IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
98            CALL lim_rhg( i_j1, i_jpj )
99         ELSE                                 ! optimization of the computational area
100            !
101            DO jj = 1, jpj
102               zswitch(jj) = SUM( 1.0 - at_i(:,jj) )   ! = REAL(jpj) if ocean everywhere on a j-line
103               zmsk(jj) = SUM( tmask(:,jj,1)    )   ! = 0         if land  everywhere on a j-line
104            END DO
105
106            IF( l_jeq ) THEN                     ! local domain include both hemisphere
107               !                                 ! Rheology is computed in each hemisphere
108               !                                 ! only over the ice cover latitude strip
109               ! Northern hemisphere
110               i_j1  = njeq
111               i_jpj = jpj
112               DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
113                  i_j1 = i_j1 + 1
114               END DO
115               i_j1 = MAX( 1, i_j1-2 )
116               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : NH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
117               CALL lim_rhg( i_j1, i_jpj )
118               !
119               ! Southern hemisphere
120               i_j1  =  1
121               i_jpj = njeq
122               DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
123                  i_jpj = i_jpj - 1
124               END DO
125               i_jpj = MIN( jpj, i_jpj+1 )
126               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : SH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
127               !
128               CALL lim_rhg( i_j1, i_jpj )
129               !
130            ELSE                                 ! local domain extends over one hemisphere only
131               !                                 ! Rheology is computed only over the ice cover
132               !                                 ! latitude strip
133               i_j1  = 1
134               DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
135                  i_j1 = i_j1 + 1
136               END DO
137               i_j1 = MAX( 1, i_j1-2 )
138
139               i_jpj  = jpj
140               DO WHILE ( i_jpj >= 1  .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
141                  i_jpj = i_jpj - 1
142               END DO
143               i_jpj = MIN( jpj, i_jpj+1)
144               !
145               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : one hemisphere:  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
146               !
147               CALL lim_rhg( i_j1, i_jpj )
148               !
149            ENDIF
150            !
151         ENDIF
152
153         ! computation of friction velocity
154         ! --------------------------------
155         ! ice-ocean velocity at U & V-points (u_ice v_ice at U- & V-points ; ssu_m, ssv_m at U- & V-points)
156         zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
157         zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
158         ! frictional velocity at T-point
159         zcoef = 0.5_wp * cw
160         DO jj = 2, jpjm1 
161            DO ji = fs_2, fs_jpim1   ! vector opt.
162               ust2s(ji,jj) = zcoef * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
163                  &                    + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj)
164            END DO
165         END DO
166         !
167         ! conservation test
168         IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
169         !
170      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
171         !
172         zcoef = SQRT( 0.5_wp ) / rau0
173         DO jj = 2, jpjm1
174            DO ji = fs_2, fs_jpim1   ! vector opt.
175               ust2s(ji,jj) = zcoef * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
176                  &                        + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1)   ) * tms(ji,jj)
177            END DO
178         END DO
179         !
180      ENDIF
181      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
182
183      IF(ln_ctl) THEN   ! Control print
184         CALL prt_ctl_info(' ')
185         CALL prt_ctl_info(' - Cell values : ')
186         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
187         CALL prt_ctl(tab2d_1=ust2s     , clinfo1=' lim_dyn  : ust2s     :')
188         CALL prt_ctl(tab2d_1=divu_i    , clinfo1=' lim_dyn  : divu_i    :')
189         CALL prt_ctl(tab2d_1=delta_i   , clinfo1=' lim_dyn  : delta_i   :')
190         CALL prt_ctl(tab2d_1=strength  , clinfo1=' lim_dyn  : strength  :')
191         CALL prt_ctl(tab2d_1=area      , clinfo1=' lim_dyn  : cell area :')
192         CALL prt_ctl(tab2d_1=at_i      , clinfo1=' lim_dyn  : at_i      :')
193         CALL prt_ctl(tab2d_1=vt_i      , clinfo1=' lim_dyn  : vt_i      :')
194         CALL prt_ctl(tab2d_1=vt_s      , clinfo1=' lim_dyn  : vt_s      :')
195         CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' lim_dyn  : stress1_i :')
196         CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' lim_dyn  : stress2_i :')
197         CALL prt_ctl(tab2d_1=stress12_i, clinfo1=' lim_dyn  : stress12_i:')
198         DO jl = 1, jpl
199            CALL prt_ctl_info(' ')
200            CALL prt_ctl_info(' - Category : ', ivar1=jl)
201            CALL prt_ctl_info('   ~~~~~~~~~~')
202            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_dyn  : a_i      : ')
203            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_i     : ')
204            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_s     : ')
205            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_dyn  : v_i      : ')
206            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_dyn  : v_s      : ')
207            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : e_s      : ')
208            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_dyn  : t_su     : ')
209            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : t_snow   : ')
210            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_dyn  : sm_i     : ')
211            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_dyn  : smv_i    : ')
212            DO ja = 1, nlay_i
213               CALL prt_ctl_info(' ')
214               CALL prt_ctl_info(' - Layer : ', ivar1=ja)
215               CALL prt_ctl_info('   ~~~~~~~')
216               CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : t_i      : ')
217               CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : e_i      : ')
218            END DO
219         END DO
220      ENDIF
221      !
222      CALL wrk_dealloc( jpi, jpj, zu_io, zv_io )
223      CALL wrk_dealloc( jpj, zswitch, zmsk )
224      !
225      IF( nn_timing == 1 )  CALL timing_stop('limdyn')
226
227   END SUBROUTINE lim_dyn
228
229
230   SUBROUTINE lim_dyn_init
231      !!-------------------------------------------------------------------
232      !!                  ***  ROUTINE lim_dyn_init  ***
233      !!
234      !! ** Purpose : Physical constants and parameters linked to the ice
235      !!      dynamics
236      !!
237      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
238      !!       parameter values called at the first timestep (nit000)
239      !!
240      !! ** input   :   Namelist namicedyn
241      !!-------------------------------------------------------------------
242      INTEGER  ::   ios                 ! Local integer output status for namelist read
243      NAMELIST/namicedyn/ epsd, om, cw, pstar,   &
244         &                c_rhg, creepl, ecc, ahi0,     &
245         &                nevp, relast, alphaevp, hminrhg
246      !!-------------------------------------------------------------------
247
248      REWIND( numnam_ice_ref )              ! Namelist namicedyn in reference namelist : Ice dynamics
249      READ  ( numnam_ice_ref, namicedyn, IOSTAT = ios, ERR = 901)
250901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in reference namelist', lwp )
251
252      REWIND( numnam_ice_cfg )              ! Namelist namicedyn in configuration namelist : Ice dynamics
253      READ  ( numnam_ice_cfg, namicedyn, IOSTAT = ios, ERR = 902 )
254902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in configuration namelist', lwp )
255      IF(lwm) WRITE ( numoni, namicedyn )
256     
257      IF(lwp) THEN                        ! control print
258         WRITE(numout,*)
259         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
260         WRITE(numout,*) '~~~~~~~~~~~~'
261         WRITE(numout,*) '   tolerance parameter                              epsd   = ', epsd
262         WRITE(numout,*) '   relaxation constant                              om     = ', om
263         WRITE(numout,*) '   drag coefficient for oceanic stress              cw     = ', cw
264         WRITE(numout,*) '   first bulk-rheology parameter                    pstar  = ', pstar
265         WRITE(numout,*) '   second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
266         WRITE(numout,*) '   creep limit                                      creepl = ', creepl
267         WRITE(numout,*) '   eccentricity of the elliptical yield curve       ecc    = ', ecc
268         WRITE(numout,*) '   horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
269         WRITE(numout,*) '   number of iterations for subcycling              nevp   = ', nevp
270         WRITE(numout,*) '   ratio of elastic timescale over ice time step    relast = ', relast
271         WRITE(numout,*) '   coefficient for the solution of int. stresses  alphaevp = ', alphaevp
272         WRITE(numout,*) '   min ice thickness for rheology calculations     hminrhg = ', hminrhg
273      ENDIF
274      !
275      usecc2 = 1._wp / ( ecc * ecc )
276      rhoco  = rau0  * cw
277
278      ! elastic damping
279      telast = relast * rdt_ice
280
281      !  Diffusion coefficients.
282      ahiu(:,:) = ahi0 * umask(:,:,1)
283      ahiv(:,:) = ahi0 * vmask(:,:,1)
284      !
285   END SUBROUTINE lim_dyn_init
286
287#else
288   !!----------------------------------------------------------------------
289   !!   Default option          Empty module           NO LIM sea-ice model
290   !!----------------------------------------------------------------------
291CONTAINS
292   SUBROUTINE lim_dyn         ! Empty routine
293   END SUBROUTINE lim_dyn
294#endif 
295
296   !!======================================================================
297END MODULE limdyn
Note: See TracBrowser for help on using the repository browser.