Changeset 4688 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

Timestamp:

2014-06-25T01:39:59+02:00 (10 years ago)

Author:

clem

Message:

new version of LIM3 with perfect conservation of heat, see ticket #1352

File:

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

trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

@@ -67 +67 @@
 
    REAL(wp) ::   epsi10 = 1.e-10_wp   !    -       -
-   REAL(wp) ::   zzero = 0.e0        !    -       -
-   REAL(wp) ::   zone  = 1.e0        !    -       -
 
    !!----------------------------------------------------------------------
@@ -113 +111 @@
                at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
                !
-               zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi10 ) ) 
+               zinda = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
                icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda  ! ice thickness
             END DO
@@ -134 +132 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - epsi10 ) ) 
-                  zindb = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi10 ) ) 
+                  zinda = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 
+                  zindb = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
                   et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                       ! snow heat content
                   smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * zinda   ! ice salinity
@@ -205 +203 @@
                DO ji = 1, jpi
                   !                                                              ! Energy of melting q(S,T) [J.m-3]
-                  zq_i    = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp) 
                   zindb   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )     ! zindb = 0 if no ice and 1 if yes
-                  zq_i    = zq_i * unit_fac * zindb                              !convert units
+                  zq_i    = zindb * e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp) 
+                  zq_i    = zq_i * unit_fac                             !convert units
                   ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt                       ! Ice layer melt temperature
                   !
@@ -231 +229 @@
                DO ji = 1, jpi
                   !Energy of melting q(S,T) [J.m-3]
-                  zq_s  = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp)
                   zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) )     ! zindb = 0 if no ice and 1 if yes
-                  zq_s  = zq_s * unit_fac * zindb                                    ! convert units
+                  zq_s  = zindb * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp)
+                  zq_s  = zq_s * unit_fac                                    ! convert units
                   !
                   t_s(ji,jj,jk,jl) = rtt + zindb * ( - zfac1 * zq_s + zfac2 )
@@ -320 +318 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  z_slope_s(ji,jj,jl) = 2._wp * sm_i(ji,jj,jl) / MAX( 0.01 , ht_i(ji,jj,jl) )
+                  z_slope_s(ji,jj,jl) = 2._wp * sm_i(ji,jj,jl) / MAX( epsi10 , ht_i(ji,jj,jl) )
                END DO
             END DO
@@ -475 +473 @@
          !
          DO ji = kideb, kiut          ! Slope of the linear profile zs_zero
-            z_slope_s(ji) = 2._wp * sm_i_b(ji) / MAX( 0.01 , ht_i_b(ji) )
+            z_slope_s(ji) = 2._wp * sm_i_b(ji) / MAX( epsi10 , ht_i_b(ji) )
          END DO