Changeset 4921 for branches/2014/dev_r4650_UKMO13_CICE_changes_take2/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90

Timestamp:

2014-11-28T14:59:01+01:00 (9 years ago)

Author:

timgraham

Message:

merged with revision 4879 of trunk

File:

• branches/2014/dev_r4650_UKMO13_CICE_changes_take2/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90 (modified) (8 diffs)

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

@@ -30 +30 @@
    USE lib_fortran      ! glob_sum
    USE timing          ! Timing
+   USE limcons        ! conservation tests
 
    IMPLICIT NONE
@@ -66 +67 @@
       REAL(wp), POINTER, DIMENSION(:)   ::   zmsk           ! i-averaged of tmask
       REAL(wp), POINTER, DIMENSION(:,:) ::   zu_io, zv_io   ! ice-ocean velocity
-      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
-      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
+      !
+      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
      !!---------------------------------------------------------------------
 
@@ -75 +76 @@
       CALL wrk_alloc( jpj, zind, zmsk )
 
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF (ln_limdiahsb) THEN
-         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_fw_b  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) )
-         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) )
-      ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
-
       IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only)
 
       IF( ln_limdyn ) THEN
          !
-         old_u_ice(:,:) = u_ice(:,:) * tmu(:,:)
-         old_v_ice(:,:) = v_ice(:,:) * tmv(:,:)
+         ! conservation test
+         IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
+         u_ice_b(:,:) = u_ice(:,:) * tmu(:,:)
+         v_ice_b(:,:) = v_ice(:,:) * tmv(:,:)
 
          ! Rheology (ice dynamics)
@@ -171 +164 @@
             END DO
          END DO
+         !
+         ! conservation test
+         IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
          !
       ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
@@ -224 +220 @@
       ENDIF
       !
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF (ln_limdiahsb) THEN
-         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
-         zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
- 
-         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice
-         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )
-
-         zchk_vmin = glob_min(v_i)
-         zchk_amax = glob_max(SUM(a_i,dim=3))
-         zchk_amin = glob_min(a_i)
-
-         IF(lwp) THEN
-            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limdyn) = ',(zchk_v_i * rday)
-            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limdyn) = ',(zchk_smv * rday)
-            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limdyn) = ',(zchk_vmin * 1.e-3)
-            !IF ( zchk_amax >  amax+1.e-10   ) WRITE(numout,*) 'violation a_i>amax            (limdyn) = ',zchk_amax
-            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limdyn) = ',zchk_amin
-         ENDIF
-      ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
-
       CALL wrk_dealloc( jpi, jpj, zu_io, zv_io )
       CALL wrk_dealloc( jpj, zind, zmsk )
@@ -269 +241 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicedyn/ epsd, alpha,     &
-         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
-         &                c_rhg, etamn, creepl, ecc, ahi0, &
-         &                nevp, telast, alphaevp, hminrhg
+      NAMELIST/namicedyn/ epsd, om, cw, angvg, pstar,   &
+         &                c_rhg, creepl, ecc, ahi0,     &
+         &                nevp, relast, alphaevp, hminrhg
       !!-------------------------------------------------------------------
 
@@ -289 +260 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   tolerance parameter                              epsd   = ', epsd
-         WRITE(numout,*) '   coefficient for semi-implicit coriolis           alpha  = ', alpha
-         WRITE(numout,*) '   diffusion constant for dynamics                  dm     = ', dm
-         WRITE(numout,*) '   number of sub-time steps for relaxation          nbiter = ', nbiter
-         WRITE(numout,*) '   maximum number of iterations for relaxation      nbitdr = ', nbitdr
          WRITE(numout,*) '   relaxation constant                              om     = ', om
-         WRITE(numout,*) '   maximum value for the residual of relaxation     resl   = ', resl
          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
-         WRITE(numout,*) '   minimun value for viscosity                      etamn  = ', etamn
          WRITE(numout,*) '   creep limit                                      creepl = ', creepl
          WRITE(numout,*) '   eccentricity of the elliptical yield curve       ecc    = ', ecc
          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
@@ -322 +287 @@
       pstarh = pstar * 0.5_wp
 
+      ! elastic damping
+      telast = relast * rdt_ice
+
       !  Diffusion coefficients.
       ahiu(:,:) = ahi0 * umask(:,:,1)