Changeset 5208 for branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90

Timestamp:

2015-04-13T15:08:59+02:00 (9 years ago)

Author:

davestorkey

Message:

Merge in changes from trunk up to 5021.

File:

• branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90 (modified) (12 diffs)

branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90

@@ -64 +64 @@
       INTEGER  ::   i_j1, i_jpj       ! Starting/ending j-indices for rheology
       REAL(wp) ::   zcoef             ! local scalar
-      REAL(wp), POINTER, DIMENSION(:)   ::   zind           ! i-averaged indicator of sea-ice
+      REAL(wp), POINTER, DIMENSION(:)   ::   zswitch        ! i-averaged indicator of sea-ice
       REAL(wp), POINTER, DIMENSION(:)   ::   zmsk           ! i-averaged of tmask
       REAL(wp), POINTER, DIMENSION(:,:) ::   zu_io, zv_io   ! ice-ocean velocity
@@ -74 +74 @@
 
       CALL wrk_alloc( jpi, jpj, zu_io, zv_io )
-      CALL wrk_alloc( jpj, zind, zmsk )
+      CALL wrk_alloc( jpj, zswitch, zmsk )
 
       IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only)
@@ -83 +83 @@
          IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
-         old_u_ice(:,:) = u_ice(:,:) * tmu(:,:)
-         old_v_ice(:,:) = v_ice(:,:) * tmv(:,:)
+         u_ice_b(:,:) = u_ice(:,:) * tmu(:,:)
+         v_ice_b(:,:) = v_ice(:,:) * tmv(:,:)
 
          ! Rheology (ice dynamics)
@@ -100 +100 @@
             !
             DO jj = 1, jpj
-               zind(jj) = SUM( 1.0 - at_i(:,jj) )   ! = REAL(jpj) if ocean everywhere on a j-line
+               zswitch(jj) = SUM( 1.0 - at_i(:,jj) )   ! = REAL(jpj) if ocean everywhere on a j-line
                zmsk(jj) = SUM( tmask(:,jj,1)    )   ! = 0         if land  everywhere on a j-line
             END DO
@@ -110 +110 @@
                i_j1  = njeq
                i_jpj = jpj
-               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
+               DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
                   i_j1 = i_j1 + 1
                END DO
@@ -120 +120 @@
                i_j1  =  1
                i_jpj = njeq
-               DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
+               DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
                   i_jpj = i_jpj - 1
                END DO
@@ -132 +132 @@
                !                                 ! latitude strip
                i_j1  = 1
-               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
+               DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
                   i_j1 = i_j1 + 1
                END DO
@@ -138 +138 @@
 
                i_jpj  = jpj
-               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
+               DO WHILE ( i_jpj >= 1  .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
                   i_jpj = i_jpj - 1
                END DO
@@ -221 +221 @@
       !
       CALL wrk_dealloc( jpi, jpj, zu_io, zv_io )
-      CALL wrk_dealloc( jpj, zind, zmsk )
+      CALL wrk_dealloc( jpj, zswitch, zmsk )
       !
       IF( nn_timing == 1 )  CALL timing_stop('limdyn')
@@ -241 +241 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicedyn/ epsd, om, cw, angvg, pstar,   &
+      NAMELIST/namicedyn/ epsd, om, cw, pstar,   &
          &                c_rhg, creepl, ecc, ahi0,     &
-         &                nevp, telast, alphaevp, hminrhg
+         &                nevp, relast, alphaevp, hminrhg
       !!-------------------------------------------------------------------
 
@@ -262 +262 @@
          WRITE(numout,*) '   relaxation constant                              om     = ', om
          WRITE(numout,*) '   drag coefficient for oceanic stress              cw     = ', cw
-         WRITE(numout,*) '   turning angle for oceanic stress                 angvg  = ', angvg
          WRITE(numout,*) '   first bulk-rheology parameter                    pstar  = ', pstar
          WRITE(numout,*) '   second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
@@ -269 +268 @@
          WRITE(numout,*) '   horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
          WRITE(numout,*) '   number of iterations for subcycling              nevp   = ', nevp
-         WRITE(numout,*) '   timescale for elastic waves                      telast = ', telast
+         WRITE(numout,*) '   ratio of elastic timescale over ice time step    relast = ', relast
          WRITE(numout,*) '   coefficient for the solution of int. stresses  alphaevp = ', alphaevp
          WRITE(numout,*) '   min ice thickness for rheology calculations     hminrhg = ', hminrhg
       ENDIF
       !
-      IF( angvg /= 0._wp ) THEN
-         CALL ctl_warn( 'lim_dyn_init: turning angle for oceanic stress not properly coded for EVP ',   &
-            &           '(see limsbc module). We force  angvg = 0._wp'  )
-         angvg = 0._wp
-      ENDIF
-      
       usecc2 = 1._wp / ( ecc * ecc )
       rhoco  = rau0  * cw
-      angvg  = angvg * rad
-      sangvg = SIN( angvg )
-      cangvg = COS( angvg )
-      pstarh = pstar * 0.5_wp
+
+      ! elastic damping
+      telast = relast * rdt_ice
 
       !  Diffusion coefficients.