Changeset 1858 for branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_zdf_2.F90

Timestamp:

2010-05-04T10:39:48+02:00 (14 years ago)

Author:

gm

Message:

ticket:#665 : step 1 - heat content of freezing-melting ice

File:

• branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_zdf_2.F90 (modified) (8 diffs)

branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_zdf_2.F90

@@ -4 +4 @@
    !!                thermodynamic growth and decay of the ice 
    !!======================================================================
-   !! History :  1.0  !  01-04 (LIM) Original code
-   !!            2.0  !  02-08 (C. Ethe, G. Madec) F90
+   !! History : LIM  ! 2001-04 (UCL) Original code
+   !!           2.0  ! 2002-08 (C. Ethe, G. Madec) re-writing F90
+   !!           2.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -158 +159 @@
           , zdrfrl1, zdrfrl2 &  ! tempory scalars
           , zihsn, zidhb, zihic, zihe, zihq, ziexp, ziqf, zihnf, zibmlt, ziqr, zihgnew, zind
+       REAL(wp) ::   ztmp   ! temporary scalar
        !!----------------------------------------------------------------------
 
@@ -168 +170 @@
        
        DO ji = kideb , kiut
+          ! do nothing if the snow (ice) thickness falls below its minimum thickness
           zihsn = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) )
-          zihic = MAX( zzero , SIGN( zone , hicdif - h_ice_1d(ji) ) )
-          !--computation of energy due to surface melting
-          zqcmlt(ji,1) = ( MAX ( zzero ,  &
-             &                   rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn )
-          !--computation of energy due to bottom melting
-          zqcmlt(ji,2) = ( MAX( zzero , &
-             &                  rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) &
-             &           + MAX( zzero , &
-             &                  rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) &
+          zihic = MAX( zzero , SIGN( zone , hicdif - h_ice_1d (ji) ) )
+          !--energy required to bring snow to its melting point (rt0_snow)
+          zqcmlt(ji,1) = ( MAX ( zzero , rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn )
+          !--energy required to bring ice to its melting point (rt0_ice)
+          zqcmlt(ji,2) = ( MAX( zzero , rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) )   &
+             &           + MAX( zzero , rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) )   &
              &           ) * ( 1.0 - zihic  )
-          !--limitation of  snow/ice system internal temperature
+          !--limitation of snow/ice system internal temperature
           tbif_1d(ji,1)   = MIN( rt0_snow, tbif_1d(ji,1) )
           tbif_1d(ji,2)   = MIN( rt0_ice , tbif_1d(ji,2) )
@@ -478 +478 @@
           dvsbq_1d(ji) =  ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnw_old(ji) - zsprecip(ji) )
           dvsbq_1d(ji) =  MIN( zzero , dvsbq_1d(ji) )
-          rdmsnif_1d(ji) =  rhosn * dvsbq_1d(ji)
+          ztmp = rhosn * dvsbq_1d(ji)
+          rdm_snw_1d(ji) =  ztmp
+          !--heat content of the water provided to the ocean (referenced to rt0)
+          rdq_snw_1d(ji) =  cpic * ztmp * ( rt0_snow - rt0 )
           !-- If the snow is completely melted the remaining heat is used to melt ice
           zqsn_mlt_rem  = MAX( zzero , -zhsn ) * xlsn
@@ -621 +624 @@
           !---updating new ice thickness and computing the newly formed ice mass
           zhicnew   =  zihgnew * zhicnew
-          rdmicif_1d(ji) =  rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) * rhoic
+          ztmp    =  ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) * rhoic
+          rdm_ice_1d(ji) =  rdm_ice_1d(ji) + ztmp
+          !---heat content of the water provided to the ocean (referenced to rt0)
+          ! use of rt0_ice is OK for melting ice, in case of freezing tfu_1d should be used. This is done in 9.5 section (see below)
+          rdq_ice_1d(ji) =  cpic * ztmp * ( rt0_ice - rt0 )
           !---updating new snow thickness and computing the newly formed snow mass
           zhsnfi   = zhsn + zdhsnm
           h_snow_1d(ji) = MAX( zzero , zhsnfi )
-          rdmsnif_1d(ji) =  rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsn ) * rhosn
+          ztmp = ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsn ) * rhosn
+          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp
+          !---updating the heat content of the water provided to the ocean (referenced to rt0)
+          rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( rt0_snow - rt0 )
           !--remaining energy in case of total ablation
           zqocea(ji) = - ( zihsn * xlic * zdhicm + xlsn * ( zhsnfi - h_snow_1d(ji) ) ) * ( 1.0 - frld_1d(ji) )
@@ -657 +667 @@
           tbif_1d(ji,3) =  zihgnew * ztb3 + ( 1.0 - zihgnew ) * tfu_1d(ji)
           h_ice_1d(ji)  =  zhicnew
-       END DO
+          ! update the ice heat content given to the ocean in freezing case (part from rt0_ice to tfu_1d)
+         ztmp = ( 1. - zidhb ) * rhoic * dvbbq_1d(ji)
+         rdqicif_1d(ji) = rdqicif_1d(ji) + cpic * ztmp * ( tfu_1d(ji) - rt0_ice )
+      END DO
 
 
@@ -698 +711 @@
           dmgwi_1d(ji) = dmgwi_1d(ji) + ( 1.0 -frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnnew ) * rhosn
           !---  volume change of ice and snow (used for ocean-ice freshwater flux computation)
-          rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhicnew - h_ice_1d (ji) ) * rhoic
-          rdmsnif_1d(ji) = rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn
+          ztmp = ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d (ji) ) * rhoic
+          rdm_ice_1d(ji) = rdm_ice_1d(ji) + ztmp
+          rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * ( tfu_1d(ji) - rt0 )
+!!gm BUG ??   snow ==>  only needed for nn_ice_embd == 0  (standard levitating sea-ice)
+          ztmp = ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn
+         rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp
+         rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( rt0_snow - rt0 )
 
           !---  Actualize new snow and ice thickness.
@@ -746 +764 @@
           !--variation of ice volume and ice mass 
           dvlbq_1d(ji)   = zihic * ( zfrl_old(ji) - frld_1d(ji) ) * h_ice_1d(ji)
-          rdmicif_1d(ji) = rdmicif_1d(ji) + dvlbq_1d(ji) * rhoic
+          ztmp = dvlbq_1d(ji) * rhoic
+          rdm_ice_1d(ji) = rdm_ice_1d(ji) + ztmp
+!!gm
+!!gm   This should be split in two parts:
+!!gm         1-  heat required to bring sea-ice at tfu  : this part should be added to the heat flux taken from the ocean
+!!gm                 cpic * ztmp * 0.5 * ( tbif_1d(ji,2) + tbif_1d(ji,3) - 2.* rt0_ice )
+!!gm         2-  heat content of lateral ablation referenced to rt0 : this part only put in rdq_ice_1d
+!!gm                 cpic * ztmp * ( rt0_ice - rt0 )
+!!gm   Currently we put all the heat in rdq_ice_1d
+          rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * 0.5 * ( tbif_1d(ji,2) + tbif_1d(ji,3) - 2.* rt0 )
+          !
           !--variation of snow volume and snow mass 
-          zdvsnvol    = zihsn * ( zfrl_old(ji) - frld_1d(ji) ) * h_snow_1d(ji)
-          rdmsnif_1d(ji) = rdmsnif_1d(ji) + zdvsnvol * rhosn
+          zdvsnvol = zihsn * ( zfrl_old(ji) - frld_1d(ji) ) * h_snow_1d(ji)
+          ztmp     = zdvsnvol * rhosn
+          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp
+!!gm
+!!gm   This should be split in two parts:
+!!gm         1-  heat required to bring snow at tfu  : this part should be added to the heat flux taken from the ocean
+!!gm                 cpic * ztmp * ( tbif_1d(ji,1) - rt0_snow )
+!!gm         2-  heat content of lateral ablation referenced to rt0 : this part only put in rdqicif_1d
+!!gm                 cpic * ztmp * ( rt0_snow - rt0 )
+!!gm   Currently we put all the heat in rdqicif_1d
+          rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( tbif_1d(ji,1) - rt0 )
+
           h_snow_1d(ji)  = ziqf * h_snow_1d(ji)