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

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

Add control prints

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   !! * Module variables
36   REAL(wp)  ::  rone    = 1.e0   ! constant value
37
38   !!----------------------------------------------------------------------
39   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)
40   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limdyn.F90,v 1.5 2005/03/27 18:34:41 opalod Exp $
41   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
42   !!----------------------------------------------------------------------
43
44CONTAINS
45
46   SUBROUTINE lim_dyn
47      !!-------------------------------------------------------------------
48      !!               ***  ROUTINE lim_dyn  ***
49      !!               
50      !! ** Purpose :   compute ice velocity and ocean-ice stress
51      !!               
52      !! ** Method  :
53      !!
54      !! ** Action  : - Initialisation
55      !!              - Call of the dynamic routine for each hemisphere
56      !!              - computation of the stress at the ocean surface         
57      !!              - treatment of the case if no ice dynamic
58      !! History :
59      !!   1.0  !  01-04   (LIM)  Original code
60      !!   2.0  !  02-08   (C. Ethe, G. Madec)  F90, mpp
61      !!   3.0  !  2007-03 (M.A. Morales Maqueda, S. Bouillon, M. Vancoppenolle)
62      !!                   LIM3, EVP, C-grid
63      !!------------------------------------------------------------------------------------
64      !! * Local variables
65      INTEGER ::   ji, jj, jl, ja    ! dummy loop indices
66      INTEGER :: i_j1, i_jpj         ! Starting/ending j-indices for rheology
67
68      REAL(wp) ::   &
69         ztairx, ztairy,          &  ! tempory scalars
70         zsang , zmod,            &
71         ztglx , ztgly ,          &
72         zt11, zt12, zt21, zt22 , &
73         zustm,                   &
74         zsfrldmx2, zsfrldmy2,    &
75         zu_ice, zv_ice, ztair2
76
77      REAL(wp),DIMENSION(jpj) ::   &
78           zind,                     &  ! i-averaged indicator of sea-ice
79           zmsk                         ! i-averaged of tmask
80      !!---------------------------------------------------------------------
81
82      WRITE(numout,*) ' lim_dyn : Ice dynamics '
83      WRITE(numout,*) ' ~~~~~~~ '
84
85      IF( numit == nstart  )   CALL lim_dyn_init   ! Initialization (first time-step only)
86     
87      IF ( ln_limdyn ) THEN
88
89         ! ocean velocity
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         ELSE                                 ! optimization of the computational area
108
109            DO jj = 1, jpj
110               zind(jj) = SUM( 1.0 - at_i (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line
111               zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line
112            END DO
113
114            IF( l_jeq ) THEN                     ! local domain include both hemisphere
115               !                                 ! Rheology is computed in each hemisphere
116               !                                 ! only over the ice cover latitude strip
117               ! Northern hemisphere
118               i_j1  = njeq
119               i_jpj = jpj
120               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
121                  i_j1 = i_j1 + 1
122               END DO
123               i_j1 = MAX( 1, i_j1-1 )
124               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : NH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
125               CALL lim_rhg( i_j1, i_jpj )
126
127               ! Southern hemisphere
128               i_j1  =  1
129               i_jpj = njeq
130               DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
131                  i_jpj = i_jpj - 1
132               END DO
133               i_jpj = MIN( jpj, i_jpj+2 )
134               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : SH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
135
136               CALL lim_rhg( i_j1, i_jpj )
137
138            ELSE                                 ! local domain extends over one hemisphere only
139               !                                 ! Rheology is computed only over the ice cover
140               !                                 ! latitude strip
141               i_j1  = 1
142               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
143                  i_j1 = i_j1 + 1
144               END DO
145               i_j1 = MAX( 1, i_j1-1 )
146
147               i_jpj  = jpj
148               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
149                  i_jpj = i_jpj - 1
150               END DO
151               i_jpj = MIN( jpj, i_jpj+2)
152
153               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : one hemisphere:  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
154
155               CALL lim_rhg( i_j1, i_jpj )
156
157            ENDIF
158
159         ENDIF
160
161         u_ice(:,1) = 0.0       !ibug  est-ce vraiment necessaire?
162         v_ice(:,1) = 0.0
163
164         IF(ln_ctl) THEN
165            CALL prt_ctl(tab2d_1=u_oce , clinfo1=' lim_dyn  : u_oce :', tab2d_2=v_oce , clinfo2=' v_oce :')
166            CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :')
167         ENDIF
168
169         ! Ice-Ocean stress
170         ! ================
171         DO jj = 2, jpjm1
172            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg
173
174            DO ji = 2, jpim1
175               ! computation of wind stress over ocean in X and Y direction
176#if defined key_coupled && defined key_lim_cp1
177!              ztairx =  ( 1.0 - at_i(ji-1,jj)   ) * gtaux(ji-1,jj)   + &
178!                        ( 1.0 - at_i(ji,jj)     ) * gtaux(ji,jj  )   + &
179!                        ( 1.0 - at_i(ji-1,jj-1) ) * gtaux(ji-1,jj-1) + &
180!                        ( 1.0 - at_i(ji,jj-1)   ) * gtaux(ji,jj-1)
181
182!              ztairy =  ( 1.0 - at_i(ji-1,jj)   ) * gtauy(ji-1,jj  ) + &
183!                        ( 1.0 - at_i(ji,jj  )   ) * gtauy(ji,jj    ) + &
184!                        ( 1.0 - at_i(ji-1,jj-1) ) * gtauy(ji-1,jj-1) + &
185!                        ( 1.0 - at_i(ji,jj-1)   ) * gtauy(ji,jj-1)
186#else
187               ztairx =  ( 2.0 - at_i(ji,jj) - at_i(ji+1,jj) ) * gtaux(ji,jj) / cai * cao
188               ztairy =  ( 2.0 - at_i(ji,jj) - at_i(ji,jj+1) ) * gtauy(ji,jj) / cai * cao
189
190               zsfrldmx2 = at_i(ji,jj) + at_i(ji+1,jj)
191               zsfrldmy2 = at_i(ji,jj) + at_i(ji,jj+1)
192
193#endif
194               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
195               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
196               zmod     = SQRT( zu_ice * zu_ice + zv_ice * zv_ice ) 
197
198               ! quadratic drag formulation
199               ztglx   = zsfrldmx2 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice ) 
200               ztgly   = zsfrldmy2 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice ) 
201!
202!              ! IMPORTANT
203!              ! these lignes are bound to prevent numerical oscillations
204!              ! in the ice-ocean stress
205!              ! They are physically ill-based. There is a cleaner solution
206!              ! to try (remember discussion in Paris Gurvan)
207!
208               ztglx   = ztglx * exp( - zmod / 0.5 )
209               ztgly   = ztglx * exp( - zmod / 0.5 )
210
211               tio_u(ji,jj) = - ( ztairx + 1.0 * ztglx ) / ( 2. * rau0 )
212               tio_v(ji,jj) = - ( ztairy + 1.0 * ztgly ) / ( 2. * rau0 )
213            END DO
214         END DO
215         
216         ! computation of friction velocity
217         DO jj = 2, jpjm1
218            DO ji = 2, jpim1
219
220               zu_ice   = u_ice(ji,jj) - u_io(ji,jj)
221               zt11  = rhoco * zu_ice * zu_ice
222
223               zu_ice   = u_ice(ji-1,jj) - u_io(ji-1,jj)
224               zt12  = rhoco * zu_ice * zu_ice
225
226               zv_ice   = v_ice(ji,jj) - v_io(ji,jj)
227               zt21  = rhoco * zv_ice * zv_ice
228
229               zv_ice   = v_ice(ji,jj-1) - v_io(ji,jj-1)
230               zt22  = rhoco * zv_ice * zv_ice
231               ztair2 = ( ( gtaux(ji,jj) + gtaux(ji-1,jj) ) / 2. )**2 + &
232                        ( ( gtauy(ji,jj) + gtauy(ji,jj-1) ) / 2. )**2
233
234               ! should not be weighted
235               zustm =  ( at_i(ji,jj)       ) * 0.5 * ( zt11 + zt12 + zt21 + zt22 )        &
236                  &  +  ( 1.0 - at_i(ji,jj) ) * SQRT( ztair2 )
237
238               ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj)
239
240            END DO
241         END DO
242
243       ELSE      ! If no ice dynamics 
244                   
245               ! virer ca (key_lim_cp1)
246          DO jj = 2, jpjm1
247             DO ji = 2, jpim1
248#if defined key_coupled && defined key_lim_cp1
249                tio_u(ji,jj) = - (  gtaux(ji  ,jj  ) + gtaux(ji-1,jj  )       &
250                   &              + gtaux(ji-1,jj-1) + gtaux(ji  ,jj-1) ) / ( 4 * rau0 )
251
252                tio_v(ji,jj) = - (  gtauy(ji  ,jj )  + gtauy(ji-1,jj  )       &
253                   &              + gtauy(ji-1,jj-1) + gtauy(ji  ,jj-1) ) / ( 4 * rau0 )
254#else
255                tio_u(ji,jj) = - gtaux(ji,jj) / cai * cao / rau0
256                tio_v(ji,jj) = - gtauy(ji,jj) / cai * cao / rau0 
257#endif
258                ztair2 = ( ( gtaux(ji,jj) + gtaux(ji-1,jj) ) / 2. )**2 + &
259                         ( ( gtauy(ji,jj) + gtauy(ji,jj-1) ) / 2. )**2
260                zustm        = SQRT( ztair2  )
261
262                ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj)
263            END DO
264         END DO
265
266      ENDIF
267
268      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
269      CALL lbc_lnk( tio_u, 'U', -1. )   ! I-point (i.e. ice U-V point)
270      CALL lbc_lnk( tio_v, 'V', -1. )   ! I-point (i.e. ice U-V point)
271
272      IF(ln_ctl) THEN   ! Control print
273         CALL prt_ctl(tab2d_1=tio_u     , clinfo1=' lim_dyn  : tio_u     :', tab2d_2=tio_v , clinfo2=' tio_v :')
274         CALL prt_ctl(tab2d_1=ust2s     , clinfo1=' lim_dyn  : ust2s     :')
275         CALL prt_ctl(tab2d_1=divu_i    , clinfo1=' lim_dyn  : divu_i    :')
276         CALL prt_ctl(tab2d_1=delta_i   , clinfo1=' lim_dyn  : delta_i   :')
277         CALL prt_ctl(tab2d_1=strength  , clinfo1=' lim_dyn  : strength  :')
278         CALL prt_ctl(tab2d_1=area      , clinfo1=' lim_dyn  : cell area :')
279         CALL prt_ctl(tab2d_1=at_i      , clinfo1=' lim_dyn  : at_i      :')
280         CALL prt_ctl(tab2d_1=vt_i      , clinfo1=' lim_dyn  : vt_i      :')
281         CALL prt_ctl(tab2d_1=vt_s      , clinfo1=' lim_dyn  : vt_s      :')
282         CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' lim_dyn  : stress1_i :')
283         CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' lim_dyn  : stress2_i :')
284         CALL prt_ctl(tab2d_1=stress12_i, clinfo1=' lim_dyn  : stress12_i:')
285         DO jl = 1, jpl
286            CALL prt_ctl_info(' - Category : ', ivar1=jl)
287            CALL prt_ctl_info('   ~~~~~~~~~~')
288            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_dyn  : a_i      : ')
289            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_i     : ')
290            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_s     : ')
291            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_dyn  : v_i      : ')
292            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_dyn  : v_s      : ')
293            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : e_s      : ')
294            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_dyn  : t_su     : ')
295            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : t_snow   : ')
296            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_dyn  : sm_i     : ')
297            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_dyn  : smv_i    : ')
298            DO ja = 1, nlay_i
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.