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 @ 862

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

Clean comments and useless lines, see ticket:#72

File size: 15.0 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 taumod
19   USE ice
20   USE ice_oce
21   USE iceini
22   USE limistate
23   USE limrhg          ! ice rheology
24   USE lbclnk
25   USE lib_mpp
26   USE prtctl          ! Print control
27
28   IMPLICIT NONE
29   PRIVATE
30
31   !! * Accessibility
32   PUBLIC lim_dyn  ! routine called by ice_step
33
34   !! * Module variables
35   REAL(wp)  ::  rone    = 1.e0   ! constant value
36
37   !!----------------------------------------------------------------------
38   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)
39   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limdyn.F90,v 1.5 2005/03/27 18:34:41 opalod Exp $
40   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
41   !!----------------------------------------------------------------------
42
43CONTAINS
44
45   SUBROUTINE lim_dyn
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      !! History :
58      !!   1.0  !  01-04   (LIM)  Original code
59      !!   2.0  !  02-08   (C. Ethe, G. Madec)  F90, mpp
60      !!   3.0  !  2007-03 (M.A. Morales Maqueda, S. Bouillon, M. Vancoppenolle)
61      !!                   LIM3, EVP, C-grid
62      !!------------------------------------------------------------------------------------
63      !! * Local variables
64      INTEGER ::   ji, jj            ! dummy loop indices
65      INTEGER :: i_j1, i_jpj         ! Starting/ending j-indices for rheology
66
67      REAL(wp) ::   &
68         ztairx, ztairy,          &  ! tempory scalars
69         zsang , zmod,            &
70         ztglx , ztgly ,          &
71         zt11, zt12, zt21, zt22 , &
72         zustm,                   &
73         zsfrldmx2, zsfrldmy2,    &
74         zu_ice, zv_ice, ztair2
75
76      REAL(wp),DIMENSION(jpj) ::   &
77           zind,                     &  ! i-averaged indicator of sea-ice
78           zmsk                         ! i-averaged of tmask
79      !!---------------------------------------------------------------------
80
81      WRITE(numout,*) ' lim_dyn : Ice dynamics '
82      WRITE(numout,*) ' ~~~~~~~ '
83
84      IF( numit == nstart  )   CALL lim_dyn_init   ! Initialization (first time-step only)
85     
86      IF ( ln_limdyn ) THEN
87
88         ! ocean velocity
89         u_oce(:,:)  = u_io(:,:) * tmu(:,:)
90         v_oce(:,:)  = v_io(:,:) * tmv(:,:)
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 ) 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         u_ice(:,1) = 0.0       !ibug  est-ce vraiment necessaire?
161         v_ice(:,1) = 0.0
162
163         IF(ln_ctl) THEN
164            CALL prt_ctl(tab2d_1=u_oce , clinfo1=' lim_dyn  : u_oce :', tab2d_2=v_oce , clinfo2=' v_oce :')
165            CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :')
166         ENDIF
167
168         ! Ice-Ocean stress
169         ! ================
170         DO jj = 2, jpjm1
171            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg
172
173            DO ji = 2, jpim1
174               ! computation of wind stress over ocean in X and Y direction
175#if defined key_coupled && defined key_lim_cp1
176!              ztairx =  ( 1.0 - at_i(ji-1,jj)   ) * gtaux(ji-1,jj)   + &
177!                        ( 1.0 - at_i(ji,jj)     ) * gtaux(ji,jj  )   + &
178!                        ( 1.0 - at_i(ji-1,jj-1) ) * gtaux(ji-1,jj-1) + &
179!                        ( 1.0 - at_i(ji,jj-1)   ) * gtaux(ji,jj-1)
180
181!              ztairy =  ( 1.0 - at_i(ji-1,jj)   ) * gtauy(ji-1,jj  ) + &
182!                        ( 1.0 - at_i(ji,jj  )   ) * gtauy(ji,jj    ) + &
183!                        ( 1.0 - at_i(ji-1,jj-1) ) * gtauy(ji-1,jj-1) + &
184!                        ( 1.0 - at_i(ji,jj-1)   ) * gtauy(ji,jj-1)
185#else
186               ztairx =  ( 2.0 - at_i(ji,jj) - at_i(ji+1,jj) ) * gtaux(ji,jj) / cai * cao
187               ztairy =  ( 2.0 - at_i(ji,jj) - at_i(ji,jj+1) ) * gtauy(ji,jj) / cai * cao
188
189               zsfrldmx2 = at_i(ji,jj) + at_i(ji+1,jj)
190               zsfrldmy2 = at_i(ji,jj) + at_i(ji,jj+1)
191
192#endif
193               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
194               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
195               zmod     = SQRT( zu_ice * zu_ice + zv_ice * zv_ice ) 
196
197               ! quadratic drag formulation
198               ztglx   = zsfrldmx2 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice ) 
199               ztgly   = zsfrldmy2 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice ) 
200!
201!              ! IMPORTANT
202!              ! these lignes are bound to prevent numerical oscillations
203!              ! in the ice-ocean stress
204!              ! They are physically ill-based. There is a cleaner solution
205!              ! to try (remember discussion in Paris Gurvan)
206!
207               ztglx   = ztglx * exp( - zmod / 0.5 )
208               ztgly   = ztglx * exp( - zmod / 0.5 )
209
210               tio_u(ji,jj) = - ( ztairx + 1.0 * ztglx ) / ( 2. * rau0 )
211               tio_v(ji,jj) = - ( ztairy + 1.0 * ztgly ) / ( 2. * rau0 )
212            END DO
213         END DO
214         
215         ! computation of friction velocity
216         DO jj = 2, jpjm1
217            DO ji = 2, jpim1
218
219               zu_ice   = u_ice(ji,jj) - u_io(ji,jj)
220               zt11  = rhoco * zu_ice * zu_ice
221
222               zu_ice   = u_ice(ji-1,jj) - u_io(ji-1,jj)
223               zt12  = rhoco * zu_ice * zu_ice
224
225               zv_ice   = v_ice(ji,jj) - v_io(ji,jj)
226               zt21  = rhoco * zv_ice * zv_ice
227
228               zv_ice   = v_ice(ji,jj-1) - v_io(ji,jj-1)
229               zt22  = rhoco * zv_ice * zv_ice
230               ztair2 = ( ( gtaux(ji,jj) + gtaux(ji-1,jj) ) / 2. )**2 + &
231                        ( ( gtauy(ji,jj) + gtauy(ji,jj-1) ) / 2. )**2
232
233               ! should not be weighted
234               zustm =  ( at_i(ji,jj)       ) * 0.5 * ( zt11 + zt12 + zt21 + zt22 )        &
235                  &  +  ( 1.0 - at_i(ji,jj) ) * SQRT( ztair2 )
236
237               ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj)
238
239            END DO
240         END DO
241
242       ELSE      ! If no ice dynamics 
243                   
244               ! virer ca (key_lim_cp1)
245          DO jj = 2, jpjm1
246             DO ji = 2, jpim1
247#if defined key_coupled && defined key_lim_cp1
248                tio_u(ji,jj) = - (  gtaux(ji  ,jj  ) + gtaux(ji-1,jj  )       &
249                   &              + gtaux(ji-1,jj-1) + gtaux(ji  ,jj-1) ) / ( 4 * rau0 )
250
251                tio_v(ji,jj) = - (  gtauy(ji  ,jj )  + gtauy(ji-1,jj  )       &
252                   &              + gtauy(ji-1,jj-1) + gtauy(ji  ,jj-1) ) / ( 4 * rau0 )
253#else
254                tio_u(ji,jj) = - gtaux(ji,jj) / cai * cao / rau0
255                tio_v(ji,jj) = - gtauy(ji,jj) / cai * cao / rau0 
256#endif
257                ztair2 = ( ( gtaux(ji,jj) + gtaux(ji-1,jj) ) / 2. )**2 + &
258                         ( ( gtauy(ji,jj) + gtauy(ji,jj-1) ) / 2. )**2
259                zustm        = SQRT( ztair2  )
260
261                ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj)
262            END DO
263         END DO
264
265      ENDIF
266
267      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
268      CALL lbc_lnk( tio_u, 'U', -1. )   ! I-point (i.e. ice U-V point)
269      CALL lbc_lnk( tio_v, 'V', -1. )   ! I-point (i.e. ice U-V point)
270
271      IF(ln_ctl) THEN
272        CALL prt_ctl(tab2d_1=tio_u , clinfo1=' lim_dyn  : tio_u :', tab2d_2=tio_v , clinfo2=' tio_v :')
273        CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :')
274        CALL prt_ctl(tab2d_1=at_i  , clinfo1=' lim_dyn  : at_i  :')
275      ENDIF
276
277   END SUBROUTINE lim_dyn
278
279    SUBROUTINE lim_dyn_init
280      !!-------------------------------------------------------------------
281      !!                  ***  ROUTINE lim_dyn_init  ***
282      !!
283      !! ** Purpose : Physical constants and parameters linked to the ice
284      !!      dynamics
285      !!
286      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
287      !!       parameter values called at the first timestep (nit000)
288      !!
289      !! ** input   :   Namelist namicedyn
290      !!
291      !! history :
292      !!  8.5  ! 03-08 (C. Ethe) original code
293      !!  9.0  ! 07-03 (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle)
294      !!               EVP-Cgrid-LIM3
295      !!-------------------------------------------------------------------
296      NAMELIST/namicedyn/ epsd, alpha,     &
297         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
298         &                c_rhg, etamn, creepl, ecc, ahi0, &
299         &                nevp, telast, alphaevp
300      !!-------------------------------------------------------------------
301
302      ! Define the initial parameters
303      ! -------------------------
304
305      ! Read Namelist namicedyn
306      REWIND ( numnam_ice )
307      READ   ( numnam_ice  , namicedyn )
308      IF(lwp) THEN
309         WRITE(numout,*)
310         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
311         WRITE(numout,*) '~~~~~~~~~~~~'
312         WRITE(numout,*) '   tolerance parameter                              epsd   = ', epsd
313         WRITE(numout,*) '   coefficient for semi-implicit coriolis           alpha  = ', alpha
314         WRITE(numout,*) '   diffusion constant for dynamics                  dm     = ', dm
315         WRITE(numout,*) '   number of sub-time steps for relaxation          nbiter = ', nbiter
316         WRITE(numout,*) '   maximum number of iterations for relaxation      nbitdr = ', nbitdr
317         WRITE(numout,*) '   relaxation constant                              om     = ', om
318         WRITE(numout,*) '   maximum value for the residual of relaxation     resl   = ', resl
319         WRITE(numout,*) '   drag coefficient for oceanic stress              cw     = ', cw
320         WRITE(numout,*) '   turning angle for oceanic stress                 angvg  = ', angvg
321         WRITE(numout,*) '   first bulk-rheology parameter                    pstar  = ', pstar
322         WRITE(numout,*) '   second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
323         WRITE(numout,*) '   minimun value for viscosity                      etamn  = ', etamn
324         WRITE(numout,*) '   creep limit                                      creepl = ', creepl
325         WRITE(numout,*) '   eccentricity of the elliptical yield curve       ecc    = ', ecc
326         WRITE(numout,*) '   horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
327         WRITE(numout,*) '   number of iterations for subcycling              nevp   = ', nevp
328         WRITE(numout,*) '   timescale for elastic waves                      telast = ', telast
329         WRITE(numout,*) '   coefficient for the solution of int. stresses  alphaevp = ', alphaevp
330
331      ENDIF
332
333      usecc2 = 1.0 / ( ecc * ecc )
334      rhoco  = rau0 * cw
335      angvg  = angvg * rad
336      sangvg = SIN( angvg )
337      cangvg = COS( angvg )
338      pstarh = pstar / 2.0
339
340      !  Diffusion coefficients.
341      ahiu(:,:) = ahi0 * umask(:,:,1)
342      ahiv(:,:) = ahi0 * vmask(:,:,1)
343
344   END SUBROUTINE lim_dyn_init
345
346#else
347   !!----------------------------------------------------------------------
348   !!   Default option          Empty module           NO LIM sea-ice model
349   !!----------------------------------------------------------------------
350CONTAINS
351   SUBROUTINE lim_dyn         ! Empty routine
352   END SUBROUTINE lim_dyn
353#endif 
354
355   !!======================================================================
356END MODULE limdyn
Note: See TracBrowser for help on using the repository browser.