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_2.F90 in branches/devmercator2010/NEMO/LIM_SRC_2 – NEMO

source: branches/devmercator2010/NEMO/LIM_SRC_2/limdyn_2.F90 @ 2076

Last change on this file since 2076 was 2076, checked in by cbricaud, 14 years ago

add change dev_1784_EVP
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 13.2 KB
Line 
1MODULE limdyn_2
2   !!======================================================================
3   !!                     ***  MODULE  limdyn_2  ***
4   !!   Sea-Ice dynamics : 
5   !!======================================================================
6   !! History :   1.0  !  01-04  (LIM)  Original code
7   !!             2.0  !  02-08  (C. Ethe, G. Madec)  F90, mpp
8   !!             2.0  !  03-08  (C. Ethe) add lim_dyn_init
9   !!             2.0  !  06-07  (G. Madec)  Surface module
10   !!             3.3  !  09-05  (G.Garric) addition of the lim2_evp cas
11   !!---------------------------------------------------------------------
12#if defined key_lim2
13   !!----------------------------------------------------------------------
14   !!   'key_lim2' :                                  LIM 2.0 sea-ice model
15   !!----------------------------------------------------------------------
16   !!    lim_dyn_2      : computes ice velocities
17   !!    lim_dyn_init_2 : initialization and namelist read
18   !!----------------------------------------------------------------------
19   USE dom_oce        ! ocean space and time domain
20   USE sbc_oce        !
21   USE phycst         !
22   USE ice_2          !
23   USE dom_ice_2      !
24   USE limistate_2    !
25#if defined key_lim2_vp
26   USE limrhg_2       ! ice rheology
27#else
28   USE limrhg         ! ice rheology
29#endif
30   USE lbclnk         !
31   USE lib_mpp        !
32   USE in_out_manager ! I/O manager
33   USE prtctl         ! Print control
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   lim_dyn_2 ! routine called by sbc_ice_lim
39
40   !! * Module variables
41   REAL(wp)  ::  rone    = 1.e0   ! constant value
42
43#  include "vectopt_loop_substitute.h90"
44   !!----------------------------------------------------------------------
45   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)
46   !! $Id$
47   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
48   !!----------------------------------------------------------------------
49
50CONTAINS
51
52   SUBROUTINE lim_dyn_2( kt )
53      !!-------------------------------------------------------------------
54      !!               ***  ROUTINE lim_dyn_2  ***
55      !!               
56      !! ** Purpose :   compute ice velocity and ocean-ice friction velocity
57      !!               
58      !! ** Method  :
59      !!
60      !! ** Action  : - Initialisation
61      !!              - Call of the dynamic routine for each hemisphere
62      !!              - computation of the friction velocity at the sea-ice base
63      !!              - treatment of the case if no ice dynamic
64      !!---------------------------------------------------------------------
65      INTEGER, INTENT(in) ::   kt     ! number of iteration
66      !!
67      INTEGER  ::   ji, jj             ! dummy loop indices
68      INTEGER  ::   i_j1, i_jpj        ! Starting/ending j-indices for rheology
69      REAL(wp) ::   zcoef              ! temporary scalar
70      REAL(wp), DIMENSION(jpj)     ::   zind           ! i-averaged indicator of sea-ice
71      REAL(wp), DIMENSION(jpj)     ::   zmsk           ! i-averaged of tmask
72      REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io   ! ice-ocean velocity
73      !!---------------------------------------------------------------------
74
75      IF( kt == nit000 )   CALL lim_dyn_init_2   ! Initialization (first time-step only)
76     
77      IF( ln_limdyn ) THEN
78         !
79         ! Mean ice and snow thicknesses.         
80         hsnm(:,:)  = ( 1.0 - frld(:,:) ) * hsnif(:,:)
81         hicm(:,:)  = ( 1.0 - frld(:,:) ) * hicif(:,:)
82         !
83         !                                     ! Rheology (ice dynamics)
84         !                                     ! ========
85         
86         !  Define the j-limits where ice rheology is computed
87         ! ---------------------------------------------------
88         
89         IF( lk_mpp .OR. nbit_cmp == 1 ) THEN                    ! mpp: compute over the whole domain
90            i_j1 = 1   
91            i_jpj = jpj
92            IF(ln_ctl)   CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
93#if defined key_lim2_vp
94            CALL lim_rhg_2( i_j1, i_jpj )
95#else
96            CALL lim_rhg( i_j1, i_jpj ) !!!!cbr CALL lim_rhg( i_j1, i_jpj, kt )
97#endif
98            !
99         ELSE                                 ! optimization of the computational area
100            !
101            DO jj = 1, jpj
102               zind(jj) = SUM( frld (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line
103               zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line
104            END DO
105            !
106            IF( l_jeq ) THEN                     ! local domain include both hemisphere
107               !                                 ! Rheology is computed in each hemisphere
108               !                                 ! only over the ice cover latitude strip
109               ! Northern hemisphere
110               i_j1  = njeq
111               i_jpj = jpj
112               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
113                  i_j1 = i_j1 + 1
114               END DO
115#if defined key_lim2_vp
116               i_j1 = MAX( 1, i_j1-1 )
117               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
118               !
119               CALL lim_rhg_2( i_j1, i_jpj )
120#else
121               i_j1 = MAX( 1, i_j1-2 )
122               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
123               CALL lim_rhg( i_j1, i_jpj )
124#endif
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#if defined key_lim2_vp
133               i_jpj = MIN( jpj, i_jpj+2 )
134               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
135               !
136               CALL lim_rhg_2( i_j1, i_jpj )
137#else
138               i_jpj = MIN( jpj, i_jpj+1 )
139               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
140               CALL lim_rhg( i_j1, i_jpj ) !!!!cbr CALL lim_rhg( i_j1, i_jpj, kt )
141#endif
142               !
143            ELSE                                 ! local domain extends over one hemisphere only
144               !                                 ! Rheology is computed only over the ice cover
145               !                                 ! latitude strip
146               i_j1  = 1
147               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
148                  i_j1 = i_j1 + 1
149               END DO
150               i_j1 = MAX( 1, i_j1-1 )
151   
152               i_jpj  = jpj
153               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
154                  i_jpj = i_jpj - 1
155               END DO
156#if defined key_lim2_vp
157               i_jpj = MIN( jpj, i_jpj+2)
158
159               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
160               !
161               CALL lim_rhg_2( i_j1, i_jpj )
162#else
163               i_j1 = MAX( 1, i_j1-2 )
164               i_jpj = MIN( jpj, i_jpj+1)
165               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
166               CALL lim_rhg( i_j1, i_jpj ) !!!!cbr CALL lim_rhg( i_j1, i_jpj, kt )
167#endif
168               !
169            ENDIF
170            !
171         ENDIF
172
173         IF(ln_ctl)   CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :')
174         
175         ! computation of friction velocity
176         ! --------------------------------
177
178     SELECT CASE( cl_grid )
179
180      CASE( 'C' )                          ! C-grid ice dynamics
181         !?????????????????????????????????
182         ! ice-ocean velocity at U & V-points (u_ice vi_ice at U- & V-points ; ssu_m, ssv_m at U- & V-points)
183         zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
184         zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
185
186
187      CASE( 'B' )                          ! B-grid ice dynamics
188         ! ice-ocean velocity at U & V-points (u_ice v_ice at I-point ; ssu_m, ssv_m at U- & V-points)
189         
190         DO jj = 1, jpjm1
191            DO ji = 1, jpim1   ! NO vector opt.
192               zu_io(ji,jj) = 0.5 * ( u_ice(ji+1,jj+1) + u_ice(ji+1,jj  ) ) - ssu_m(ji,jj)
193               zv_io(ji,jj) = 0.5 * ( v_ice(ji+1,jj+1) + v_ice(ji  ,jj+1) ) - ssv_m(ji,jj)
194            END DO
195         END DO
196
197      END SELECT
198
199         ! frictional velocity at T-point
200         DO jj = 2, jpjm1
201            DO ji = 2, jpim1   ! NO vector opt. because of zu_io
202               ust2s(ji,jj) = 0.5 * cw                                                          &
203                  &         * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
204                  &            + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj)
205            END DO
206         END DO
207         !
208      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
209         !
210         zcoef = SQRT( 0.5 ) / rau0
211         DO jj = 2, jpjm1
212            DO ji = fs_2, fs_jpim1   ! vector opt.
213               ust2s(ji,jj) = zcoef * tms(ji,jj) * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
214                  &                                     + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) )
215            END DO
216         END DO
217         !
218      ENDIF
219      !
220      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
221      !
222      IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :')
223
224   END SUBROUTINE lim_dyn_2
225
226
227   SUBROUTINE lim_dyn_init_2
228      !!-------------------------------------------------------------------
229      !!                  ***  ROUTINE lim_dyn_init_2  ***
230      !!
231      !! ** Purpose :   Physical constants and parameters linked to the ice
232      !!              dynamics
233      !!
234      !! ** Method  :   Read the namicedyn namelist and check the ice-dynamic
235      !!              parameter values
236      !!
237      !! ** input   :   Namelist namicedyn
238      !!-------------------------------------------------------------------
239      NAMELIST/namicedyn/ epsd, alpha,     &
240         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
241         &                c_rhg, etamn, creepl, ecc, ahi0, &
242         &                nevp, telast, alphaevp
243      !!-------------------------------------------------------------------
244
245      REWIND ( numnam_ice )                       ! Read Namelist namicedyn
246      READ   ( numnam_ice  , namicedyn )
247
248      IF(lwp) THEN                                ! Control print
249         WRITE(numout,*)
250         WRITE(numout,*) 'lim_dyn_init_2: ice parameters for ice dynamics '
251         WRITE(numout,*) '~~~~~~~~~~~~~~'
252         WRITE(numout,*) '       tolerance parameter                              epsd   = ', epsd
253         WRITE(numout,*) '       coefficient for semi-implicit coriolis           alpha  = ', alpha
254         WRITE(numout,*) '       diffusion constant for dynamics                  dm     = ', dm
255         WRITE(numout,*) '       number of sub-time steps for relaxation          nbiter = ', nbiter
256         WRITE(numout,*) '       maximum number of iterations for relaxation      nbitdr = ', nbitdr
257         WRITE(numout,*) '       relaxation constant                              om     = ', om
258         WRITE(numout,*) '       maximum value for the residual of relaxation     resl   = ', resl
259         WRITE(numout,*) '       drag coefficient for oceanic stress              cw     = ', cw
260         WRITE(numout,*) '       turning angle for oceanic stress                 angvg  = ', angvg, ' degrees'
261         WRITE(numout,*) '       first bulk-rheology parameter                    pstar  = ', pstar
262         WRITE(numout,*) '       second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
263         WRITE(numout,*) '       minimun value for viscosity                      etamn  = ', etamn
264         WRITE(numout,*) '       creep limit                                      creepl = ', creepl
265         WRITE(numout,*) '       eccentricity of the elliptical yield curve       ecc    = ', ecc
266         WRITE(numout,*) '       horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
267         WRITE(numout,*) '       number of iterations for subcycling              nevp   = ', nevp
268         WRITE(numout,*) '       timescale for elastic waves                      telast = ', telast
269         WRITE(numout,*) '       coefficient for the solution of int. stresses  alphaevp = ', alphaevp
270      ENDIF
271
272      !  Initialization
273      usecc2 = 1.0 / ( ecc * ecc )
274      rhoco  = rau0 * cw
275      angvg  = angvg * rad      ! convert angvg from degree to radian
276      sangvg = SIN( angvg )
277      cangvg = COS( angvg )
278      pstarh = pstar / 2.0
279      !
280      ahiu(:,:) = ahi0 * umask(:,:,1)            ! Ice eddy Diffusivity coefficients.
281      ahiv(:,:) = ahi0 * vmask(:,:,1)
282      !
283   END SUBROUTINE lim_dyn_init_2
284
285#else
286   !!----------------------------------------------------------------------
287   !!   Default option          Empty module       NO LIM 2.0 sea-ice model
288   !!----------------------------------------------------------------------
289CONTAINS
290   SUBROUTINE lim_dyn_2         ! Empty routine
291   END SUBROUTINE lim_dyn_2
292#endif 
293
294   !!======================================================================
295END MODULE limdyn_2
Note: See TracBrowser for help on using the repository browser.