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

source: branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limdyn_2.F90 @ 1857

Last change on this file since 1857 was 1857, checked in by gm, 14 years ago

ticket:#665 Reverting previous commit and going back to revision 1850

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