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

source: branches/dev_1784_EVP/NEMO/LIM_SRC_2/limdyn_2.F90 @ 2119

Last change on this file since 2119 was 2095, checked in by cbricaud, 14 years ago

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