MODULE limthd_lab #if defined key_lim3 !!====================================================================== !! *** MODULE limthd_lab *** !! thermodynamic lateral ablation of the ice !!====================================================================== !!---------------------------------------------------------------------- !! lim_thd_lab : vertical accr./abl. and lateral ablation of sea ice !! * Modules used USE par_oce ! ocean parameters USE phycst ! physical constants (OCE directory) USE ice_oce ! ice variables USE thd_ice USE iceini USE limistate USE in_out_manager USE ice USE par_ice USE limtab USE limicepoints IMPLICIT NONE PRIVATE !! * Routine accessibility PUBLIC lim_thd_lab ! called by lim_thd !! * Module variables REAL(wp) :: & ! constant values epsi20 = 1e-20 , & epsi13 = 1e-13 , & zzero = 0.e0 , & zone = 1.e0 !!---------------------------------------------------------------------- !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2005) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE lim_thd_lab !!------------------------------------------------------------------ !! *** ROUTINE lim_thd_lab *** !! !! ** Purpose : This routine determines the time evolution of !! ice concentration, snow thickness, ice thickness !! ice temperature after lateral ablation !! !! ** Method : Thermal lateral ablation !! - Heat comes from excessive bottom or total ablation !! Automatic lateral ablation !! - Concentration is "automatically reduced" !! when it melts at the bottom !! !! ** Action : 0) switches !! 1) Thermal lateral ablation !! Lead budget -> fscbq_1d, qfvbq_1d !! Ice heat content, new ice fraction !! 2) Automatic lateral ablation !! 3) Constraints on lead fraction !! 4) Corrections on surface and inner temperatures !! 5) Variations in ice mass and volume !! !! References : !! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646 !! Fichefet T. and M. Maqueda 1999, Clim. Dyn, 15(4), 251-268 !! !! History : !! These lines were formerly included into the vertical growth/decay routine !! original : 01-04 (M.A. Morales Maqueda, H. Goosse, T. Fichefet) !! addition : 02-08 (C. Ethe, G. Madec) !! addition : 05-12 (M. Vancoppenolle ) New routine, Multilayer code, !! Salinity variation in the ice !!------------------------------------------------------------------ !! * Arguments !! * Local variables INTEGER :: ji,jj,jk,jl,index, & ! dummy loop indices nbpb REAL(wp) :: & zihic, & zihsn, & zihicsn, & zqice, & zqicetot, & ziqf, & zdfrl, & zdvsnvol, & zdrfrl1, & zdfseqv INTEGER :: & zji, & zjj REAL(wp), DIMENSION(jpij) :: & zfrld_old, zfrld_1d WRITE(numout,*) 'lim_thd_lab' WRITE(numout,*) '~~~~~~~~~~~' DO jl = 1, jpl !--------------------------------------------------------------------------------------------------! ! 1) Select icy points !--------------------------------------------------------------------------------------------------! nbpb = 0 DO jj = 1, jpj DO ji = 1, jpi ! IF ( a_i(ji,jj,jl) .gt. zareamin ) THEN IF ( frld(ji,jj).lt. 1.0 ) THEN nbpb = nbpb + 1 npb(nbpb) = (jj - 1) * jpi + ji ENDIF END DO END DO !--------------------------------------------------------------------------------------------------! ! 2) Convert vectors !--------------------------------------------------------------------------------------------------! ! If there is no ice, do nothing. Otherwise, compute lateral ablation IF ( nbpb > 0) THEN CALL tab_2d_1d( nbpb, frld_1d (1:nbpb) , frld , jpi, jpj, npb(1:nbpb) ) CALL tab_2d_1d( nbpb, ht_i_b (1:nbpb) , ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) ) CALL tab_2d_1d( nbpb, ht_s_b (1:nbpb) , ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) ) ! no use ! CALL tab_2d_1d( nbpb, old_ht_i_b (1:nbpb) , old_ht_i(:,:,jl), jpi, jpj, npb(1:nbpb) ) ! CALL tab_2d_1d( nbpb, old_ht_s_b (1:nbpb) , old_ht_s(:,:,jl), jpi, jpj, npb(1:nbpb) ) CALL tab_2d_1d( nbpb, t_su_b (1:nbpb) , t_su(:,:,jl), jpi, jpj, npb(1:nbpb) ) CALL tab_2d_1d( nbpb, sm_i_b (1:nbpb) , sm_i(:,:,jl), jpi, jpj, npb(1:nbpb) ) DO jk = 1, nlay_s CALL tab_2d_1d( nbpb, t_s_b(1:nbpb,jk) , t_s(:,:,jk,jl), jpi, jpj, npb(1:nbpb) ) END DO DO jk = 1, nlay_i CALL tab_2d_1d( nbpb, t_i_b(1:nbpb,jk) , t_i(:,:,jk,jl), jpi, jpj, npb(1:nbpb) ) CALL tab_2d_1d( nbpb, q_i_b(1:nbpb,jk) , e_i(:,:,jk,jl), jpi, jpj, npb(1:nbpb) ) CALL tab_2d_1d( nbpb, s_i_b(1:nbpb,jk) , s_i(:,:,jk,jl), jpi, jpj, npb(1:nbpb) ) END DO CALL tab_2d_1d( nbpb, t_bo_b (1:nbpb) , t_bo , jpi, jpj, npb(1:nbpb) ) ! CALL tab_2d_1d( nbpb, thcm_1d (1:nbpb) , thcm , jpi, jpj, npb(1:nbpb) ) CALL tab_2d_1d( nbpb, qldif_1d (1:nbpb) , qldif , jpi, jpj, npb(1:nbpb) ) CALL tab_2d_1d( nbpb, rdmicif_1d (1:nbpb) , rdmicif , jpi, jpj, npb(1:nbpb) ) CALL tab_2d_1d( nbpb, fseqv_1d (1:nbpb) , fseqv , jpi, jpj, npb(1:nbpb) ) !dummy for temporary reasons CALL tab_2d_1d( nbpb, dh_i_surf (1:nbpb) , dh_i_surf2D, jpi, jpj, npb(1:nbpb) ) CALL tab_2d_1d( nbpb, dh_i_bott (1:nbpb) , dh_i_bott2D, jpi, jpj, npb(1:nbpb) ) CALL tab_2d_1d( nbpb, fstbif_1d (1:nbpb) , fstbif , jpi, jpj, npb(1:nbpb) ) CALL tab_2d_1d( nbpb, fsup (1:nbpb) , fsup2D , jpi, jpj, npb(1:nbpb) ) CALL tab_2d_1d( nbpb, focea (1:nbpb) , focea2D , jpi, jpj, npb(1:nbpb) ) !- Ice Sample points One-dimensional Grids (Useful only for debugging the code) !----------------------------------------------------------------------------------- IF (ln_nicep) THEN DO ji = 1, nbpb zji = MOD( npb(ji) - 1, jpi ) + 1 zjj = ( npb(ji) - 1 ) / jpi + 1 mask_sp_1d( ji ) = mask_sp( zji , zjj ) DO index = 1, arc_sp_num IF ( (zji .eq.arc_sp_grid(index,1)).and.(zjj .eq.arc_sp_grid(index,2)) ) THEN arc_sp_grid_1d(index) = ji*mask_sp_1d(ji) ENDIF IF ( (zji .eq.ant_sp_grid(index,1)).and.(zjj .eq.ant_sp_grid(index,2)) ) THEN ant_sp_grid_1d(index) = ji*mask_sp_1d(ji) ENDIF END DO END DO ENDIF !----------------------------------------------------------------------------------- ! Lateral ablation starts !----------------------------------------------------------------------------------- DO ji = 1, nbpb zfrld_old(ji) = frld_1d(ji) zihic = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) ) zihsn = 1.0 - MAX( zzero , SIGN( zone , -ht_s_b(ji) ) ) !--0) Case of total vertical ablation !------------------------------------- ! In the case of total ablation (all the ice ice has melted) there are only leads ! frld = 1 frld_1d(ji) = ( 1.0 - zihic ) + zihic * zfrld_old(ji) !--1) Lateral ablation !---------------------- !-- Part of solar radiation transmitted through the ice and going to the ocean !-- fscbq_1d or fscmbq in 3d routines ! fscbq_1d(ji) = ( 1.0 - zfrld_old(ji) ) * ( 1.0 - thcm_1d(ji) ) * fstbif_1d(ji) fscbq_1d(ji) = ( 1.0 - zfrld_old(ji) ) * fstbif_1d(ji) !--fstbif_1d(ji) = part of the solar radiation transmitted throught the ice !--thcm = portion of the solar radiation used for the lead budget... !--Energy remaining after excessive bottom melt or total ablation qfvbq_1d(ji) = fsup(ji) + ( 1.0 - zihic ) * focea(ji) !!fsup latent heat remaining if bottom ablation is too strong !!focea latent heat remaining in case of total ablation !--Total lead budget !!--fscbq_1d(ji) is the part of solar radiation transmitted through the ice to the ocean !! put in the lead !!--qfvbq_1d is the remaining energy after excessive bottom melt or total ablation !!--qldif contains -> formation of ice /snow-ice etc ... !!--no heat from the ocean !!--residuals from excessive bottom melt or total ablation qldif_1d(ji) = qldif_1d(ji) + qfvbq_1d(ji) + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice !-- Total heat content per unit cell area per unit mass in the snow-ice system !-- [J.kg-1.m-2] ! q_s_b(ji,1) = rcpsn*( rtt-t_s_b(ji,1) ) + xlic q_s_b(ji,1) = rhosn*cpic*( rtt-t_s_b(ji,1) ) + rhosn*lfus zqice = q_s_b(ji,1)*ht_s_b(ji) DO jk = 1, nlay_i zqice = zqice + q_i_b(ji,jk)*ht_i_b(ji)/FLOAT(nlay_i) END DO zqicetot = ( 1.0 - frld_1d(ji) ) * zqice !--The concentration of ice is reduced (frld increases) if the heat !--exchange between ice and ocean is positive ziqf = MAX( zzero , SIGN( zone , zqicetot - qldif_1d(ji) ) ) zdfrl = qldif_1d(ji) / MAX( epsi20 , zqice ) frld_1d(ji) = ( 1.0 - ziqf ) & & + ziqf * ( frld_1d(ji) + MAX( zzero , zdfrl ) ) !-- ice total energy per unit cell area per unit mass Q/delta_t !-- fltbif_1d or ffltbif in 3d Routintes fltbif_1d(ji) = ( - ( 1.0 - zfrld_old(ji) ) * zqicetot ) / rdt_ice !-->used for computing heat fluxes !-- 2) Automatic induced lateral ablation !---------------------------------------------------------------------------- !-- Hakkinen & Mellor, 1992. zdfrl = - ( dh_i_surf(ji) + dh_i_bott(ji) ) * hakspl * ( 1.0 - zfrld_old(ji) ) & / MAX( epsi13 , ht_i_b(ji) + ht_s_b(ji) * rhosn/rhoic ) zfrld_1d(ji) = frld_1d(ji) + MAX( zzero , zdfrl ) !-- 3) Constraints on lead fraction !------------------------------------------------ ! Sometimes, it happens that all the ice melts, but that all the heat content is not ! consumed. Then, we fix the lead fraction to 0.99 (the ice concentration to 0.01) ! if there is no more heat content in the ice (ziqf = 0) -> frld = 1 (only leads) ! otherwise there is still some ice (ziqf = 1) -> frld <= 0.99 zji = mod(npb(ji)-1,182)+1 zjj = ( npb(ji)-1 )/ 182 + 1 !-- 4) Update surface and internal temperature and snow/ice thicknesses !----------------------------------------------------------------------- t_su_b(ji) = t_su_b(ji) + ( 1.0 - ziqf ) * ( rtt - t_su_b(ji) ) DO jk = 1, nlay_i q_i_b(ji,jk) = q_i_b(ji,jk)*ziqf END DO zfrld_1d(ji) = ziqf * MIN( 0.99 * zone , zfrld_1d(ji) ) + ( 1 - ziqf ) zihicsn = 1.0 - MAX( zzero , SIGN( zone , - ht_i_b(ji) - ht_s_b(ji) ) ) zfrld_1d(ji) = ( 1.0 - zihicsn ) + zihicsn * zfrld_1d(ji) !-- 5) Variation of ice/snow volume and mass !------------------------------------------- dvlbq_1d(ji) = zihic * ( zfrld_old(ji) - frld_1d(ji) ) * ht_i_b(ji) rdmicif_1d(ji) = rdmicif_1d(ji) + dvlbq_1d(ji) * rhoic zdvsnvol = zihsn * ( zfrld_old(ji) - frld_1d(ji) ) * ht_s_b(ji) rdmsnif_1d(ji) = rdmsnif_1d(ji) + zdvsnvol * rhosn ! Equivalent salt flux (3) Lateral Ice melt ! ----------------------------------------- fseqv_1d(ji) = fseqv_1d(ji) + & ( soce - sm_i_b(ji) ) * & MIN( dvlbq_1d(ji) , 0.0) * rhoic / rdt_ice !-- 6) Final update !------------------ ht_s_b(ji) = ziqf * ht_s_b(ji) zdrfrl1 = ziqf * ( 1.0 - frld_1d(ji) ) / MAX( epsi20 , 1.0 - zfrld_1d(ji) ) ht_s_b(ji) = zdrfrl1 * ht_s_b(ji) ht_i_b(ji) = zdrfrl1 * ht_i_b(ji) frld_1d(ji) = zfrld_1d(ji) END DO !-4.4) Back to the geographic grid !---------------------------------------------- CALL tab_1d_2d( nbpb, frld , npb, frld_1d(1:nbpb), jpi, jpj ) CALL tab_1d_2d( nbpb, ht_i(:,:,jl), npb, ht_i_b(1:nbpb), jpi, jpj ) CALL tab_1d_2d( nbpb, ht_s(:,:,jl), npb, ht_s_b(1:nbpb), jpi, jpj ) CALL tab_1d_2d( nbpb, t_su(:,:,jl), npb, t_su_b(1:nbpb), jpi, jpj ) CALL tab_1d_2d( nbpb, sm_i(:,:,jl), npb, sm_i_b(1:nbpb), jpi, jpj ) DO jk = 1, nlay_s CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_b(1:nbpb,jk), jpi, jpj) END DO DO jk = 1, nlay_i CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_b(1:nbpb,jk), jpi, jpj) CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_b(1:nbpb,jk), jpi, jpj) CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_b(1:nbpb,jk), jpi, jpj) END DO CALL tab_1d_2d( nbpb, qldif , npb, qldif_1d (1:nbpb) , jpi, jpj ) CALL tab_1d_2d( nbpb, qfvbq , npb, qfvbq_1d (1:nbpb) , jpi, jpj ) CALL tab_1d_2d( nbpb, rdmicif , npb, rdmicif_1d(1:nbpb) , jpi, jpj ) CALL tab_1d_2d( nbpb, dmgwi , npb, dmgwi_1d (1:nbpb) , jpi, jpj ) CALL tab_1d_2d( nbpb, rdmsnif , npb, rdmsnif_1d(1:nbpb) , jpi, jpj ) ENDIF !-4.5) Update sea ice thickness !---------------------------------------------------------- DO jj = 1, jpj DO ji = 1, jpi phicif(ji,jj) = ht_i(ji,jj,1) ! ht_i(ji,jj,1) = ht_i(ji,jj,1) * ( zone - MAX( zzero, SIGN( zone, - at_i(ji,jj) ) ) ) ht_i(ji,jj,1) = ht_i(ji,jj,1) * ( zone - MAX( zzero, SIGN( zone, - ( 1.0 - frld(ji,jj) ) ) ) ) END DO END DO END DO END SUBROUTINE lim_thd_lab #else !!====================================================================== !! *** MODULE limthd_lab *** !! no sea ice model !!====================================================================== CONTAINS SUBROUTINE lim_thd_lab ! Empty routine END SUBROUTINE lim_thd_lab #endif END MODULE limthd_lab