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

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90 @ 4624

Last change on this file since 4624 was 4624, checked in by acc, 10 years ago

#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.
  • Property svn:keywords set to Id
File size: 17.0 KB
Line 
1MODULE limdyn
2   !!======================================================================
3   !!                     ***  MODULE  limdyn  ***
4   !!   Sea-Ice dynamics : 
5   !!======================================================================
6   !! history :  1.0  ! 2002-08  (C. Ethe, G. Madec)  original VP code
7   !!            3.0  ! 2007-03  (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle)  LIM3: EVP-Cgrid
8   !!            4.0  ! 2011-02  (G. Madec) dynamical allocation
9   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3' :                                 LIM3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!    lim_dyn      : computes ice velocities
15   !!    lim_dyn_init : initialization and namelist read
16   !!----------------------------------------------------------------------
17   USE phycst           ! physical constants
18   USE dom_oce          ! ocean space and time domain
19   USE sbc_oce          ! Surface boundary condition: ocean fields
20   USE sbc_ice          ! Surface boundary condition: ice   fields
21   USE ice              ! LIM-3 variables
22   USE par_ice          ! LIM-3 parameters
23   USE dom_ice          ! LIM-3 domain
24   USE limrhg           ! LIM-3 rheology
25   USE lbclnk           ! lateral boundary conditions - MPP exchanges
26   USE lib_mpp          ! MPP library
27   USE wrk_nemo         ! work arrays
28   USE in_out_manager   ! I/O manager
29   USE prtctl           ! Print control
30   USE lib_fortran      ! glob_sum
31   USE timing          ! Timing
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   lim_dyn   ! routine called by ice_step
37
38   !! * Substitutions
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
42   !! $Id$
43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE lim_dyn( kt )
48      !!-------------------------------------------------------------------
49      !!               ***  ROUTINE lim_dyn  ***
50      !!               
51      !! ** Purpose :   compute ice velocity and ocean-ice stress
52      !!               
53      !! ** Method  :
54      !!
55      !! ** Action  : - Initialisation
56      !!              - Call of the dynamic routine for each hemisphere
57      !!              - computation of the stress at the ocean surface         
58      !!              - treatment of the case if no ice dynamic
59      !!------------------------------------------------------------------------------------
60      INTEGER, INTENT(in) ::   kt     ! number of iteration
61      !!
62      INTEGER  ::   ji, jj, jl, ja    ! dummy loop indices
63      INTEGER  ::   i_j1, i_jpj       ! Starting/ending j-indices for rheology
64      REAL(wp) ::   zcoef             ! local scalar
65      REAL(wp), POINTER, DIMENSION(:)   ::   zind           ! i-averaged indicator of sea-ice
66      REAL(wp), POINTER, DIMENSION(:)   ::   zmsk           ! i-averaged of tmask
67      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_io, zv_io   ! ice-ocean velocity
68      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
69      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
70     !!---------------------------------------------------------------------
71
72      IF( nn_timing == 1 )  CALL timing_start('limdyn')
73
74      CALL wrk_alloc( jpi, jpj, zu_io, zv_io )
75      CALL wrk_alloc( jpj, zind, zmsk )
76
77      ! -------------------------------
78      !- check conservation (C Rousset)
79      IF (ln_limdiahsb) THEN
80         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
81         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
82         zchk_fw_b  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) )
83         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) )
84      ENDIF
85      !- check conservation (C Rousset)
86      ! -------------------------------
87
88      IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only)
89
90      IF( ln_limdyn ) THEN
91         !
92         old_u_ice(:,:) = u_ice(:,:) * tmu(:,:)
93         old_v_ice(:,:) = v_ice(:,:) * tmv(:,:)
94
95         ! Rheology (ice dynamics)
96         ! ========
97
98         !  Define the j-limits where ice rheology is computed
99         ! ---------------------------------------------------
100
101         IF( lk_mpp .OR. lk_mpp_rep ) THEN                    ! mpp: compute over the whole domain
102            i_j1 = 1
103            i_jpj = jpj
104            IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
105            CALL lim_rhg( i_j1, i_jpj )
106         ELSE                                 ! optimization of the computational area
107            !
108            DO jj = 1, jpj
109               zind(jj) = SUM( 1.0 - at_i(:,jj) )   ! = REAL(jpj) if ocean everywhere on a j-line
110               zmsk(jj) = SUM( tmask(:,jj,1)    )   ! = 0         if land  everywhere on a j-line
111            END DO
112
113            IF( l_jeq ) THEN                     ! local domain include both hemisphere
114               !                                 ! Rheology is computed in each hemisphere
115               !                                 ! only over the ice cover latitude strip
116               ! Northern hemisphere
117               i_j1  = njeq
118               i_jpj = jpj
119               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
120                  i_j1 = i_j1 + 1
121               END DO
122               i_j1 = MAX( 1, i_j1-2 )
123               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : NH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
124               CALL lim_rhg( i_j1, i_jpj )
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               i_jpj = MIN( jpj, i_jpj+1 )
133               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : SH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
134               !
135               CALL lim_rhg( i_j1, i_jpj )
136               !
137            ELSE                                 ! local domain extends over one hemisphere only
138               !                                 ! Rheology is computed only over the ice cover
139               !                                 ! latitude strip
140               i_j1  = 1
141               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
142                  i_j1 = i_j1 + 1
143               END DO
144               i_j1 = MAX( 1, i_j1-2 )
145
146               i_jpj  = jpj
147               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
148                  i_jpj = i_jpj - 1
149               END DO
150               i_jpj = MIN( jpj, i_jpj+1)
151               !
152               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : one hemisphere:  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
153               !
154               CALL lim_rhg( i_j1, i_jpj )
155               !
156            ENDIF
157            !
158         ENDIF
159
160         ! computation of friction velocity
161         ! --------------------------------
162         ! ice-ocean velocity at U & V-points (u_ice v_ice at U- & V-points ; ssu_m, ssv_m at U- & V-points)
163         zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
164         zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
165         ! frictional velocity at T-point
166         zcoef = 0.5_wp * cw
167         DO jj = 2, jpjm1 
168            DO ji = fs_2, fs_jpim1   ! vector opt.
169               ust2s(ji,jj) = zcoef * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
170                  &                    + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj)
171            END DO
172         END DO
173         !
174      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
175         !
176         zcoef = SQRT( 0.5_wp ) / rau0
177         DO jj = 2, jpjm1
178            DO ji = fs_2, fs_jpim1   ! vector opt.
179               ust2s(ji,jj) = zcoef * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
180                  &                        + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1)   ) * tms(ji,jj)
181            END DO
182         END DO
183         !
184      ENDIF
185      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
186
187      IF(ln_ctl) THEN   ! Control print
188         CALL prt_ctl_info(' ')
189         CALL prt_ctl_info(' - Cell values : ')
190         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
191         CALL prt_ctl(tab2d_1=ust2s     , clinfo1=' lim_dyn  : ust2s     :')
192         CALL prt_ctl(tab2d_1=divu_i    , clinfo1=' lim_dyn  : divu_i    :')
193         CALL prt_ctl(tab2d_1=delta_i   , clinfo1=' lim_dyn  : delta_i   :')
194         CALL prt_ctl(tab2d_1=strength  , clinfo1=' lim_dyn  : strength  :')
195         CALL prt_ctl(tab2d_1=area      , clinfo1=' lim_dyn  : cell area :')
196         CALL prt_ctl(tab2d_1=at_i      , clinfo1=' lim_dyn  : at_i      :')
197         CALL prt_ctl(tab2d_1=vt_i      , clinfo1=' lim_dyn  : vt_i      :')
198         CALL prt_ctl(tab2d_1=vt_s      , clinfo1=' lim_dyn  : vt_s      :')
199         CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' lim_dyn  : stress1_i :')
200         CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' lim_dyn  : stress2_i :')
201         CALL prt_ctl(tab2d_1=stress12_i, clinfo1=' lim_dyn  : stress12_i:')
202         DO jl = 1, jpl
203            CALL prt_ctl_info(' ')
204            CALL prt_ctl_info(' - Category : ', ivar1=jl)
205            CALL prt_ctl_info('   ~~~~~~~~~~')
206            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_dyn  : a_i      : ')
207            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_i     : ')
208            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_s     : ')
209            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_dyn  : v_i      : ')
210            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_dyn  : v_s      : ')
211            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : e_s      : ')
212            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_dyn  : t_su     : ')
213            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : t_snow   : ')
214            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_dyn  : sm_i     : ')
215            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_dyn  : smv_i    : ')
216            DO ja = 1, nlay_i
217               CALL prt_ctl_info(' ')
218               CALL prt_ctl_info(' - Layer : ', ivar1=ja)
219               CALL prt_ctl_info('   ~~~~~~~')
220               CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : t_i      : ')
221               CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : e_i      : ')
222            END DO
223         END DO
224      ENDIF
225      !
226      ! -------------------------------
227      !- check conservation (C Rousset)
228      IF (ln_limdiahsb) THEN
229         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
230         zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
231 
232         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice
233         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )
234
235         zchk_vmin = glob_min(v_i)
236         zchk_amax = glob_max(SUM(a_i,dim=3))
237         zchk_amin = glob_min(a_i)
238
239         IF(lwp) THEN
240            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limdyn) = ',(zchk_v_i * rday)
241            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limdyn) = ',(zchk_smv * rday)
242            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limdyn) = ',(zchk_vmin * 1.e-3)
243            !IF ( zchk_amax >  amax+1.e-10   ) WRITE(numout,*) 'violation a_i>amax            (limdyn) = ',zchk_amax
244            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limdyn) = ',zchk_amin
245         ENDIF
246      ENDIF
247      !- check conservation (C Rousset)
248      ! -------------------------------
249
250      CALL wrk_dealloc( jpi, jpj, zu_io, zv_io )
251      CALL wrk_dealloc( jpj, zind, zmsk )
252      !
253      IF( nn_timing == 1 )  CALL timing_stop('limdyn')
254
255   END SUBROUTINE lim_dyn
256
257
258   SUBROUTINE lim_dyn_init
259      !!-------------------------------------------------------------------
260      !!                  ***  ROUTINE lim_dyn_init  ***
261      !!
262      !! ** Purpose : Physical constants and parameters linked to the ice
263      !!      dynamics
264      !!
265      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
266      !!       parameter values called at the first timestep (nit000)
267      !!
268      !! ** input   :   Namelist namicedyn
269      !!-------------------------------------------------------------------
270      INTEGER  ::   ios                 ! Local integer output status for namelist read
271      NAMELIST/namicedyn/ epsd, alpha,     &
272         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
273         &                c_rhg, etamn, creepl, ecc, ahi0, &
274         &                nevp, telast, alphaevp, hminrhg
275      !!-------------------------------------------------------------------
276
277      REWIND( numnam_ice_ref )              ! Namelist namicedyn in reference namelist : Ice dynamics
278      READ  ( numnam_ice_ref, namicedyn, IOSTAT = ios, ERR = 901)
279901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in reference namelist', lwp )
280
281      REWIND( numnam_ice_cfg )              ! Namelist namicedyn in configuration namelist : Ice dynamics
282      READ  ( numnam_ice_cfg, namicedyn, IOSTAT = ios, ERR = 902 )
283902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in configuration namelist', lwp )
284      IF(lwm) WRITE ( numoni, namicedyn )
285     
286      IF(lwp) THEN                        ! control print
287         WRITE(numout,*)
288         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
289         WRITE(numout,*) '~~~~~~~~~~~~'
290         WRITE(numout,*) '   tolerance parameter                              epsd   = ', epsd
291         WRITE(numout,*) '   coefficient for semi-implicit coriolis           alpha  = ', alpha
292         WRITE(numout,*) '   diffusion constant for dynamics                  dm     = ', dm
293         WRITE(numout,*) '   number of sub-time steps for relaxation          nbiter = ', nbiter
294         WRITE(numout,*) '   maximum number of iterations for relaxation      nbitdr = ', nbitdr
295         WRITE(numout,*) '   relaxation constant                              om     = ', om
296         WRITE(numout,*) '   maximum value for the residual of relaxation     resl   = ', resl
297         WRITE(numout,*) '   drag coefficient for oceanic stress              cw     = ', cw
298         WRITE(numout,*) '   turning angle for oceanic stress                 angvg  = ', angvg
299         WRITE(numout,*) '   first bulk-rheology parameter                    pstar  = ', pstar
300         WRITE(numout,*) '   second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
301         WRITE(numout,*) '   minimun value for viscosity                      etamn  = ', etamn
302         WRITE(numout,*) '   creep limit                                      creepl = ', creepl
303         WRITE(numout,*) '   eccentricity of the elliptical yield curve       ecc    = ', ecc
304         WRITE(numout,*) '   horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
305         WRITE(numout,*) '   number of iterations for subcycling              nevp   = ', nevp
306         WRITE(numout,*) '   timescale for elastic waves                      telast = ', telast
307         WRITE(numout,*) '   coefficient for the solution of int. stresses  alphaevp = ', alphaevp
308         WRITE(numout,*) '   min ice thickness for rheology calculations     hminrhg = ', hminrhg
309      ENDIF
310      !
311      IF( angvg /= 0._wp ) THEN
312         CALL ctl_warn( 'lim_dyn_init: turning angle for oceanic stress not properly coded for EVP ',   &
313            &           '(see limsbc module). We force  angvg = 0._wp'  )
314         angvg = 0._wp
315      ENDIF
316     
317      usecc2 = 1._wp / ( ecc * ecc )
318      rhoco  = rau0  * cw
319      angvg  = angvg * rad
320      sangvg = SIN( angvg )
321      cangvg = COS( angvg )
322      pstarh = pstar * 0.5_wp
323
324      !  Diffusion coefficients.
325      ahiu(:,:) = ahi0 * umask(:,:,1)
326      ahiv(:,:) = ahi0 * vmask(:,:,1)
327      !
328   END SUBROUTINE lim_dyn_init
329
330#else
331   !!----------------------------------------------------------------------
332   !!   Default option          Empty module           NO LIM sea-ice model
333   !!----------------------------------------------------------------------
334CONTAINS
335   SUBROUTINE lim_dyn         ! Empty routine
336   END SUBROUTINE lim_dyn
337#endif 
338
339   !!======================================================================
340END MODULE limdyn
Note: See TracBrowser for help on using the repository browser.