Changeset 4649 for branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

Timestamp:

2014-05-27T11:28:12+02:00 (10 years ago)

Author:

clem

Message:

finalizing LIM3 heat budget conservation + multiple minor bugs corrections

File:

• branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 (modified) (9 diffs)

branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

@@ -30 +30 @@
    USE limvar          ! clem for ice thickness correction
    USE timing          ! Timing
+   USE limcons        ! conservation tests
 
    IMPLICIT NONE
@@ -72 +73 @@
       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi
       REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   zs0e
-      REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b ! Check conservation (C Rousset)
-      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax, zchk_umax ! Check errors (C Rousset)
       ! mass and salt flux (clem)
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold   ! old ice volume...
@@ -79 +78 @@
       REAL(wp), POINTER, DIMENSION(:,:)   ::   zeiold, zesold   ! old enthalpies
       REAL(wp) :: zdv, zda, zvi, zvs, zsmv, zes, zei
+      !
+      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
       !!---------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('limtrp')
@@ -87 +88 @@
 
       CALL wrk_alloc( jpi, jpj, jpl, zaiold, zhimax, zviold, zvsold )   ! clem
-
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF( ln_limdiahsb ) THEN
-         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) )
-         zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) )
-         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) )
-         zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) )
-      ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
 
       IF( numit == nstart .AND. lwp ) THEN
@@ -114 +102 @@
       IF( ln_limdyn ) THEN          !   Advection of sea ice properties   !
          !                          !-------------------------------------!
+
+         ! conservation test
+         IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
          ! mass and salt flux init (clem)
          zviold(:,:,:)  = v_i(:,:,:)
@@ -368 +360 @@
 
                  ! Update fluxes
-                  wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 
-                  wfx_snw(ji,jj) = wfx_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
+                  wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 
+                  wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
                   sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 
                   hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
@@ -425 +417 @@
 
                      ! Update mass fluxes
-                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
-                     wfx_snw(ji,jj) = wfx_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
+                     wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
+                     wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
                      sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 
                      hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
@@ -474 +466 @@
             END DO
          END DO      
+
+         ! conservation test
+         IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
       ENDIF
@@ -508 +503 @@
          END DO
       ENDIF
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF( ln_limdiahsb ) THEN
-         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
-         zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b
-         zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b
- 
-         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw 
-         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )
-         zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft
-
-         zchk_vmin = glob_min(v_i)
-         zchk_amax = glob_max(SUM(a_i,dim=3))
-         zchk_amin = glob_min(a_i)
-         zchk_umax = glob_max(SQRT(u_ice**2 + v_ice**2))
-
-         IF(lwp) THEN
-            IF ( ABS( zchk_v_i   ) >  1.e-4 ) THEN
-               WRITE(numout,*) 'violation volume [kg/day]     (limtrp) = ',(zchk_v_i * rday)
-               WRITE(numout,*) 'u_ice max [m/s]               (limtrp) = ',zchk_umax
-               WRITE(numout,*) 'number of time steps          (limtrp) =',kt
-            ENDIF
-            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limtrp) = ',(zchk_smv * rday)
-            IF ( ABS( zchk_e_i   ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limtrp) = ',(zchk_e_i)
-            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limtrp) = ',(zchk_vmin * 1.e-3)
-            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limtrp) = ',zchk_amin
-         ENDIF
-      ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
       !
       CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold )