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

source: branches/dev_002_LIM/NEMO/LIM_SRC_2/limdyn_2.F90 @ 823

Last change on this file since 823 was 823, checked in by rblod, 16 years ago

Final step to rename LIM_SRC in LIM_SRC_2
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.0 KB
Line 
1MODULE limdyn_2
2   !!======================================================================
3   !!                     ***  MODULE  limdyn_2  ***
4   !!   Sea-Ice dynamics : 
5   !!======================================================================
6#if defined key_lim2
7   !!----------------------------------------------------------------------
8   !!   'key_lim2' :                                  LIM 2.0 sea-ice model
9   !!----------------------------------------------------------------------
10   !!    lim_dyn_2      : computes ice velocities
11   !!    lim_dyn_init_2 : initialization and namelist read
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE phycst
15   USE in_out_manager  ! I/O manager
16   USE dom_ice_2
17   USE dom_oce         ! ocean space and time domain
18   USE ice_2
19   USE ice_oce
20   USE iceini_2
21   USE limistate_2
22   USE limrhg_2        ! ice rheology
23   USE lbclnk
24   USE lib_mpp
25   USE prtctl          ! Print control
26
27   IMPLICIT NONE
28   PRIVATE
29
30   !! * Accessibility
31   PUBLIC lim_dyn_2  ! routine called by ice_step
32
33   !! * Module variables
34   REAL(wp)  ::  rone    = 1.e0   ! constant value
35
36   !!----------------------------------------------------------------------
37   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)
38   !! $Header$
39   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
40   !!----------------------------------------------------------------------
41
42CONTAINS
43
44   SUBROUTINE lim_dyn_2( kt )
45      !!-------------------------------------------------------------------
46      !!               ***  ROUTINE lim_dyn_2  ***
47      !!               
48      !! ** Purpose :   compute ice velocity and ocean-ice stress
49      !!               
50      !! ** Method  :
51      !!
52      !! ** Action  : - Initialisation
53      !!              - Call of the dynamic routine for each hemisphere
54      !!              - computation of the stress at the ocean surface         
55      !!              - treatment of the case if no ice dynamic
56      !! History :
57      !!   1.0  !  01-04  (LIM)  Original code
58      !!   2.0  !  02-08  (C. Ethe, G. Madec)  F90, mpp
59      !!---------------------------------------------------------------------
60      INTEGER, INTENT(in) ::   kt     ! number of iteration
61
62      INTEGER ::   ji, jj             ! dummy loop indices
63      INTEGER ::   i_j1, i_jpj        ! Starting/ending j-indices for rheology
64      REAL(wp) ::   &
65         ztairx, ztairy,           &  ! tempory scalars
66         zsang , zmod,             &
67         ztglx , ztgly ,           &
68         zt11, zt12, zt21, zt22 ,  &
69         zustm, zsfrld, zsfrldm4,  &
70         zu_ice, zv_ice, ztair2
71      REAL(wp),DIMENSION(jpj) ::   &
72         zind,                     &  ! i-averaged indicator of sea-ice
73         zmsk                         ! i-averaged of tmask
74      !!---------------------------------------------------------------------
75
76      IF( kt == nit000  )   CALL lim_dyn_init_2   ! Initialization (first time-step only)
77     
78      IF ( ln_limdyn ) THEN
79
80         ! Mean ice and snow thicknesses.         
81         hsnm(:,:)  = ( 1.0 - frld(:,:) ) * hsnif(:,:)
82         hicm(:,:)  = ( 1.0 - frld(:,:) ) * hicif(:,:)
83
84         u_oce(:,:)  = u_io(:,:) * tmu(:,:)
85         v_oce(:,:)  = v_io(:,:) * tmu(:,:)
86       
87         !                                         ! Rheology (ice dynamics)
88         !                                         ! ========
89         
90         !  Define the j-limits where ice rheology is computed
91         ! ---------------------------------------------------
92         
93         IF( lk_mpp .OR. nbit_cmp == 1 ) THEN                    ! mpp: compute over the whole domain
94            i_j1 = 1   
95            i_jpj = jpj
96            IF(ln_ctl)    THEN
97               CALL prt_ctl_info('lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj)
98            ENDIF
99            CALL lim_rhg_2( i_j1, i_jpj )
100
101         ELSE                                 ! optimization of the computational area
102
103            DO jj = 1, jpj
104               zind(jj) = SUM( frld (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line
105               zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line
106   !!i         write(numout,*) narea, 'limdyn' , jj, zind(jj), zmsk(jj)
107            END DO
108
109            IF( l_jeq ) THEN                     ! local domain include both hemisphere
110               !                                 ! Rheology is computed in each hemisphere
111               !                                 ! only over the ice cover latitude strip
112               ! Northern hemisphere
113               i_j1  = njeq
114               i_jpj = jpj
115               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
116                  i_j1 = i_j1 + 1
117               END DO
118               i_j1 = MAX( 1, i_j1-1 )
119               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
120   
121               CALL lim_rhg_2( i_j1, i_jpj )
122   
123               ! Southern hemisphere
124               i_j1  =  1 
125               i_jpj = njeq
126               DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
127                  i_jpj = i_jpj - 1
128               END DO
129               i_jpj = MIN( jpj, i_jpj+2 )
130               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
131   
132               CALL lim_rhg_2( i_j1, i_jpj )
133   
134            ELSE                                 ! local domain extends over one hemisphere only
135               !                                 ! Rheology is computed only over the ice cover
136               !                                 ! latitude strip
137               i_j1  = 1
138               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
139                  i_j1 = i_j1 + 1
140               END DO
141               i_j1 = MAX( 1, i_j1-1 )
142   
143               i_jpj  = jpj
144               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
145                  i_jpj = i_jpj - 1
146               END DO
147               i_jpj = MIN( jpj, i_jpj+2)
148   
149               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
150   
151               CALL lim_rhg_2( i_j1, i_jpj )
152
153            ENDIF
154
155         ENDIF
156
157         IF(ln_ctl)   THEN
158            CALL prt_ctl(tab2d_1=u_oce , clinfo1=' lim_dyn  : u_oce :', tab2d_2=v_oce , clinfo2=' v_oce :')
159            CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :')
160         ENDIF
161         
162         !                                         ! Ice-Ocean stress
163         !                                         ! ================
164         DO jj = 2, jpjm1
165            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg
166            DO ji = 2, jpim1
167               ! computation of wind stress over ocean in X and Y direction
168#if defined key_coupled && defined key_lim_cp1
169               ztairx =  frld(ji-1,jj  ) * gtaux(ji-1,jj  ) + frld(ji,jj  ) * gtaux(ji,jj  )      &
170                  &    + frld(ji-1,jj-1) * gtaux(ji-1,jj-1) + frld(ji,jj-1) * gtaux(ji,jj-1)
171
172               ztairy =  frld(ji-1,jj  ) * gtauy(ji-1,jj  ) + frld(ji,jj  ) * gtauy(ji,jj  )      &
173                  &    + frld(ji-1,jj-1) * gtauy(ji-1,jj-1) + frld(ji,jj-1) * gtauy(ji,jj-1)
174#else
175               zsfrld  = frld(ji,jj) + frld(ji-1,jj) + frld(ji-1,jj-1) + frld(ji,jj-1)
176               ztairx  = zsfrld * gtaux(ji,jj)
177               ztairy  = zsfrld * gtauy(ji,jj)
178#endif
179               zsfrldm4 = 4 - frld(ji,jj) - frld(ji-1,jj) - frld(ji-1,jj-1) - frld(ji,jj-1)
180               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
181               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
182               zmod     = SQRT( zu_ice * zu_ice + zv_ice * zv_ice ) 
183               ztglx   = zsfrldm4 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice ) 
184               ztgly   = zsfrldm4 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice ) 
185
186               tio_u(ji,jj) = - ( ztairx + 1.0 * ztglx ) / ( 4 * rau0 )
187               tio_v(ji,jj) = - ( ztairy + 1.0 * ztgly ) / ( 4 * rau0 )
188            END DO
189         END DO
190         
191         ! computation of friction velocity
192         DO jj = 2, jpjm1
193            DO ji = 2, jpim1
194
195               zu_ice   = u_ice(ji-1,jj-1) - u_oce(ji-1,jj-1)
196               zv_ice   = v_ice(ji-1,jj-1) - v_oce(ji-1,jj-1)
197               zt11  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice )
198
199               zu_ice   = u_ice(ji-1,jj) - u_oce(ji-1,jj)
200               zv_ice   = v_ice(ji-1,jj) - v_oce(ji-1,jj)
201               zt12  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 
202
203               zu_ice   = u_ice(ji,jj-1) - u_oce(ji,jj-1)
204               zv_ice   = v_ice(ji,jj-1) - v_oce(ji,jj-1)
205               zt21  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 
206
207               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
208               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
209               zt22  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 
210
211               ztair2 = gtaux(ji,jj) * gtaux(ji,jj) + gtauy(ji,jj) * gtauy(ji,jj)
212
213               zustm =  ( 1 - frld(ji,jj) ) * 0.25 * ( zt11 + zt12 + zt21 + zt22 )        &
214                  &  +        frld(ji,jj)   * SQRT( ztair2 )
215
216               ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj)
217            END DO
218         END DO
219
220       ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
221                   
222          DO jj = 2, jpjm1
223             DO ji = 2, jpim1
224#if defined key_coupled && defined key_lim_cp1
225                tio_u(ji,jj) = - (  gtaux(ji  ,jj  ) + gtaux(ji-1,jj  )       &
226                   &              + gtaux(ji-1,jj-1) + gtaux(ji  ,jj-1) ) / ( 4 * rau0 )
227
228                tio_v(ji,jj) = - (  gtauy(ji  ,jj )  + gtauy(ji-1,jj  )       &
229                   &              + gtauy(ji-1,jj-1) + gtauy(ji  ,jj-1) ) / ( 4 * rau0 )
230#else
231                tio_u(ji,jj) = - gtaux(ji,jj) / rau0
232                tio_v(ji,jj) = - gtauy(ji,jj) / rau0 
233#endif
234                ztair2       = gtaux(ji,jj) * gtaux(ji,jj) + gtauy(ji,jj) * gtauy(ji,jj)
235                zustm        = SQRT( ztair2  )
236
237                ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj)
238            END DO
239         END DO
240
241      ENDIF
242
243      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
244      CALL lbc_lnk( tio_u, 'I', -1. )   ! I-point (i.e. ice U-V point)
245      CALL lbc_lnk( tio_v, 'I', -1. )   ! I-point (i.e. ice U-V point)
246
247      IF(ln_ctl) THEN
248            CALL prt_ctl(tab2d_1=tio_u , clinfo1=' lim_dyn  : tio_u :', tab2d_2=tio_v , clinfo2=' tio_v :')
249            CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :')
250      ENDIF
251
252   END SUBROUTINE lim_dyn_2
253
254
255   SUBROUTINE lim_dyn_init_2
256      !!-------------------------------------------------------------------
257      !!                  ***  ROUTINE lim_dyn_init_2  ***
258      !!
259      !! ** Purpose : Physical constants and parameters linked to the ice
260      !!      dynamics
261      !!
262      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
263      !!       parameter values called at the first timestep (nit000)
264      !!
265      !! ** input   :   Namelist namicedyn
266      !!
267      !! history :
268      !!  8.5  ! 03-08 (C. Ethe) original code
269      !!-------------------------------------------------------------------
270      NAMELIST/namicedyn/ epsd, alpha,     &
271         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
272         &                c_rhg, etamn, creepl, ecc, ahi0
273      !!-------------------------------------------------------------------
274
275      ! Define the initial parameters
276      ! -------------------------
277
278      ! Read Namelist namicedyn
279      REWIND ( numnam_ice )
280      READ   ( numnam_ice  , namicedyn )
281      IF(lwp) THEN
282         WRITE(numout,*)
283         WRITE(numout,*) 'lim_dyn_init_2: ice parameters for ice dynamics '
284         WRITE(numout,*) '~~~~~~~~~~~~~~'
285         WRITE(numout,*) '       tolerance parameter                              epsd   = ', epsd
286         WRITE(numout,*) '       coefficient for semi-implicit coriolis           alpha  = ', alpha
287         WRITE(numout,*) '       diffusion constant for dynamics                  dm     = ', dm
288         WRITE(numout,*) '       number of sub-time steps for relaxation          nbiter = ', nbiter
289         WRITE(numout,*) '       maximum number of iterations for relaxation      nbitdr = ', nbitdr
290         WRITE(numout,*) '       relaxation constant                              om     = ', om
291         WRITE(numout,*) '       maximum value for the residual of relaxation     resl   = ', resl
292         WRITE(numout,*) '       drag coefficient for oceanic stress              cw     = ', cw
293         WRITE(numout,*) '       turning angle for oceanic stress                 angvg  = ', angvg
294         WRITE(numout,*) '       first bulk-rheology parameter                    pstar  = ', pstar
295         WRITE(numout,*) '       second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
296         WRITE(numout,*) '       minimun value for viscosity                      etamn  = ', etamn
297         WRITE(numout,*) '       creep limit                                      creepl = ', creepl
298         WRITE(numout,*) '       eccentricity of the elliptical yield curve       ecc    = ', ecc
299         WRITE(numout,*) '       horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
300      ENDIF
301
302      !  Initialization
303      usecc2 = 1.0 / ( ecc * ecc )
304      rhoco  = rau0 * cw
305      angvg  = angvg * rad
306      sangvg = SIN( angvg )
307      cangvg = COS( angvg )
308      pstarh = pstar / 2.0
309      sdvt(:,:) = 0.e0
310
311      !  Diffusion coefficients.
312      ahiu(:,:) = ahi0 * umask(:,:,1)
313      ahiv(:,:) = ahi0 * vmask(:,:,1)
314
315   END SUBROUTINE lim_dyn_init_2
316
317#else
318   !!----------------------------------------------------------------------
319   !!   Default option          Empty module       NO LIM 2.0 sea-ice model
320   !!----------------------------------------------------------------------
321CONTAINS
322   SUBROUTINE lim_dyn_2         ! Empty routine
323   END SUBROUTINE lim_dyn_2
324#endif 
325
326   !!======================================================================
327END MODULE limdyn_2
Note: See TracBrowser for help on using the repository browser.