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

source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90 @ 3963

Last change on this file since 3963 was 3963, checked in by clem, 11 years ago

bugs correction + creation of glob_max and glob_min in lib_fortran.F90, see ticket:#1116
  • Property svn:keywords set to Id
File size: 17.4 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
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   lim_dyn   ! routine called by ice_step
36
37   !! * Substitutions
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
41   !! $Id$
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE lim_dyn( kt )
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      !!------------------------------------------------------------------------------------
59      INTEGER, INTENT(in) ::   kt     ! number of iteration
60      !!
61      INTEGER  ::   ji, jj, jl, ja    ! dummy loop indices
62      INTEGER  ::   i_j1, i_jpj       ! Starting/ending j-indices for rheology
63      REAL(wp) ::   zcoef             ! local scalar
64      REAL(wp), POINTER, DIMENSION(:)   ::   zind           ! i-averaged indicator of sea-ice
65      REAL(wp), POINTER, DIMENSION(:)   ::   zmsk           ! i-averaged of tmask
66      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_io, zv_io   ! ice-ocean velocity
67      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)
68      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
69     !!---------------------------------------------------------------------
70
71      CALL wrk_alloc( jpi, jpj, zu_io, zv_io )
72      CALL wrk_alloc( jpj, zind, zmsk )
73
74      ! -------------------------------
75      !- check conservation (C Rousset)
76      IF (ln_limdiahsb) THEN
77         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
78         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
79         zchk_fw_b  = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) )
80         zchk_fs_b  = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) )
81      ENDIF
82      !- check conservation (C Rousset)
83      ! -------------------------------
84
85      IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only)
86
87      IF( ln_limdyn ) THEN
88         !
89         old_u_ice(:,:) = u_ice(:,:) * tmu(:,:)
90         old_v_ice(:,:) = v_ice(:,:) * tmv(:,:)
91
92         ! Rheology (ice dynamics)
93         ! ========
94
95         !  Define the j-limits where ice rheology is computed
96         ! ---------------------------------------------------
97
98         IF( lk_mpp .OR. lk_mpp_rep ) THEN                    ! mpp: compute over the whole domain
99            i_j1 = 1
100            i_jpj = jpj
101            IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
102            CALL lim_rhg( i_j1, i_jpj )
103         ELSE                                 ! optimization of the computational area
104            !
105            DO jj = 1, jpj
106               zind(jj) = SUM( 1.0 - at_i(:,jj) )   ! = REAL(jpj) if ocean everywhere on a j-line
107               zmsk(jj) = SUM( tmask(:,jj,1)    )   ! = 0         if land  everywhere on a j-line
108            END DO
109
110            IF( l_jeq ) THEN                     ! local domain include both hemisphere
111               !                                 ! Rheology is computed in each hemisphere
112               !                                 ! only over the ice cover latitude strip
113               ! Northern hemisphere
114               i_j1  = njeq
115               i_jpj = jpj
116               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
117                  i_j1 = i_j1 + 1
118               END DO
119               i_j1 = MAX( 1, i_j1-2 )
120               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : NH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
121               CALL lim_rhg( 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+1 )
130               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : SH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
131               !
132               CALL lim_rhg( 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-2 )
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+1)
148               !
149               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : one hemisphere:  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
150               !
151               CALL lim_rhg( i_j1, i_jpj )
152               !
153            ENDIF
154            !
155         ENDIF
156
157         ! computation of friction velocity
158         ! --------------------------------
159         ! ice-ocean velocity at U & V-points (u_ice v_ice at U- & V-points ; ssu_m, ssv_m at U- & V-points)
160         zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
161         zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
162         ! frictional velocity at T-point
163         zcoef = 0.5_wp * cw
164         DO jj = 2, jpjm1 
165            DO ji = fs_2, fs_jpim1   ! vector opt.
166               ust2s(ji,jj) = zcoef * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
167                  &                    + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj)
168            END DO
169         END DO
170         !
171      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
172         !
173         zcoef = SQRT( 0.5_wp ) / rau0
174         DO jj = 2, jpjm1
175            DO ji = fs_2, fs_jpim1   ! vector opt.
176               ust2s(ji,jj) = zcoef * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
177                  &                        + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1)   ) * tms(ji,jj)
178            END DO
179         END DO
180         !
181      ENDIF
182      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
183
184      IF(ln_ctl) THEN   ! Control print
185         CALL prt_ctl_info(' ')
186         CALL prt_ctl_info(' - Cell values : ')
187         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
188         CALL prt_ctl(tab2d_1=ust2s     , clinfo1=' lim_dyn  : ust2s     :')
189         CALL prt_ctl(tab2d_1=divu_i    , clinfo1=' lim_dyn  : divu_i    :')
190         CALL prt_ctl(tab2d_1=delta_i   , clinfo1=' lim_dyn  : delta_i   :')
191         CALL prt_ctl(tab2d_1=strength  , clinfo1=' lim_dyn  : strength  :')
192         CALL prt_ctl(tab2d_1=area      , clinfo1=' lim_dyn  : cell area :')
193         CALL prt_ctl(tab2d_1=at_i      , clinfo1=' lim_dyn  : at_i      :')
194         CALL prt_ctl(tab2d_1=vt_i      , clinfo1=' lim_dyn  : vt_i      :')
195         CALL prt_ctl(tab2d_1=vt_s      , clinfo1=' lim_dyn  : vt_s      :')
196         CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' lim_dyn  : stress1_i :')
197         CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' lim_dyn  : stress2_i :')
198         CALL prt_ctl(tab2d_1=stress12_i, clinfo1=' lim_dyn  : stress12_i:')
199         DO jl = 1, jpl
200            CALL prt_ctl_info(' ')
201            CALL prt_ctl_info(' - Category : ', ivar1=jl)
202            CALL prt_ctl_info('   ~~~~~~~~~~')
203            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_dyn  : a_i      : ')
204            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_i     : ')
205            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_s     : ')
206            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_dyn  : v_i      : ')
207            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_dyn  : v_s      : ')
208            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : e_s      : ')
209            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_dyn  : t_su     : ')
210            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : t_snow   : ')
211            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_dyn  : sm_i     : ')
212            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_dyn  : smv_i    : ')
213            DO ja = 1, nlay_i
214               CALL prt_ctl_info(' ')
215               CALL prt_ctl_info(' - Layer : ', ivar1=ja)
216               CALL prt_ctl_info('   ~~~~~~~')
217               CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : t_i      : ')
218               CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : e_i      : ')
219            END DO
220         END DO
221      ENDIF
222      !
223
224      ! -------------------------------
225      !- check conservation (C Rousset)
226      IF (ln_limdiahsb) THEN
227         !INTEGER                                 ::   numhsb
228         !CHARACTER (len=32) ::   cl_name  ! output file name
229         !cl_name    = 'heat_salt_volume_budgets.txt'                         ! name of output file
230
231         zchk_fs  = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
232         zchk_fw  = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
233 
234         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice
235         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )
236
237         zchk_vmin = glob_min(v_i)
238         zchk_amax = glob_max(SUM(a_i,dim=3))
239         zchk_amin = glob_min(a_i)
240
241         IF(lwp) THEN
242            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limdyn) = ',(zchk_v_i * 86400.)
243            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limdyn) = ',(zchk_smv * 86400.)
244            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limdyn) = ',(zchk_vmin * 1.e-3)
245            !IF ( zchk_amax >  amax+1.e-10   ) WRITE(numout,*) 'violation a_i>amax            (limdyn) = ',zchk_amax
246            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limdyn) = ',zchk_amin
247         ENDIF
248         !CALL ctl_opn( numhsb , cl_name , 'UNKNOWN' , 'FORMATTED' , 'SEQUENTIAL' , 1 , numout , lwp , 1 )
249         !
250         !WRITE( numhsb, 9010 ) "kt   |     heat content budget     |            salt content budget             ",   &
251         !     &                                                  "|            volume budget (ssh)             ",   &
252         !     &                                                  "|            volume budget (e3t)             "
253         !WRITE( numhsb, 9010 ) "     |      [C]         [W/m2]     |     [psu]        [mmm/s]          [SV]     ",   &
254         !     &                                                  "|     [m3]         [mmm/s]          [SV]     ",   &
255         !     &                                                  "|     [m3]         [mmm/s]          [SV]     "
256         !IF ( kt == nitend ) CLOSE( numhsb )
257         
258!9010     FORMAT(A80,A45,A45)
259      ENDIF
260      !- check conservation (C Rousset)
261      ! -------------------------------
262
263      CALL wrk_dealloc( jpi, jpj, zu_io, zv_io )
264      CALL wrk_dealloc( jpj, zind, zmsk )
265      !
266   END SUBROUTINE lim_dyn
267
268
269   SUBROUTINE lim_dyn_init
270      !!-------------------------------------------------------------------
271      !!                  ***  ROUTINE lim_dyn_init  ***
272      !!
273      !! ** Purpose : Physical constants and parameters linked to the ice
274      !!      dynamics
275      !!
276      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
277      !!       parameter values called at the first timestep (nit000)
278      !!
279      !! ** input   :   Namelist namicedyn
280      !!-------------------------------------------------------------------
281      NAMELIST/namicedyn/ epsd, alpha,     &
282         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
283         &                c_rhg, etamn, creepl, ecc, ahi0, &
284         &                nevp, telast, alphaevp, hminrhg
285      !!-------------------------------------------------------------------
286
287      REWIND( numnam_ice )                ! Read Namelist namicedyn
288      READ  ( numnam_ice  , namicedyn )
289     
290      IF(lwp) THEN                        ! control print
291         WRITE(numout,*)
292         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
293         WRITE(numout,*) '~~~~~~~~~~~~'
294         WRITE(numout,*) '   tolerance parameter                              epsd   = ', epsd
295         WRITE(numout,*) '   coefficient for semi-implicit coriolis           alpha  = ', alpha
296         WRITE(numout,*) '   diffusion constant for dynamics                  dm     = ', dm
297         WRITE(numout,*) '   number of sub-time steps for relaxation          nbiter = ', nbiter
298         WRITE(numout,*) '   maximum number of iterations for relaxation      nbitdr = ', nbitdr
299         WRITE(numout,*) '   relaxation constant                              om     = ', om
300         WRITE(numout,*) '   maximum value for the residual of relaxation     resl   = ', resl
301         WRITE(numout,*) '   drag coefficient for oceanic stress              cw     = ', cw
302         WRITE(numout,*) '   turning angle for oceanic stress                 angvg  = ', angvg
303         WRITE(numout,*) '   first bulk-rheology parameter                    pstar  = ', pstar
304         WRITE(numout,*) '   second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
305         WRITE(numout,*) '   minimun value for viscosity                      etamn  = ', etamn
306         WRITE(numout,*) '   creep limit                                      creepl = ', creepl
307         WRITE(numout,*) '   eccentricity of the elliptical yield curve       ecc    = ', ecc
308         WRITE(numout,*) '   horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
309         WRITE(numout,*) '   number of iterations for subcycling              nevp   = ', nevp
310         WRITE(numout,*) '   timescale for elastic waves                      telast = ', telast
311         WRITE(numout,*) '   coefficient for the solution of int. stresses  alphaevp = ', alphaevp
312         WRITE(numout,*) '   min ice thickness for rheology calculations     hminrhg = ', hminrhg
313      ENDIF
314      !
315      IF( angvg /= 0._wp ) THEN
316         CALL ctl_warn( 'lim_dyn_init: turning angle for oceanic stress not properly coded for EVP ',   &
317            &           '(see limsbc module). We force  angvg = 0._wp'  )
318         angvg = 0._wp
319      ENDIF
320     
321      usecc2 = 1._wp / ( ecc * ecc )
322      rhoco  = rau0  * cw
323      angvg  = angvg * rad
324      sangvg = SIN( angvg )
325      cangvg = COS( angvg )
326      pstarh = pstar * 0.5_wp
327
328      !  Diffusion coefficients.
329      ahiu(:,:) = ahi0 * umask(:,:,1)
330      ahiv(:,:) = ahi0 * vmask(:,:,1)
331      !
332   END SUBROUTINE lim_dyn_init
333
334#else
335   !!----------------------------------------------------------------------
336   !!   Default option          Empty module           NO LIM sea-ice model
337   !!----------------------------------------------------------------------
338CONTAINS
339   SUBROUTINE lim_dyn         ! Empty routine
340   END SUBROUTINE lim_dyn
341#endif 
342
343   !!======================================================================
344END MODULE limdyn
Note: See TracBrowser for help on using the repository browser.