Changeset 5167 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

Timestamp:

2015-03-24T18:35:00+01:00 (9 years ago)

Author:

clem

Message:

LIM3 bug fix. see ticket #1492 (forthcoming update) which also solve ticket #1497

File:

• trunk/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90 (modified) (7 diffs)

trunk/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

@@ -22 +22 @@
    USE lib_mpp        ! MPP library
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE sbc_oce , ONLY : sfx  ! Surface boundary condition: ocean fields
 
    IMPLICIT NONE
@@ -30 +31 @@
    PUBLIC   lim_cons_check
    PUBLIC   lim_cons_hsm
+   PUBLIC   lim_cons_final
 
    !!----------------------------------------------------------------------
@@ -167 +169 @@
       REAL(wp)                        :: zvi,   zsmv,   zei,   zfs,   zfw,   zft
       REAL(wp)                        :: zvmin, zamin, zamax 
-      REAL(wp)                        :: zconv
-
-      zconv = 1.e-9
+      REAL(wp)                        :: zvtrp, zetrp
+      REAL(wp), PARAMETER             :: zconv = 1.e-9
 
       IF( icount == 0 ) THEN
@@ -187 +188 @@
          zvi_b  = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tmask(:,:,1) )
 
-         zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) )
+         zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) * rhoic )
 
          zei_b  = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  &
@@ -210 +211 @@
             &                    * e12t(:,:) * tmask(:,:,1) ) - zvi_b ) * r1_rdtice - zfw 
 
-         zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) - zsmv_b ) * r1_rdtice + ( zfs * r1_rhoic )
+         zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) * rhoic ) - zsmv_b ) * r1_rdtice + zfs
 
          zei  =   glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  &
@@ -216 +217 @@
             &                ) * e12t(:,:) * tmask(:,:,1) * zconv ) * r1_rdtice - zei_b * r1_rdtice + zft
 
+         zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e12t(:,:) * tmask(:,:,1) ) 
+         zetrp = glob_sum( ( diag_trp_ei + diag_trp_es ) * e12t(:,:) * tmask(:,:,1) * zconv ) 
          zvmin = glob_min( v_i )
          zamax = glob_max( SUM( a_i, dim=3 ) )
          zamin = glob_min( a_i )
+
        
          IF(lwp) THEN
-            IF ( ABS( zvi    ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (',cd_routine,') = ',(zvi * rday)
-            IF ( ABS( zsmv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (',cd_routine,') = ',(zsmv * rday)
-            IF ( ABS( zei    ) >  1.e-4 ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',(zei)
-            IF ( zvmin <  -epsi10       ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',(zvmin)
+            IF ( ABS( zvi  * rday ) >  0.5 * 1.e9 ) WRITE(numout,*) 'violation volume [kg/day]     (',cd_routine,') = ',(zvi * rday)
+            IF ( ABS( zsmv * rday ) >  5.  * 1.e9 ) WRITE(numout,*) 'violation saline [psu*kg/day] (',cd_routine,') = ',(zsmv * rday)
+            IF ( ABS( zei         ) >  2.  * 1.e9 ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',(zei)
+            IF ( zvmin <  -epsi10          ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',(zvmin)
             IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > rn_amax+epsi10 ) THEN
-                                          WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax
+                                             WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax
             ENDIF
-            IF ( zamin <  -epsi10       ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
+            IF ( zamin <  -epsi10          ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
+            IF( cd_routine == 'limtrp' .AND. ABS( zvtrp * rday ) > 0.5*1.e9 ) THEN
+                                             WRITE(numout,*) 'violation vtrp [kg/day]        (',cd_routine,') = ',(zvtrp * rday)
+                                             WRITE(numout,*) 'violation etrp [GW]            (',cd_routine,') = ',(zetrp )
+            ENDIF
          ENDIF
 
@@ -234 +242 @@
 
    END SUBROUTINE lim_cons_hsm
+
+   SUBROUTINE lim_cons_final( cd_routine )
+      CHARACTER(len=*), INTENT(in)    :: cd_routine  ! name of the routine
+      REAL(wp)                        :: zhfx, zsfx, zvfx
+      REAL(wp), PARAMETER             :: zconv = 1.e-9
+
+      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - hfx_sub ) * e12t(:,:) * tmask(:,:,1) * zconv ) 
+      zsfx  = glob_sum( ( sfx + diag_smvi ) * e12t(:,:) * tmask(:,:,1) ) * rday
+      zvfx  = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e12t(:,:) * tmask(:,:,1) ) * rday 
+
+      ! if error > 1 mm / 100 years over the Arctic Basin
+      IF( ABS( zvfx ) > 0.5 * 1.e9    ) WRITE(numout,*) 'violation vfx [kg/day]       (',cd_routine,') = ',(zvfx)
+      ! if error > 1 mm / 100 years over the Arctic Basin (ice with latent heat = 3e6 J/kg) 
+      IF( ABS( zhfx ) > 2.  * 1.e9   ) WRITE(numout,*) 'violation hfx [GW]           (',cd_routine,') = ',(zhfx)
+      ! if error > 1 mm / 100 years over the Arctic Basin (ice of salinity = 10 pss)
+      IF( ABS( zsfx ) > 5.  * 1.e9   ) WRITE(numout,*) 'violation sfx [psu*kg/day]   (',cd_routine,') = ',(zsfx)
+
+   END SUBROUTINE lim_cons_final
 
 #else