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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90 @ 2370

Last change on this file since 2370 was 2370, checked in by gm, 13 years ago

v3.3beta: ice-ocean stress at kt with VP & EVP (LIM-2 and -3)

  • Property svn:keywords set to Id
File size: 13.5 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   !!----------------------------------------------------------------------
9#if defined key_lim3
10   !!----------------------------------------------------------------------
11   !!   'key_lim3' :                                 LIM3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!    lim_dyn      : computes ice velocities
14   !!    lim_dyn_init : initialization and namelist read
15   !!----------------------------------------------------------------------
16   USE phycst
17   USE in_out_manager  ! I/O manager
18   USE dom_ice
19   USE dom_oce         ! ocean space and time domain
20   USE ice
21   USE par_ice
22   USE sbc_oce         ! Surface boundary condition: ocean fields
23   USE sbc_ice         ! Surface boundary condition: ice fields
24   USE iceini
25   USE limistate
26   USE limrhg          ! ice rheology
27   USE lbclnk
28   USE lib_mpp
29   USE prtctl          ! Print control
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   lim_dyn   ! routine called by ice_step
35
36   !! * Substitutions
37#  include "vectopt_loop_substitute.h90"
38   !!----------------------------------------------------------------------
39   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
40   !! $Id$
41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE lim_dyn( kt )
46      !!-------------------------------------------------------------------
47      !!               ***  ROUTINE lim_dyn  ***
48      !!               
49      !! ** Purpose :   compute ice velocity and ocean-ice stress
50      !!               
51      !! ** Method  :
52      !!
53      !! ** Action  : - Initialisation
54      !!              - Call of the dynamic routine for each hemisphere
55      !!              - computation of the stress at the ocean surface         
56      !!              - treatment of the case if no ice dynamic
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) ::   zcoef             ! local scalar
63      REAL(wp), DIMENSION(jpj)     ::   zind           ! i-averaged indicator of sea-ice
64      REAL(wp), DIMENSION(jpj)     ::   zmsk           ! i-averaged of tmask
65      REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io   ! ice-ocean velocity
66      !!---------------------------------------------------------------------
67
68      IF( kt == nit000 .AND. lwp ) THEN
69         WRITE(numout,*) ' lim_dyn : Ice dynamics '
70         WRITE(numout,*) ' ~~~~~~~ '
71      ENDIF
72
73      IF( numit == nstart  )   CALL lim_dyn_init   ! Initialization (first time-step only)
74
75      IF ( ln_limdyn ) THEN
76
77         old_u_ice(:,:) = u_ice(:,:) * tmu(:,:)
78         old_v_ice(:,:) = v_ice(:,:) * tmv(:,:)
79
80         ! Rheology (ice dynamics)
81         ! ========
82
83         !  Define the j-limits where ice rheology is computed
84         ! ---------------------------------------------------
85
86         IF( lk_mpp .OR. lk_mpp_rep ) THEN                    ! mpp: compute over the whole domain
87            i_j1 = 1
88            i_jpj = jpj
89            IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
90            CALL lim_rhg( i_j1, i_jpj )
91         ELSE                                 ! optimization of the computational area
92
93            DO jj = 1, jpj
94               zind(jj) = SUM( 1.0 - at_i (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line
95               zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line
96            END DO
97
98            IF( l_jeq ) THEN                     ! local domain include both hemisphere
99               !                                 ! Rheology is computed in each hemisphere
100               !                                 ! only over the ice cover latitude strip
101               ! Northern hemisphere
102               i_j1  = njeq
103               i_jpj = jpj
104               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
105                  i_j1 = i_j1 + 1
106               END DO
107               i_j1 = MAX( 1, i_j1-2 )
108               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : NH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
109               CALL lim_rhg( i_j1, i_jpj )
110
111               ! Southern hemisphere
112               i_j1  =  1
113               i_jpj = njeq
114               DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
115                  i_jpj = i_jpj - 1
116               END DO
117               i_jpj = MIN( jpj, i_jpj+1 )
118               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : SH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
119
120       CALL lim_rhg( i_j1, i_jpj )
121
122    ELSE                                 ! local domain extends over one hemisphere only
123       !                                 ! Rheology is computed only over the ice cover
124       !                                 ! latitude strip
125       i_j1  = 1
126               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
127                  i_j1 = i_j1 + 1
128               END DO
129               i_j1 = MAX( 1, i_j1-2 )
130
131               i_jpj  = jpj
132               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
133                  i_jpj = i_jpj - 1
134               END DO
135               i_jpj = MIN( jpj, i_jpj+1)
136
137               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : one hemisphere:  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
138
139               CALL lim_rhg( i_j1, i_jpj )
140
141            ENDIF
142
143         ENDIF
144
145         ! computation of friction velocity
146         ! --------------------------------
147         ! ice-ocean velocity at U & V-points (u_ice v_ice at U- & V-points ; ssu_m, ssv_m at U- & V-points)
148         zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
149         zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
150         ! frictional velocity at T-point
151         zcoef = 0.5 * cw
152         DO jj = 2, jpjm1 
153            DO ji = fs_2, fs_jpim1   ! vector opt.
154               ust2s(ji,jj) = zcoef * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
155                  &                    + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj)
156            END DO
157         END DO
158         !
159      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
160         !
161         zcoef = SQRT( 0.5 ) / rau0
162         DO jj = 2, jpjm1
163            DO ji = fs_2, fs_jpim1   ! vector opt.
164               ust2s(ji,jj) = zcoef * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
165                  &                        + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1)   ) * tms(ji,jj)
166            END DO
167         END DO
168         !
169      ENDIF
170      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
171
172      IF(ln_ctl) THEN   ! Control print
173         CALL prt_ctl_info(' ')
174         CALL prt_ctl_info(' - Cell values : ')
175         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
176         CALL prt_ctl(tab2d_1=ust2s     , clinfo1=' lim_dyn  : ust2s     :')
177         CALL prt_ctl(tab2d_1=divu_i    , clinfo1=' lim_dyn  : divu_i    :')
178         CALL prt_ctl(tab2d_1=delta_i   , clinfo1=' lim_dyn  : delta_i   :')
179         CALL prt_ctl(tab2d_1=strength  , clinfo1=' lim_dyn  : strength  :')
180         CALL prt_ctl(tab2d_1=area      , clinfo1=' lim_dyn  : cell area :')
181         CALL prt_ctl(tab2d_1=at_i      , clinfo1=' lim_dyn  : at_i      :')
182         CALL prt_ctl(tab2d_1=vt_i      , clinfo1=' lim_dyn  : vt_i      :')
183         CALL prt_ctl(tab2d_1=vt_s      , clinfo1=' lim_dyn  : vt_s      :')
184         CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' lim_dyn  : stress1_i :')
185         CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' lim_dyn  : stress2_i :')
186         CALL prt_ctl(tab2d_1=stress12_i, clinfo1=' lim_dyn  : stress12_i:')
187         DO jl = 1, jpl
188            CALL prt_ctl_info(' ')
189            CALL prt_ctl_info(' - Category : ', ivar1=jl)
190            CALL prt_ctl_info('   ~~~~~~~~~~')
191            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_dyn  : a_i      : ')
192            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_i     : ')
193            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_s     : ')
194            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_dyn  : v_i      : ')
195            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_dyn  : v_s      : ')
196            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : e_s      : ')
197            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_dyn  : t_su     : ')
198            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : t_snow   : ')
199            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_dyn  : sm_i     : ')
200            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_dyn  : smv_i    : ')
201            DO ja = 1, nlay_i
202               CALL prt_ctl_info(' ')
203               CALL prt_ctl_info(' - Layer : ', ivar1=ja)
204               CALL prt_ctl_info('   ~~~~~~~')
205               CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : t_i      : ')
206               CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : e_i      : ')
207            END DO
208         END DO
209      ENDIF
210      !
211   END SUBROUTINE lim_dyn
212
213
214   SUBROUTINE lim_dyn_init
215      !!-------------------------------------------------------------------
216      !!                  ***  ROUTINE lim_dyn_init  ***
217      !!
218      !! ** Purpose : Physical constants and parameters linked to the ice
219      !!      dynamics
220      !!
221      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
222      !!       parameter values called at the first timestep (nit000)
223      !!
224      !! ** input   :   Namelist namicedyn
225      !!-------------------------------------------------------------------
226      NAMELIST/namicedyn/ epsd, alpha,     &
227         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
228         &                c_rhg, etamn, creepl, ecc, ahi0, &
229         &                nevp, telast, alphaevp
230      !!-------------------------------------------------------------------
231
232      REWIND( numnam_ice )                ! Read Namelist namicedyn
233      READ  ( numnam_ice  , namicedyn )
234     
235      IF(lwp) THEN                        ! control print
236         WRITE(numout,*)
237         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
238         WRITE(numout,*) '~~~~~~~~~~~~'
239         WRITE(numout,*) '   tolerance parameter                              epsd   = ', epsd
240         WRITE(numout,*) '   coefficient for semi-implicit coriolis           alpha  = ', alpha
241         WRITE(numout,*) '   diffusion constant for dynamics                  dm     = ', dm
242         WRITE(numout,*) '   number of sub-time steps for relaxation          nbiter = ', nbiter
243         WRITE(numout,*) '   maximum number of iterations for relaxation      nbitdr = ', nbitdr
244         WRITE(numout,*) '   relaxation constant                              om     = ', om
245         WRITE(numout,*) '   maximum value for the residual of relaxation     resl   = ', resl
246         WRITE(numout,*) '   drag coefficient for oceanic stress              cw     = ', cw
247         WRITE(numout,*) '   turning angle for oceanic stress                 angvg  = ', angvg
248         WRITE(numout,*) '   first bulk-rheology parameter                    pstar  = ', pstar
249         WRITE(numout,*) '   second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
250         WRITE(numout,*) '   minimun value for viscosity                      etamn  = ', etamn
251         WRITE(numout,*) '   creep limit                                      creepl = ', creepl
252         WRITE(numout,*) '   eccentricity of the elliptical yield curve       ecc    = ', ecc
253         WRITE(numout,*) '   horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
254         WRITE(numout,*) '   number of iterations for subcycling              nevp   = ', nevp
255         WRITE(numout,*) '   timescale for elastic waves                      telast = ', telast
256         WRITE(numout,*) '   coefficient for the solution of int. stresses  alphaevp = ', alphaevp
257      ENDIF
258      !
259      IF( angvg /= 0._wp ) THEN
260         CALL ctl_warn( 'lim_dyn_init: turning angle for oceanic stress not properly coded for EVP ',   &
261            &           '(see limsbc module). We force  angvg = 0._wp'  )
262         angvg = 0._wp
263      ENDIF
264     
265      usecc2 = 1._wp / ( ecc * ecc )
266      rhoco  = rau0  * cw
267      angvg  = angvg * rad
268      sangvg = SIN( angvg )
269      cangvg = COS( angvg )
270      pstarh = pstar * 0.5_wp
271
272      !  Diffusion coefficients.
273      ahiu(:,:) = ahi0 * umask(:,:,1)
274      ahiv(:,:) = ahi0 * vmask(:,:,1)
275
276   END SUBROUTINE lim_dyn_init
277
278#else
279   !!----------------------------------------------------------------------
280   !!   Default option          Empty module           NO LIM sea-ice model
281   !!----------------------------------------------------------------------
282CONTAINS
283   SUBROUTINE lim_dyn         ! Empty routine
284   END SUBROUTINE lim_dyn
285#endif 
286
287   !!======================================================================
288END MODULE limdyn
Note: See TracBrowser for help on using the repository browser.