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

source: trunk/NEMO/LIM_SRC_3/limdyn.F90 @ 904

Last change on this file since 904 was 904, checked in by ctlod, 16 years ago

remove the sdvt(:,:) variable which is not initialized at all but is used in the computation of the friction velocity ust2s(:,:), see ticket: #121
File size: 13.9 KB
Line 
1MODULE limdyn
2   !!======================================================================
3   !!                     ***  MODULE  limdyn  ***
4   !!   Sea-Ice dynamics : 
5   !!======================================================================
6#if defined key_lim3
7   !!----------------------------------------------------------------------
8   !!   'key_lim3' :                                 LIM3 sea-ice model
9   !!----------------------------------------------------------------------
10   !!    lim_dyn      : computes ice velocities
11   !!    lim_dyn_init : initialization and namelist read
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE phycst
15   USE in_out_manager  ! I/O manager
16   USE dom_ice
17   USE dom_oce         ! ocean space and time domain
18   USE ice
19   USE par_ice
20   USE sbc_ice         ! Surface boundary condition: ice fields
21   USE ice_oce
22   USE iceini
23   USE limistate
24   USE limrhg          ! ice rheology
25   USE lbclnk
26   USE lib_mpp
27   USE prtctl          ! Print control
28
29   IMPLICIT NONE
30   PRIVATE
31
32   !! * Accessibility
33   PUBLIC lim_dyn  ! routine called by ice_step
34
35   !! * Substitutions
36#  include "vectopt_loop_substitute.h90"
37
38   !! * Module variables
39   REAL(wp)  ::  rone    = 1.e0   ! constant value
40
41   !!----------------------------------------------------------------------
42   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)
43   !! $ Id: $
44   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46
47CONTAINS
48
49   SUBROUTINE lim_dyn
50      !!-------------------------------------------------------------------
51      !!               ***  ROUTINE lim_dyn  ***
52      !!               
53      !! ** Purpose :   compute ice velocity and ocean-ice stress
54      !!               
55      !! ** Method  :
56      !!
57      !! ** Action  : - Initialisation
58      !!              - Call of the dynamic routine for each hemisphere
59      !!              - computation of the stress at the ocean surface         
60      !!              - treatment of the case if no ice dynamic
61      !! History :
62      !!   1.0  !  01-04   (LIM)  Original code
63      !!   2.0  !  02-08   (C. Ethe, G. Madec)  F90, mpp
64      !!   3.0  !  2007-03 (M.A. Morales Maqueda, S. Bouillon, M. Vancoppenolle)
65      !!                   LIM3, EVP, C-grid
66      !!------------------------------------------------------------------------------------
67      !! * Local variables
68      INTEGER ::   ji, jj, jl, ja    ! dummy loop indices
69      INTEGER :: i_j1, i_jpj         ! Starting/ending j-indices for rheology
70
71      REAL(wp) ::   &
72         ztairx, ztairy,          &  ! tempory scalars
73         zsang , zmod,            &
74         ztglx , ztgly ,          &
75         zt11, zt12, zt21, zt22 , &
76         zustm,                   &
77         zsfrldmx2, zsfrldmy2,    &
78         zu_ice, zv_ice, ztair2
79
80      REAL(wp),DIMENSION(jpj) ::   &
81           zind,                     &  ! i-averaged indicator of sea-ice
82           zmsk                         ! i-averaged of tmask
83      !!---------------------------------------------------------------------
84
85      WRITE(numout,*) ' lim_dyn : Ice dynamics '
86      WRITE(numout,*) ' ~~~~~~~ '
87
88      IF( numit == nstart  )   CALL lim_dyn_init   ! Initialization (first time-step only)
89     
90      IF ( ln_limdyn ) THEN
91
92         old_u_ice(:,:) = u_ice(:,:) * tmu(:,:)
93         old_v_ice(:,:) = v_ice(:,:) * tmv(:,:)
94
95         ! Rheology (ice dynamics)
96         ! ========
97
98         !  Define the j-limits where ice rheology is computed
99         ! ---------------------------------------------------
100
101         IF( lk_mpp .OR. nbit_cmp == 1 ) THEN                    ! mpp: compute over the whole domain
102            i_j1 = 1
103            i_jpj = jpj
104            IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
105            CALL lim_rhg( i_j1, i_jpj )
106         ELSE                                 ! optimization of the computational area
107
108            DO jj = 1, jpj
109               zind(jj) = SUM( 1.0 - at_i (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line
110               zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line
111            END DO
112
113            IF( l_jeq ) THEN                     ! local domain include both hemisphere
114               !                                 ! Rheology is computed in each hemisphere
115               !                                 ! only over the ice cover latitude strip
116               ! Northern hemisphere
117               i_j1  = njeq
118               i_jpj = jpj
119               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
120                  i_j1 = i_j1 + 1
121               END DO
122               i_j1 = MAX( 1, i_j1-1 )
123               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : NH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
124               CALL lim_rhg( i_j1, i_jpj )
125
126               ! Southern hemisphere
127               i_j1  =  1
128               i_jpj = njeq
129               DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
130                  i_jpj = i_jpj - 1
131               END DO
132               i_jpj = MIN( jpj, i_jpj+2 )
133               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : SH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
134
135               CALL lim_rhg( i_j1, i_jpj )
136
137            ELSE                                 ! local domain extends over one hemisphere only
138               !                                 ! Rheology is computed only over the ice cover
139               !                                 ! latitude strip
140               i_j1  = 1
141               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
142                  i_j1 = i_j1 + 1
143               END DO
144               i_j1 = MAX( 1, i_j1-1 )
145
146               i_jpj  = jpj
147               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
148                  i_jpj = i_jpj - 1
149               END DO
150               i_jpj = MIN( jpj, i_jpj+2)
151
152               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : one hemisphere:  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
153
154               CALL lim_rhg( i_j1, i_jpj )
155
156            ENDIF
157
158         ENDIF
159
160         ! computation of friction velocity
161         DO jj = 2, jpjm1
162            DO ji = fs_2, fs_jpim1
163
164               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
165               zt11  = rhoco * zu_ice * zu_ice
166
167               zu_ice   = u_ice(ji-1,jj) - u_oce(ji-1,jj)
168               zt12  = rhoco * zu_ice * zu_ice
169
170               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
171               zt21  = rhoco * zv_ice * zv_ice
172
173               zv_ice   = v_ice(ji,jj-1) - v_oce(ji,jj-1)
174               zt22  = rhoco * zv_ice * zv_ice
175               ztair2 = ( ( utaui_ice(ji,jj) + utaui_ice(ji-1,jj) ) / 2. )**2 + &
176                        ( ( vtaui_ice(ji,jj) + vtaui_ice(ji,jj-1) ) / 2. )**2
177
178               ! should not be weighted
179               zustm =  ( at_i(ji,jj)       ) * 0.5 * ( zt11 + zt12 + zt21 + zt22 )        &
180                  &  +  ( 1.0 - at_i(ji,jj) ) * SQRT( ztair2 )
181
182               ust2s(ji,jj) = ( zustm / rau0 ) * tms(ji,jj)
183
184            END DO
185         END DO
186
187       ELSE      ! If no ice dynamics 
188                   
189               ! virer ca (key_lim_cp1)
190          DO jj = 2, jpjm1
191             DO ji = fs_2, fs_jpim1
192                ztair2 = ( ( utaui_ice(ji,jj) + utaui_ice(ji-1,jj) ) / 2. )**2 + &
193                         ( ( vtaui_ice(ji,jj) + vtaui_ice(ji,jj-1) ) / 2. )**2
194                zustm        = SQRT( ztair2  )
195
196                ust2s(ji,jj) = ( zustm / rau0 ) * tms(ji,jj)
197            END DO
198         END DO
199
200      ENDIF
201
202      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
203
204      IF(ln_ctl) THEN   ! Control print
205         CALL prt_ctl_info(' ')
206         CALL prt_ctl_info(' - Cell values : ')
207         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
208         CALL prt_ctl(tab2d_1=ust2s     , clinfo1=' lim_dyn  : ust2s     :')
209         CALL prt_ctl(tab2d_1=divu_i    , clinfo1=' lim_dyn  : divu_i    :')
210         CALL prt_ctl(tab2d_1=delta_i   , clinfo1=' lim_dyn  : delta_i   :')
211         CALL prt_ctl(tab2d_1=strength  , clinfo1=' lim_dyn  : strength  :')
212         CALL prt_ctl(tab2d_1=area      , clinfo1=' lim_dyn  : cell area :')
213         CALL prt_ctl(tab2d_1=at_i      , clinfo1=' lim_dyn  : at_i      :')
214         CALL prt_ctl(tab2d_1=vt_i      , clinfo1=' lim_dyn  : vt_i      :')
215         CALL prt_ctl(tab2d_1=vt_s      , clinfo1=' lim_dyn  : vt_s      :')
216         CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' lim_dyn  : stress1_i :')
217         CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' lim_dyn  : stress2_i :')
218         CALL prt_ctl(tab2d_1=stress12_i, clinfo1=' lim_dyn  : stress12_i:')
219         DO jl = 1, jpl
220            CALL prt_ctl_info(' ')
221            CALL prt_ctl_info(' - Category : ', ivar1=jl)
222            CALL prt_ctl_info('   ~~~~~~~~~~')
223            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_dyn  : a_i      : ')
224            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_i     : ')
225            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_s     : ')
226            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_dyn  : v_i      : ')
227            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_dyn  : v_s      : ')
228            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : e_s      : ')
229            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_dyn  : t_su     : ')
230            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : t_snow   : ')
231            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_dyn  : sm_i     : ')
232            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_dyn  : smv_i    : ')
233            DO ja = 1, nlay_i
234               CALL prt_ctl_info(' ')
235               CALL prt_ctl_info(' - Layer : ', ivar1=ja)
236               CALL prt_ctl_info('   ~~~~~~~')
237               CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : t_i      : ')
238               CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : e_i      : ')
239            END DO
240         END DO
241      ENDIF
242
243   END SUBROUTINE lim_dyn
244
245    SUBROUTINE lim_dyn_init
246      !!-------------------------------------------------------------------
247      !!                  ***  ROUTINE lim_dyn_init  ***
248      !!
249      !! ** Purpose : Physical constants and parameters linked to the ice
250      !!      dynamics
251      !!
252      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
253      !!       parameter values called at the first timestep (nit000)
254      !!
255      !! ** input   :   Namelist namicedyn
256      !!
257      !! history :
258      !!  8.5  ! 03-08 (C. Ethe) original code
259      !!  9.0  ! 07-03 (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle)
260      !!               EVP-Cgrid-LIM3
261      !!-------------------------------------------------------------------
262      NAMELIST/namicedyn/ epsd, alpha,     &
263         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
264         &                c_rhg, etamn, creepl, ecc, ahi0, &
265         &                nevp, telast, alphaevp
266      !!-------------------------------------------------------------------
267
268      ! Define the initial parameters
269      ! -------------------------
270
271      ! Read Namelist namicedyn
272      REWIND ( numnam_ice )
273      READ   ( numnam_ice  , namicedyn )
274      IF(lwp) THEN
275         WRITE(numout,*)
276         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
277         WRITE(numout,*) '~~~~~~~~~~~~'
278         WRITE(numout,*) '   tolerance parameter                              epsd   = ', epsd
279         WRITE(numout,*) '   coefficient for semi-implicit coriolis           alpha  = ', alpha
280         WRITE(numout,*) '   diffusion constant for dynamics                  dm     = ', dm
281         WRITE(numout,*) '   number of sub-time steps for relaxation          nbiter = ', nbiter
282         WRITE(numout,*) '   maximum number of iterations for relaxation      nbitdr = ', nbitdr
283         WRITE(numout,*) '   relaxation constant                              om     = ', om
284         WRITE(numout,*) '   maximum value for the residual of relaxation     resl   = ', resl
285         WRITE(numout,*) '   drag coefficient for oceanic stress              cw     = ', cw
286         WRITE(numout,*) '   turning angle for oceanic stress                 angvg  = ', angvg
287         WRITE(numout,*) '   first bulk-rheology parameter                    pstar  = ', pstar
288         WRITE(numout,*) '   second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
289         WRITE(numout,*) '   minimun value for viscosity                      etamn  = ', etamn
290         WRITE(numout,*) '   creep limit                                      creepl = ', creepl
291         WRITE(numout,*) '   eccentricity of the elliptical yield curve       ecc    = ', ecc
292         WRITE(numout,*) '   horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
293         WRITE(numout,*) '   number of iterations for subcycling              nevp   = ', nevp
294         WRITE(numout,*) '   timescale for elastic waves                      telast = ', telast
295         WRITE(numout,*) '   coefficient for the solution of int. stresses  alphaevp = ', alphaevp
296
297      ENDIF
298
299      usecc2 = 1.0 / ( ecc * ecc )
300      rhoco  = rau0 * cw
301      angvg  = angvg * rad
302      sangvg = SIN( angvg )
303      cangvg = COS( angvg )
304      pstarh = pstar / 2.0
305
306      !  Diffusion coefficients.
307      ahiu(:,:) = ahi0 * umask(:,:,1)
308      ahiv(:,:) = ahi0 * vmask(:,:,1)
309
310   END SUBROUTINE lim_dyn_init
311
312#else
313   !!----------------------------------------------------------------------
314   !!   Default option          Empty module           NO LIM sea-ice model
315   !!----------------------------------------------------------------------
316CONTAINS
317   SUBROUTINE lim_dyn         ! Empty routine
318   END SUBROUTINE lim_dyn
319#endif 
320
321   !!======================================================================
322END MODULE limdyn
Note: See TracBrowser for help on using the repository browser.