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

source: branches/DEV_1879_mpp_rep/NEMO/LIM_SRC_3/limdyn.F90 @ 1920

Last change on this file since 1920 was 1920, checked in by rblod, 14 years ago

Add modifications for mpp reproducibility, see ticket #677

  • Property svn:keywords set to Id
File size: 13.7 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 ice
19   USE par_ice
20   USE sbc_oce         ! Surface boundary condition: ocean fields
21   USE sbc_ice         ! Surface boundary condition: ice fields
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   !! $Id$
44   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46
47CONTAINS
48
49   SUBROUTINE lim_dyn( kt )
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      INTEGER, INTENT(in) ::   kt     ! number of iteration
68      !! * Local variables
69      INTEGER  ::   ji, jj, jl, ja    ! dummy loop indices
70      INTEGER  ::   i_j1, i_jpj       ! Starting/ending j-indices for rheology
71      REAL(wp) ::   zcoef             ! temporary scalar
72      REAL(wp), DIMENSION(jpj)     ::   zind           ! i-averaged indicator of sea-ice
73      REAL(wp), DIMENSION(jpj)     ::   zmsk           ! i-averaged of tmask
74      REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io   ! ice-ocean velocity
75      !!---------------------------------------------------------------------
76
77      IF( kt == nit000 .AND. lwp ) THEN
78         WRITE(numout,*) ' lim_dyn : Ice dynamics '
79         WRITE(numout,*) ' ~~~~~~~ '
80      ENDIF
81
82      IF( numit == nstart  )   CALL lim_dyn_init   ! Initialization (first time-step only)
83
84      IF ( ln_limdyn ) THEN
85
86         old_u_ice(:,:) = u_ice(:,:) * tmu(:,:)
87         old_v_ice(:,:) = v_ice(:,:) * tmv(:,:)
88
89         ! Rheology (ice dynamics)
90         ! ========
91
92         !  Define the j-limits where ice rheology is computed
93         ! ---------------------------------------------------
94
95         IF( lk_mpp .OR. lk_mpp_rep ) THEN                    ! mpp: compute over the whole domain
96            i_j1 = 1
97            i_jpj = jpj
98            IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
99            CALL lim_rhg( i_j1, i_jpj )
100         ELSE                                 ! optimization of the computational area
101
102            DO jj = 1, jpj
103               zind(jj) = SUM( 1.0 - at_i (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line
104               zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line
105            END DO
106
107            IF( l_jeq ) THEN                     ! local domain include both hemisphere
108               !                                 ! Rheology is computed in each hemisphere
109               !                                 ! only over the ice cover latitude strip
110               ! Northern hemisphere
111               i_j1  = njeq
112               i_jpj = jpj
113               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
114                  i_j1 = i_j1 + 1
115               END DO
116               i_j1 = MAX( 1, i_j1-2 )
117               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : NH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
118               CALL lim_rhg( i_j1, i_jpj )
119
120               ! Southern hemisphere
121               i_j1  =  1
122               i_jpj = njeq
123               DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
124                  i_jpj = i_jpj - 1
125               END DO
126               i_jpj = MIN( jpj, i_jpj+1 )
127               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : SH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
128
129       CALL lim_rhg( i_j1, i_jpj )
130
131    ELSE                                 ! local domain extends over one hemisphere only
132       !                                 ! Rheology is computed only over the ice cover
133       !                                 ! latitude strip
134       i_j1  = 1
135               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
136                  i_j1 = i_j1 + 1
137               END DO
138               i_j1 = MAX( 1, i_j1-2 )
139
140               i_jpj  = jpj
141               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
142                  i_jpj = i_jpj - 1
143               END DO
144               i_jpj = MIN( jpj, i_jpj+1)
145
146               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : one hemisphere:  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
147
148               CALL lim_rhg( i_j1, i_jpj )
149
150            ENDIF
151
152         ENDIF
153
154         ! computation of friction velocity
155         ! --------------------------------
156         ! ice-ocean velocity at U & V-points (u_ice v_ice at U- & V-points ; ssu_m, ssv_m at U- & V-points)
157         zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
158         zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
159         ! frictional velocity at T-point
160         DO jj = 2, jpjm1 
161            DO ji = fs_2, fs_jpim1   ! vector opt.
162               ust2s(ji,jj) = 0.5 * cw                                                          &
163                  &         * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
164                  &            + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj)
165            END DO
166         END DO
167         !
168      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
169         !
170         zcoef = SQRT( 0.5 ) / rau0
171         DO jj = 2, jpjm1
172            DO ji = fs_2, fs_jpim1   ! vector opt.
173               ust2s(ji,jj) = zcoef * tms(ji,jj) * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
174                  &                                     + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) )
175            END DO
176         END DO
177         !
178      ENDIF
179
180      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
181
182      IF(ln_ctl) THEN   ! Control print
183         CALL prt_ctl_info(' ')
184         CALL prt_ctl_info(' - Cell values : ')
185         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
186         CALL prt_ctl(tab2d_1=ust2s     , clinfo1=' lim_dyn  : ust2s     :')
187         CALL prt_ctl(tab2d_1=divu_i    , clinfo1=' lim_dyn  : divu_i    :')
188         CALL prt_ctl(tab2d_1=delta_i   , clinfo1=' lim_dyn  : delta_i   :')
189         CALL prt_ctl(tab2d_1=strength  , clinfo1=' lim_dyn  : strength  :')
190         CALL prt_ctl(tab2d_1=area      , clinfo1=' lim_dyn  : cell area :')
191         CALL prt_ctl(tab2d_1=at_i      , clinfo1=' lim_dyn  : at_i      :')
192         CALL prt_ctl(tab2d_1=vt_i      , clinfo1=' lim_dyn  : vt_i      :')
193         CALL prt_ctl(tab2d_1=vt_s      , clinfo1=' lim_dyn  : vt_s      :')
194         CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' lim_dyn  : stress1_i :')
195         CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' lim_dyn  : stress2_i :')
196         CALL prt_ctl(tab2d_1=stress12_i, clinfo1=' lim_dyn  : stress12_i:')
197         DO jl = 1, jpl
198            CALL prt_ctl_info(' ')
199            CALL prt_ctl_info(' - Category : ', ivar1=jl)
200            CALL prt_ctl_info('   ~~~~~~~~~~')
201            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_dyn  : a_i      : ')
202            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_i     : ')
203            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_s     : ')
204            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_dyn  : v_i      : ')
205            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_dyn  : v_s      : ')
206            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : e_s      : ')
207            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_dyn  : t_su     : ')
208            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : t_snow   : ')
209            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_dyn  : sm_i     : ')
210            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_dyn  : smv_i    : ')
211            DO ja = 1, nlay_i
212               CALL prt_ctl_info(' ')
213               CALL prt_ctl_info(' - Layer : ', ivar1=ja)
214               CALL prt_ctl_info('   ~~~~~~~')
215               CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : t_i      : ')
216               CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : e_i      : ')
217            END DO
218         END DO
219      ENDIF
220
221   END SUBROUTINE lim_dyn
222
223   SUBROUTINE lim_dyn_init
224      !!-------------------------------------------------------------------
225      !!                  ***  ROUTINE lim_dyn_init  ***
226      !!
227      !! ** Purpose : Physical constants and parameters linked to the ice
228      !!      dynamics
229      !!
230      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
231      !!       parameter values called at the first timestep (nit000)
232      !!
233      !! ** input   :   Namelist namicedyn
234      !!
235      !! history :
236      !!  8.5  ! 03-08 (C. Ethe) original code
237      !!  9.0  ! 07-03 (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle)
238      !!               EVP-Cgrid-LIM3
239      !!-------------------------------------------------------------------
240      NAMELIST/namicedyn/ epsd, alpha,     &
241         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
242         &                c_rhg, etamn, creepl, ecc, ahi0, &
243         &                nevp, telast, alphaevp
244      !!-------------------------------------------------------------------
245
246      ! Define the initial parameters
247      ! -------------------------
248
249      ! Read Namelist namicedyn
250      REWIND ( numnam_ice )
251      READ   ( numnam_ice  , namicedyn )
252      IF(lwp) THEN
253         WRITE(numout,*)
254         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
255         WRITE(numout,*) '~~~~~~~~~~~~'
256         WRITE(numout,*) '   tolerance parameter                              epsd   = ', epsd
257         WRITE(numout,*) '   coefficient for semi-implicit coriolis           alpha  = ', alpha
258         WRITE(numout,*) '   diffusion constant for dynamics                  dm     = ', dm
259         WRITE(numout,*) '   number of sub-time steps for relaxation          nbiter = ', nbiter
260         WRITE(numout,*) '   maximum number of iterations for relaxation      nbitdr = ', nbitdr
261         WRITE(numout,*) '   relaxation constant                              om     = ', om
262         WRITE(numout,*) '   maximum value for the residual of relaxation     resl   = ', resl
263         WRITE(numout,*) '   drag coefficient for oceanic stress              cw     = ', cw
264         WRITE(numout,*) '   turning angle for oceanic stress                 angvg  = ', angvg
265         WRITE(numout,*) '   first bulk-rheology parameter                    pstar  = ', pstar
266         WRITE(numout,*) '   second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
267         WRITE(numout,*) '   minimun value for viscosity                      etamn  = ', etamn
268         WRITE(numout,*) '   creep limit                                      creepl = ', creepl
269         WRITE(numout,*) '   eccentricity of the elliptical yield curve       ecc    = ', ecc
270         WRITE(numout,*) '   horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
271         WRITE(numout,*) '   number of iterations for subcycling              nevp   = ', nevp
272         WRITE(numout,*) '   timescale for elastic waves                      telast = ', telast
273         WRITE(numout,*) '   coefficient for the solution of int. stresses  alphaevp = ', alphaevp
274
275      ENDIF
276
277      usecc2 = 1.0 / ( ecc * ecc )
278      rhoco  = rau0 * cw
279      angvg  = angvg * rad
280      sangvg = SIN( angvg )
281      cangvg = COS( angvg )
282      pstarh = pstar / 2.0
283
284      !  Diffusion coefficients.
285      ahiu(:,:) = ahi0 * umask(:,:,1)
286      ahiv(:,:) = ahi0 * vmask(:,:,1)
287
288   END SUBROUTINE lim_dyn_init
289
290#else
291   !!----------------------------------------------------------------------
292   !!   Default option          Empty module           NO LIM sea-ice model
293   !!----------------------------------------------------------------------
294CONTAINS
295   SUBROUTINE lim_dyn         ! Empty routine
296   END SUBROUTINE lim_dyn
297#endif 
298
299   !!======================================================================
300END MODULE limdyn
Note: See TracBrowser for help on using the repository browser.