Changeset 8726
- Timestamp:
- 2017-11-17T10:22:38+01:00 (6 years ago)
- Location:
- branches/2017/dev_r8183_ICEMODEL/NEMOGCM
- Files:
-
- 11 edited
-
CONFIG/SHARED/namelist_ice_lim3_ref (modified) (2 diffs)
-
NEMO/LIM_SRC_3/ice.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/ice1d.F90 (modified) (4 diffs)
-
NEMO/LIM_SRC_3/iceforcing.F90 (modified) (4 diffs)
-
NEMO/LIM_SRC_3/icestp.F90 (modified) (2 diffs)
-
NEMO/LIM_SRC_3/icethd.F90 (modified) (4 diffs)
-
NEMO/LIM_SRC_3/icethd_dh.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/icethd_zdf.F90 (modified) (20 diffs)
-
NEMO/OPA_SRC/SBC/sbc_ice.F90 (modified) (4 diffs)
-
NEMO/OPA_SRC/SBC/sbcblk.F90 (modified) (9 diffs)
-
NEMO/OPA_SRC/SBC/sbccpl.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
r8680 r8726 122 122 ! = 2 Redistribute a single flux over categories 123 123 ! ==> coupled mode only 124 nice_jules = 1 ! Jules coupling (0=OFF, 1=EMULATED, 2=ACTIVE) 124 125 / 125 126 !------------------------------------------------------------------------------ … … 140 141 ! Obs: 0.1-0.5 (Lecomte et al, JAMES 2013) 141 142 rn_kappa_i = 1.0 ! radiation attenuation coefficient in sea ice [1/m] 142 ln_dqns_i = .true. ! change the surface non-solar flux with surface temperature (T) or not (F)143 143 / 144 144 !------------------------------------------------------------------------------ -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r8678 r8726 196 196 ! ! = 2 Redistribute a single flux over categories 197 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) 203 198 204 ! !!** ice-salinity namelist (namthd_sal) ** 199 205 INTEGER , PUBLIC :: nn_icesal !: salinity configuration used in the model -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1d.F90
r8691 r8726 33 33 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ftr_ice_1d 34 34 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qsr_ice_1d 35 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fr1_i0_1d36 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fr2_i0_1d37 35 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qns_ice_1d 38 36 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: t_bo_1d 39 37 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: rn_amax_1d 38 39 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qml_ice_1d !: heat available for snow / ice surface melting [W/m2] 40 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qcn_ice_1d !: heat available for snow / ice surface sublimation¬ [W/m2] 41 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qsr_ice_tr_1d !: solar flux transmitted below the ice surface¬ [W/m2] 40 42 41 43 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_sum_1d … … 96 98 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: evap_ice_1d !: <==> the 2D evap_ice 97 99 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qprec_ice_1d !: <==> the 2D qprec_ice 98 ! ! to reintegrate longwave flux inside the ice thermodynamics99 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: i0 !: fraction of radiation transmitted to the ice100 100 101 101 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: t_su_1d !: <==> the 2D t_su … … 178 178 ALLOCATE( nptidx (jpij) , & 179 179 & qlead_1d (jpij) , ftr_ice_1d(jpij) , qsr_ice_1d (jpij) , & 180 & fr1_i0_1d(jpij) , fr2_i0_1d (jpij) , qns_ice_1d(jpij) , & 180 & qns_ice_1d(jpij) , & 181 & qml_ice_1d(jpij), qcn_ice_1d(jpij) , qsr_ice_tr_1d(jpij) , & 181 182 & t_bo_1d (jpij) , & 182 183 & hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) , & … … 194 195 & wfx_snw_sub_1d(jpij), wfx_snw_dyn_1d(jpij), wfx_ice_sub_1d(jpij), wfx_err_sub_1d(jpij) , & 195 196 & wfx_lam_1d(jpij) , wfx_dyn_1d(jpij), wfx_pnd_1d(jpij),dqns_ice_1d(jpij) , evap_ice_1d (jpij), & 196 & qprec_ice_1d(jpij), i0 (jpij) ,&197 & qprec_ice_1d(jpij), & 197 198 & sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij), & 198 199 & sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij), & -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90
r8597 r8726 99 99 !! 100 100 !! ** Action : It provides the following fields used in sea ice model: 101 !! fr1_i0 , fr2_i0 = 1sr & 2nd fraction of qsr penetration in ice [%]102 101 !! emp_oce , emp_ice = E-P over ocean and sea ice [Kg/m2/s] 103 102 !! sprecip = solid precipitation [Kg/m2/s] … … 142 141 ! 143 142 CASE( jp_blk ) !--- bulk formulation 144 CALL blk_ice_flx( t_su, alb_ice ) !145 IF( ln_mixcpl ) CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su )143 CALL blk_ice_flx( t_su, h_s, h_i, alb_ice ) ! 144 IF( ln_mixcpl ) CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i ) 146 145 IF( nn_iceflx /= 2 ) CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_iceflx ) 147 146 ! 148 147 CASE ( jp_purecpl ) !--- coupled formulation 149 CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su )148 CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i ) 150 149 IF( nn_iceflx == 2 ) CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_iceflx ) 151 150 ! … … 268 267 INTEGER :: ios, ioptio ! Local integer output status for namelist read 269 268 !! 270 NAMELIST/namforcing/ rn_cio, rn_blow_s, nn_iceflx 269 NAMELIST/namforcing/ rn_cio, rn_blow_s, nn_iceflx, nice_jules 271 270 !!------------------------------------------------------------------- 272 271 ! … … 285 284 WRITE(numout,*) '~~~~~~~~~~~~~~~' 286 285 WRITE(numout,*) ' Namelist namforcing:' 287 WRITE(numout,*) ' drag coefficient for oceanic stress rn_cio = ', rn_cio 288 WRITE(numout,*) ' coefficient for ice-lead partition of snowfall rn_blow_s = ', rn_blow_s 289 WRITE(numout,*) ' Multicategory heat flux formulation nn_iceflx = ', nn_iceflx 286 WRITE(numout,*) ' drag coefficient for oceanic stress rn_cio = ', rn_cio 287 WRITE(numout,*) ' coefficient for ice-lead partition of snowfall rn_blow_s = ', rn_blow_s 288 WRITE(numout,*) ' Multicategory heat flux formulation nn_iceflx = ', nn_iceflx 289 WRITE(numout,*) ' Jules coupling (0=no, 1=emulated, 2=active) nice_jules = ', nice_jules 290 290 ENDIF 291 291 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
r8598 r8726 153 153 !------------------------------------------------------! 154 154 ! It provides the following fields used in sea ice model: 155 ! fr1_i0 , fr2_i0 = 1sr & 2nd fraction of qsr penetration in ice [%]156 155 ! emp_oce , emp_ice = E-P over ocean and sea ice [Kg/m2/s] 157 156 ! sprecip = solid precipitation [Kg/m2/s] … … 320 319 IF(lwp) WRITE(numout,*) ' nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen' 321 320 ENDIF 322 323 324 321 ! IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN 322 ! CALL ctl_stop( 'STOP', 'par_init : if jpl=1 then nn_monocat should be between 1 and 4' ) 323 ! ENDIF 325 324 ! 326 325 IF( ln_bdy .AND. ln_icediachk ) CALL ctl_warn('par_init: online conservation check does not work with BDY') -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
r8626 r8726 26 26 USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot, sprecip, ln_cpl 27 27 USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, & 28 & fr1_i0, fr2_i028 & qml_ice, qcn_ice, qsr_ice_tr 29 29 USE ice1D ! sea-ice: thermodynamics variables 30 30 USE icethd_zdf ! sea-ice: vertical heat diffusion … … 383 383 CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d(1:npti), qprec_ice ) 384 384 CALL tab_2d_1d( npti, nptidx(1:npti), qsr_ice_1d (1:npti), qsr_ice (:,:,kl) ) 385 CALL tab_2d_1d( npti, nptidx(1:npti), fr1_i0_1d (1:npti), fr1_i0 )386 CALL tab_2d_1d( npti, nptidx(1:npti), fr2_i0_1d (1:npti), fr2_i0 )387 385 CALL tab_2d_1d( npti, nptidx(1:npti), qns_ice_1d (1:npti), qns_ice (:,:,kl) ) 388 386 CALL tab_2d_1d( npti, nptidx(1:npti), ftr_ice_1d (1:npti), ftr_ice (:,:,kl) ) … … 393 391 CALL tab_2d_1d( npti, nptidx(1:npti), fhtur_1d (1:npti), fhtur ) 394 392 CALL tab_2d_1d( npti, nptidx(1:npti), fhld_1d (1:npti), fhld ) 393 394 CALL tab_2d_1d( npti, nptidx(1:npti), qml_ice_1d (1:npti), qml_ice (:,:,kl) ) 395 CALL tab_2d_1d( npti, nptidx(1:npti), qcn_ice_1d (1:npti), qcn_ice (:,:,kl) ) 396 CALL tab_2d_1d( npti, nptidx(1:npti), qsr_ice_tr_1d(1:npti), qsr_ice_tr (:,:,kl) ) 395 397 ! 396 398 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) … … 498 500 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr ) 499 501 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam ) 500 CALL tab_ 2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd )502 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd ) 501 503 ! 502 504 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog ) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90
r8626 r8726 130 130 !------------------------------------------------------------------------------! 131 131 ! 132 DO ji = 1, npti 133 zdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 132 SELECT CASE ( nice_jules ) 133 134 CASE ( np_jules_ACTIVE ) 135 136 DO ji = 1, npti 137 zq_su(ji) = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice ) 138 END DO 139 140 CASE ( np_jules_EMULE ) 141 142 DO ji = 1, npti 143 zdum = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) ) 144 qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 145 zq_su(ji) = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice ) 146 END DO 147 148 CASE ( np_jules_OFF ) 149 150 DO ji = 1, npti 151 zdum = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) ) 152 zdum = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 153 zq_su(ji) = MAX( 0._wp, zdum * rdt_ice ) 154 END DO 155 156 END SELECT 157 158 DO ji = 1, npti 134 159 zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 135 136 zq_su (ji) = MAX( 0._wp, zdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )137 160 zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 138 161 END DO -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
r8585 r8726 34 34 LOGICAL :: ln_cndi_U64 ! thermal conductivity: Untersteiner (1964) 35 35 LOGICAL :: ln_cndi_P07 ! thermal conductivity: Pringle et al (2007) 36 REAL(wp) :: rn_cnd_s ! thermal conductivity of the snow [W/m/K]37 36 REAL(wp) :: rn_kappa_i ! coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 38 LOGICAL :: ln_dqns_i ! change non-solar surface flux with changing surface temperature (T) or not (F) 39 37 REAL(wp), PUBLIC :: rn_cnd_s ! thermal conductivity of the snow [W/m/K] 38 39 INTEGER :: nice_zdf ! Choice of the type of vertical heat diffusion formulation 40 ! ! associated indices: 41 INTEGER, PARAMETER :: np_BL99 = 1 ! Bitz and Lipscomb (1999) 42 43 INTEGER , PARAMETER :: np_zdf_jules_OFF = 0 ! compute all temperatures from qsr and qns 44 INTEGER , PARAMETER :: np_zdf_jules_SND = 1 ! compute conductive heat flux and surface temperature from qsr and qns 45 INTEGER , PARAMETER :: np_zdf_jules_RCV = 2 ! compute snow and inner ice temperatures from qcnd 46 40 47 !!---------------------------------------------------------------------- 41 48 !! NEMO/ICE 4.0 , NEMO Consortium (2017) … … 45 52 CONTAINS 46 53 47 SUBROUTINE ice_thd_zdf 54 SUBROUTINE ice_thd_zdf 55 56 !! 48 57 !!------------------------------------------------------------------- 49 58 !! *** ROUTINE ice_thd_zdf *** 50 59 !! ** Purpose : 60 !! This chooses between the appropriate routine for the 61 !! computation of vertical diffusion 62 !! 63 !!------------------------------------------------------------------- 64 !! 65 66 SELECT CASE ( nice_zdf ) ! Choose the vertical heat diffusion solver 67 68 !------------- 69 CASE( np_BL99 ) ! BL99 solver 70 !------------- 71 72 IF ( nice_jules == np_jules_OFF ) THEN ! No Jules coupler 73 74 CALL ice_thd_zdf_BL99 ( np_zdf_jules_OFF ) 75 76 ENDIF 77 78 IF ( nice_jules == np_jules_EMULE ) THEN ! Jules coupler is emulated 79 80 CALL ice_thd_zdf_BL99 ( np_zdf_jules_SND ) 81 CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV ) 82 83 ENDIF 84 85 IF ( nice_jules == np_jules_ACTIVE ) THEN ! Jules coupler is emulated 86 87 CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV ) 88 89 ENDIF 90 91 END SELECT 92 93 END SUBROUTINE ice_thd_zdf 94 95 96 97 SUBROUTINE ice_thd_zdf_BL99(k_jules) 98 !!------------------------------------------------------------------- 99 !! *** ROUTINE ice_thd_zdf_BL99 *** 100 !! ** Purpose : 51 101 !! This routine determines the time evolution of snow and sea-ice 52 !! temperature profiles. 102 !! temperature profiles, using the original Bitz and Lipscomb (1999) algorithm 103 !! 53 104 !! ** Method : 54 105 !! This is done by solving the heat equation diffusion with … … 83 134 !! total ice/snow thickness : h_i_1d, h_s_1d 84 135 !!------------------------------------------------------------------- 136 INTEGER, INTENT(in) :: k_jules ! Jules coupling (0=OFF, 1=RECEIVE, 2=SEND) 137 85 138 INTEGER :: ji, jk ! spatial loop index 86 139 INTEGER :: jm ! current reference number of equation … … 101 154 REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature 102 155 REAL(wp) :: ztmelt_i ! ice melting temperature 103 REAL(wp) :: z1_hsu104 156 REAL(wp) :: zdti_max ! current maximal error on temperature 105 157 REAL(wp) :: zcpi ! Ice specific heat … … 111 163 REAL(wp), DIMENSION(jpij) :: zh_i, z1_h_i ! ice layer thickness 112 164 REAL(wp), DIMENSION(jpij) :: zh_s, z1_h_s ! snow layer thickness 113 REAL(wp), DIMENSION(jpij) :: zfsw ! solar radiation absorbed at the surface114 165 REAL(wp), DIMENSION(jpij) :: zqns_ice_b ! solar radiation absorbed at the surface 115 166 REAL(wp), DIMENSION(jpij) :: zfnet ! surface flux function 116 167 REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function 117 REAL(wp), DIMENSION(jpij) :: zftrice ! solar radiation transmitted through the ice118 168 169 REAL(wp), DIMENSION(jpij ) :: ztsuold ! Old surface temperature in the ice 119 170 REAL(wp), DIMENSION(jpij,nlay_i) :: ztiold ! Old temperature in the ice 120 171 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsold ! Old temperature in the snow … … 136 187 REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat 137 188 REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat 189 190 REAL(wp) :: zfr1, zfr2, zfrqsr_tr_i ! dummy factor 138 191 139 192 ! Mono-category … … 166 219 ELSEWHERE ; z1_h_s(1:npti) = 0._wp 167 220 END WHERE 168 ! 169 ! temperatures 170 ztsub (1:npti) = t_su_1d(1:npti) ! temperature at the previous iteration 221 222 ! Store initial temperatures and non solar heat fluxes 223 IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND ) THEN ! OFF or SND mode 224 225 ztsub (1:npti) = t_su_1d(1:npti) ! surface temperature at iteration n-1 226 ztsuold(1:npti) = t_su_1d(1:npti) ! surface temperature initial value 227 zdqns_ice_b(1:npti) = dqns_ice_1d(1:npti) ! derivative of incoming nonsolar flux 228 zqns_ice_b (1:npti) = qns_ice_1d(1:npti) ! store previous qns_ice_1d value 229 230 t_su_1d(1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err ) ! required to leave the choice between melting or not 231 232 ENDIF 233 171 234 ztsold (1:npti,:) = t_s_1d(1:npti,:) ! Old snow temperature 172 235 ztiold (1:npti,:) = t_i_1d(1:npti,:) ! Old ice temperature 173 t_su_1d(1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err ) ! necessary 174 ! 236 175 237 !------------- 176 238 ! 2) Radiation 177 239 !------------- 178 z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0179 DO ji = 1, npti180 ! --- Computation of i0 --- !181 ! i0 describes the fraction of solar radiation which does not contribute182 ! to the surface energy budget but rather penetrates inside the ice.183 ! We assume that no radiation is transmitted through the snow184 ! If there is no no snow185 ! zfsw = (1-i0).qsr_ice is absorbed at the surface186 ! zftrice = io.qsr_ice is below the surface187 ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice188 ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover189 zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * z1_hsu ) )190 i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zfac * fr2_i0_1d(ji) )191 192 ! --- Solar radiation absorbed / transmitted at the surface --- !193 ! Derivative of the non solar flux194 zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface195 zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer196 zdqns_ice_b(ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux197 zqns_ice_b (ji) = qns_ice_1d(ji) ! store previous qns_ice_1d value198 END DO199 200 240 ! --- Transmission/absorption of solar radiation in the ice --- ! 201 zradtr_s(1:npti,0) = zftrice(1:npti) 241 ! zfr1 = ( 0.18 * ( 1.0 - 0.81 ) + 0.35 * 0.81 ) ! standard value 242 ! zfr2 = ( 0.82 * ( 1.0 - 0.81 ) + 0.65 * 0.81 ) ! zfr2 such that zfr1 + zfr2 to equal 1 243 244 ! DO ji = 1, npti 245 246 ! zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * 10._wp ) ) 247 248 ! zfrqsr_tr_i = zfr1 + zfac * zfr2 ! below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness 249 ! IF ( h_s_1d(ji) >= 0.0_wp ) zfrqsr_tr_i = 0._wp ! snow fully opaque 250 251 ! qsr_ice_tr_1d(ji) = zfrqsr_tr_i * qsr_ice_1d(ji) ! transmitted solar radiation 252 253 ! zfsw(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) 254 ! zftrice(ji) = qsr_ice_tr_1d(ji) 255 ! i0(ji) = zfrqsr_tr_i 256 257 ! END DO 258 259 zradtr_s(1:npti,0) = qsr_ice_tr_1d(1:npti) 202 260 DO jk = 1, nlay_s 203 261 DO ji = 1, npti … … 209 267 END DO 210 268 211 zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + zftrice(1:npti) * ( 1._wp - isnow(1:npti) )269 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) ) 212 270 DO jk = 1, nlay_i 213 271 DO ji = 1, npti … … 330 388 END DO 331 389 END DO 332 ! 333 !---------------------------- 334 ! 6) surface flux computation 335 !---------------------------- 336 IF ( ln_dqns_i ) THEN 337 DO ji = 1, npti 338 ! update of the non solar flux according to the update in T_su 390 391 ! 392 !----------------------------------------! 393 ! ! 394 ! JULES IF (OFF or SND MODE) ! 395 ! ! 396 !----------------------------------------! 397 ! 398 399 IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND ) THEN ! OFF or SND mode 400 401 ! 402 ! In OFF mode the original BL99 temperature computation is used 403 ! (with qsr_ice, qns_ice and dqns_ice as inputs) 404 ! 405 ! In SND mode, the computation is required to compute the conduction fluxes 406 ! 407 408 !---------------------------- 409 ! 6) surface flux computation 410 !---------------------------- 411 412 DO ji = 1, npti 413 ! update of the non solar flux according to the update in T_su 339 414 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 340 415 END DO 341 ENDIF 342 343 DO ji = 1, npti 344 zfnet(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH) 345 END DO 346 ! 347 !---------------------------- 348 ! 7) tridiagonal system terms 349 !---------------------------- 350 !!layer denotes the number of the layer in the snow or in the ice 351 !!jm denotes the reference number of the equation in the tridiagonal 352 !!system, terms of tridiagonal system are indexed as following : 353 !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 354 355 !!ice interior terms (top equation has the same form as the others) 356 ztrid (1:npti,:,:) = 0._wp 357 zindterm(1:npti,:) = 0._wp 358 zindtbis(1:npti,:) = 0._wp 359 zdiagbis(1:npti,:) = 0._wp 360 361 DO jm = nlay_s + 2, nlay_s + nlay_i 416 417 DO ji = 1, npti 418 zfnet(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) + qns_ice_1d(ji) ! net heat flux = net solar - transmitted solar + non solar 419 END DO 420 ! 421 !---------------------------- 422 ! 7) tridiagonal system terms 423 !---------------------------- 424 !!layer denotes the number of the layer in the snow or in the ice 425 !!jm denotes the reference number of the equation in the tridiagonal 426 !!system, terms of tridiagonal system are indexed as following : 427 !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 428 429 !!ice interior terms (top equation has the same form as the others) 430 ztrid (1:npti,:,:) = 0._wp 431 zindterm(1:npti,:) = 0._wp 432 zindtbis(1:npti,:) = 0._wp 433 zdiagbis(1:npti,:) = 0._wp 434 435 DO jm = nlay_s + 2, nlay_s + nlay_i 362 436 DO ji = 1, npti 363 437 jk = jm - nlay_s - 1 … … 367 441 zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 368 442 END DO 369 ENDDO370 371 jm = nlay_s + nlay_i + 1372 DO ji = 1, npti443 ENDDO 444 445 jm = nlay_s + nlay_i + 1 446 DO ji = 1, npti 373 447 !!ice bottom term 374 448 ztrid(ji,jm,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) … … 377 451 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & 378 452 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 379 ENDDO380 381 382 DO ji = 1, npti453 ENDDO 454 455 456 DO ji = 1, npti 383 457 ! !---------------------! 384 IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells !458 IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells ! 385 459 ! !---------------------! 386 460 ! snow interior terms (bottom equation has the same form as the others) 387 461 DO jm = 3, nlay_s + 1 388 389 390 391 392 462 jk = jm - 1 463 ztrid(ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 464 ztrid(ji,jm,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 465 ztrid(ji,jm,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 466 zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 393 467 END DO 394 468 395 469 ! case of only one layer in the ice (ice equation is altered) 396 470 IF ( nlay_i == 1 ) THEN 397 398 471 ztrid(ji,nlay_s+2,3) = 0.0 472 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 399 473 ENDIF 400 474 401 475 IF ( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting 402 476 403 404 405 406 407 408 409 410 411 412 413 414 415 416 477 jm_min(ji) = 1 478 jm_max(ji) = nlay_i + nlay_s + 1 479 480 ! surface equation 481 ztrid(ji,1,1) = 0.0 482 ztrid(ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0) 483 ztrid(ji,1,3) = zg1s * zkappa_s(ji,0) 484 zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 485 486 ! first layer of snow equation 487 ztrid(ji,2,1) = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1) 488 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 489 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) 490 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 417 491 418 492 ELSE !-- case 2 : surface is melting 419 420 421 422 423 424 425 426 427 428 493 ! 494 jm_min(ji) = 2 495 jm_max(ji) = nlay_i + nlay_s + 1 496 497 ! first layer of snow equation 498 ztrid(ji,2,1) = 0.0 499 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 500 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 501 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * & 502 & ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 429 503 ENDIF 430 504 ! !---------------------! … … 433 507 ! 434 508 IF ( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 509 ! 510 jm_min(ji) = nlay_s + 1 511 jm_max(ji) = nlay_i + nlay_s + 1 512 513 ! surface equation 514 ztrid(ji,jm_min(ji),1) = 0.0 515 ztrid(ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1 516 ztrid(ji,jm_min(ji),3) = zkappa_i(ji,0)*zg1 517 zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji) 518 519 ! first layer of ice equation 520 ztrid(ji,jm_min(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) 521 ztrid(ji,jm_min(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 522 ztrid(ji,jm_min(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1) 523 zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 524 525 ! case of only one layer in the ice (surface & ice equations are altered) 526 IF ( nlay_i == 1 ) THEN 527 ztrid(ji,jm_min(ji),1) = 0.0 528 ztrid(ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 529 ztrid(ji,jm_min(ji),3) = zkappa_i(ji,0) * 2.0 530 ztrid(ji,jm_min(ji)+1,1) = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 531 ztrid(ji,jm_min(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 532 ztrid(ji,jm_min(ji)+1,3) = 0.0 533 zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * & 534 & ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 535 ENDIF 462 536 463 537 ELSE !-- case 2 : surface is melting 464 538 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 539 jm_min(ji) = nlay_s + 2 540 jm_max(ji) = nlay_i + nlay_s + 1 541 542 ! first layer of ice equation 543 ztrid(ji,jm_min(ji),1) = 0.0 544 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 545 ztrid(ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 546 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * & 547 & ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 548 549 ! case of only one layer in the ice (surface & ice equations are altered) 550 IF ( nlay_i == 1 ) THEN 551 ztrid(ji,jm_min(ji),1) = 0.0 552 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 553 ztrid(ji,jm_min(ji),3) = 0.0 554 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) & 555 & + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0 556 ENDIF 483 557 484 558 ENDIF … … 488 562 zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 489 563 ! 490 END DO491 !492 !------------------------------493 ! 8) tridiagonal system solving494 !------------------------------495 ! Solve the tridiagonal system with Gauss elimination method.496 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984497 jm_maxt = 0498 jm_mint = nlay_i+5499 DO ji = 1, npti564 END DO 565 ! 566 !------------------------------ 567 ! 8) tridiagonal system solving 568 !------------------------------ 569 ! Solve the tridiagonal system with Gauss elimination method. 570 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 571 jm_maxt = 0 572 jm_mint = nlay_i+5 573 DO ji = 1, npti 500 574 jm_mint = MIN(jm_min(ji),jm_mint) 501 575 jm_maxt = MAX(jm_max(ji),jm_maxt) 502 END DO503 504 DO jk = jm_mint+1, jm_maxt576 END DO 577 578 DO jk = jm_mint+1, jm_maxt 505 579 DO ji = 1, npti 506 580 jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) … … 508 582 zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 509 583 END DO 510 END DO511 512 DO ji = 1, npti584 END DO 585 586 DO ji = 1, npti 513 587 ! ice temperatures 514 588 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 515 END DO516 517 DO jm = nlay_i + nlay_s, nlay_s + 2, -1589 END DO 590 591 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 518 592 DO ji = 1, npti 519 593 jk = jm - nlay_s - 1 520 594 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 521 595 END DO 522 END DO523 524 DO ji = 1, npti596 END DO 597 598 DO ji = 1, npti 525 599 ! snow temperatures 526 600 IF( h_s_1d(ji) > 0._wp ) THEN 527 601 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) & 528 602 & / zdiagbis(ji,nlay_s+1) 529 603 ENDIF 530 604 ! surface temperature … … 532 606 IF( t_su_1d(ji) < rt0 ) THEN 533 607 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 534 535 ENDIF 536 END DO537 !538 !--------------------------------------------------------------539 ! 9) Has the scheme converged ?, end of the iterative procedure540 !--------------------------------------------------------------541 ! check that nowhere it has started to melt542 ! zdti_max is a measure of error, it has to be under zdti_bnd543 zdti_max = 0._wp544 DO ji = 1, npti608 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 609 ENDIF 610 END DO 611 ! 612 !-------------------------------------------------------------- 613 ! 9) Has the scheme converged ?, end of the iterative procedure 614 !-------------------------------------------------------------- 615 ! check that nowhere it has started to melt 616 ! zdti_max is a measure of error, it has to be under zdti_bnd 617 zdti_max = 0._wp 618 DO ji = 1, npti 545 619 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 546 620 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 547 END DO548 549 DO jk = 1, nlay_s621 END DO 622 623 DO jk = 1, nlay_s 550 624 DO ji = 1, npti 551 625 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 552 626 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 553 627 END DO 554 END DO555 556 DO jk = 1, nlay_i628 END DO 629 630 DO jk = 1, nlay_i 557 631 DO ji = 1, npti 558 632 ztmelt_i = -tmut * sz_i_1d(ji,jk) + rt0 … … 560 634 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 561 635 END DO 562 END DO 563 564 ! Compute spatial maximum over all errors 565 ! note that this could be optimized substantially by iterating only the non-converging points 566 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 567 636 END DO 637 638 ! Compute spatial maximum over all errors 639 ! note that this could be optimized substantially by iterating only the non-converging points 640 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 641 ! 642 !----------------------------------------! 643 ! ! 644 ! JULES IF (RCV MODE) ! 645 ! ! 646 !----------------------------------------! 647 ! 648 649 ELSE IF ( k_jules == np_zdf_jules_RCV ) THEN ! RCV mode 650 651 ! 652 ! In RCV mode, we use a modified BL99 solver 653 ! with conduction flux (qcn_ice) as forcing term 654 ! 655 !---------------------------- 656 ! 7) tridiagonal system terms 657 !---------------------------- 658 !!layer denotes the number of the layer in the snow or in the ice 659 !!jm denotes the reference number of the equation in the tridiagonal 660 !!system, terms of tridiagonal system are indexed as following : 661 !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 662 663 !!ice interior terms (top equation has the same form as the others) 664 ztrid (1:npti,:,:) = 0._wp 665 zindterm(1:npti,:) = 0._wp 666 zindtbis(1:npti,:) = 0._wp 667 zdiagbis(1:npti,:) = 0._wp 668 669 DO jm = nlay_s + 2, nlay_s + nlay_i 670 DO ji = 1, npti 671 jk = jm - nlay_s - 1 672 ztrid(ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 673 ztrid(ji,jm,2) = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 674 ztrid(ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk) 675 zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 676 END DO 677 ENDDO 678 679 jm = nlay_s + nlay_i + 1 680 DO ji = 1, npti 681 !!ice bottom term 682 ztrid(ji,jm,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) 683 ztrid(ji,jm,2) = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 684 ztrid(ji,jm,3) = 0.0 685 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & 686 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 687 ENDDO 688 689 690 DO ji = 1, npti 691 ! !---------------------! 692 IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells ! 693 ! !---------------------! 694 ! snow interior terms (bottom equation has the same form as the others) 695 DO jm = 3, nlay_s + 1 696 jk = jm - 1 697 ztrid(ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 698 ztrid(ji,jm,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 699 ztrid(ji,jm,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 700 zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 701 END DO 702 703 ! case of only one layer in the ice (ice equation is altered) 704 IF ( nlay_i == 1 ) THEN 705 ztrid(ji,nlay_s+2,3) = 0.0 706 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 707 ENDIF 708 709 jm_min(ji) = 2 710 jm_max(ji) = nlay_i + nlay_s + 1 711 712 ! first layer of snow equation 713 ztrid(ji,2,1) = 0.0 714 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * zkappa_s(ji,1) 715 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 716 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * & 717 & ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 718 719 ! !---------------------! 720 ELSE ! cells without snow ! 721 ! !---------------------! 722 723 jm_min(ji) = nlay_s + 2 724 jm_max(ji) = nlay_i + nlay_s + 1 725 726 ! first layer of ice equation 727 ztrid(ji,jm_min(ji),1) = 0.0 728 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) 729 ztrid(ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 730 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * & 731 & ( zradab_i(ji,1) + qcn_ice_1d(ji) ) 732 733 ! case of only one layer in the ice (surface & ice equations are altered) 734 IF ( nlay_i == 1 ) THEN 735 ztrid(ji,jm_min(ji),1) = 0.0 736 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) 737 ztrid(ji,jm_min(ji),3) = 0.0 738 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) & 739 & + qcn_ice_1d(ji) ) 740 741 ENDIF 742 743 ENDIF 744 ! 745 zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji)) 746 zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 747 ! 748 END DO 749 ! 750 !------------------------------ 751 ! 8) tridiagonal system solving 752 !------------------------------ 753 ! Solve the tridiagonal system with Gauss elimination method. 754 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 755 jm_maxt = 0 756 jm_mint = nlay_i+5 757 DO ji = 1, npti 758 jm_mint = MIN(jm_min(ji),jm_mint) 759 jm_maxt = MAX(jm_max(ji),jm_maxt) 760 END DO 761 762 DO jk = jm_mint+1, jm_maxt 763 DO ji = 1, npti 764 jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) 765 zdiagbis(ji,jm) = ztrid(ji,jm,2) - ztrid(ji,jm,1) * ztrid(ji,jm-1,3) / zdiagbis(ji,jm-1) 766 zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 767 END DO 768 END DO 769 770 DO ji = 1, npti 771 ! ice temperatures 772 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 773 END DO 774 775 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 776 DO ji = 1, npti 777 jk = jm - nlay_s - 1 778 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 779 END DO 780 END DO 781 782 DO ji = 1, npti 783 ! snow temperatures 784 IF( h_s_1d(ji) > 0._wp ) THEN 785 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) & 786 & / zdiagbis(ji,nlay_s+1) 787 ENDIF 788 END DO 789 ! 790 !-------------------------------------------------------------- 791 ! 9) Has the scheme converged ?, end of the iterative procedure 792 !-------------------------------------------------------------- 793 ! check that nowhere it has started to melt 794 ! zdti_max is a measure of error, it has to be under zdti_bnd 795 zdti_max = 0._wp 796 797 DO jk = 1, nlay_s 798 DO ji = 1, npti 799 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 800 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 801 END DO 802 END DO 803 804 DO jk = 1, nlay_i 805 DO ji = 1, npti 806 ztmelt_i = -tmut * sz_i_1d(ji,jk) + rt0 807 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) 808 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 809 END DO 810 END DO 811 812 ! Compute spatial maximum over all errors 813 ! note that this could be optimized substantially by iterating only the non-converging points 814 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 815 816 ENDIF ! k_jules 817 568 818 END DO ! End of the do while iterative procedure 569 819 570 820 IF( ln_icectl .AND. lwp ) THEN 571 821 WRITE(numout,*) ' zdti_max : ', zdti_max 572 822 WRITE(numout,*) ' iconv : ', iconv 573 823 ENDIF 824 574 825 ! 575 826 !----------------------------- 576 827 ! 10) Fluxes at the interfaces 577 828 !----------------------------- 829 ! 830 ! --- update conduction fluxes 831 ! 578 832 DO ji = 1, npti 579 833 ! ! surface ice conduction flux 580 834 fc_su(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) & 581 835 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji)) 582 836 ! ! bottom ice conduction flux 583 837 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 584 838 END DO 585 586 ! --- computes sea ice energy of melting compulsory for icethd_dh --- ! 587 CALL ice_thd_enmelt 588 589 ! --- diagnose the change in non-solar flux due to surface temperature change --- ! 590 IF ( ln_dqns_i ) THEN 839 840 ! 841 ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- ! 842 ! 843 DO ji = 1, npti 844 IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND ) THEN 845 ! OFF or SND mode 846 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 847 ELSE ! RCV mode 848 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji) 849 ENDIF 850 END DO 851 852 ! 853 ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- ! 854 ! 855 856 IF ( ( k_jules == np_zdf_jules_OFF ) .OR. ( k_jules == np_zdf_jules_RCV ) ) THEN ! OFF 857 858 CALL ice_thd_enmelt 859 860 ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 591 861 DO ji = 1, npti 592 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 862 zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + & 863 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 864 865 IF ( ( k_jules == np_zdf_jules_OFF ) ) THEN 866 867 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 868 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) 869 ELSE ! case T_su = 0degC 870 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) 871 ENDIF 872 873 ELSE ! RCV CASE 874 875 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) 876 877 ENDIF 878 879 ! total heat sink to be sent to the ocean 880 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 881 882 ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2 883 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 884 885 END DO 886 887 ! 888 ! --- SIMIP diagnostics 889 ! 890 891 DO ji = 1, npti 892 !--- Conduction fluxes (positive downwards) 893 diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) 894 diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji) * a_i_1d(ji) / at_i_1d(ji) 895 896 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 897 zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 898 IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN 899 t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + & 900 & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac 901 ELSE 902 t_si_1d(ji) = t_su_1d(ji) 903 ENDIF 593 904 END DO 594 END IF 595 596 ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 597 ! hfx_dif = Heat flux used to warm/cool ice in W.m-2 598 ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 599 DO ji = 1, npti 600 zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + & 601 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 602 603 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 604 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) 605 ELSE ! case T_su = 0degC 606 zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 607 ENDIF 608 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 609 610 ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation) 611 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 612 613 END DO 614 615 ! --- Diagnostics SIMIP --- ! 616 DO ji = 1, npti 617 !--- Conduction fluxes (positive downwards) 618 diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) 619 diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji) * a_i_1d(ji) / at_i_1d(ji) 620 621 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 622 zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 623 IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN 624 t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + & 625 & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac 626 ELSE 627 t_si_1d(ji) = t_su_1d(ji) 628 ENDIF 629 END DO 630 ! 631 END SUBROUTINE ice_thd_zdf 905 906 ENDIF 907 ! 908 !--------------------------------------------------------------------------------------- 909 ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes 910 !--------------------------------------------------------------------------------------- 911 ! 912 IF ( k_jules == np_zdf_jules_SND ) THEN ! --- Jules coupling in "SND" mode 913 914 ! Restore temperatures to their initial values 915 t_s_1d(1:npti,:) = ztsold (1:npti,:) 916 t_i_1d(1:npti,:) = ztiold (1:npti,:) 917 qcn_ice_1d(1:npti) = fc_su(1:npti) 918 919 ENDIF 920 921 END SUBROUTINE ice_thd_zdf_BL99 922 632 923 633 924 … … 663 954 664 955 956 665 957 SUBROUTINE ice_thd_zdf_init 666 958 !!----------------------------------------------------------------------- … … 675 967 !! ** input : Namelist namthd_zdf 676 968 !!------------------------------------------------------------------- 677 INTEGER :: ios ! Local integer output status for namelist read969 INTEGER :: ios, ioptio ! Local integer output status for namelist read 678 970 !! 679 NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i , ln_dqns_i971 NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i 680 972 !!------------------------------------------------------------------- 681 973 ! … … 694 986 WRITE(numout,*) '~~~~~~~~~~~~~~~~' 695 987 WRITE(numout,*) ' Namelist namthd_zdf:' 696 WRITE(numout,*) ' Diffusion follows a Bitz and Lipscomb (1999)ln_zdf_BL99 = ', ln_zdf_BL99988 WRITE(numout,*) ' Bitz and Lipscomb (1999) formulation ln_zdf_BL99 = ', ln_zdf_BL99 697 989 WRITE(numout,*) ' thermal conductivity in the ice (Untersteiner 1964) ln_cndi_U64 = ', ln_cndi_U64 698 990 WRITE(numout,*) ' thermal conductivity in the ice (Pringle et al 2007) ln_cndi_P07 = ', ln_cndi_P07 699 991 WRITE(numout,*) ' thermal conductivity in the snow rn_cnd_s = ', rn_cnd_s 700 992 WRITE(numout,*) ' extinction radiation parameter in sea ice rn_kappa_i = ', rn_kappa_i 701 WRITE(numout,*) ' change the surface non-solar flux with Tsu or not ln_dqns_i = ', ln_dqns_i702 993 ENDIF 994 703 995 ! 704 996 IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN 705 997 CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' ) 706 998 ENDIF 999 1000 ! !== set the choice of ice vertical thermodynamic formulation ==! 1001 ioptio = 0 1002 ! !--- BL99 thermo dynamics (linear liquidus + constant thermal properties) 1003 IF( ln_zdf_BL99 ) THEN ; ioptio = ioptio + 1 ; nice_zdf = np_BL99 ; ENDIF 1004 IF( ioptio /= 1 ) CALL ctl_stop( 'ice_thd_init: one and only one ice thermo option has to be defined ' ) 707 1005 ! 708 1006 END SUBROUTINE ice_thd_zdf_init … … 711 1009 !!---------------------------------------------------------------------- 712 1010 !! Default option Dummy Module No ESIM sea-ice model 713 !!--------------------------------------------------------------------- -1011 !!--------------------------------------------------------------------- 714 1012 #endif 715 1013 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90
r8563 r8726 48 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: alb_ice !: ice albedo [-] 49 49 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qml_ice !: heat available for snow / ice surface melting [W/m2] 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qcn_ice !: heat conduction flux in the layer below surface [W/m2] 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qsr_ice_tr !: solar flux transmitted below the ice surface [W/m2] 53 50 54 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau_ice !: atmos-ice u-stress. VP: I-pt ; EVP: U,V-pts [N/m2] 51 55 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtau_ice !: atmos-ice v-stress. VP: I-pt ; EVP: U,V-pts [N/m2] 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fr1_i0 !: Solar surface transmission parameter, thick ice [-]53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fr2_i0 !: Solar surface transmission parameter, thin ice [-]54 56 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: emp_ice !: sublimation - precip over sea ice [kg/m2/s] 55 57 … … 123 125 ALLOCATE( qns_ice (jpi,jpj,jpl) , qsr_ice (jpi,jpj,jpl) , & 124 126 & qla_ice (jpi,jpj,jpl) , dqla_ice(jpi,jpj,jpl) , & 125 & dqns_ice(jpi,jpj,jpl) , tn_ice (jpi,jpj,jpl) , alb_ice (jpi,jpj,jpl) , & 127 & dqns_ice(jpi,jpj,jpl) , tn_ice (jpi,jpj,jpl) , alb_ice (jpi,jpj,jpl) , & 128 & qml_ice(jpi,jpj,jpl) , qcn_ice(jpi,jpj,jpl) , qsr_ice_tr(jpi,jpj,jpl), & 126 129 & utau_ice(jpi,jpj) , vtau_ice(jpi,jpj) , wndm_ice(jpi,jpj) , & 127 & fr1_i0 (jpi,jpj) , fr2_i0 (jpi,jpj) , &128 130 & evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) , & 129 131 & qemp_ice(jpi,jpj) , qevap_ice(jpi,jpj,jpl) , qemp_oce (jpi,jpj) , & … … 139 141 a_i(jpi,jpj,ncat) , topmelt(jpi,jpj,ncat) , botmelt(jpi,jpj,ncat) , & 140 142 STAT= ierr(2) ) 141 IF( ln_cpl ) ALLOCATE( u_ice(jpi,jpj) , fr1_i0(jpi,jpj) ,tn_ice (jpi,jpj,1) , &142 & v_ice(jpi,jpj) , fr2_i0(jpi,jpj) ,alb_ice(jpi,jpj,1) , &143 IF( ln_cpl ) ALLOCATE( u_ice(jpi,jpj) , tn_ice (jpi,jpj,1) , & 144 & v_ice(jpi,jpj) , alb_ice(jpi,jpj,1) , & 143 145 & emp_ice(jpi,jpj) , qns_ice(jpi,jpj,1) , dqns_ice(jpi,jpj,1) , & 144 146 & STAT= ierr(3) ) … … 168 170 REAL , PUBLIC, PARAMETER :: cldf_ice = 0.81 !: cloud fraction over sea ice, summer CLIO value [-] 169 171 INTEGER , PUBLIC, PARAMETER :: jpl = 1 170 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_ice, v_ice ,fr1_i0,fr2_i0! jpi, jpj172 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_ice, v_ice ! jpi, jpj 171 173 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tn_ice, alb_ice, qns_ice, dqns_ice ! (jpi,jpj,jpl) 172 174 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
r8589 r8726 17 17 !! ==> based on AeroBulk (http://aerobulk.sourceforge.net/) 18 18 !! 4.0 ! 2016-10 (G. Madec) introduce a sbc_blk_init routine 19 !! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce Jules emulator (M. Vancoppenolle) 19 20 !!---------------------------------------------------------------------- 20 21 … … 42 43 USE ice , ONLY : u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su 43 44 USE icethd_dh ! for CALL ice_thd_snwblow 45 USE icethd_zdf, ONLY: rn_cnd_s ! for blk_ice_qcn 44 46 #endif 45 47 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) … … 64 66 PUBLIC blk_ice_tau ! routine called in icestp module 65 67 PUBLIC blk_ice_flx ! routine called in icestp module 66 #endif 68 PUBLIC blk_ice_qcn ! routine called in icestp module 69 #endif 67 70 68 71 !!Lolo: should ultimately be moved in the module with all physical constants ? … … 697 700 698 701 699 SUBROUTINE blk_ice_flx( ptsu, p alb )702 SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb ) 700 703 !!--------------------------------------------------------------------- 701 704 !! *** ROUTINE blk_ice_flx *** … … 710 713 !!--------------------------------------------------------------------- 711 714 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature 715 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness 716 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness 712 717 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies) 713 718 !! … … 716 721 REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - - 717 722 REAL(wp) :: zztmp, z1_lsub ! - - 723 REAL(wp) :: zfrqsr_tr_i ! sea ice shortwave fraction transmitted below through the ice 724 REAL(wp) :: zfr1, zfr2, zfac ! local variables 718 725 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qlw ! long wave heat flux over ice 719 726 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qsb ! sensible heat flux over ice … … 722 729 REAL(wp), DIMENSION(:,:) , POINTER :: zevap, zsnw ! evaporation and snw distribution after wind blowing (LIM3) 723 730 REAL(wp), DIMENSION(:,:) , POINTER :: zrhoa 731 724 732 !!--------------------------------------------------------------------- 725 733 ! … … 830 838 CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 831 839 832 !-------------------------------------------------------------------- 833 ! FRACTIONs of net shortwave radiation which is not absorbed in the 834 ! thin surface layer and penetrates inside the ice cover 835 ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 836 ! 837 fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 838 fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 839 ! 840 ! 840 ! --- absorbed and transmitted shortwave radiation (W/m2) --- ! 841 ! 842 ! former coding was 843 ! fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 844 ! fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 845 846 ! --- surface transmission parameter (i0, Grenfell Maykut 77) --- ! 847 zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) ! standard value 848 zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) ! zfr2 such that zfr1 + zfr2 to equal 1 849 850 qsr_ice_tr(:,:,:) = 0._wp 851 852 DO jl = 1, jpl 853 DO jj = 1, jpj 854 DO ji = 1, jpi 855 856 zfac = MAX( 0._wp , 1._wp - phi(ji,jj,jl) * 10._wp ) ! linear weighting factor: =0 for phi=0, =1 for phi = 0.1 857 zfrqsr_tr_i = zfr1 + zfac * zfr2 ! below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness 858 859 IF ( phs(ji,jj,jl) <= 0.0_wp ) THEN 860 zfrqsr_tr_i = zfr1 + zfac * zfr2 861 ELSE 862 zfrqsr_tr_i = 0._wp ! snow fully opaque 863 ENDIF 864 865 qsr_ice_tr(ji,jj,jl) = zfrqsr_tr_i * qsr_ice(ji,jj,jl) ! transmitted solar radiation 866 867 END DO 868 END DO 869 END DO 870 871 841 872 IF(ln_ctl) THEN 842 873 CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=jpl) … … 854 885 855 886 END SUBROUTINE blk_ice_flx 887 888 889 890 SUBROUTINE blk_ice_qcn( k_monocat, ptsu, ptb, phs, phi ) 891 892 !!--------------------------------------------------------------------- 893 !! *** ROUTINE blk_ice_qcn *** 894 !! 895 !! ** Purpose : Compute surface temperature and snow/ice conduction flux 896 !! to force sea ice / snow thermodynamics 897 !! in the case JULES coupler is emulated 898 !! 899 !! ** Method : compute surface energy balance assuming neglecting heat storage 900 !! following the 0-layer Semtner (1976) approach 901 !! 902 !! ** Outputs : - ptsu : sea-ice / snow surface temperature (K) 903 !! - qcn_ice : surface inner conduction flux (W/m2) 904 !! 905 !!--------------------------------------------------------------------- 906 !! 907 INTEGER, INTENT(in) :: k_monocat ! single-category option 908 909 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptsu ! sea ice / snow surface temperature 910 911 REAL(wp), DIMENSION(:,:) , INTENT(in) :: ptb ! sea ice base temperature 912 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness 913 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! sea ice thickness 914 915 INTEGER :: ji, jj, jl ! dummy loop indices 916 INTEGER :: iter ! 917 REAL(wp) :: zfac, zfac2, zfac3 ! dummy factors 918 REAL(wp) :: zkeff_h, ztsu ! 919 REAL(wp) :: zqc, zqnet ! 920 REAL(wp) :: zhe, zqa0 ! 921 922 INTEGER , PARAMETER :: nit = 10 ! number of iterations 923 REAL(wp), PARAMETER :: zepsilon = 0.1_wp ! characteristic thickness for enhanced conduction 924 925 REAL(wp), DIMENSION(:,:,:), POINTER :: zgfac ! enhanced conduction factor 926 927 !!--------------------------------------------------------------------- 928 929 IF( nn_timing == 1 ) CALL timing_start('blk_ice_qcn') 930 ! 931 CALL wrk_alloc( jpi,jpj,jpl, zgfac ) 932 933 ! -------------------------------------! 934 ! I Enhanced conduction factor ! 935 ! -------------------------------------! 936 ! 937 ! Emulates the enhancement of conduction by unresolved thin ice (k_monocat = 1/3) 938 ! Fichefet and Morales Maqueda, JGR 1997 939 ! 940 zgfac(:,:,:) = 1._wp 941 942 SELECT CASE ( k_monocat ) 943 944 CASE ( 1 , 3 ) 945 946 zfac = 1._wp / ( rn_cnd_s + rcdic ) 947 zfac2 = EXP(1._wp) * 0.5_wp * zepsilon 948 zfac3 = 2._wp / zepsilon 949 950 DO jl = 1, jpl 951 DO jj = 1 , jpj 952 DO ji = 1, jpi 953 ! ! Effective thickness 954 zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac 955 ! ! Enhanced conduction factor 956 IF( zhe >= zfac2 ) & 957 zgfac(ji,jj,jl) = MIN( 2._wp, ( 0.5_wp + 0.5 * LOG( zhe * zfac3 ) ) ) 958 END DO 959 END DO 960 END DO 961 962 END SELECT 963 964 ! -------------------------------------------------------------! 965 ! II Surface temperature and conduction flux ! 966 ! -------------------------------------------------------------! 967 968 zfac = rcdic * rn_cnd_s 969 ! ========================== ! 970 DO jl = 1, jpl ! Loop over ice categories ! 971 ! ! ========================== ! 972 DO jj = 1 , jpj 973 DO ji = 1, jpi 974 ! ! Effective conductivity of the snow-ice system divided by thickness 975 zkeff_h = zfac * zgfac(ji,jj,jl) / ( rcdic * phs(ji,jj,jl) + rn_cnd_s * phi(ji,jj,jl) ) 976 ! Store initial surface temperature 977 ztsu = ptsu(ji,jj,jl) 978 ! Net initial atmospheric heat flux 979 zqa0 = qsr_ice(ji,jj,jl) - qsr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl) 980 981 DO iter = 1, nit ! --- Iteration loop 982 983 ! ! Conduction heat flux through snow-ice system (>0 downwards) 984 zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) 985 ! ! Surface energy budget 986 zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc 987 ! ! Temperature update 988 ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) 989 990 END DO 991 992 ptsu(ji,jj,jl) = MIN( rt0, ztsu ) 993 994 qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 995 996 END DO 997 END DO 998 999 END DO 1000 1001 CALL wrk_dealloc( jpi,jpj,jpl, zgfac ) 1002 1003 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_qcn') 1004 1005 END SUBROUTINE blk_ice_qcn 856 1006 857 1007 #endif -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
r8592 r8726 1524 1524 1525 1525 1526 SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist )1526 SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist, phs, phi ) 1527 1527 !!---------------------------------------------------------------------- 1528 1528 !! *** ROUTINE sbc_cpl_ice_flx *** … … 1579 1579 REAL(wp), INTENT(in), DIMENSION(:,: ), OPTIONAL :: psst ! sea surface temperature [Celsius] 1580 1580 REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL :: pist ! ice surface temperature [Kelvin] 1581 ! 1582 INTEGER :: jl ! dummy loop index 1581 REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL :: phs ! snow depth [m] 1582 REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL :: phi ! ice thickness [m] 1583 ! 1584 INTEGER :: ji,jj,jl ! dummy loop index 1583 1585 REAL(wp), POINTER, DIMENSION(:,: ) :: zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw 1584 1586 REAL(wp), POINTER, DIMENSION(:,: ) :: zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice 1585 1587 REAL(wp), POINTER, DIMENSION(:,: ) :: zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice 1586 REAL(wp), POINTER, DIMENSION(:,:,:) :: zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice 1588 REAL(wp), POINTER, DIMENSION(:,:,:) :: zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zfrqsr_tr_i 1587 1589 !!---------------------------------------------------------------------- 1588 1590 ! … … 1592 1594 CALL wrk_alloc( jpi,jpj, zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice ) 1593 1595 CALL wrk_alloc( jpi,jpj, zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice ) 1594 CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )1596 CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zfrqsr_tr_i ) 1595 1597 1596 1598 IF( ln_mixcpl ) zmsk(:,:) = 1. - xcplmask(:,:,0) … … 1945 1947 END SELECT 1946 1948 1947 ! Surface transimission parameter io (Maykut Untersteiner , 1971 ; Ebert and Curry, 1993 ) 1948 ! Used for LIM3 1949 ! Coupled case: since cloud cover is not received from atmosphere 1950 ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 1951 fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 1952 fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 1949 ! --- Transmitted shortwave radiation (W/m2) --- ! 1950 1951 IF ( nice_jules == 0 ) THEN 1952 1953 zfrqsr_tr_i(:,:,:) = 0._wp ! surface transmission parameter 1954 1955 ! former coding was 1956 ! Coupled case: since cloud cover is not received from atmosphere 1957 ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 1958 ! fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 1959 ! fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 1960 1961 ! to retrieve that coding, we needed to access h_i & h_s from here 1962 ! we could even retrieve cloud fraction from the coupler 1963 1964 DO jl = 1, jpl 1965 DO jj = 1 , jpj 1966 DO ji = 1, jpi 1967 1968 !--- surface transmission parameter (Grenfell Maykut 77) --- ! 1969 zfrqsr_tr_i(ji,jj,jl) = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice 1970 1971 ! --- influence of snow and thin ice --- ! 1972 IF ( phs(ji,jj,jl) >= 0.0_wp ) zfrqsr_tr_i(ji,jj,jl) = 0._wp ! snow fully opaque 1973 IF ( phi(ji,jj,jl) <= 0.1_wp ) zfrqsr_tr_i(ji,jj,jl) = 1._wp ! thin ice transmits all solar radiation 1974 END DO 1975 END DO 1976 END DO 1977 1978 qsr_ice_tr(:,:,:) = zfrqsr_tr_i(:,:,:) * qsr_ice(:,:,:) ! transmitted solar radiation 1979 1980 ENDIF 1981 1982 IF ( nice_jules == 2 ) THEN 1983 1984 ! here we must receive the qsr_ice_tr array from the coupler 1985 ! for now just assume zero 1986 1987 qsr_ice_tr(:,:,:) = 0.0_wp 1988 1989 ENDIF 1990 1991 1953 1992 1954 1993 CALL wrk_dealloc( jpi,jpj, zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw ) 1955 1994 CALL wrk_dealloc( jpi,jpj, zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice ) 1956 1995 CALL wrk_dealloc( jpi,jpj, zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice ) 1957 CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )1996 CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zfrqsr_tr_i ) 1958 1997 ! 1959 1998 IF( nn_timing == 1 ) CALL timing_stop('sbc_cpl_ice_flx')
Note: See TracChangeset
for help on using the changeset viewer.
