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

source: trunk/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90 @ 2855

Last change on this file since 2855 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk

  • Property svn:keywords set to Id
File size: 14.3 KB
Line 
1MODULE limdyn_2
2   !!======================================================================
3   !!                     ***  MODULE  limdyn_2  ***
4   !!   Sea-Ice dynamics : 
5   !!======================================================================
6   !! History :  1.0  ! 2001-04  (LIM)  Original code
7   !!            2.0  ! 2002-08  (C. Ethe, G. Madec)  F90, mpp
8   !!            2.0  ! 2003-08  (C. Ethe) add lim_dyn_init
9   !!            2.0  ! 2006-07  (G. Madec)  Surface module
10   !!            3.3  ! 2009-05 (G. Garric, C. Bricaud) addition of the lim2_evp case
11   !!---------------------------------------------------------------------
12#if defined key_lim2
13   !!----------------------------------------------------------------------
14   !!   'key_lim2' :                                  LIM 2.0 sea-ice model
15   !!----------------------------------------------------------------------
16   !!    lim_dyn_2      : computes ice velocities
17   !!    lim_dyn_init_2 : initialization and namelist read
18   !!----------------------------------------------------------------------
19   USE dom_oce          ! ocean space and time domain
20   USE sbc_oce          ! ocean surface boundary condition
21   USE phycst           ! physical constant
22   USE ice_2            ! LIM-2: ice variables
23   USE sbc_ice          ! Surface boundary condition: sea-ice fields
24   USE dom_ice_2        ! LIM-2: ice domain
25   USE limistate_2      ! LIM-2: initial state
26   USE limrhg_2         ! LIM-2: VP  ice rheology
27   USE limrhg           ! LIM  : EVP ice rheology
28   USE lbclnk           ! lateral boundary condition - MPP link
29   USE lib_mpp          ! MPP library
30   USE in_out_manager   ! I/O manager
31   USE prtctl           ! Print control
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   lim_dyn_2   ! routine called by sbc_ice_lim
37
38   !! * Substitutions
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
42   !! $Id$
43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE lim_dyn_2( kt )
48      !!-------------------------------------------------------------------
49      !!               ***  ROUTINE lim_dyn_2  ***
50      !!               
51      !! ** Purpose :   compute ice velocity and ocean-ice friction velocity
52      !!               
53      !! ** Method  :
54      !!
55      !! ** Action  : - Initialisation
56      !!              - Call of the dynamic routine for each hemisphere
57      !!              - computation of the friction velocity at the sea-ice base
58      !!              - treatment of the case if no ice dynamic
59      !!---------------------------------------------------------------------
60      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
61      USE wrk_nemo, ONLY:   wrk_1d_1, wrk_1d_2
62      USE wrk_nemo, ONLY:   zu_io => wrk_2d_1, zv_io => wrk_2d_2  ! ice-ocean velocity
63      !
64      INTEGER, INTENT(in) ::   kt     ! number of iteration
65      !!
66      INTEGER  ::   ji, jj             ! dummy loop indices
67      INTEGER  ::   i_j1, i_jpj        ! Starting/ending j-indices for rheology
68      REAL(wp) ::   zcoef              ! temporary scalar
69      REAL(wp), POINTER, DIMENSION(:) ::   zind     ! i-averaged indicator of sea-ice
70      REAL(wp), POINTER, DIMENSION(:) ::   zmsk     ! i-averaged of tmask
71      !!---------------------------------------------------------------------
72
73      IF(  wrk_in_use(1, 1,2)  .OR.  wrk_in_use(2, 1,2)  ) THEN
74         CALL ctl_stop( 'lim_dyn_2 : requested workspace arrays unavailable' )   ;   RETURN
75      ENDIF
76      zind => wrk_1d_1(1:jpj)      ! Set-up pointers to sub-arrays of workspaces
77      zmsk => wrk_1d_2(1:jpj)
78
79      IF( kt == nit000 )   CALL lim_dyn_init_2   ! Initialization (first time-step only)
80     
81      IF( ln_limdyn ) THEN
82         !
83         ! Mean ice and snow thicknesses.         
84         hsnm(:,:)  = ( 1.0 - frld(:,:) ) * hsnif(:,:)
85         hicm(:,:)  = ( 1.0 - frld(:,:) ) * hicif(:,:)
86         !
87         !                                     ! Rheology (ice dynamics)
88         !                                     ! ========
89         
90         !  Define the j-limits where ice rheology is computed
91         ! ---------------------------------------------------
92         
93         IF( lk_mpp .OR. lk_mpp_rep ) THEN                    ! mpp: compute over the whole domain
94            i_j1 = 1   
95            i_jpj = jpj
96            IF(ln_ctl)   CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
97            IF( lk_lim2_vp )   THEN   ;   CALL lim_rhg_2( i_j1, i_jpj )             !  VP rheology
98            ELSE                      ;   CALL lim_rhg  ( i_j1, i_jpj )             ! EVP rheology
99            ENDIF
100            !
101         ELSE                                 ! optimization of the computational area
102            !
103            DO jj = 1, jpj
104               zind(jj) = SUM( frld (:,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. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
115                  i_j1 = i_j1 + 1
116               END DO
117               IF( lk_lim2_vp )   THEN             ! VP  rheology
118                  i_j1 = MAX( 1, i_j1-1 )
119                  CALL lim_rhg_2( i_j1, i_jpj )
120               ELSE                                ! EVP rheology
121                  i_j1 = MAX( 1, i_j1-2 )
122                  CALL lim_rhg( i_j1, i_jpj )
123               ENDIF
124               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, 'ij_jpj = ', 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               IF( lk_lim2_vp )   THEN             ! VP  rheology
133                  i_jpj = MIN( jpj, i_jpj+2 )
134                  CALL lim_rhg_2( i_j1, i_jpj )
135               ELSE                                ! EVP rheology
136                  i_jpj = MIN( jpj, i_jpj+1 )
137                  CALL lim_rhg( i_j1, i_jpj )
138               ENDIF
139               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, 'ij_jpj = ', i_jpj
140               !
141            ELSE                                 ! local domain extends over one hemisphere only
142               !                                 ! Rheology is computed only over the ice cover
143               !                                 ! latitude strip
144               i_j1  = 1
145               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
146                  i_j1 = i_j1 + 1
147               END DO
148               i_j1 = MAX( 1, i_j1-1 )
149   
150               i_jpj  = jpj
151               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
152                  i_jpj = i_jpj - 1
153               END DO
154               i_jpj = MIN( jpj, i_jpj+2 )
155               !
156               IF( lk_lim2_vp )   THEN             ! VP  rheology
157                  i_jpj = MIN( jpj, i_jpj+2 )
158                  CALL lim_rhg_2( i_j1, i_jpj )                !  VP rheology
159               ELSE                                ! EVP rheology
160                  i_j1  = MAX( 1  , i_j1-2  )
161                  i_jpj = MIN( jpj, i_jpj+1 )
162                  CALL lim_rhg  ( i_j1, i_jpj )                ! EVP rheology
163               ENDIF
164               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
165               !
166            ENDIF
167            !
168         ENDIF
169
170         IF(ln_ctl)   CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :')
171         
172         ! computation of friction velocity
173         ! --------------------------------
174         SELECT CASE( cp_ice_msh )           ! ice-ocean relative velocity at u- & v-pts
175         CASE( 'C' )                               ! EVP : C-grid ice dynamics
176            zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)           ! ice-ocean & ice velocity at ocean velocity points
177            zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
178         CASE( 'I' )                               ! VP  : B-grid ice dynamics (I-point)
179            DO jj = 1, jpjm1                               ! u_ice v_ice at I-point ; ssu_m, ssv_m at U- & V-points
180               DO ji = 1, jpim1   ! NO vector opt.         !
181                  zu_io(ji,jj) = 0.5_wp * ( u_ice(ji+1,jj+1) + u_ice(ji+1,jj  ) ) - ssu_m(ji,jj)
182                  zv_io(ji,jj) = 0.5_wp * ( v_ice(ji+1,jj+1) + v_ice(ji  ,jj+1) ) - ssv_m(ji,jj)
183               END DO
184            END DO
185         END SELECT
186
187         ! frictional velocity at T-point
188         zcoef = 0.5_wp * cw
189         DO jj = 2, jpjm1
190            DO ji = 2, jpim1   ! NO vector opt. because of zu_io
191               ust2s(ji,jj) = zcoef * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
192                  &                    + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj)
193            END DO
194         END DO
195         !
196      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
197         !
198         zcoef = SQRT( 0.5 ) / rau0
199         DO jj = 2, jpjm1
200            DO ji = fs_2, fs_jpim1   ! vector opt.
201               ust2s(ji,jj) = zcoef * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
202                  &                        + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1)   ) * tms(ji,jj)
203            END DO
204         END DO
205         !
206      ENDIF
207      !
208      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
209      !
210      IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :')
211      !
212      IF( wrk_not_released(1, 1,2)     .OR.   &
213          wrk_not_released(2, 1,2) )   CALL ctl_stop('lim_dyn_2 : failed to release workspace arrays')
214      !
215   END SUBROUTINE lim_dyn_2
216
217
218   SUBROUTINE lim_dyn_init_2
219      !!-------------------------------------------------------------------
220      !!                  ***  ROUTINE lim_dyn_init_2  ***
221      !!
222      !! ** Purpose :   Physical constants and parameters linked to the ice
223      !!              dynamics
224      !!
225      !! ** Method  :   Read the namicedyn namelist and check the ice-dynamic
226      !!              parameter values
227      !!
228      !! ** input   :   Namelist namicedyn
229      !!-------------------------------------------------------------------
230      NAMELIST/namicedyn/ epsd, alpha,     &
231         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
232         &                c_rhg, etamn, creepl, ecc, ahi0,                  &
233         &                nevp, telast,alphaevp
234      !!-------------------------------------------------------------------
235
236      REWIND ( numnam_ice )                       ! Read Namelist namicedyn
237      READ   ( numnam_ice  , namicedyn )
238
239      IF(lwp) THEN                                ! Control print
240         WRITE(numout,*)
241         WRITE(numout,*) 'lim_dyn_init_2: ice parameters for ice dynamics '
242         WRITE(numout,*) '~~~~~~~~~~~~~~'
243         WRITE(numout,*) '       tolerance parameter                              epsd   = ', epsd
244         WRITE(numout,*) '       coefficient for semi-implicit coriolis           alpha  = ', alpha
245         WRITE(numout,*) '       diffusion constant for dynamics                  dm     = ', dm
246         WRITE(numout,*) '       number of sub-time steps for relaxation          nbiter = ', nbiter
247         WRITE(numout,*) '       maximum number of iterations for relaxation      nbitdr = ', nbitdr
248         WRITE(numout,*) '       relaxation constant                              om     = ', om
249         WRITE(numout,*) '       maximum value for the residual of relaxation     resl   = ', resl
250         WRITE(numout,*) '       drag coefficient for oceanic stress              cw     = ', cw
251         WRITE(numout,*) '       turning angle for oceanic stress                 angvg  = ', angvg, ' degrees'
252         WRITE(numout,*) '       first bulk-rheology parameter                    pstar  = ', pstar
253         WRITE(numout,*) '       second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
254         WRITE(numout,*) '       minimun value for viscosity                      etamn  = ', etamn
255         WRITE(numout,*) '       creep limit                                      creepl = ', creepl
256         WRITE(numout,*) '       eccentricity of the elliptical yield curve       ecc    = ', ecc
257         WRITE(numout,*) '       horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
258         WRITE(numout,*) '       number of iterations for subcycling nevp   = ', nevp
259         WRITE(numout,*) '       timescale for elastic waves telast = ', telast
260         WRITE(numout,*) '       coefficient for the solution of int. stresses alphaevp = ', alphaevp
261      ENDIF
262      !
263      IF( angvg /= 0._wp .AND. .NOT.lk_lim2_vp ) THEN
264         CALL ctl_warn( 'lim_dyn_init_2: turning angle for oceanic stress not properly coded for EVP ',   &
265            &           '(see limsbc_2 module). We force  angvg = 0._wp'  )
266         angvg = 0._wp
267      ENDIF
268
269      !  Initialization
270      usecc2 = 1.0 / ( ecc * ecc )
271      rhoco  = rau0 * cw
272      angvg  = angvg * rad      ! convert angvg from degree to radian
273      sangvg = SIN( angvg )
274      cangvg = COS( angvg )
275      pstarh = pstar / 2.0
276      !
277      ahiu(:,:) = ahi0 * umask(:,:,1)            ! Ice eddy Diffusivity coefficients.
278      ahiv(:,:) = ahi0 * vmask(:,:,1)
279      !
280   END SUBROUTINE lim_dyn_init_2
281
282#else
283   !!----------------------------------------------------------------------
284   !!   Default option          Empty module       NO LIM 2.0 sea-ice model
285   !!----------------------------------------------------------------------
286CONTAINS
287   SUBROUTINE lim_dyn_2         ! Empty routine
288   END SUBROUTINE lim_dyn_2
289#endif 
290
291   !!======================================================================
292END MODULE limdyn_2
Note: See TracBrowser for help on using the repository browser.