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

source: branches/dev_002_LIM/NEMO/LIM_SRC_3/limdyn.F90 @ 830

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

dev_002_LIM : add the LIM 3.0 component, see ticketr: #71

File size: 15.3 KB
Line 
1MODULE limdyn
2   !!======================================================================
3   !!                     ***  MODULE  limdyn  ***
4   !!   Sea-Ice dynamics : 
5   !!======================================================================
6#if defined key_lim3
7   !!----------------------------------------------------------------------
8   !!   'key_lim3' :                                   LIM 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 2.0,  UCL-LOCEAN-IPSL (2005)
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) LIM3, EVP, C-grid
61      !!------------------------------------------------------------------------------------
62      !! * Local variables
63      INTEGER ::   ji, jj, jl        ! dummy loop indices
64      INTEGER :: i_j1, i_jpj         ! Starting/ending j-indices for rheology
65! nemo modif
66!        jhemis                      ! jhemis = 1 (NH) ; jhemis = -1 (SH)
67      REAL(wp) ::   &
68         ztairx, ztairy,          &  ! tempory scalars
69         zsang , zmod,            &
70         ztglx , ztgly ,          &
71         zt11, zt12, zt21, zt22 , &
72         zustm, zsfrld, zsfrldm4, &
73         zsfrldmx2, zsfrldmy2,    &
74         zu_ice, zv_ice, ztair2
75! nemo modif
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         ! Mean ice and snow thicknesses.         
89
90         u_oce(:,:)  = u_io(:,:) * tmu(:,:)
91         v_oce(:,:)  = v_io(:,:) * tmv(:,:)
92         
93         old_u_ice(:,:) = u_ice(:,:) * tmu(:,:)
94         old_v_ice(:,:) = v_ice(:,:) * tmv(:,:)
95
96         !                                         ! Rheology (ice dynamics)
97         !                                         ! ========
98
99         !  Define the j-limits where ice rheology is computed
100         ! ---------------------------------------------------
101
102         IF( lk_mpp ) THEN                    ! mpp: compute over the whole domain
103            i_j1 = 1
104            i_jpj = jpj
105            IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
106            CALL lim_rhg( i_j1, i_jpj )
107
108         ELSE                                 ! optimization of the computational area
109
110            DO jj = 1, jpj
111               zind(jj) = SUM( 1.0 - at_i (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line
112               zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line
113            END DO
114
115            IF( l_jeq ) THEN                     ! local domain include both hemisphere
116               !                                 ! Rheology is computed in each hemisphere
117               !                                 ! only over the ice cover latitude strip
118               ! Northern hemisphere
119               i_j1  = njeq
120               i_jpj = jpj
121               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
122                  i_j1 = i_j1 + 1
123               END DO
124               i_j1 = MAX( 1, i_j1-1 )
125               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : NH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
126
127               CALL lim_rhg( i_j1, i_jpj )
128
129               ! Southern hemisphere
130               i_j1  =  1
131               i_jpj = njeq
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+2 )
136               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : SH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
137
138               CALL lim_rhg( i_j1, i_jpj )
139
140            ELSE                                 ! local domain extends over one hemisphere only
141               !                                 ! Rheology is computed only over the ice cover
142               !                                 ! latitude strip
143               i_j1  = 1
144               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
145                  i_j1 = i_j1 + 1
146               END DO
147               i_j1 = MAX( 1, i_j1-1 )
148
149               i_jpj  = jpj
150               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
151                  i_jpj = i_jpj - 1
152               END DO
153               i_jpj = MIN( jpj, i_jpj+2)
154
155               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : one hemisphere:  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
156
157               CALL lim_rhg( i_j1, i_jpj )
158
159            ENDIF
160
161         ENDIF
162
163         u_ice(:,1) = 0.0       !ibug  est-ce vraiment necessaire?
164         v_ice(:,1) = 0.0
165
166         IF(ln_ctl) THEN
167            CALL prt_ctl(tab2d_1=u_oce , clinfo1=' lim_dyn  : u_oce :', tab2d_2=v_oce , clinfo2=' v_oce :')
168            CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :')
169         ENDIF
170
171         !                                         ! Ice-Ocean stress
172         !                                         ! ================
173         DO jj = 2, jpjm1
174!           jhemis = SIGN(1, jj - jeq )
175!           zsang  = REAL(jhemis) * sangvg
176!nemo new version modif
177            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg
178
179            DO ji = 2, jpim1
180               ! computation of wind stress over ocean in X and Y direction
181#if defined key_coupled && defined key_lim_cp1
182!              ztairx =  ( 1.0 - at_i(ji-1,jj)   ) * gtaux(ji-1,jj)   + &
183!                        ( 1.0 - at_i(ji,jj)     ) * gtaux(ji,jj  )   + &
184!                        ( 1.0 - at_i(ji-1,jj-1) ) * gtaux(ji-1,jj-1) + &
185!                        ( 1.0 - at_i(ji,jj-1)   ) * gtaux(ji,jj-1)
186
187!              ztairy =  ( 1.0 - at_i(ji-1,jj)   ) * gtauy(ji-1,jj  ) + &
188!                        ( 1.0 - at_i(ji,jj  )   ) * gtauy(ji,jj    ) + &
189!                        ( 1.0 - at_i(ji-1,jj-1) ) * gtauy(ji-1,jj-1) + &
190!                        ( 1.0 - at_i(ji,jj-1)   ) * gtauy(ji,jj-1)
191#else
192               ztairx =  ( 2.0 - at_i(ji,jj) - at_i(ji+1,jj) ) * gtaux(ji,jj) / cai * cao
193               ztairy =  ( 2.0 - at_i(ji,jj) - at_i(ji,jj+1) ) * gtauy(ji,jj) / cai * cao
194
195!              ! test on oscillations
196!              ztairx = 0.0
197!              ztairy = 0.0
198
199               zsfrldmx2 = at_i(ji,jj) + at_i(ji+1,jj)
200               zsfrldmy2 = at_i(ji,jj) + at_i(ji,jj+1)
201
202#endif
203               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
204               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
205               zmod     = SQRT( zu_ice * zu_ice + zv_ice * zv_ice ) 
206
207               ! quadratic drag formulation
208               ztglx   = zsfrldmx2 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice ) 
209               ztgly   = zsfrldmy2 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice ) 
210
211!              ! test on oscillations
212!              ztglx  = 0.0
213!              ztgly  = 0.0
214               ztglx   = ztglx * exp( - zmod / 0.5 )
215               ztgly   = ztglx * exp( - zmod / 0.5 )
216
217               tio_u(ji,jj) = - ( ztairx + 1.0 * ztglx ) / ( 2. * rau0 )
218               tio_v(ji,jj) = - ( ztairy + 1.0 * ztgly ) / ( 2. * rau0 )
219            END DO
220         END DO
221         
222         ! computation of friction velocity
223         DO jj = 2, jpjm1
224            DO ji = 2, jpim1
225
226               zu_ice   = u_ice(ji,jj) - u_io(ji,jj)
227               zt11  = rhoco * zu_ice * zu_ice
228
229               zu_ice   = u_ice(ji-1,jj) - u_io(ji-1,jj)
230               zt12  = rhoco * zu_ice * zu_ice
231
232               zv_ice   = v_ice(ji,jj) - v_io(ji,jj)
233               zt21  = rhoco * zv_ice * zv_ice
234
235               zv_ice   = v_ice(ji,jj-1) - v_io(ji,jj-1)
236               zt22  = rhoco * zv_ice * zv_ice
237               ztair2 = ( ( gtaux(ji,jj) + gtaux(ji-1,jj) ) / 2. )**2 + &
238                        ( ( gtauy(ji,jj) + gtauy(ji,jj-1) ) / 2. )**2
239
240               ! should not be weighted
241               zustm =  ( at_i(ji,jj)       ) * 0.5 * ( zt11 + zt12 + zt21 + zt22 )        &
242                  &  +  ( 1.0 - at_i(ji,jj) ) * SQRT( ztair2 )
243
244               ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj)
245
246            END DO
247         END DO
248
249       ELSE      ! If no ice dynamics 
250                   
251               ! virer ca (key_lim_cp1)
252          DO jj = 2, jpjm1
253             DO ji = 2, jpim1
254#if defined key_coupled && defined key_lim_cp1
255                tio_u(ji,jj) = - (  gtaux(ji  ,jj  ) + gtaux(ji-1,jj  )       &
256                   &              + gtaux(ji-1,jj-1) + gtaux(ji  ,jj-1) ) / ( 4 * rau0 )
257
258                tio_v(ji,jj) = - (  gtauy(ji  ,jj )  + gtauy(ji-1,jj  )       &
259                   &              + gtauy(ji-1,jj-1) + gtauy(ji  ,jj-1) ) / ( 4 * rau0 )
260#else
261                tio_u(ji,jj) = - gtaux(ji,jj) / cai * cao / rau0
262                tio_v(ji,jj) = - gtauy(ji,jj) / cai * cao / rau0 
263#endif
264                ztair2 = ( ( gtaux(ji,jj) + gtaux(ji-1,jj) ) / 2. )**2 + &
265                         ( ( gtauy(ji,jj) + gtauy(ji,jj-1) ) / 2. )**2
266                zustm        = SQRT( ztair2  )
267
268                ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj)
269            END DO
270         END DO
271
272      ENDIF
273
274      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
275      CALL lbc_lnk( tio_u, 'U', -1. )   ! I-point (i.e. ice U-V point)
276      CALL lbc_lnk( tio_v, 'V', -1. )   ! I-point (i.e. ice U-V point)
277
278      IF(ln_ctl) THEN
279        CALL prt_ctl(tab2d_1=tio_u , clinfo1=' lim_dyn  : tio_u :', tab2d_2=tio_v , clinfo2=' tio_v :')
280        CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :')
281        CALL prt_ctl(tab2d_1=at_i  , clinfo1=' lim_dyn  : at_i  :')
282      ENDIF
283
284!OPA9_DYN LIM2.1 2005
285!END OPA9_DYN LIM2.1 2005
286
287   END SUBROUTINE lim_dyn
288
289    SUBROUTINE lim_dyn_init
290      !!-------------------------------------------------------------------
291      !!                  ***  ROUTINE lim_dyn_init  ***
292      !!
293      !! ** Purpose : Physical constants and parameters linked to the ice
294      !!      dynamics
295      !!
296      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
297      !!       parameter values called at the first timestep (nit000)
298      !!
299      !! ** input   :   Namelist namicedyn
300      !!
301      !! history :
302      !!  8.5  ! 03-08 (C. Ethe) original code
303      !!  9.0  ! 07-03 (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle)
304      !!               EVP-Cgrid-LIM3
305      !!-------------------------------------------------------------------
306      NAMELIST/namicedyn/ epsd, alpha,     &
307         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
308         &                c_rhg, etamn, creepl, ecc, ahi0, &
309         &                nevp, telast, alphaevp
310      !!-------------------------------------------------------------------
311
312      ! Define the initial parameters
313      ! -------------------------
314
315      ! Read Namelist namicedyn
316      REWIND ( numnam_ice )
317      READ   ( numnam_ice  , namicedyn )
318      IF(lwp) THEN
319         WRITE(numout,*)
320         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
321         WRITE(numout,*) '~~~~~~~~~~~~'
322         WRITE(numout,*) '   tolerance parameter                              epsd   = ', epsd
323         WRITE(numout,*) '   coefficient for semi-implicit coriolis           alpha  = ', alpha
324         WRITE(numout,*) '   diffusion constant for dynamics                  dm     = ', dm
325         WRITE(numout,*) '   number of sub-time steps for relaxation          nbiter = ', nbiter
326         WRITE(numout,*) '   maximum number of iterations for relaxation      nbitdr = ', nbitdr
327         WRITE(numout,*) '   relaxation constant                              om     = ', om
328         WRITE(numout,*) '   maximum value for the residual of relaxation     resl   = ', resl
329         WRITE(numout,*) '   drag coefficient for oceanic stress              cw     = ', cw
330         WRITE(numout,*) '   turning angle for oceanic stress                 angvg  = ', angvg
331         WRITE(numout,*) '   first bulk-rheology parameter                    pstar  = ', pstar
332         WRITE(numout,*) '   second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
333         WRITE(numout,*) '   minimun value for viscosity                      etamn  = ', etamn
334         WRITE(numout,*) '   creep limit                                      creepl = ', creepl
335         WRITE(numout,*) '   eccentricity of the elliptical yield curve       ecc    = ', ecc
336         WRITE(numout,*) '   horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
337         WRITE(numout,*) '   number of iterations for subcycling              nevp   = ', nevp
338         WRITE(numout,*) '   timescale for elastic waves                      telast = ', telast
339         WRITE(numout,*) '   coefficient for the solution of int. stresses  alphaevp = ', alphaevp
340
341      ENDIF
342
343      usecc2 = 1.0 / ( ecc * ecc )
344      rhoco  = rau0 * cw
345      angvg  = angvg * rad
346      sangvg = SIN( angvg )
347      cangvg = COS( angvg )
348      pstarh = pstar / 2.0
349
350      !  Diffusion coefficients.
351      ahiu(:,:) = ahi0 * umask(:,:,1)
352      ahiv(:,:) = ahi0 * vmask(:,:,1)
353
354   END SUBROUTINE lim_dyn_init
355
356#else
357   !!----------------------------------------------------------------------
358   !!   Default option          Empty module           NO LIM sea-ice model
359   !!----------------------------------------------------------------------
360CONTAINS
361   SUBROUTINE lim_dyn         ! Empty routine
362   END SUBROUTINE lim_dyn
363#endif 
364
365   !!======================================================================
366END MODULE limdyn
Note: See TracBrowser for help on using the repository browser.