source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.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: 14.7 KB
RevLine 
[821]1MODULE limdyn_2
[3]2   !!======================================================================
[821]3   !!                     ***  MODULE  limdyn_2  ***
[3]4   !!   Sea-Ice dynamics : 
5   !!======================================================================
[2528]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
[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   !!----------------------------------------------------------------------
[2528]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
[3294]30   USE wrk_nemo         ! work arrays
[2528]31   USE in_out_manager   ! I/O manager
32   USE prtctl           ! Print control
[3625]33   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
[3]34
35   IMPLICIT NONE
36   PRIVATE
37
[2528]38   PUBLIC   lim_dyn_2   ! routine called by sbc_ice_lim
[3]39
[2528]40   !! * Substitutions
[888]41#  include "vectopt_loop_substitute.h90"
[3]42   !!----------------------------------------------------------------------
[2528]43   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
[1156]44   !! $Id$
[2528]45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]46   !!----------------------------------------------------------------------
47CONTAINS
48
[821]49   SUBROUTINE lim_dyn_2( kt )
[3]50      !!-------------------------------------------------------------------
[821]51      !!               ***  ROUTINE lim_dyn_2  ***
[3]52      !!               
[888]53      !! ** Purpose :   compute ice velocity and ocean-ice friction velocity
[3]54      !!               
55      !! ** Method  :
56      !!
57      !! ** Action  : - Initialisation
58      !!              - Call of the dynamic routine for each hemisphere
[888]59      !!              - computation of the friction velocity at the sea-ice base
[3]60      !!              - treatment of the case if no ice dynamic
61      !!---------------------------------------------------------------------
[508]62      INTEGER, INTENT(in) ::   kt     ! number of iteration
[888]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
[3294]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
[3]70      !!---------------------------------------------------------------------
71
[3294]72      CALL wrk_alloc( jpi, jpj, zu_io, zv_io )
73      CALL wrk_alloc(      jpj, zind , zmsk  )
[2715]74
[888]75      IF( kt == nit000 )   CALL lim_dyn_init_2   ! Initialization (first time-step only)
[3]76     
[888]77      IF( ln_limdyn ) THEN
78         !
[3]79         ! Mean ice and snow thicknesses.         
80         hsnm(:,:)  = ( 1.0 - frld(:,:) ) * hsnif(:,:)
81         hicm(:,:)  = ( 1.0 - frld(:,:) ) * hicif(:,:)
[888]82         !
83         !                                     ! Rheology (ice dynamics)
84         !                                     ! ========
[76]85         
86         !  Define the j-limits where ice rheology is computed
87         ! ---------------------------------------------------
88         
[2528]89         IF( lk_mpp .OR. lk_mpp_rep ) THEN                    ! mpp: compute over the whole domain
[76]90            i_j1 = 1   
91            i_jpj = jpj
[888]92            IF(ln_ctl)   CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
[2528]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
[888]96            !
[76]97         ELSE                                 ! optimization of the computational area
[888]98            !
[76]99            DO jj = 1, jpj
[2715]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
[76]102            END DO
[888]103            !
[76]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
[2528]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               !
[76]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
[2528]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               !
[76]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
[2528]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
[258]160               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
[888]161               !
[76]162            ENDIF
[888]163            !
[76]164         ENDIF
165
[1470]166         IF(ln_ctl)   CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :')
[3]167         
[888]168         ! computation of friction velocity
169         ! --------------------------------
[2528]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
[3]180            END DO
[2528]181         END SELECT
182
[888]183         ! frictional velocity at T-point
[2528]184         zcoef = 0.5_wp * cw
[3]185         DO jj = 2, jpjm1
[1694]186            DO ji = 2, jpim1   ! NO vector opt. because of zu_io
[2528]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)
[3]189            END DO
190         END DO
[888]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.
[2528]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)
[3]199            END DO
200         END DO
[888]201         !
[3]202      ENDIF
[888]203      !
[3]204      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
[888]205      !
206      IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :')
[2528]207      !
[3294]208      CALL wrk_dealloc( jpi, jpj, zu_io, zv_io )
209      CALL wrk_dealloc(      jpj, zind , zmsk  )
[2715]210      !
[821]211   END SUBROUTINE lim_dyn_2
[3]212
[76]213
[821]214   SUBROUTINE lim_dyn_init_2
[3]215      !!-------------------------------------------------------------------
[821]216      !!                  ***  ROUTINE lim_dyn_init_2  ***
[3]217      !!
[888]218      !! ** Purpose :   Physical constants and parameters linked to the ice
219      !!              dynamics
[3]220      !!
[888]221      !! ** Method  :   Read the namicedyn namelist and check the ice-dynamic
222      !!              parameter values
[3]223      !!
224      !! ** input   :   Namelist namicedyn
225      !!-------------------------------------------------------------------
[4147]226      INTEGER  ::   ios                 ! Local integer output status for namelist read
[12]227      NAMELIST/namicedyn/ epsd, alpha,     &
[3]228         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
[5600]229         &                c_rhg, etamn, rn_creepl, rn_ecc, ahi0,                  &
230         &                nn_nevp, telast, alphaevp
[3]231      !!-------------------------------------------------------------------
[4147]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 )
[3]236
[4147]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 )
[4624]240      IF(lwm) WRITE ( numoni, namicedyn )
[719]241
[888]242      IF(lwp) THEN                                ! Control print
[3]243         WRITE(numout,*)
[821]244         WRITE(numout,*) 'lim_dyn_init_2: ice parameters for ice dynamics '
245         WRITE(numout,*) '~~~~~~~~~~~~~~'
[76]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
[888]254         WRITE(numout,*) '       turning angle for oceanic stress                 angvg  = ', angvg, ' degrees'
[76]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
[5600]258         WRITE(numout,*) '       creep limit                                      rn_creepl = ', rn_creepl
259         WRITE(numout,*) '       eccentricity of the elliptical yield curve       rn_ecc = ', rn_ecc
[76]260         WRITE(numout,*) '       horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
[5600]261         WRITE(numout,*) '       number of iterations for subcycling              nn_nevp= ', nn_nevp
[2528]262         WRITE(numout,*) '       timescale for elastic waves telast = ', telast
263         WRITE(numout,*) '       coefficient for the solution of int. stresses alphaevp = ', alphaevp
[3]264      ENDIF
[2528]265      !
266      IF( angvg /= 0._wp .AND. .NOT.lk_lim2_vp ) THEN
267         CALL ctl_warn( 'lim_dyn_init_2: turning angle for oceanic stress not properly coded for EVP ',   &
268            &           '(see limsbc_2 module). We force  angvg = 0._wp'  )
269         angvg = 0._wp
270      ENDIF
[3]271
[719]272      !  Initialization
[5600]273      usecc2 = 1.0 / ( rn_ecc * rn_ecc )
[3]274      rhoco  = rau0 * cw
[888]275      angvg  = angvg * rad      ! convert angvg from degree to radian
[3]276      sangvg = SIN( angvg )
277      cangvg = COS( angvg )
278      pstarh = pstar / 2.0
[888]279      !
280      ahiu(:,:) = ahi0 * umask(:,:,1)            ! Ice eddy Diffusivity coefficients.
[3]281      ahiv(:,:) = ahi0 * vmask(:,:,1)
[888]282      !
[821]283   END SUBROUTINE lim_dyn_init_2
[3]284
285#else
286   !!----------------------------------------------------------------------
[821]287   !!   Default option          Empty module       NO LIM 2.0 sea-ice model
[3]288   !!----------------------------------------------------------------------
289CONTAINS
[821]290   SUBROUTINE lim_dyn_2         ! Empty routine
291   END SUBROUTINE lim_dyn_2
[3]292#endif 
293
294   !!======================================================================
[821]295END MODULE limdyn_2
Note: See TracBrowser for help on using the repository browser.