source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 5 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

  • Property svn:keywords set to Id
File size: 16.2 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   !!            3.5  ! 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 dom_ice          ! LIM-3 domain
23   USE limrhg           ! LIM-3 rheology
24   USE lbclnk           ! lateral boundary conditions - MPP exchanges
25   USE lib_mpp          ! MPP library
26   USE wrk_nemo         ! work arrays
27   USE in_out_manager   ! I/O manager
28   USE prtctl           ! Print control
29   USE lib_fortran      ! glob_sum
30   USE timing          ! Timing
31   USE limcons        ! conservation tests
32   USE limvar
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   lim_dyn   ! routine called by ice_step
38
39   !! * Substitutions
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
43   !! $Id$
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE lim_dyn( kt )
49      !!-------------------------------------------------------------------
50      !!               ***  ROUTINE lim_dyn  ***
51      !!               
52      !! ** Purpose :   compute ice velocity and ocean-ice stress
53      !!               
54      !! ** Method  :
55      !!
56      !! ** Action  : - Initialisation
57      !!              - Call of the dynamic routine for each hemisphere
58      !!              - computation of the stress at the ocean surface         
59      !!              - treatment of the case if no ice dynamic
60      !!------------------------------------------------------------------------------------
61      INTEGER, INTENT(in) ::   kt     ! number of iteration
62      !!
63      INTEGER  ::   ji, jj, jl, ja    ! dummy loop indices
64      INTEGER  ::   i_j1, i_jpj       ! Starting/ending j-indices for rheology
65      REAL(wp) ::   zcoef             ! local scalar
66      REAL(wp), POINTER, DIMENSION(:)   ::   zswitch        ! i-averaged indicator of sea-ice
67      REAL(wp), POINTER, DIMENSION(:)   ::   zmsk           ! i-averaged of tmask
68      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_io, zv_io   ! ice-ocean velocity
69      !
70      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
71     !!---------------------------------------------------------------------
72
73      IF( nn_timing == 1 )  CALL timing_start('limdyn')
74
75      CALL wrk_alloc( jpi, jpj, zu_io, zv_io )
76      CALL wrk_alloc( jpj, zswitch, zmsk )
77
78      CALL lim_var_agg(1)             ! aggregate ice categories
79
80      IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only)
81
82      IF( ln_limdyn ) THEN
83         !
84         ! conservation test
85         IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
86
87         u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1)
88         v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1)
89
90         ! Rheology (ice dynamics)
91         ! ========
92
93         !  Define the j-limits where ice rheology is computed
94         ! ---------------------------------------------------
95
96         IF( lk_mpp .OR. lk_mpp_rep ) THEN                    ! mpp: compute over the whole domain
97            i_j1 = 1
98            i_jpj = jpj
99            IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
100            CALL lim_rhg( i_j1, i_jpj )
101         ELSE                                 ! optimization of the computational area
102            !
103            DO jj = 1, jpj
104               zswitch(jj) = SUM( 1.0 - at_i(:,jj) )   ! = REAL(jpj) if ocean everywhere on a j-line
105               zmsk   (jj) = SUM( tmask(:,jj,1)    )   ! = 0         if land  everywhere on a j-line
106            END DO
107
108            IF( l_jeq ) THEN                     ! local domain include both hemisphere
109               !                                 ! Rheology is computed in each hemisphere
110               !                                 ! only over the ice cover latitude strip
111               ! Northern hemisphere
112               i_j1  = njeq
113               i_jpj = jpj
114               DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
115                  i_j1 = i_j1 + 1
116               END DO
117               i_j1 = MAX( 1, i_j1-2 )
118               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : NH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
119               CALL lim_rhg( i_j1, i_jpj )
120               !
121               ! Southern hemisphere
122               i_j1  =  1
123               i_jpj = njeq
124               DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
125                  i_jpj = i_jpj - 1
126               END DO
127               i_jpj = MIN( jpj, i_jpj+1 )
128               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : SH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
129               !
130               CALL lim_rhg( i_j1, i_jpj )
131               !
132            ELSE                                 ! local domain extends over one hemisphere only
133               !                                 ! Rheology is computed only over the ice cover
134               !                                 ! latitude strip
135               i_j1  = 1
136               DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
137                  i_j1 = i_j1 + 1
138               END DO
139               i_j1 = MAX( 1, i_j1-2 )
140
141               i_jpj  = jpj
142               DO WHILE ( i_jpj >= 1  .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
143                  i_jpj = i_jpj - 1
144               END DO
145               i_jpj = MIN( jpj, i_jpj+1)
146               !
147               IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : one hemisphere:  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
148               !
149               CALL lim_rhg( i_j1, i_jpj )
150               !
151            ENDIF
152            !
153         ENDIF
154
155         ! computation of friction velocity
156         ! --------------------------------
157         ! ice-ocean velocity at U & V-points (u_ice v_ice at U- & V-points ; ssu_m, ssv_m at U- & V-points)
158         zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
159         zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
160         ! frictional velocity at T-point
161         zcoef = 0.5_wp * rn_cio
162         DO jj = 2, jpjm1 
163            DO ji = fs_2, fs_jpim1   ! vector opt.
164               ust2s(ji,jj) = zcoef * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
165                  &                    + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) * tmask(ji,jj,1)
166            END DO
167         END DO
168         !
169         ! conservation test
170         IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
171         !
172      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
173         !
174         zcoef = SQRT( 0.5_wp ) * r1_rau0
175         DO jj = 2, jpjm1
176            DO ji = fs_2, fs_jpim1   ! vector opt.
177               ust2s(ji,jj) = zcoef * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
178                  &                        + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) * tmask(ji,jj,1)
179            END DO
180         END DO
181         !
182      ENDIF
183      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
184
185      IF(ln_ctl) THEN   ! Control print
186         CALL prt_ctl_info(' ')
187         CALL prt_ctl_info(' - Cell values : ')
188         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
189         CALL prt_ctl(tab2d_1=ust2s     , clinfo1=' lim_dyn  : ust2s     :')
190         CALL prt_ctl(tab2d_1=divu_i    , clinfo1=' lim_dyn  : divu_i    :')
191         CALL prt_ctl(tab2d_1=delta_i   , clinfo1=' lim_dyn  : delta_i   :')
192         CALL prt_ctl(tab2d_1=strength  , clinfo1=' lim_dyn  : strength  :')
193         CALL prt_ctl(tab2d_1=e12t      , clinfo1=' lim_dyn  : cell area :')
194         CALL prt_ctl(tab2d_1=at_i      , clinfo1=' lim_dyn  : at_i      :')
195         CALL prt_ctl(tab2d_1=vt_i      , clinfo1=' lim_dyn  : vt_i      :')
196         CALL prt_ctl(tab2d_1=vt_s      , clinfo1=' lim_dyn  : vt_s      :')
197         CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' lim_dyn  : stress1_i :')
198         CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' lim_dyn  : stress2_i :')
199         CALL prt_ctl(tab2d_1=stress12_i, clinfo1=' lim_dyn  : stress12_i:')
200         DO jl = 1, jpl
201            CALL prt_ctl_info(' ')
202            CALL prt_ctl_info(' - Category : ', ivar1=jl)
203            CALL prt_ctl_info('   ~~~~~~~~~~')
204            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_dyn  : a_i      : ')
205            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_i     : ')
206            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_dyn  : ht_s     : ')
207            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_dyn  : v_i      : ')
208            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_dyn  : v_s      : ')
209            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : e_s      : ')
210            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_dyn  : t_su     : ')
211            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_dyn  : t_snow   : ')
212            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_dyn  : sm_i     : ')
213            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_dyn  : smv_i    : ')
214            DO ja = 1, nlay_i
215               CALL prt_ctl_info(' ')
216               CALL prt_ctl_info(' - Layer : ', ivar1=ja)
217               CALL prt_ctl_info('   ~~~~~~~')
218               CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : t_i      : ')
219               CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : e_i      : ')
220            END DO
221         END DO
222      ENDIF
223      !
224      CALL wrk_dealloc( jpi, jpj, zu_io, zv_io )
225      CALL wrk_dealloc( jpj, zswitch, zmsk )
226      !
227      IF( nn_timing == 1 )  CALL timing_stop('limdyn')
228
229   END SUBROUTINE lim_dyn
230
231
232   SUBROUTINE lim_dyn_init
233      !!-------------------------------------------------------------------
234      !!                  ***  ROUTINE lim_dyn_init  ***
235      !!
236      !! ** Purpose : Physical constants and parameters linked to the ice
237      !!      dynamics
238      !!
239      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
240      !!       parameter values called at the first timestep (nit000)
241      !!
242      !! ** input   :   Namelist namicedyn
243      !!-------------------------------------------------------------------
244      INTEGER  ::   ios                 ! Local integer output status for namelist read
245      NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, &
246         &                nn_nevp, rn_relast, nn_ahi0, rn_ahi0_ref
247      INTEGER  ::   ji, jj
248      REAL(wp) ::   za00, zd_max
249      !!-------------------------------------------------------------------
250
251      REWIND( numnam_ice_ref )              ! Namelist namicedyn in reference namelist : Ice dynamics
252      READ  ( numnam_ice_ref, namicedyn, IOSTAT = ios, ERR = 901)
253901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in reference namelist', lwp )
254
255      REWIND( numnam_ice_cfg )              ! Namelist namicedyn in configuration namelist : Ice dynamics
256      READ  ( numnam_ice_cfg, namicedyn, IOSTAT = ios, ERR = 902 )
257902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in configuration namelist', lwp )
258      IF(lwm) WRITE ( numoni, namicedyn )
259     
260      IF(lwp) THEN                        ! control print
261         WRITE(numout,*)
262         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
263         WRITE(numout,*) '~~~~~~~~~~~~'
264         WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)  nn_icestr     = ', nn_icestr 
265         WRITE(numout,*)'    Including brine volume in ice strength comp.         ln_icestr_bvf = ', ln_icestr_bvf
266         WRITE(numout,*)'   Ratio of ridging work to PotEner change in ridging    rn_pe_rdg     = ', rn_pe_rdg 
267         WRITE(numout,*) '   drag coefficient for oceanic stress                  rn_cio        = ', rn_cio
268         WRITE(numout,*) '   first bulk-rheology parameter                        rn_pstar      = ', rn_pstar
269         WRITE(numout,*) '   second bulk-rhelogy parameter                        rn_crhg       = ', rn_crhg
270         WRITE(numout,*) '   creep limit                                          rn_creepl     = ', rn_creepl
271         WRITE(numout,*) '   eccentricity of the elliptical yield curve           rn_ecc        = ', rn_ecc
272         WRITE(numout,*) '   number of iterations for subcycling                  nn_nevp       = ', nn_nevp
273         WRITE(numout,*) '   ratio of elastic timescale over ice time step        rn_relast     = ', rn_relast
274         WRITE(numout,*) '   horizontal diffusivity calculation                   nn_ahi0       = ', nn_ahi0
275         WRITE(numout,*) '   horizontal diffusivity coeff. (orca2 grid)           rn_ahi0_ref   = ', rn_ahi0_ref
276      ENDIF
277      !
278      usecc2 = 1._wp / ( rn_ecc * rn_ecc )
279      rhoco  = rau0  * rn_cio
280      !
281      !  Diffusion coefficients
282      SELECT CASE( nn_ahi0 )
283
284      CASE( 0 )
285         ahiu(:,:) = rn_ahi0_ref
286         ahiv(:,:) = rn_ahi0_ref
287
288         IF(lwp) WRITE(numout,*) ''
289         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim constant = rn_ahi0_ref'
290
291      CASE( 1 ) 
292
293         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
294         IF( lk_mpp )   CALL mpp_max( zd_max )          ! max over the global domain
295         
296         ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp   ! 1.e05 = 100km = max grid space at 60° latitude in orca2
297                                                        !                    (60° = min latitude for ice cover) 
298         ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp
299
300         IF(lwp) WRITE(numout,*) ''
301         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')'
302         IF(lwp) WRITE(numout,*) '   value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp 
303         
304      CASE( 2 ) 
305
306         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
307         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
308         
309         za00 = rn_ahi0_ref * 1.e-05_wp          ! 1.e05 = 100km = max grid space at 60° latitude in orca2
310                                                 !                    (60° = min latitude for ice cover) 
311         DO jj = 1, jpj
312            DO ji = 1, jpi
313               ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1)
314               ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1)
315            END DO
316         END DO
317         !
318         IF(lwp) WRITE(numout,*) ''
319         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to e1'
320         IF(lwp) WRITE(numout,*) '   maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max
321         
322      END SELECT
323
324   END SUBROUTINE lim_dyn_init
325
326#else
327   !!----------------------------------------------------------------------
328   !!   Default option          Empty module           NO LIM sea-ice model
329   !!----------------------------------------------------------------------
330CONTAINS
331   SUBROUTINE lim_dyn         ! Empty routine
332   END SUBROUTINE lim_dyn
333#endif 
334
335   !!======================================================================
336END MODULE limdyn
Note: See TracBrowser for help on using the repository browser.