Changeset 8984 for branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3
- Timestamp:
- 2017-12-12T00:32:40+01:00 (6 years ago)
- Location:
- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 1 added
- 4 edited
-
ice.F90 (modified) (9 diffs)
-
icedyn_adv_pra.F90 (modified) (1 diff)
-
icethd_zdf.F90 (modified) (5 diffs)
-
icethd_zdf_bl99.F90 (added)
-
icevar.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r8933 r8984 154 154 !!---------------------------------------------------------------------- 155 155 ! !!** ice-generic parameters namelist (nampar) ** 156 INTEGER , PUBLIC :: jpl !: number of ice categories157 INTEGER , PUBLIC :: nlay_i !: number of ice layers158 INTEGER , PUBLIC :: nlay_s !: number of snow layers159 INTEGER , PUBLIC :: nn_monocat !: virtual ITD mono-category parameterizations (1-4) or not (0)160 LOGICAL , PUBLIC :: ln_icedyn !: flag for ice dynamics (T) or not (F)161 LOGICAL , PUBLIC :: ln_icethd !: flag for ice thermo (T) or not (F)162 REAL(wp) , PUBLIC :: rn_amax_n !: maximum ice concentration Northern hemisphere163 REAL(wp) , PUBLIC :: rn_amax_s !: maximum ice concentration Southern hemisphere164 CHARACTER(len=256), PUBLIC :: cn_icerst_in !: suffix of ice restart name (input)165 CHARACTER(len=256), PUBLIC :: cn_icerst_out !: suffix of ice restart name (output)166 CHARACTER(len=256), PUBLIC :: cn_icerst_indir !: ice restart input directory167 CHARACTER(len=256), PUBLIC :: cn_icerst_outdir !: ice restart output directory156 INTEGER , PUBLIC :: jpl !: number of ice categories 157 INTEGER , PUBLIC :: nlay_i !: number of ice layers 158 INTEGER , PUBLIC :: nlay_s !: number of snow layers 159 INTEGER , PUBLIC :: nn_monocat !: virtual ITD mono-category parameterizations (1-4) or not (0) 160 LOGICAL , PUBLIC :: ln_icedyn !: flag for ice dynamics (T) or not (F) 161 LOGICAL , PUBLIC :: ln_icethd !: flag for ice thermo (T) or not (F) 162 REAL(wp) , PUBLIC :: rn_amax_n !: maximum ice concentration Northern hemisphere 163 REAL(wp) , PUBLIC :: rn_amax_s !: maximum ice concentration Southern hemisphere 164 CHARACTER(len=256), PUBLIC :: cn_icerst_in !: suffix of ice restart name (input) 165 CHARACTER(len=256), PUBLIC :: cn_icerst_out !: suffix of ice restart name (output) 166 CHARACTER(len=256), PUBLIC :: cn_icerst_indir !: ice restart input directory 167 CHARACTER(len=256), PUBLIC :: cn_icerst_outdir !: ice restart output directory 168 168 169 169 ! !!** ice-itd namelist (namitd) ** … … 192 192 ! ! =-1 Do nothing (needs N(cat) fluxes) 193 193 ! ! = 0 Average N(cat) fluxes then apply the average over the N(cat) ice 194 ! ! = 1 Average N(cat) fluxes then redistribute over the N(cat) ice 195 ! ! using T-ice and albedo sensitivity 194 ! ! = 1 Average N(cat) fluxes then redistribute over the N(cat) ice using T-ice and albedo sensitivity 196 195 ! ! = 2 Redistribute a single flux over categories 197 198 INTEGER , PUBLIC :: nice_jules !: Choice of jules coupling 199 ! ! Associated indices 200 INTEGER , PUBLIC, PARAMETER :: np_jules_OFF = 0 !: no Jules coupling (ice thermodynamics forced via qsr and qns) 201 INTEGER , PUBLIC, PARAMETER :: np_jules_EMULE = 1 !: emulated Jules coupling via icethd_zdf.F90 (BL99) (1st round compute qcn and qsr_tr, 2nd round use it) 202 INTEGER , PUBLIC, PARAMETER :: np_jules_ACTIVE = 2 !: active Jules coupling (SM0L) (compute qcn and qsr_tr via sbcblk.F90 or sbccpl.F90) 196 INTEGER , PUBLIC :: nice_jules !: Choice of jules coupling 197 INTEGER , PUBLIC, PARAMETER :: np_jules_OFF = 0 !: no Jules coupling (ice thermodynamics forced via qsr and qns) 198 INTEGER , PUBLIC, PARAMETER :: np_jules_EMULE = 1 !: emulated Jules coupling via icethd_zdf.F90 (BL99) (1st round compute qcn and qsr_tr, 2nd round use it) 199 INTEGER , PUBLIC, PARAMETER :: np_jules_ACTIVE = 2 !: active Jules coupling (SM0L) (compute qcn and qsr_tr via sbcblk.F90 or sbccpl.F90) 200 201 ! !!** ice-vertical diffusion namelist (namthd_zdf) ** 202 LOGICAL , PUBLIC :: ln_cndi_U64 !: thermal conductivity: Untersteiner (1964) 203 LOGICAL , PUBLIC :: ln_cndi_P07 !: thermal conductivity: Pringle et al (2007) 204 REAL(wp), PUBLIC :: rn_kappa_i !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 205 REAL(wp), PUBLIC :: rn_cnd_s !: thermal conductivity of the snow [W/m/K] 203 206 204 207 ! !!** ice-salinity namelist (namthd_sal) ** … … 211 214 REAL(wp), PUBLIC :: rn_simin !: minimum ice salinity [PSU] 212 215 213 ! !!** namelist namthd_pnd216 ! !!** ice-ponds namelist (namthd_pnd) 214 217 LOGICAL , PUBLIC :: ln_pnd_H12 !: Melt ponds scheme from Holland et al 2012 215 218 LOGICAL , PUBLIC :: ln_pnd_fwb !: melt ponds store freshwater … … 255 258 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_snw !: snow-ocean mass exchange [kg.m-2.s-1] 256 259 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_snw_sni !: snow ice growth component of wfx_snw [kg.m-2.s-1] 257 ! MV MP 2016258 260 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_snw_sum !: surface melt component of wfx_snw [kg.m-2.s-1] 259 261 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_pnd !: melt pond-ocean mass exchange [kg.m-2.s-1] 260 ! END MV MP 2016261 262 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_spr !: snow precipitation on ice [kg.m-2.s-1] 262 263 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_sub !: sublimation of snow/ice [kg.m-2.s-1] … … 353 354 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sz_i !: ice salinity [PSU] 354 355 355 ! MV MP 2016356 356 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_ip !: melt pond fraction per grid cell area 357 357 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_ip !: melt pond volume per grid cell area [m] … … 361 361 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_ip !: total melt pond fraction 362 362 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vt_ip !: total melt pond volume per unit area [m] 363 ! END MV MP 2016364 363 365 364 !!---------------------------------------------------------------------- … … 376 375 !! * Ice thickness distribution variables 377 376 !!---------------------------------------------------------------------- 378 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) 379 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) 377 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hi_max !: Boundary of ice thickness categories in thickness space 378 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hi_mean !: Mean ice thickness in catgories 380 379 ! 381 380 !!---------------------------------------------------------------------- … … 385 384 ! trp '' '' '' advection (transport of ice) 386 385 ! 387 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) 388 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) 389 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) 390 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) 391 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) 392 ! 393 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) 394 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) 395 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) 396 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) 386 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_vi !: transport of ice volume 387 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_vs !: transport of snw volume 388 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_ei !: transport of ice enthalpy (W/m2) 389 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_es !: transport of snw enthalpy (W/m2) 390 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_sv !: transport of salt content 391 ! 392 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_heat !: snw/ice heat content variation [W/m2] 393 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_sice !: ice salt content variation [] 394 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vice !: ice volume variation [m/s] 395 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vsnw !: snw volume variation [m/s] 397 396 398 397 ! … … 401 400 !!---------------------------------------------------------------------- 402 401 ! Extra sea ice diagnostics to address the data request 403 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) 404 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) 405 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) 406 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) 402 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: t_si !: Temperature at Snow-ice interface (K) 403 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tm_si !: mean temperature at the snow-ice interface (K) 404 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_fc_bo !: Bottom conduction flux (W/m2) 405 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_fc_su !: Surface conduction flux (W/m2) 407 406 408 407 ! -
branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90
r8882 r8984 251 251 pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) 252 252 END DO 253 ! MV MP 2016254 253 IF ( ln_pnd_H12 ) THEN 255 pa_ip (:,:,jl) 256 pv_ip (:,:,jl) 254 pa_ip (:,:,jl) = z0ap (:,:,jl) * r1_e1e2t(:,:) 255 pv_ip (:,:,jl) = z0vp (:,:,jl) * r1_e1e2t(:,:) 257 256 ENDIF 258 ! END MV MP 2016259 257 END DO 260 258 ! -
branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
r8939 r8984 17 17 !! ice_thd_zdf : select the appropriate routine for vertical heat diffusion calculation 18 18 !! ice_thd_zdf_BL99 : 19 !! ice_thd_enmelt :20 19 !! ice_thd_zdf_init : 21 20 !!---------------------------------------------------------------------- 22 USE dom_oce ! ocean space and time domain23 USE phycst ! physical constants (ocean directory)24 USE ice ! sea-ice: variables25 USE ice 1D ! sea-ice: thermodynamics variables21 USE dom_oce ! ocean space and time domain 22 USE phycst ! physical constants (ocean directory) 23 USE ice ! sea-ice: variables 24 USE icethd_zdf_BL99 ! sea-ice: vertical diffusion (Bitz and Lipscomb, 1999) 26 25 ! 27 USE in_out_manager ! I/O manager28 USE lib_mpp ! MPP library29 USE lib_fortran ! fortran utilities (glob_sum + no signed zero)26 USE in_out_manager ! I/O manager 27 USE lib_mpp ! MPP library 28 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) 30 29 31 30 IMPLICIT NONE … … 35 34 PUBLIC ice_thd_zdf_init ! called by icestp 36 35 36 INTEGER :: nice_zdf ! Choice of the type of vertical heat diffusion formulation 37 ! ! associated indices: 38 INTEGER, PARAMETER :: np_BL99 = 1 ! Bitz and Lipscomb (1999) 39 !! INTEGER, PARAMETER :: np_XXXX = 2 40 37 41 !!** namelist (namthd_zdf) ** 38 LOGICAL :: ln_zdf_BL99 ! Heat diffusion follows Bitz and Lipscomb (1999) 39 LOGICAL :: ln_cndi_U64 ! thermal conductivity: Untersteiner (1964) 40 LOGICAL :: ln_cndi_P07 ! thermal conductivity: Pringle et al (2007) 41 REAL(wp) :: rn_kappa_i ! coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 42 REAL(wp), PUBLIC :: rn_cnd_s ! thermal conductivity of the snow [W/m/K] 42 LOGICAL :: ln_zdf_BL99 ! Heat diffusion follows Bitz and Lipscomb (1999) 43 43 44 INTEGER :: nice_zdf ! Choice of the type of vertical heat diffusion formulation45 ! ! associated indices:46 INTEGER, PARAMETER :: np_BL99 = 1 ! Bitz and Lipscomb (1999)47 48 INTEGER, PARAMETER :: np_zdf_jules_OFF = 0 ! compute all temperatures from qsr and qns49 INTEGER, PARAMETER :: np_zdf_jules_SND = 1 ! compute conductive heat flux and surface temperature from qsr and qns50 INTEGER, PARAMETER :: np_zdf_jules_RCV = 2 ! compute snow and inner ice temperatures from qcnd51 52 44 !!---------------------------------------------------------------------- 53 45 !! NEMO/ICE 4.0 , NEMO Consortium (2017) … … 67 59 SELECT CASE ( nice_zdf ) ! Choose the vertical heat diffusion solver 68 60 ! 69 ! !------------- 70 CASE( np_BL99 ) ! BL99 solver 71 ! !------------- 61 ! !-------------! 62 CASE( np_BL99 ) ! BL99 solver ! 63 ! !-------------! 72 64 SELECT CASE( nice_jules ) 73 CASE( np_jules_OFF ) ; CALL ice_thd_zdf_BL99 ( np_zdf_jules_OFF ) ! No Jules coupler 74 CASE( np_jules_EMULE ) ; CALL ice_thd_zdf_BL99 ( np_zdf_jules_SND ) ! Jules coupler is emulated (send/recieve) 75 CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV ) 76 CASE( np_jules_ACTIVE ) ; CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV ) ! Jules coupler is active (receive only) 65 ! ! No Jules coupler ==> default option 66 CASE( np_jules_OFF ) ; CALL ice_thd_zdf_BL99 ( np_jules_OFF ) 67 ! 68 ! ! Jules coupler is emulated => 1st call to get the needed fields (conduction...) 69 ! 2nd call to use these fields to calculate heat diffusion 70 CASE( np_jules_EMULE ) ; CALL ice_thd_zdf_BL99 ( np_jules_EMULE ) 71 CALL ice_thd_zdf_BL99 ( np_jules_ACTIVE ) 72 ! 73 ! ! Jules coupler is active ==> Met Office default option 74 CASE( np_jules_ACTIVE ) ; CALL ice_thd_zdf_BL99 ( np_jules_ACTIVE ) 75 ! 77 76 END SELECT 77 ! 78 78 END SELECT 79 79 … … 81 81 82 82 83 SUBROUTINE ice_thd_zdf_BL99( k_jules )84 !!-------------------------------------------------------------------85 !! *** ROUTINE ice_thd_zdf_BL99 ***86 !!87 !! ** Purpose : computes the time evolution of snow and sea-ice temperature88 !! profiles, using the original Bitz and Lipscomb (1999) algorithm89 !!90 !! ** Method : solves the heat equation diffusion with a Neumann boundary91 !! condition at the surface and a Dirichlet one at the bottom.92 !! Solar radiation is partially absorbed into the ice.93 !! The specific heat and thermal conductivities depend on ice94 !! salinity and temperature to take into account brine pocket95 !! melting. The numerical scheme is an iterative Crank-Nicolson96 !! on a non-uniform multilayer grid in the ice and snow system.97 !!98 !! The successive steps of this routine are99 !! 1. initialization of ice-snow layers thicknesses100 !! 2. Internal absorbed and transmitted radiation101 !! Then iterative procedure begins102 !! 3. Thermal conductivity103 !! 4. Kappa factors104 !! 5. specific heat in the ice105 !! 6. eta factors106 !! 7. surface flux computation107 !! 8. tridiagonal system terms108 !! 9. solving the tridiagonal system with Gauss elimination109 !! Iterative procedure ends according to a criterion on evolution110 !! of temperature111 !! 10. Fluxes at the interfaces112 !!113 !! ** Inputs / Ouputs : (global commons)114 !! surface temperature : t_su_1d115 !! ice/snow temperatures : t_i_1d, t_s_1d116 !! ice salinities : sz_i_1d117 !! number of layers in the ice/snow: nlay_i, nlay_s118 !! total ice/snow thickness : h_i_1d, h_s_1d119 !!-------------------------------------------------------------------120 INTEGER, INTENT(in) :: k_jules ! Jules coupling (0=OFF, 1=RECEIVE, 2=SEND)121 !122 INTEGER :: ji, jk ! spatial loop index123 INTEGER :: jm ! current reference number of equation124 INTEGER :: jm_mint, jm_maxt125 INTEGER :: iconv ! number of iterations in iterative procedure126 INTEGER :: iconv_max = 50 ! max number of iterations in iterative procedure127 !128 INTEGER, DIMENSION(jpij) :: jm_min ! reference number of top equation129 INTEGER, DIMENSION(jpij) :: jm_max ! reference number of bottom equation130 !131 REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system132 REAL(wp) :: zg1 = 2._wp !133 REAL(wp) :: zgamma = 18009._wp ! for specific heat134 REAL(wp) :: zbeta = 0.117_wp ! for thermal conductivity (could be 0.13)135 REAL(wp) :: zraext_s = 10._wp ! extinction coefficient of radiation in the snow136 REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity137 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered at 0C138 REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature139 REAL(wp) :: ztmelt_i ! ice melting temperature140 REAL(wp) :: zdti_max ! current maximal error on temperature141 REAL(wp) :: zcpi ! Ice specific heat142 REAL(wp) :: zhfx_err, zdq ! diag errors on heat143 REAL(wp) :: zfac ! dummy factor144 !145 REAL(wp), DIMENSION(jpij) :: isnow ! switch for presence (1) or absence (0) of snow146 REAL(wp), DIMENSION(jpij) :: ztsub ! surface temperature at previous iteration147 REAL(wp), DIMENSION(jpij) :: zh_i, z1_h_i ! ice layer thickness148 REAL(wp), DIMENSION(jpij) :: zh_s, z1_h_s ! snow layer thickness149 REAL(wp), DIMENSION(jpij) :: zqns_ice_b ! solar radiation absorbed at the surface150 REAL(wp), DIMENSION(jpij) :: zfnet ! surface flux function151 REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function152 !153 REAL(wp), DIMENSION(jpij ) :: ztsuold ! Old surface temperature in the ice154 REAL(wp), DIMENSION(jpij,nlay_i) :: ztiold ! Old temperature in the ice155 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsold ! Old temperature in the snow156 REAL(wp), DIMENSION(jpij,nlay_i) :: ztib ! Temporary temperature in the ice to check the convergence157 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsb ! Temporary temperature in the snow to check the convergence158 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity159 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice160 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice161 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zkappa_i ! Kappa factor in the ice162 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zeta_i ! Eta factor in the ice163 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow164 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow165 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zkappa_s ! Kappa factor in the snow166 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeta_s ! Eta factor in the snow167 REAL(wp), DIMENSION(jpij,nlay_i+3) :: zindterm ! 'Ind'ependent term168 REAL(wp), DIMENSION(jpij,nlay_i+3) :: zindtbis ! Temporary 'ind'ependent term169 REAL(wp), DIMENSION(jpij,nlay_i+3) :: zdiagbis ! Temporary 'dia'gonal term170 REAL(wp), DIMENSION(jpij,nlay_i+3,3) :: ztrid ! Tridiagonal system terms171 REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat172 REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat173 !174 REAL(wp) :: zfr1, zfr2, zfrqsr_tr_i ! dummy factor175 !176 ! Mono-category177 REAL(wp) :: zepsilon ! determines thres. above which computation of G(h) is done178 REAL(wp) :: zhe ! dummy factor179 REAL(wp) :: zcnd_i ! mean sea ice thermal conductivity180 !!------------------------------------------------------------------181 182 ! --- diag error on heat diffusion - PART 1 --- !183 DO ji = 1, npti184 zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + &185 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )186 END DO187 188 !------------------189 ! 1) Initialization190 !------------------191 DO ji = 1, npti192 isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - h_s_1d(ji) ) ) ! is there snow or not193 ! layer thickness194 zh_i(ji) = h_i_1d(ji) * r1_nlay_i195 zh_s(ji) = h_s_1d(ji) * r1_nlay_s196 END DO197 !198 WHERE( zh_i(1:npti) >= epsi10 ) ; z1_h_i(1:npti) = 1._wp / zh_i(1:npti)199 ELSEWHERE ; z1_h_i(1:npti) = 0._wp200 END WHERE201 !202 WHERE( zh_s(1:npti) >= epsi10 ) ; z1_h_s(1:npti) = 1._wp / zh_s(1:npti)203 ELSEWHERE ; z1_h_s(1:npti) = 0._wp204 END WHERE205 !206 ! Store initial temperatures and non solar heat fluxes207 IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND ) THEN ! OFF or SND mode208 !209 ztsub (1:npti) = t_su_1d(1:npti) ! surface temperature at iteration n-1210 ztsuold(1:npti) = t_su_1d(1:npti) ! surface temperature initial value211 zdqns_ice_b(1:npti) = dqns_ice_1d(1:npti) ! derivative of incoming nonsolar flux212 zqns_ice_b (1:npti) = qns_ice_1d(1:npti) ! store previous qns_ice_1d value213 !214 t_su_1d(1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err ) ! required to leave the choice between melting or not215 !216 ENDIF217 !218 ztsold (1:npti,:) = t_s_1d(1:npti,:) ! Old snow temperature219 ztiold (1:npti,:) = t_i_1d(1:npti,:) ! Old ice temperature220 221 !-------------222 ! 2) Radiation223 !-------------224 ! --- Transmission/absorption of solar radiation in the ice --- !225 ! zfr1 = ( 0.18 * ( 1.0 - 0.81 ) + 0.35 * 0.81 ) ! standard value226 ! zfr2 = ( 0.82 * ( 1.0 - 0.81 ) + 0.65 * 0.81 ) ! zfr2 such that zfr1 + zfr2 to equal 1227 !228 ! DO ji = 1, npti229 !230 ! zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * 10._wp ) )231 !232 ! zfrqsr_tr_i = zfr1 + zfac * zfr2 ! below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness233 ! IF ( h_s_1d(ji) >= 0.0_wp ) zfrqsr_tr_i = 0._wp ! snow fully opaque234 !235 ! qsr_ice_tr_1d(ji) = zfrqsr_tr_i * qsr_ice_1d(ji) ! transmitted solar radiation236 !237 ! zfsw(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji)238 ! zftrice(ji) = qsr_ice_tr_1d(ji)239 ! i0(ji) = zfrqsr_tr_i240 !241 ! END DO242 243 zradtr_s(1:npti,0) = qsr_ice_tr_1d(1:npti)244 DO jk = 1, nlay_s245 DO ji = 1, npti246 ! ! radiation transmitted below the layer-th snow layer247 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * zh_s(ji) * REAL(jk) )248 ! ! radiation absorbed by the layer-th snow layer249 zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)250 END DO251 END DO252 !253 zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qsr_ice_tr_1d(1:npti) * ( 1._wp - isnow(1:npti) )254 DO jk = 1, nlay_i255 DO ji = 1, npti256 ! ! radiation transmitted below the layer-th ice layer257 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )258 ! ! radiation absorbed by the layer-th ice layer259 zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)260 END DO261 END DO262 !263 ftr_ice_1d(1:npti) = zradtr_i(1:npti,nlay_i) ! record radiation transmitted below the ice264 !265 iconv = 0 ! number of iterations266 zdti_max = 1000._wp ! maximal value of error on all points267 !268 ! !============================!269 DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max ) ! Iterative procedure begins !270 ! !============================!271 iconv = iconv + 1272 !273 ztib(1:npti,:) = t_i_1d(1:npti,:)274 ztsb(1:npti,:) = t_s_1d(1:npti,:)275 !276 !--------------------------------277 ! 3) Sea ice thermal conductivity278 !--------------------------------279 IF( ln_cndi_U64 ) THEN !-- Untersteiner (1964) formula: k = k0 + beta.S/T280 !281 DO ji = 1, npti282 ztcond_i(ji,0) = rcdic + zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )283 ztcond_i(ji,nlay_i) = rcdic + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )284 END DO285 DO jk = 1, nlay_i-1286 DO ji = 1, npti287 ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / &288 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )289 END DO290 END DO291 !292 ELSEIF( ln_cndi_P07 ) THEN !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T293 !294 DO ji = 1, npti295 ztcond_i(ji,0) = rcdic + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) &296 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 )297 ztcond_i(ji,nlay_i) = rcdic + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) &298 & - 0.011_wp * ( t_bo_1d(ji) - rt0 )299 END DO300 DO jk = 1, nlay_i-1301 DO ji = 1, npti302 ztcond_i(ji,jk) = rcdic + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / &303 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) &304 & - 0.011_wp * ( 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )305 END DO306 END DO307 !308 ENDIF309 ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) )310 !311 !--- G(he) : enhancement of thermal conductivity in mono-category case312 ! Computation of effective thermal conductivity G(h)313 ! Used in mono-category case only to simulate an ITD implicitly314 ! Fichefet and Morales Maqueda, JGR 1997315 zghe(1:npti) = 1._wp316 !317 SELECT CASE ( nn_monocat )318 !319 CASE ( 1 , 3 )320 !321 zepsilon = 0.1_wp322 DO ji = 1, npti323 zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp ) ! Mean sea ice thermal conductivity324 zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i ) ! Effective thickness he (zhe)325 IF( zhe >= zepsilon * 0.5_wp * EXP(1._wp) ) THEN326 zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) ) ! G(he)327 ENDIF328 END DO329 !330 END SELECT331 !332 !-----------------333 ! 4) kappa factors334 !-----------------335 !--- Snow336 DO jk = 0, nlay_s-1337 DO ji = 1, npti338 zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji)339 END DO340 END DO341 DO ji = 1, npti ! Snow-ice interface342 zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )343 IF( zfac > epsi10 ) THEN344 zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac345 ELSE346 zkappa_s(ji,nlay_s) = 0._wp347 ENDIF348 END DO349 350 !--- Ice351 DO jk = 0, nlay_i352 DO ji = 1, npti353 zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)354 END DO355 END DO356 DO ji = 1, npti ! Snow-ice interface357 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )358 END DO359 !360 !--------------------------------------361 ! 5) Sea ice specific heat, eta factors362 !--------------------------------------363 DO jk = 1, nlay_i364 DO ji = 1, npti365 zcpi = cpic + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )366 zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi )367 END DO368 END DO369 370 DO jk = 1, nlay_s371 DO ji = 1, npti372 zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji)373 END DO374 END DO375 376 !377 !----------------------------------------!378 ! !379 ! JULES IF (OFF or SND MODE) !380 ! !381 !----------------------------------------!382 !383 384 IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND ) THEN ! OFF or SND mode385 386 !387 ! In OFF mode the original BL99 temperature computation is used388 ! (with qsr_ice, qns_ice and dqns_ice as inputs)389 !390 ! In SND mode, the computation is required to compute the conduction fluxes391 !392 393 !----------------------------394 ! 6) surface flux computation395 !----------------------------396 397 ! update of the non solar flux according to the update in T_su398 DO ji = 1, npti399 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )400 END DO401 402 DO ji = 1, npti403 zfnet(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) + qns_ice_1d(ji) ! net heat flux = net solar - transmitted solar + non solar404 END DO405 !406 !----------------------------407 ! 7) tridiagonal system terms408 !----------------------------409 !!layer denotes the number of the layer in the snow or in the ice410 !!jm denotes the reference number of the equation in the tridiagonal411 !!system, terms of tridiagonal system are indexed as following :412 !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one413 414 !!ice interior terms (top equation has the same form as the others)415 ztrid (1:npti,:,:) = 0._wp416 zindterm(1:npti,:) = 0._wp417 zindtbis(1:npti,:) = 0._wp418 zdiagbis(1:npti,:) = 0._wp419 420 DO jm = nlay_s + 2, nlay_s + nlay_i421 DO ji = 1, npti422 jk = jm - nlay_s - 1423 ztrid(ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)424 ztrid(ji,jm,2) = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )425 ztrid(ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk)426 zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)427 END DO428 END DO429 430 jm = nlay_s + nlay_i + 1431 DO ji = 1, npti432 !!ice bottom term433 ztrid(ji,jm,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)434 ztrid(ji,jm,2) = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )435 ztrid(ji,jm,3) = 0.0436 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * &437 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) )438 END DO439 440 441 DO ji = 1, npti442 ! !---------------------!443 IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells !444 ! !---------------------!445 ! snow interior terms (bottom equation has the same form as the others)446 DO jm = 3, nlay_s + 1447 jk = jm - 1448 ztrid(ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)449 ztrid(ji,jm,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )450 ztrid(ji,jm,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk)451 zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)452 END DO453 454 ! case of only one layer in the ice (ice equation is altered)455 IF ( nlay_i == 1 ) THEN456 ztrid(ji,nlay_s+2,3) = 0.0457 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)458 ENDIF459 460 IF ( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting461 462 jm_min(ji) = 1463 jm_max(ji) = nlay_i + nlay_s + 1464 465 ! surface equation466 ztrid(ji,1,1) = 0.0467 ztrid(ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)468 ztrid(ji,1,3) = zg1s * zkappa_s(ji,0)469 zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)470 471 ! first layer of snow equation472 ztrid(ji,2,1) = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)473 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )474 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1)475 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)476 477 ELSE !-- case 2 : surface is melting478 !479 jm_min(ji) = 2480 jm_max(ji) = nlay_i + nlay_s + 1481 482 ! first layer of snow equation483 ztrid(ji,2,1) = 0.0484 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )485 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1)486 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )487 ENDIF488 ! !---------------------!489 ELSE ! cells without snow !490 ! !---------------------!491 !492 IF ( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting493 !494 jm_min(ji) = nlay_s + 1495 jm_max(ji) = nlay_i + nlay_s + 1496 497 ! surface equation498 ztrid(ji,jm_min(ji),1) = 0.0499 ztrid(ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1500 ztrid(ji,jm_min(ji),3) = zkappa_i(ji,0)*zg1501 zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji)502 503 ! first layer of ice equation504 ztrid(ji,jm_min(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)505 ztrid(ji,jm_min(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )506 ztrid(ji,jm_min(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1)507 zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)508 509 ! case of only one layer in the ice (surface & ice equations are altered)510 IF ( nlay_i == 1 ) THEN511 ztrid(ji,jm_min(ji),1) = 0.0512 ztrid(ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0513 ztrid(ji,jm_min(ji),3) = zkappa_i(ji,0) * 2.0514 ztrid(ji,jm_min(ji)+1,1) = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)515 ztrid(ji,jm_min(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )516 ztrid(ji,jm_min(ji)+1,3) = 0.0517 zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )518 ENDIF519 520 ELSE !-- case 2 : surface is melting521 522 jm_min(ji) = nlay_s + 2523 jm_max(ji) = nlay_i + nlay_s + 1524 525 ! first layer of ice equation526 ztrid(ji,jm_min(ji),1) = 0.0527 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )528 ztrid(ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1)529 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * &530 & ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) )531 532 ! case of only one layer in the ice (surface & ice equations are altered)533 IF ( nlay_i == 1 ) THEN534 ztrid(ji,jm_min(ji),1) = 0.0535 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )536 ztrid(ji,jm_min(ji),3) = 0.0537 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) &538 & + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0539 ENDIF540 541 ENDIF542 ENDIF543 !544 zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))545 zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2)546 !547 END DO548 !549 !------------------------------550 ! 8) tridiagonal system solving551 !------------------------------552 ! Solve the tridiagonal system with Gauss elimination method.553 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984554 jm_maxt = 0555 jm_mint = nlay_i+5556 DO ji = 1, npti557 jm_mint = MIN(jm_min(ji),jm_mint)558 jm_maxt = MAX(jm_max(ji),jm_maxt)559 END DO560 561 DO jk = jm_mint+1, jm_maxt562 DO ji = 1, npti563 jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))564 zdiagbis(ji,jm) = ztrid(ji,jm,2) - ztrid(ji,jm,1) * ztrid(ji,jm-1,3) / zdiagbis(ji,jm-1)565 zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)566 END DO567 END DO568 569 ! ice temperatures570 DO ji = 1, npti571 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))572 END DO573 574 DO jm = nlay_i + nlay_s, nlay_s + 2, -1575 DO ji = 1, npti576 jk = jm - nlay_s - 1577 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)578 END DO579 END DO580 581 DO ji = 1, npti582 ! snow temperatures583 IF( h_s_1d(ji) > 0._wp ) THEN584 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)585 ENDIF586 ! surface temperature587 ztsub(ji) = t_su_1d(ji)588 IF( t_su_1d(ji) < rt0 ) THEN589 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * &590 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))591 ENDIF592 END DO593 !594 !--------------------------------------------------------------595 ! 9) Has the scheme converged ?, end of the iterative procedure596 !--------------------------------------------------------------597 ! check that nowhere it has started to melt598 ! zdti_max is a measure of error, it has to be under zdti_bnd599 zdti_max = 0._wp600 DO ji = 1, npti601 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )602 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )603 END DO604 605 DO jk = 1, nlay_s606 DO ji = 1, npti607 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )608 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )609 END DO610 END DO611 612 DO jk = 1, nlay_i613 DO ji = 1, npti614 ztmelt_i = -tmut * sz_i_1d(ji,jk) + rt0615 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )616 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )617 END DO618 END DO619 620 ! Compute spatial maximum over all errors621 ! note that this could be optimized substantially by iterating only the non-converging points622 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )623 !624 !----------------------------------------!625 ! !626 ! JULES IF (RCV MODE) !627 ! !628 !----------------------------------------!629 !630 631 ELSE IF ( k_jules == np_zdf_jules_RCV ) THEN ! RCV mode632 633 !634 ! In RCV mode, we use a modified BL99 solver635 ! with conduction flux (qcn_ice) as forcing term636 !637 !----------------------------638 ! 7) tridiagonal system terms639 !----------------------------640 !!layer denotes the number of the layer in the snow or in the ice641 !!jm denotes the reference number of the equation in the tridiagonal642 !!system, terms of tridiagonal system are indexed as following :643 !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one644 645 !!ice interior terms (top equation has the same form as the others)646 ztrid (1:npti,:,:) = 0._wp647 zindterm(1:npti,:) = 0._wp648 zindtbis(1:npti,:) = 0._wp649 zdiagbis(1:npti,:) = 0._wp650 651 DO jm = nlay_s + 2, nlay_s + nlay_i652 DO ji = 1, npti653 jk = jm - nlay_s - 1654 ztrid(ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)655 ztrid(ji,jm,2) = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )656 ztrid(ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk)657 zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)658 END DO659 ENDDO660 661 jm = nlay_s + nlay_i + 1662 DO ji = 1, npti663 !!ice bottom term664 ztrid(ji,jm,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)665 ztrid(ji,jm,2) = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )666 ztrid(ji,jm,3) = 0.0667 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * &668 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) )669 ENDDO670 671 672 DO ji = 1, npti673 ! !---------------------!674 IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells !675 ! !---------------------!676 ! snow interior terms (bottom equation has the same form as the others)677 DO jm = 3, nlay_s + 1678 jk = jm - 1679 ztrid(ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)680 ztrid(ji,jm,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )681 ztrid(ji,jm,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk)682 zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)683 END DO684 685 ! case of only one layer in the ice (ice equation is altered)686 IF ( nlay_i == 1 ) THEN687 ztrid(ji,nlay_s+2,3) = 0.0688 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)689 ENDIF690 691 jm_min(ji) = 2692 jm_max(ji) = nlay_i + nlay_s + 1693 694 ! first layer of snow equation695 ztrid(ji,2,1) = 0.0696 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * zkappa_s(ji,1)697 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1)698 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) )699 700 ! !---------------------!701 ELSE ! cells without snow !702 ! !---------------------!703 jm_min(ji) = nlay_s + 2704 jm_max(ji) = nlay_i + nlay_s + 1705 706 ! first layer of ice equation707 ztrid(ji,jm_min(ji),1) = 0.0708 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)709 ztrid(ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1)710 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) )711 712 ! case of only one layer in the ice (surface & ice equations are altered)713 IF ( nlay_i == 1 ) THEN714 ztrid(ji,jm_min(ji),1) = 0.0715 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)716 ztrid(ji,jm_min(ji),3) = 0.0717 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) &718 & + qcn_ice_1d(ji) )719 ENDIF720 721 ENDIF722 !723 zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))724 zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2)725 !726 END DO727 !728 !------------------------------729 ! 8) tridiagonal system solving730 !------------------------------731 ! Solve the tridiagonal system with Gauss elimination method.732 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984733 jm_maxt = 0734 jm_mint = nlay_i+5735 DO ji = 1, npti736 jm_mint = MIN(jm_min(ji),jm_mint)737 jm_maxt = MAX(jm_max(ji),jm_maxt)738 END DO739 740 DO jk = jm_mint+1, jm_maxt741 DO ji = 1, npti742 jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))743 zdiagbis(ji,jm) = ztrid(ji,jm,2) - ztrid(ji,jm,1) * ztrid(ji,jm-1,3) / zdiagbis(ji,jm-1)744 zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)745 END DO746 END DO747 748 ! ice temperatures749 DO ji = 1, npti750 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))751 END DO752 753 DO jm = nlay_i + nlay_s, nlay_s + 2, -1754 DO ji = 1, npti755 jk = jm - nlay_s - 1756 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)757 END DO758 END DO759 760 ! snow temperatures761 DO ji = 1, npti762 IF( h_s_1d(ji) > 0._wp ) THEN763 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)764 ENDIF765 END DO766 !767 !--------------------------------------------------------------768 ! 9) Has the scheme converged ?, end of the iterative procedure769 !--------------------------------------------------------------770 ! check that nowhere it has started to melt771 ! zdti_max is a measure of error, it has to be under zdti_bnd772 zdti_max = 0._wp773 774 DO jk = 1, nlay_s775 DO ji = 1, npti776 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )777 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )778 END DO779 END DO780 781 DO jk = 1, nlay_i782 DO ji = 1, npti783 ztmelt_i = -tmut * sz_i_1d(ji,jk) + rt0784 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )785 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )786 END DO787 END DO788 789 ! Compute spatial maximum over all errors790 ! note that this could be optimized substantially by iterating only the non-converging points791 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )792 793 ENDIF ! k_jules794 795 END DO ! End of the do while iterative procedure796 797 IF( ln_icectl .AND. lwp ) THEN798 WRITE(numout,*) ' zdti_max : ', zdti_max799 WRITE(numout,*) ' iconv : ', iconv800 ENDIF801 802 !803 !-----------------------------804 ! 10) Fluxes at the interfaces805 !-----------------------------806 !807 ! --- update conduction fluxes808 !809 DO ji = 1, npti810 ! ! surface ice conduction flux811 fc_su(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) &812 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji))813 ! ! bottom ice conduction flux814 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )815 END DO816 817 !818 ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- !819 !820 IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND ) THEN ! OFF or SND mode821 DO ji = 1, npti822 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji)823 END DO824 ELSE ! RCV mode825 DO ji = 1, npti826 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji)827 END DO828 ENDIF829 830 !831 ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- !832 !833 834 IF ( ( k_jules == np_zdf_jules_OFF ) .OR. ( k_jules == np_zdf_jules_RCV ) ) THEN ! OFF835 836 CALL ice_thd_enmelt837 838 ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation839 DO ji = 1, npti840 zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + &841 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )842 843 IF ( ( k_jules == np_zdf_jules_OFF ) ) THEN844 845 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC846 zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice )*a_i_1d(ji)847 ELSE ! case T_su = 0degC848 zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice )*a_i_1d(ji)849 ENDIF850 851 ELSE ! RCV CASE852 853 zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)854 855 ENDIF856 !857 ! total heat sink to be sent to the ocean858 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err859 !860 ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2861 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)862 !863 END DO864 !865 ! --- SIMIP diagnostics866 !867 DO ji = 1, npti868 !--- Conduction fluxes (positive downwards)869 diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)870 diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji) * a_i_1d(ji) / at_i_1d(ji)871 872 !--- Snow-ice interfacial temperature (diagnostic SIMIP)873 zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)874 IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN875 t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + &876 & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac877 ELSE878 t_si_1d(ji) = t_su_1d(ji)879 ENDIF880 END DO881 !882 ENDIF883 !884 !---------------------------------------------------------------------------------------885 ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes886 !---------------------------------------------------------------------------------------887 ! effective conductivity and 1st layer temperature (needed by Met Office)888 DO ji = 1, npti889 IF( h_s_1d(ji) > 0.1_wp ) THEN890 cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0)891 ELSE892 IF( h_i_1d(ji) > 0.1_wp ) THEN893 cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)894 ELSE895 cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / 0.1_wp896 ENDIF897 ENDIF898 t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)899 END DO900 !901 IF ( k_jules == np_zdf_jules_SND ) THEN ! --- Jules coupling in "SND" mode902 ! Restore temperatures to their initial values903 t_s_1d (1:npti,:) = ztsold(1:npti,:)904 t_i_1d (1:npti,:) = ztiold(1:npti,:)905 qcn_ice_1d(1:npti) = fc_su (1:npti)906 ENDIF907 !908 END SUBROUTINE ice_thd_zdf_BL99909 910 911 SUBROUTINE ice_thd_enmelt912 !!-------------------------------------------------------------------913 !! *** ROUTINE ice_thd_enmelt ***914 !!915 !! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature916 !!917 !! ** Method : Formula (Bitz and Lipscomb, 1999)918 !!-------------------------------------------------------------------919 INTEGER :: ji, jk ! dummy loop indices920 REAL(wp) :: ztmelts ! local scalar921 !!-------------------------------------------------------------------922 !923 DO jk = 1, nlay_i ! Sea ice energy of melting924 DO ji = 1, npti925 ztmelts = - tmut * sz_i_1d(ji,jk)926 t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point927 ! (sometimes dif scheme produces abnormally high temperatures)928 e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) &929 & + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) &930 & - rcp * ztmelts )931 END DO932 END DO933 DO jk = 1, nlay_s ! Snow energy of melting934 DO ji = 1, npti935 e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )936 END DO937 END DO938 !939 END SUBROUTINE ice_thd_enmelt940 941 942 83 SUBROUTINE ice_thd_zdf_init 943 84 !!----------------------------------------------------------------------- … … 976 117 WRITE(numout,*) ' thermal conductivity in the snow rn_cnd_s = ', rn_cnd_s 977 118 WRITE(numout,*) ' extinction radiation parameter in sea ice rn_kappa_i = ', rn_kappa_i 978 ENDIF 979 119 ENDIF 980 120 ! 981 121 IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN 982 122 CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' ) 983 123 ENDIF 984 985 124 ! !== set the choice of ice vertical thermodynamic formulation ==! 986 125 ioptio = 0 987 ! !--- BL99 thermo dynamics(linear liquidus + constant thermal properties)988 IF( ln_zdf_BL99 ) THEN ; ioptio = ioptio + 1 ; nice_zdf = np_BL99; ENDIF989 IF( ioptio /= 1 ) 126 IF( ln_zdf_BL99 ) THEN ; ioptio = ioptio + 1 ; nice_zdf = np_BL99 ; ENDIF ! BL99 thermodynamics (linear liquidus + constant thermal properties) 127 !! IF( ln_zdf_XXXX ) THEN ; ioptio = ioptio + 1 ; nice_zdf = np_XXXX ; ENDIF 128 IF( ioptio /= 1 ) CALL ctl_stop( 'ice_thd_init: one and only one ice thermo option has to be defined ' ) 990 129 ! 991 130 END SUBROUTINE ice_thd_zdf_init -
branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
r8906 r8984 47 47 !! ice_var_itd2 : convert N-cat to jpl-cat 48 48 !! ice_var_bv : brine volume 49 !! ice_var_enthalpy : compute ice and snow enthalpies from temperature 49 50 !!---------------------------------------------------------------------- 50 51 USE dom_oce ! ocean space and time domain … … 70 71 PUBLIC ice_var_itd2 71 72 PUBLIC ice_var_bv 73 PUBLIC ice_var_enthalpy 72 74 73 75 !!---------------------------------------------------------------------- … … 396 398 ! ! Slope of the linear profile 397 399 WHERE( h_i_1d(1:npti) > epsi20 ) ; z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti) 398 ELSEWHERE 400 ELSEWHERE ; z_slope_s(1:npti) = 0._wp 399 401 END WHERE 400 402 … … 822 824 823 825 826 SUBROUTINE ice_var_enthalpy 827 !!------------------------------------------------------------------- 828 !! *** ROUTINE ice_var_enthalpy *** 829 !! 830 !! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature 831 !! 832 !! ** Method : Formula (Bitz and Lipscomb, 1999) 833 !!------------------------------------------------------------------- 834 INTEGER :: ji, jk ! dummy loop indices 835 REAL(wp) :: ztmelts ! local scalar 836 !!------------------------------------------------------------------- 837 ! 838 DO jk = 1, nlay_i ! Sea ice energy of melting 839 DO ji = 1, npti 840 ztmelts = - tmut * sz_i_1d(ji,jk) 841 t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point 842 ! (sometimes zdf scheme produces abnormally high temperatures) 843 e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) & 844 & + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) & 845 & - rcp * ztmelts ) 846 END DO 847 END DO 848 DO jk = 1, nlay_s ! Snow energy of melting 849 DO ji = 1, npti 850 e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) 851 END DO 852 END DO 853 ! 854 END SUBROUTINE ice_var_enthalpy 855 824 856 #else 825 857 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
