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 branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90 @ 4810

Last change on this file since 4810 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: 14.7 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 wrk_nemo         ! work arrays
31   USE in_out_manager   ! I/O manager
32   USE prtctl           ! Print control
33   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   lim_dyn_2   ! routine called by sbc_ice_lim
39
40   !! * Substitutions
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
44   !! $Id$
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE lim_dyn_2( kt )
50      !!-------------------------------------------------------------------
51      !!               ***  ROUTINE lim_dyn_2  ***
52      !!               
53      !! ** Purpose :   compute ice velocity and ocean-ice friction velocity
54      !!               
55      !! ** Method  :
56      !!
57      !! ** Action  : - Initialisation
58      !!              - Call of the dynamic routine for each hemisphere
59      !!              - computation of the friction velocity at the sea-ice base
60      !!              - treatment of the case if no ice dynamic
61      !!---------------------------------------------------------------------
62      INTEGER, INTENT(in) ::   kt     ! number of iteration
63      !!
64      INTEGER  ::   ji, jj             ! dummy loop indices
65      INTEGER  ::   i_j1, i_jpj        ! Starting/ending j-indices for rheology
66      REAL(wp) ::   zcoef              ! temporary scalar
67      REAL(wp), POINTER, DIMENSION(:  ) ::   zind           ! i-averaged indicator of sea-ice
68      REAL(wp), POINTER, DIMENSION(:  ) ::   zmsk           ! i-averaged of tmask
69      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_io, zv_io   ! ice-ocean velocity
70      !!---------------------------------------------------------------------
71
72      CALL wrk_alloc( jpi, jpj, zu_io, zv_io )
73      CALL wrk_alloc(      jpj, zind , zmsk  )
74
75      IF( kt == nit000 )   CALL lim_dyn_init_2   ! Initialization (first time-step only)
76     
77      IF( ln_limdyn ) THEN
78         !
79         ! Mean ice and snow thicknesses.         
80         hsnm(:,:)  = ( 1.0 - frld(:,:) ) * hsnif(:,:)
81         hicm(:,:)  = ( 1.0 - frld(:,:) ) * hicif(:,:)
82         !
83         !                                     ! Rheology (ice dynamics)
84         !                                     ! ========
85         
86         !  Define the j-limits where ice rheology is computed
87         ! ---------------------------------------------------
88         
89         IF( lk_mpp .OR. lk_mpp_rep ) THEN                    ! mpp: compute over the whole domain
90            i_j1 = 1   
91            i_jpj = jpj
92            IF(ln_ctl)   CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
93            IF( lk_lim2_vp )   THEN   ;   CALL lim_rhg_2( i_j1, i_jpj )             !  VP rheology
94            ELSE                      ;   CALL lim_rhg  ( i_j1, i_jpj )             ! EVP rheology
95            ENDIF
96            !
97         ELSE                                 ! optimization of the computational area
98            !
99            DO jj = 1, jpj
100               zind(jj) = SUM( frld (:,jj  ) )   ! = REAL(jpj) if ocean everywhere on a j-line
101               zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0         if land  everywhere on a j-line
102            END DO
103            !
104            IF( l_jeq ) THEN                     ! local domain include both hemisphere
105               !                                 ! Rheology is computed in each hemisphere
106               !                                 ! only over the ice cover latitude strip
107               ! Northern hemisphere
108               i_j1  = njeq
109               i_jpj = jpj
110               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
111                  i_j1 = i_j1 + 1
112               END DO
113               IF( lk_lim2_vp )   THEN             ! VP  rheology
114                  i_j1 = MAX( 1, i_j1-1 )
115                  CALL lim_rhg_2( i_j1, i_jpj )
116               ELSE                                ! EVP rheology
117                  i_j1 = MAX( 1, i_j1-2 )
118                  CALL lim_rhg( i_j1, i_jpj )
119               ENDIF
120               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, 'ij_jpj = ', i_jpj
121               !
122               ! Southern hemisphere
123               i_j1  =  1 
124               i_jpj = njeq
125               DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
126                  i_jpj = i_jpj - 1
127               END DO
128               IF( lk_lim2_vp )   THEN             ! VP  rheology
129                  i_jpj = MIN( jpj, i_jpj+2 )
130                  CALL lim_rhg_2( i_j1, i_jpj )
131               ELSE                                ! EVP rheology
132                  i_jpj = MIN( jpj, i_jpj+1 )
133                  CALL lim_rhg( i_j1, i_jpj )
134               ENDIF
135               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, 'ij_jpj = ', 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-1 )
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+2 )
151               !
152               IF( lk_lim2_vp )   THEN             ! VP  rheology
153                  i_jpj = MIN( jpj, i_jpj+2 )
154                  CALL lim_rhg_2( i_j1, i_jpj )                !  VP rheology
155               ELSE                                ! EVP rheology
156                  i_j1  = MAX( 1  , i_j1-2  )
157                  i_jpj = MIN( jpj, i_jpj+1 )
158                  CALL lim_rhg  ( i_j1, i_jpj )                ! EVP rheology
159               ENDIF
160               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
161               !
162            ENDIF
163            !
164         ENDIF
165
166         IF(ln_ctl)   CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :')
167         
168         ! computation of friction velocity
169         ! --------------------------------
170         SELECT CASE( cp_ice_msh )           ! ice-ocean relative velocity at u- & v-pts
171         CASE( 'C' )                               ! EVP : C-grid ice dynamics
172            zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)           ! ice-ocean & ice velocity at ocean velocity points
173            zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
174         CASE( 'I' )                               ! VP  : B-grid ice dynamics (I-point)
175            DO jj = 1, jpjm1                               ! u_ice v_ice at I-point ; ssu_m, ssv_m at U- & V-points
176               DO ji = 1, jpim1   ! NO vector opt.         !
177                  zu_io(ji,jj) = 0.5_wp * ( u_ice(ji+1,jj+1) + u_ice(ji+1,jj  ) ) - ssu_m(ji,jj)
178                  zv_io(ji,jj) = 0.5_wp * ( v_ice(ji+1,jj+1) + v_ice(ji  ,jj+1) ) - ssv_m(ji,jj)
179               END DO
180            END DO
181         END SELECT
182
183         ! frictional velocity at T-point
184         zcoef = 0.5_wp * cw
185         DO jj = 2, jpjm1
186            DO ji = 2, jpim1   ! NO vector opt. because of zu_io
187               ust2s(ji,jj) = zcoef * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
188                  &                    + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj)
189            END DO
190         END DO
191         !
192      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
193         !
194         zcoef = SQRT( 0.5 ) / rau0
195         DO jj = 2, jpjm1
196            DO ji = fs_2, fs_jpim1   ! vector opt.
197               ust2s(ji,jj) = zcoef * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
198                  &                        + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1)   ) * tms(ji,jj)
199            END DO
200         END DO
201         !
202      ENDIF
203      !
204      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
205      !
206      IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :')
207      !
208      CALL wrk_dealloc( jpi, jpj, zu_io, zv_io )
209      CALL wrk_dealloc(      jpj, zind , zmsk  )
210      !
211   END SUBROUTINE lim_dyn_2
212
213
214   SUBROUTINE lim_dyn_init_2
215      !!-------------------------------------------------------------------
216      !!                  ***  ROUTINE lim_dyn_init_2  ***
217      !!
218      !! ** Purpose :   Physical constants and parameters linked to the ice
219      !!              dynamics
220      !!
221      !! ** Method  :   Read the namicedyn namelist and check the ice-dynamic
222      !!              parameter values
223      !!
224      !! ** input   :   Namelist namicedyn
225      !!-------------------------------------------------------------------
226      INTEGER  ::   ios                 ! Local integer output status for namelist read
227      NAMELIST/namicedyn/ epsd, alpha,     &
228         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
229         &                c_rhg, etamn, creepl, ecc, ahi0,                  &
230         &                nevp, telast, alphaevp, hminrhg
231      !!-------------------------------------------------------------------
232                   
233      REWIND( numnam_ice_ref )              ! Namelist namicedyn in reference namelist : Ice dynamics
234      READ  ( numnam_ice_ref, namicedyn, IOSTAT = ios, ERR = 901)
235901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in reference namelist', lwp )
236
237      REWIND( numnam_ice_cfg )              ! Namelist namicedyn in configuration namelist : Ice dynamics
238      READ  ( numnam_ice_cfg, namicedyn, IOSTAT = ios, ERR = 902 )
239902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in configuration namelist', lwp )
240      IF(lwm) WRITE ( numoni, namicedyn )
241
242      IF(lwp) THEN                                ! Control print
243         WRITE(numout,*)
244         WRITE(numout,*) 'lim_dyn_init_2: ice parameters for ice dynamics '
245         WRITE(numout,*) '~~~~~~~~~~~~~~'
246         WRITE(numout,*) '       tolerance parameter                              epsd   = ', epsd
247         WRITE(numout,*) '       coefficient for semi-implicit coriolis           alpha  = ', alpha
248         WRITE(numout,*) '       diffusion constant for dynamics                  dm     = ', dm
249         WRITE(numout,*) '       number of sub-time steps for relaxation          nbiter = ', nbiter
250         WRITE(numout,*) '       maximum number of iterations for relaxation      nbitdr = ', nbitdr
251         WRITE(numout,*) '       relaxation constant                              om     = ', om
252         WRITE(numout,*) '       maximum value for the residual of relaxation     resl   = ', resl
253         WRITE(numout,*) '       drag coefficient for oceanic stress              cw     = ', cw
254         WRITE(numout,*) '       turning angle for oceanic stress                 angvg  = ', angvg, ' degrees'
255         WRITE(numout,*) '       first bulk-rheology parameter                    pstar  = ', pstar
256         WRITE(numout,*) '       second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
257         WRITE(numout,*) '       minimun value for viscosity                      etamn  = ', etamn
258         WRITE(numout,*) '       creep limit                                      creepl = ', creepl
259         WRITE(numout,*) '       eccentricity of the elliptical yield curve       ecc    = ', ecc
260         WRITE(numout,*) '       horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
261         WRITE(numout,*) '       number of iterations for subcycling nevp   = ', nevp
262         WRITE(numout,*) '       timescale for elastic waves telast = ', telast
263         WRITE(numout,*) '       coefficient for the solution of int. stresses alphaevp = ', alphaevp
264         WRITE(numout,*) '       min ice thickness for rheology calculations     hminrhg = ', hminrhg
265      ENDIF
266      !
267      IF( angvg /= 0._wp .AND. .NOT.lk_lim2_vp ) THEN
268         CALL ctl_warn( 'lim_dyn_init_2: turning angle for oceanic stress not properly coded for EVP ',   &
269            &           '(see limsbc_2 module). We force  angvg = 0._wp'  )
270         angvg = 0._wp
271      ENDIF
272
273      !  Initialization
274      usecc2 = 1.0 / ( ecc * ecc )
275      rhoco  = rau0 * cw
276      angvg  = angvg * rad      ! convert angvg from degree to radian
277      sangvg = SIN( angvg )
278      cangvg = COS( angvg )
279      pstarh = pstar / 2.0
280      !
281      ahiu(:,:) = ahi0 * umask(:,:,1)            ! Ice eddy Diffusivity coefficients.
282      ahiv(:,:) = ahi0 * vmask(:,:,1)
283      !
284   END SUBROUTINE lim_dyn_init_2
285
286#else
287   !!----------------------------------------------------------------------
288   !!   Default option          Empty module       NO LIM 2.0 sea-ice model
289   !!----------------------------------------------------------------------
290CONTAINS
291   SUBROUTINE lim_dyn_2         ! Empty routine
292   END SUBROUTINE lim_dyn_2
293#endif 
294
295   !!======================================================================
296END MODULE limdyn_2
Note: See TracBrowser for help on using the repository browser.