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

source: branches/dev_001_SBC/NEMO/LIM_SRC_3/limdyn.F90 @ 884

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

dev_001_SBC: Step I: add the LIM 3.0 component into LIM_SRC_3 directory, see ticket: #111

File size: 17.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 par_ice
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   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limdyn.F90,v 1.5 2005/03/27 18:34:41 opalod Exp $
44   !! This software is governed by the CeCILL licence see 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         ! ocean velocity
93         u_oce(:,:)  = u_io(:,:) * tmu(:,:)
94         v_oce(:,:)  = v_io(:,:) * tmv(:,:)
95         
96         old_u_ice(:,:) = u_ice(:,:) * tmu(:,:)
97         old_v_ice(:,:) = v_ice(:,:) * tmv(:,:)
98
99         ! Rheology (ice dynamics)
100         ! ========
101
102         !  Define the j-limits where ice rheology is computed
103         ! ---------------------------------------------------
104
105         IF( lk_mpp .OR. nbit_cmp == 1 ) THEN                    ! mpp: compute over the whole domain
106            i_j1 = 1
107            i_jpj = jpj
108            IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
109            CALL lim_rhg( i_j1, i_jpj )
110         ELSE                                 ! optimization of the computational area
111
112            DO jj = 1, jpj
113               zind(jj) = SUM( 1.0 - at_i (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line
114               zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line
115            END DO
116
117            IF( l_jeq ) THEN                     ! local domain include both hemisphere
118               !                                 ! Rheology is computed in each hemisphere
119               !                                 ! only over the ice cover latitude strip
120               ! Northern hemisphere
121               i_j1  = njeq
122               i_jpj = jpj
123               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
124                  i_j1 = i_j1 + 1
125               END DO
126               i_j1 = MAX( 1, i_j1-1 )
127               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : NH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
128               CALL lim_rhg( i_j1, i_jpj )
129
130               ! Southern hemisphere
131               i_j1  =  1
132               i_jpj = njeq
133               DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
134                  i_jpj = i_jpj - 1
135               END DO
136               i_jpj = MIN( jpj, i_jpj+2 )
137               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : SH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
138
139               CALL lim_rhg( i_j1, i_jpj )
140
141            ELSE                                 ! local domain extends over one hemisphere only
142               !                                 ! Rheology is computed only over the ice cover
143               !                                 ! latitude strip
144               i_j1  = 1
145               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
146                  i_j1 = i_j1 + 1
147               END DO
148               i_j1 = MAX( 1, i_j1-1 )
149
150               i_jpj  = jpj
151               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
152                  i_jpj = i_jpj - 1
153               END DO
154               i_jpj = MIN( jpj, i_jpj+2)
155
156               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : one hemisphere:  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
157
158               CALL lim_rhg( i_j1, i_jpj )
159
160            ENDIF
161
162         ENDIF
163
164         ! Ice-Ocean stress
165         ! ================
166         DO jj = 2, jpjm1
167            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg
168
169            DO ji = fs_2, fs_jpim1
170               ! computation of wind stress over ocean in X and Y direction
171#if defined key_coupled && defined key_lim_cp1
172!              ztairx =  ( 1.0 - at_i(ji-1,jj)   ) * gtaux(ji-1,jj)   + &
173!                        ( 1.0 - at_i(ji,jj)     ) * gtaux(ji,jj  )   + &
174!                        ( 1.0 - at_i(ji-1,jj-1) ) * gtaux(ji-1,jj-1) + &
175!                        ( 1.0 - at_i(ji,jj-1)   ) * gtaux(ji,jj-1)
176
177!              ztairy =  ( 1.0 - at_i(ji-1,jj)   ) * gtauy(ji-1,jj  ) + &
178!                        ( 1.0 - at_i(ji,jj  )   ) * gtauy(ji,jj    ) + &
179!                        ( 1.0 - at_i(ji-1,jj-1) ) * gtauy(ji-1,jj-1) + &
180!                        ( 1.0 - at_i(ji,jj-1)   ) * gtauy(ji,jj-1)
181#else
182               ztairx =  ( 2.0 - at_i(ji,jj) - at_i(ji+1,jj) ) * gtaux(ji,jj) / cai * cao
183               ztairy =  ( 2.0 - at_i(ji,jj) - at_i(ji,jj+1) ) * gtauy(ji,jj) / cai * cao
184
185               zsfrldmx2 = at_i(ji,jj) + at_i(ji+1,jj)
186               zsfrldmy2 = at_i(ji,jj) + at_i(ji,jj+1)
187
188#endif
189               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
190               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
191               zmod     = SQRT( zu_ice * zu_ice + zv_ice * zv_ice ) 
192
193               ! quadratic drag formulation
194               ztglx   = zsfrldmx2 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice ) 
195               ztgly   = zsfrldmy2 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice ) 
196!
197!              ! IMPORTANT
198!              ! these lignes are bound to prevent numerical oscillations
199!              ! in the ice-ocean stress
200!              ! They are physically ill-based. There is a cleaner solution
201!              ! to try (remember discussion in Paris Gurvan)
202!
203               ztglx   = ztglx * exp( - zmod / 0.5 )
204               ztgly   = ztglx * exp( - zmod / 0.5 )
205
206               tio_u(ji,jj) = - ( ztairx + 1.0 * ztglx ) / ( 2. * rau0 )
207               tio_v(ji,jj) = - ( ztairy + 1.0 * ztgly ) / ( 2. * rau0 )
208            END DO
209         END DO
210         
211         ! computation of friction velocity
212         DO jj = 2, jpjm1
213            DO ji = fs_2, fs_jpim1
214
215               zu_ice   = u_ice(ji,jj) - u_io(ji,jj)
216               zt11  = rhoco * zu_ice * zu_ice
217
218               zu_ice   = u_ice(ji-1,jj) - u_io(ji-1,jj)
219               zt12  = rhoco * zu_ice * zu_ice
220
221               zv_ice   = v_ice(ji,jj) - v_io(ji,jj)
222               zt21  = rhoco * zv_ice * zv_ice
223
224               zv_ice   = v_ice(ji,jj-1) - v_io(ji,jj-1)
225               zt22  = rhoco * zv_ice * zv_ice
226               ztair2 = ( ( gtaux(ji,jj) + gtaux(ji-1,jj) ) / 2. )**2 + &
227                        ( ( gtauy(ji,jj) + gtauy(ji,jj-1) ) / 2. )**2
228
229               ! should not be weighted
230               zustm =  ( at_i(ji,jj)       ) * 0.5 * ( zt11 + zt12 + zt21 + zt22 )        &
231                  &  +  ( 1.0 - at_i(ji,jj) ) * SQRT( ztair2 )
232
233               ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj)
234
235            END DO
236         END DO
237
238       ELSE      ! If no ice dynamics 
239                   
240               ! virer ca (key_lim_cp1)
241          DO jj = 2, jpjm1
242             DO ji = fs_2, fs_jpim1
243#if defined key_coupled && defined key_lim_cp1
244                tio_u(ji,jj) = - (  gtaux(ji  ,jj  ) + gtaux(ji-1,jj  )       &
245                   &              + gtaux(ji-1,jj-1) + gtaux(ji  ,jj-1) ) / ( 4 * rau0 )
246
247                tio_v(ji,jj) = - (  gtauy(ji  ,jj )  + gtauy(ji-1,jj  )       &
248                   &              + gtauy(ji-1,jj-1) + gtauy(ji  ,jj-1) ) / ( 4 * rau0 )
249#else
250                tio_u(ji,jj) = - gtaux(ji,jj) / cai * cao / rau0
251                tio_v(ji,jj) = - gtauy(ji,jj) / cai * cao / rau0 
252#endif
253                ztair2 = ( ( gtaux(ji,jj) + gtaux(ji-1,jj) ) / 2. )**2 + &
254                         ( ( gtauy(ji,jj) + gtauy(ji,jj-1) ) / 2. )**2
255                zustm        = SQRT( ztair2  )
256
257                ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj)
258            END DO
259         END DO
260
261      ENDIF
262
263      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
264      CALL lbc_lnk( tio_u, 'U', -1. )   ! I-point (i.e. ice U-V point)
265      CALL lbc_lnk( tio_v, 'V', -1. )   ! I-point (i.e. ice U-V point)
266
267      IF(ln_ctl) THEN   ! Control print
268         CALL prt_ctl_info(' ')
269         CALL prt_ctl_info(' - Cell values : ')
270         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
271         CALL prt_ctl(tab2d_1=tio_u     , clinfo1=' lim_dyn  : tio_u     :', tab2d_2=tio_v , clinfo2=' tio_v :')
272         CALL prt_ctl(tab2d_1=ust2s     , clinfo1=' lim_dyn  : ust2s     :')
273         CALL prt_ctl(tab2d_1=divu_i    , clinfo1=' lim_dyn  : divu_i    :')
274         CALL prt_ctl(tab2d_1=delta_i   , clinfo1=' lim_dyn  : delta_i   :')
275         CALL prt_ctl(tab2d_1=strength  , clinfo1=' lim_dyn  : strength  :')
276         CALL prt_ctl(tab2d_1=area      , clinfo1=' lim_dyn  : cell area :')
277         CALL prt_ctl(tab2d_1=at_i      , clinfo1=' lim_dyn  : at_i      :')
278         CALL prt_ctl(tab2d_1=vt_i      , clinfo1=' lim_dyn  : vt_i      :')
279         CALL prt_ctl(tab2d_1=vt_s      , clinfo1=' lim_dyn  : vt_s      :')
280         CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' lim_dyn  : stress1_i :')
281         CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' lim_dyn  : stress2_i :')
282         CALL prt_ctl(tab2d_1=stress12_i, clinfo1=' lim_dyn  : stress12_i:')
283         DO jl = 1, jpl
284            CALL prt_ctl_info(' ')
285            CALL prt_ctl_info(' - Category : ', ivar1=jl)
286            CALL prt_ctl_info('   ~~~~~~~~~~')
287            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_dyn  : a_i      : ')
288            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_i     : ')
289            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_s     : ')
290            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_dyn  : v_i      : ')
291            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_dyn  : v_s      : ')
292            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : e_s      : ')
293            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_dyn  : t_su     : ')
294            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : t_snow   : ')
295            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_dyn  : sm_i     : ')
296            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_dyn  : smv_i    : ')
297            DO ja = 1, nlay_i
298               CALL prt_ctl_info(' ')
299               CALL prt_ctl_info(' - Layer : ', ivar1=ja)
300               CALL prt_ctl_info('   ~~~~~~~')
301               CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : t_i      : ')
302               CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : e_i      : ')
303            END DO
304         END DO
305      ENDIF
306
307   END SUBROUTINE lim_dyn
308
309    SUBROUTINE lim_dyn_init
310      !!-------------------------------------------------------------------
311      !!                  ***  ROUTINE lim_dyn_init  ***
312      !!
313      !! ** Purpose : Physical constants and parameters linked to the ice
314      !!      dynamics
315      !!
316      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
317      !!       parameter values called at the first timestep (nit000)
318      !!
319      !! ** input   :   Namelist namicedyn
320      !!
321      !! history :
322      !!  8.5  ! 03-08 (C. Ethe) original code
323      !!  9.0  ! 07-03 (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle)
324      !!               EVP-Cgrid-LIM3
325      !!-------------------------------------------------------------------
326      NAMELIST/namicedyn/ epsd, alpha,     &
327         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
328         &                c_rhg, etamn, creepl, ecc, ahi0, &
329         &                nevp, telast, alphaevp
330      !!-------------------------------------------------------------------
331
332      ! Define the initial parameters
333      ! -------------------------
334
335      ! Read Namelist namicedyn
336      REWIND ( numnam_ice )
337      READ   ( numnam_ice  , namicedyn )
338      IF(lwp) THEN
339         WRITE(numout,*)
340         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
341         WRITE(numout,*) '~~~~~~~~~~~~'
342         WRITE(numout,*) '   tolerance parameter                              epsd   = ', epsd
343         WRITE(numout,*) '   coefficient for semi-implicit coriolis           alpha  = ', alpha
344         WRITE(numout,*) '   diffusion constant for dynamics                  dm     = ', dm
345         WRITE(numout,*) '   number of sub-time steps for relaxation          nbiter = ', nbiter
346         WRITE(numout,*) '   maximum number of iterations for relaxation      nbitdr = ', nbitdr
347         WRITE(numout,*) '   relaxation constant                              om     = ', om
348         WRITE(numout,*) '   maximum value for the residual of relaxation     resl   = ', resl
349         WRITE(numout,*) '   drag coefficient for oceanic stress              cw     = ', cw
350         WRITE(numout,*) '   turning angle for oceanic stress                 angvg  = ', angvg
351         WRITE(numout,*) '   first bulk-rheology parameter                    pstar  = ', pstar
352         WRITE(numout,*) '   second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
353         WRITE(numout,*) '   minimun value for viscosity                      etamn  = ', etamn
354         WRITE(numout,*) '   creep limit                                      creepl = ', creepl
355         WRITE(numout,*) '   eccentricity of the elliptical yield curve       ecc    = ', ecc
356         WRITE(numout,*) '   horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
357         WRITE(numout,*) '   number of iterations for subcycling              nevp   = ', nevp
358         WRITE(numout,*) '   timescale for elastic waves                      telast = ', telast
359         WRITE(numout,*) '   coefficient for the solution of int. stresses  alphaevp = ', alphaevp
360
361      ENDIF
362
363      usecc2 = 1.0 / ( ecc * ecc )
364      rhoco  = rau0 * cw
365      angvg  = angvg * rad
366      sangvg = SIN( angvg )
367      cangvg = COS( angvg )
368      pstarh = pstar / 2.0
369
370      !  Diffusion coefficients.
371      ahiu(:,:) = ahi0 * umask(:,:,1)
372      ahiv(:,:) = ahi0 * vmask(:,:,1)
373
374   END SUBROUTINE lim_dyn_init
375
376#else
377   !!----------------------------------------------------------------------
378   !!   Default option          Empty module           NO LIM sea-ice model
379   !!----------------------------------------------------------------------
380CONTAINS
381   SUBROUTINE lim_dyn         ! Empty routine
382   END SUBROUTINE lim_dyn
383#endif 
384
385   !!======================================================================
386END MODULE limdyn
Note: See TracBrowser for help on using the repository browser.