Changeset 4672 for branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

Timestamp:

2014-06-17T17:06:59+02:00 (10 years ago)

Author:

clem

Message:

 

File:

• branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 (modified) (7 diffs)

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

@@ -40 +40 @@
 CONTAINS
 
-   SUBROUTINE lim_thd_dif( kideb , kiut , jl )
+   SUBROUTINE lim_thd_dif( kideb , kiut )
       !!------------------------------------------------------------------
       !!                ***  ROUTINE lim_thd_dif  ***
@@ -93 +93 @@
       !!------------------------------------------------------------------
       INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
-      INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number
 
       !! * Local variables
@@ -146 +145 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zdiagbis
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid   ! tridiagonal system terms
+      ! diag errors on heat
+      REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini
+      REAL(wp)                        :: zhfx_err
       !!------------------------------------------------------------------     
       ! 
@@ -155 +157 @@
       CALL wrk_alloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis  )
       CALL wrk_alloc( jpij, jkmax+2, 3, ztrid )
+
+      CALL wrk_alloc( jpij, zdq, zq_ini )
+
+      ! --- diag error on heat diffusion - PART 1 --- !
+      zdq(:) = 0._wp ; zq_ini(:) = 0._wp      
+      DO ji = kideb, kiut
+         zq_ini(ji) = ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) +  &
+            &           SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) ) 
+      END DO
 
       !------------------------------------------------------------------------------!
@@ -671 +682 @@
          DO layer  =  1, nlay_s
             DO ji = kideb , kiut
-               ii = MOD( npb(ji) - 1, jpi ) + 1
-               ij = ( npb(ji) - 1 ) / jpi + 1
                t_s_b(ji,layer) = MAX(  MIN( t_s_b(ji,layer), rtt ), 190._wp  )
                zerrit(ji)      = MAX(zerrit(ji),ABS(t_s_b(ji,layer) - ztstemp(ji,layer)))
@@ -722 +731 @@
          IF( t_su_b(ji) < rtt ) THEN  ! case T_su < 0degC
             hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji)
-         ELSE                                    ! case T_su = 0degC
+         ELSE                         ! case T_su = 0degC
             hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji)
          ENDIF
       END DO
 
+      ! --- computes sea ice energy of melting compulsory for limthd_dh --- !
+      CALL lim_thd_enmelt( kideb, kiut )
+
+      ! --- diag error on heat diffusion - PART 2 --- !
+      DO ji = kideb, kiut
+         zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) +  &
+            &                              SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) )
+         zhfx_err    = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice ) 
+         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err * a_i_b(ji)
+         ! --- correction of qns_ice and surface conduction flux --- !
+         qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err 
+         fc_su     (ji) = fc_su     (ji) - zhfx_err 
+         ! --- Heat flux at the ice surface in W.m-2 --- !
+         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
+         hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) )
+      END DO
+   
       !
       CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow )
@@ -735 +761 @@
       CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis )
       CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid )
+      CALL wrk_dealloc( jpij, zdq, zq_ini )
 
    END SUBROUTINE lim_thd_dif
+
+   SUBROUTINE lim_thd_enmelt( kideb, kiut )
+      !!-----------------------------------------------------------------------
+      !!                   ***  ROUTINE lim_thd_enmelt *** 
+      !!                 
+      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
+      !!
+      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
+      !!-------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop
+      !
+      INTEGER  ::   ji, jk   ! dummy loop indices
+      REAL(wp) ::   ztmelts, zindb  ! local scalar 
+      !!-------------------------------------------------------------------
+      !
+      DO jk = 1, nlay_i             ! Sea ice energy of melting
+         DO ji = kideb, kiut
+            ztmelts      = - tmut  * s_i_b(ji,jk) + rtt 
+            zindb        = MAX( 0._wp , SIGN( 1._wp , -(t_i_b(ji,jk) - rtt) - epsi10 ) )
+            q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) )                                             &
+               &                   + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) )   &
+               &                   - rcp  *                 ( ztmelts-rtt )  ) 
+         END DO
+      END DO
+      DO jk = 1, nlay_s             ! Snow energy of melting
+         DO ji = kideb, kiut
+            q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
+         END DO
+      END DO
+      !
+   END SUBROUTINE lim_thd_enmelt
 
 #else