SUBROUTINE ice_rad(nlay_s,nlay_i,kideb,kiut) !=============================================================================! ! ! ice_rad : ! ! Transmission / absorption of radiation ! through the snow-ice system ! ! (c) v1.0 Martouf, UCL-ASTR, June 2007. Nadal Federer 1 set partout ! v2.0 Oct 2011: clarification of snow and the SSL. ! Les diables n'y sont pas ! ! Practically, downwelling radiation decreases exponentially ! through the snow-ice surface, except near the surface ! The uppermost portion of the snow is highly scattering (Perovich, ! J. Glaciol., 2007). For sea ice, if melt has occured, a ! a highly scattering layer is formed near the surface (Light et ! al., JGR 2008). This layer is referred to as "Surface Scattering Layer" (SSL) ! ! We assume that the surface layer, both in snow and ice, absorbs ! all near infrared radiation, so that the remaining radiation is ! only in the visible part of the SW spectrum. ! ! Here, snow and sea ice are represented as a highly-scattering ! surface layer (thickness h_not_s / h_not_i) and a deeper layer ! ! The surface layer has radiative properties zrad_kappa_i_su_x, ! zrad_kappa_s_su_x, ! Depth has radiative properties zrad_kappa_s_de_x, ! zrad_kappa_i_de_x ! ! Algae and detritus also absorb radiation ! ! If sea ice is snow covered the SSL has a depth equal to zero ! ! 2 schemes for the SSL are used ! ... more ... ! ! 2 discretizations for radiative transfer are coded ! ... more ... ! !=============================================================================! ! USE lib_fortran INCLUDE 'type.com' INCLUDE 'para.com' INCLUDE 'const.com' INCLUDE 'ice.com' INCLUDE 'thermo.com' INCLUDE 'bio.com' !-----------------------------------------------------------------------------! ! Local variables REAL(8), DIMENSION(maxnlay) :: & zkappa_alg , !: extinction coefficient for algae & zkappa_det !: extinction coefficient for detritus LOGICAL :: ln_write zeps = 1.0e-20 ln_write = .FALSE. !==============================================================================! WRITE(numout,*) WRITE(numout,*) ' ** ice_rad : ' WRITE(numout,*) ' ~~~~~~~~~~~~ ' WRITE(numout,*) ' c_rad_scheme : ', c_rad_scheme WRITE(numout,*) ' c_rad_discr : ', c_rad_discr WRITE(numout,*) ' h_not_s: ', h_not_s WRITE(numout,*) ' h_not_i: ', h_not_i WRITE(numout,*) DO ji = kideb, kiut ! !------------------------------------------------------------------------------! ! 1) Surface transmission parameter ! !------------------------------------------------------------------------------! ! IF ( ln_write ) THEN WRITE(numout,*) WRITE(numout,*) ' Surface transmission parameter ' WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' ENDIF !-------------------------------- ! Snow and surface melt switches !-------------------------------- ! Is there snow or not ? isnow = INT( 1.0 - MAX( 0.0 , SIGN( 1.0d0 , - ht_s_b(ji) ) ) ) IF ( ln_write ) THEN WRITE(numout,*) ' Test isnow ', ji WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji) WRITE(numout,*) ' SIGN : ', SIGN( 1.0d0 , - ht_s_b(ji) ) WRITE(numout,*) ' MAX : ', MAX( 0.0 , & SIGN( 1.0d0 , - ht_s_b(ji) ) ) WRITE(numout,*) ' ARG : ', 1.0 - MAX( 0.0 , SIGN( 1.0d0 , & - ht_s_b(ji) ) ) ENDIF ! Melting or not imelt_su = INT( MAX( 0.0 , SIGN( 1.0d0 , t_su_b(ji)-tpw) ) ) ! 1 if surface melt imelt_sn = INT( MAX( 0.0 , SIGN( 1.0d0 , t_s_b(ji,1)-tpw) ) ) ! 1 if snow melts imelt = 0 IF ( ( imelt_su .EQ. 1 ) .OR. ( imelt_sn .EQ. 1 ) ) imelt = 1 !----------------------------------- ! Surface transmissivity i_o (inot) !----------------------------------- IF ( c_rad_scheme .EQ. 'INOT' ) THEN z_inot_vis_snow = rad_inot_s_dry * ( 1.0 - imelt ) + & rad_inot_s_wet * imelt z_inot_vis_ice = rad_inot_i_dry * ( 1.0 - imelt ) + & rad_inot_i_wet * imelt z_inot_vis = isnow * z_inot_vis_snow + & ( 1.0 - isnow ) * z_inot_vis_ice ! z_inot_vis = isnow * rad_inot_s + ( 1.0 - isnow ) * rad_inot_i z_inot = z_inot_vis * fpar_fsw ENDIF IF ( c_rad_scheme .EQ. 'EXTC' ) THEN z_inot = 1.0 ENDIF ab(ji) = 1.0 - z_inot ! used in other routines !------------------------------------------- ! Solar radiation transmitted below the SSL !------------------------------------------- ftrice = fsolgb(ji) * ( 1.0 - ab(ji) ) IF ( ln_write ) THEN WRITE(numout,*) ' isnow : ', isnow WRITE(numout,*) ' imelt : ', imelt WRITE(numout,*) ' z_inot: ', z_inot WRITE(numout,*) ' ab : ', ab(ji) WRITE(numout,*) ' ftrice : ', ftrice ENDIF ! !------------------------------------------------------------------------------! ! 2) Extinction coefficients !------------------------------------------------------------------------------! ! IF ( ln_write ) THEN WRITE(numout,*) WRITE(numout,*) ' Extinction coefficients ' WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~' WRITE(numout,*) ENDIF zmurad = 0.656 ! angle factor zchla = 0.0 ! temporary value of chla in mg/m-3 ! Chlorophyll a interpolation IF ( c_bio_model .EQ. 'KRILL' ) THEN CALL ice_bio_interp_bio2phy(kideb,kiut,nlay_i,.FALSE.) ELSE DO layer = 1, nlay_bio chla_i(layer) = 0.0 END DO ENDIF IF ( ln_write ) THEN WRITE(numout,*) ' chla_i : ', & ( chla_i(layer), layer = 1, nlay_i ) WRITE(numout,*) ENDIF ! chlorophyll and detritus extinction coefficients DO layer = 1, nlay_i zchla = chla_i(layer) zkappa_alg(layer) = zchla * astar_alg / zmurad zkappa_det(layer) = zkappa_alg(layer) * fdet_alg END DO ! snow and ice extinction coefficients IF ( c_rad_scheme .EQ. 'INOT' ) THEN zrad_kappa_s_su = 0.0 zrad_kappa_i_su = 0.0 ENDIF IF ( c_rad_scheme .EQ. 'EXTC' ) THEN zrad_kappa_s_su = ( 1 - imelt ) * rad_kappa_s_su_d + & imelt * rad_kappa_s_su_m zrad_kappa_i_su = ( 1 - imelt ) * rad_kappa_i_su_d + & imelt * rad_kappa_i_su_m ENDIF ! deep snow extinction coefficient zrad_kappa_s_de = ( 1 - imelt ) * rad_kappa_s_de_d + & imelt * rad_kappa_s_de_m ! deep ice extinction coefficient zrad_kappa_i_de = ( 1 - imelt ) * rad_kappa_i_de_d + & imelt * rad_kappa_i_de_m IF ( ln_write ) THEN WRITE(numout,*) ' zrad_kappa_s_su : ', zrad_kappa_s_su WRITE(numout,*) ' zrad_kappa_s_de : ', zrad_kappa_s_de WRITE(numout,*) ' zrad_kappa_i_de : ', zrad_kappa_i_de WRITE(numout,*) ' zrad_kappa_i_su : ', zrad_kappa_i_su WRITE(numout,*) ' zkappa_alg : ', ( zkappa_alg(layer), & layer = 1, nlay_i ) WRITE(numout,*) ' zkappa_det : ', ( zkappa_det(layer), & layer = 1, nlay_i ) ENDIF ! !------------------------------------------------------------------------------! ! 3) Radiation transmitted through / absorbed by snow ! !------------------------------------------------------------------------------! ! IF ( ln_write ) THEN WRITE(numout,*) WRITE(numout,*) ' Radiation absorption in snow ' WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' WRITE(numout,*) ENDIF radtr_s(0) = ftrice zdh1 = MIN( h_not_s, ht_s_b(ji) ) zdh2 = MIN( MAX( ht_s_b(ji) - h_not_s, 0. ), ht_s_b(ji) ) IF ( ln_write ) THEN WRITE(numout,*) ' ftrice : ', ftrice WRITE(numout,*) ' zdh1 : ', zdh1 WRITE(numout,*) ' zdh2 : ', zdh2 ENDIF ! scheme based on transmission ! other scheme based on transmission, more exact, the other one seems ! not to work for high transmission coefficient and small SSL IF ( c_rad_discr .EQ. 'TRA' ) THEN zi1 = radtr_s(0) * EXP( -zrad_kappa_s_su * zdh1 ) radtr_s(1) = zi1 * EXP( -zrad_kappa_s_de * zdh2 ) radab_s(1) = radtr_s(0) - radtr_s(1) IF ( ln_write) WRITE(numout,*) ' zi1 : ', zi1 ENDIF ! scheme based on absorption IF ( c_rad_discr .EQ. 'ABS' ) THEN zrada1 = radtr_s(0) * zdh1 * zrad_kappa_s_su * & EXP( - zrad_kappa_s_su * zdh1 ) zi1 = radtr_s(0) - zrada1 zrada2 = zi1 * zdh2 * zrad_kappa_s_de * & EXP( - zrad_kappa_s_de * zdh2 ) IF ( ln_write ) THEN WRITE(numout,*) ' zrada1 : ', zrada1 WRITE(numout,*) ' zi1 : ', zi1 WRITE(numout,*) ' zrada2 : ', zrada2 ENDIF radab_s(1) = zrada1 + zrada2 radtr_s(1) = radtr_s(0) - radab_s(1) ENDIF IF ( ln_write ) THEN WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji) WRITE(numout,*) ' radtr_s : ', & ( radtr_s(layer) , layer = 0, nlay_s ) WRITE(numout,*) ' radab_s : ', & ( radab_s(layer) , layer = 1, nlay_s ) ENDIF ! !------------------------------------------------------------------------------! ! 4) Radiation transmitted through / absorbed by ice ! !------------------------------------------------------------------------------! ! IF ( ln_write ) THEN WRITE(numout,*) WRITE(numout,*) ' Radiation absorption in ice ' WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' WRITE(numout,*) ENDIF ! transmitted at the upper ice layer radtr_i(0) = isnow * radtr_s(nlay_s) + ( 1. - isnow ) * ftrice ! thickness of surface SSL in sea ice is zero if there is no snow zh_ssl = isnow * 0. + ( 1. - isnow ) * h_not_i IF ( ln_write ) THEN WRITE(numout,*) ' zh_ssl : ', zh_ssl ENDIF ! radiation absorbed physically and biologically in each layer zz0 = 0. zz1 = 0. DO layer = 1, nlay_i zz1 = zz1 + deltaz_i_phy(layer) IF ( ln_write ) THEN WRITE(numout,*) ' zz0 : ', zz0 WRITE(numout,*) ' zz1 : ', zz1 ENDIF zdh1 = MIN ( MAX ( zh_ssl - zz0 , 0. ) , & deltaz_i_phy(layer) ) ! part of the ice layer belonging to the SSL zdh2 = MIN ( MAX ( zz1 - zh_ssl , 0. ) , & deltaz_i_phy(layer) ) ! part of the ice layer belonging to the deep ice zz0 = zz1 IF ( ln_write ) THEN WRITE(numout,*) ' zh_ssl : ', zh_ssl WRITE(numout,*) ' zz1-zh_ssl ',zz1 - zh_ssl WRITE(numout,*) ' deltaz_i_phy : ', deltaz_i_phy(layer) WRITE(numout,*) ' ' WRITE(numout,*) ' layer : ', layer WRITE(numout,*) ' zdh1 : ', zdh1 WRITE(numout,*) ' zdh2 : ', zdh2 WRITE(numout,*) ' deltaz_i_phy : ', deltaz_i_phy(layer) WRITE(numout,*) ' ' ENDIF ! Extinction coefficient in the SSL zkappa1_phy = zrad_kappa_i_su + zkappa_det(layer) ! physical extinction in SSL zkappa1 = zkappa1_phy + zkappa_alg(layer) ! total extinction ! Extinction coefficient deeper in the ice zkappa2_phy = zrad_kappa_i_de + zkappa_det(layer) ! physical extinction at depth zkappa2 = zkappa2_phy + zkappa_alg(layer) ! total extinction at depth IF ( ln_write ) THEN WRITE(numout,*) ' zkappa1_phy : ', zkappa1_phy WRITE(numout,*) ' zkappa1 : ', zkappa1 WRITE(numout,*) ' zkappa2_phy : ', zkappa2_phy WRITE(numout,*) ' zkappa2 : ', zkappa2 ENDIF !------------------------------ ! Scheme based on transmission !------------------------------ IF ( c_rad_discr .EQ. 'TRA' ) THEN zi1 = radtr_i(layer-1) * EXP ( -zkappa1*zdh1 ) WRITE(numout,*) ' Discretization : ', 'TRA' WRITE(numout,*) ' zi1 : ', zi1 WRITE(numout,*) ! transmitted radiation radtr_i(layer) = zi1 * EXP( -zkappa2 * zdh2 ) ! absorbed radiation radab_i(layer) = radtr_i(layer-1) - radtr_i(layer) ! physically absorbed radab_phy_i(layer) = radab_i(layer) * deltaz_i_phy(layer) & * ( zkappa1_phy * zdh1 / MAX( zkappa1, zeps) & + zkappa2_phy * zdh2 / MAX( zkappa2, zeps) ) IF ( ln_write ) THEN WRITE(numout,*) ' deltaz_i_phy : ', deltaz_i_phy(layer) WRITE(numout,*) ' zkappa2 : ', zkappa2 WRITE(numout,*) ' zkappa2_phy : ', zkappa2_phy WRITE(numout,*) ' zdh2 : ', zdh2 WRITE(numout,*) ' radab_i : ', radab_i(layer) WRITE(numout,*) ' zkappa1 : ', zkappa1 WRITE(numout,*) ' zkappa1_phy : ', zkappa1_phy WRITE(numout,*) ' zdh1 : ', zdh1 ENDIF ! biologically absorbed radab_alg_i(layer) = zkappa_alg(layer) * radab_i(layer) & / deltaz_i_phy(layer) & * ( zdh1 / MAX ( zkappa1, zeps ) & + zdh2 / MAX ( zkappa2, zeps ) ) ENDIF !---------------------------- ! Scheme based on absorption !---------------------------- IF ( c_rad_discr .EQ. 'ABS' ) THEN zdummy = radtr_i(layer-1) * zdh1 * & EXP( - zkappa1 * zdh1 ) zrada1_phy = zkappa1_phy * zdummy ! physically absorbed zrada1_alg = zkappa_alg(layer) * zdummy ! biologically absorbed zrada1 = zrada1_phy + zrada1_alg ! total absorption zi1 = radtr_i(layer-1) - zrada1 ! radiance available below the SSL IF ( ln_write) THEN WRITE(numout,*) ' Discretization : ', 'ABS' WRITE(numout,*) ' zi1 : ', zi1 WRITE(numout,*) WRITE(numout,*) ' zdummy : ', zdummy WRITE(numout,*) ' zrada1_phy : ', zrada1_phy WRITE(numout,*) ' zrada1_alg : ', zrada1_alg WRITE(numout,*) ' zrada1 : ', zrada1 ENDIF zdummy = zi1 * zdh2 * EXP ( - zkappa2 * zdh2 ) zrada2_phy = zkappa2_phy * zdummy ! physically absorbed zrada2_alg = zkappa_alg(layer) * zdummy ! biologically absorbed zrada2 = zrada2_phy + zrada2_alg ! total absorption IF ( ln_write) THEN WRITE(numout,*) ' zdummy : ', zdummy WRITE(numout,*) ' zrada2_phy : ', zrada2_phy WRITE(numout,*) ' zrada2_alg : ', zrada2_alg WRITE(numout,*) ' zrada2 : ', zrada2 ENDIF ! Cumulated absorption radab_phy_i(layer) = zrada1_phy + zrada2_phy radab_alg_i(layer) = zrada1_alb + zrada2_alg radab_i(layer) = radab_phy_i(layer) + radab_alg_i(layer) ! Transmission below the layer radtr_i(layer) = radtr_i(layer-1) - radab_phy_i(layer) & - radab_alg_i(layer) ENDIF IF ( ln_write ) THEN WRITE(numout,*) ' radab_phy_i : ', radab_phy_i(layer) WRITE(numout,*) ' radab_alg_i : ', radab_alg_i(layer) WRITE(numout,*) ' radab_i : ', radab_i(layer) WRITE(numout,*) ' sum of both : ', radab_phy_i(layer) & + radab_alg_i(layer) WRITE(numout,*) ' radtr_i : ', radtr_i(layer) ENDIF END DO ! radiation sent to the ocean ftroce = radtr_i(nlay_i) ! !------------------------------------------------------------------------------! ! 5) PAR (phy and bio grids) !------------------------------------------------------------------------------! ! IF ( ln_write ) THEN WRITE(numout,*) WRITE(numout,*) ' PAR ' WRITE(numout,*) ' ~~~~' WRITE(numout,*) ENDIF ! PAR on phy grid (top of each layer) ! this choice was made to minimize sensitivity to grid type at comparable resolution) IF ( ln_write ) WRITE(numout,*) ' * physical grid ... ' WRITE(numout,*) ' qpar_fsw : ', qpar_fsw WRITE(numout,*) ' fpar_fsw : ', fpar_fsw DO layer = 1, nlay_i par(layer) = qpar_fsw / fpar_fsw * radtr_i(layer-1) END DO ! ! PAR on bio grid (top of each layer) IF ( ln_write ) WRITE(numout,*) ' * biological grid ... ' DO layer_bio = 1, nlay_bio ! identify in which layer the upper boundary of the bio layer is DO layer_phy = 1, nlay_bio zzb = zb_i_bio(layer_bio-1) zzf1 = zb_i_phy(layer_phy-1) zzf2 = zb_i_phy(layer_phy) IF ( ( zzb .GE. zzf1 ) .AND. ( zzb .LT. zzf2 ) ) & zindex = layer_phy END DO zm = ( radtr_i(zindex) - radtr_i(zindex-1) ) / & ( zb_i_phy(zindex) - zb_i_phy(zindex-1) ) zdh = zb_i_bio(layer_bio-1) - zb_i_phy(zindex-1) zdrad = zm * zdh par_bio(layer_bio) = radtr_i(zindex-1) + zdrad ! at the upper interface of the layer ! IF ( ln_write ) THEN ! WRITE(numout,*) ! WRITE(numout,*) ' layer : ', layer_bio ! WRITE(numout,*) ' zindex : ', zindex ! WRITE(numout,*) ' zb_i_bio : ', zb_i_bio(layer_bio-1) ! WRITE(numout,*) ' zb_i_phy : ', zb_i_phy(zindex-1), ! & zb_i_phy(zindex) ! WRITE(numout,*) ' radtr_i : ', radtr_i(zindex-1), ! & radtr_i(zindex) ! WRITE(numout,*) ' In W/m2 ... ' ! WRITE(numout,*) ' par_bio : ', par_bio(layer_bio) ! ENDIF par_bio(layer_bio) = par_bio(layer_bio) * qpar_fsw / fpar_fsw ! IF ( ln_write ) THEN ! WRITE(numout,*) ' In microE/m2/s ... ' ! WRITE(numout,*) ' par_bio : ', par_bio(layer_bio) ! ENDIF END DO ! !------------------------------------------------------------------------------! ! 6) Conservation check !------------------------------------------------------------------------------! ! sumrad = radab_s(1) DO layer = 1, nlay_i sumrad = sumrad + radab_i(layer) END DO WRITE(numout,*) ' Conservation check ' WRITE(numout,*) ' ftrice - ftroce : ', ftrice-ftroce WRITE(numout,*) ' sumrad : ', sumrad ! Final write WRITE(numout,*) WRITE(numout,*) ' i0 : ', 1.0-ab(ji) WRITE(numout,*) ' ftrice : ', ftrice WRITE(numout,*) ' ftroce : ', ftroce WRITE(numout,*) ' radtr_i : ', & ( radtr_i(layer) , layer = 0, nlay_i ) WRITE(numout,*) ' radab_phy_i : ', & ( radab_phy_i(layer) , layer = 1, nlay_i ) WRITE(numout,*) ' radab_alg_i : ', & ( radab_alg_i(layer) , layer = 1, nlay_i ) WRITE(numout,*) ' par : ', ( par(layer), layer = 1, nlay_i ) WRITE(numout,*) ' par_bio : ', ( par_bio(layer), & layer = 1, nlay_bio ) WRITE(numout,*) WRITE(numout,*) ' End of ice_rad ' WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' END DO !ji !==============================================================================! ! end of the subroutine END SUBROUTINE