From 16572fdfacaec2449301c3fba89e0fa6a57d0d35 Mon Sep 17 00:00:00 2001 From: Nikolay Koldunov Date: Sun, 2 Mar 2025 23:11:41 +0100 Subject: [PATCH 1/8] bring changes from Frank to latest version --- src/ice_thermo_oce.F90 | 216 +++++++++++++++++++++++++++++++++-------- 1 file changed, 178 insertions(+), 38 deletions(-) diff --git a/src/ice_thermo_oce.F90 b/src/ice_thermo_oce.F90 index afc2fd501..89d942f76 100755 --- a/src/ice_thermo_oce.F90 +++ b/src/ice_thermo_oce.F90 @@ -26,18 +26,18 @@ module ice_therm_interface interface subroutine therm_ice(ithermp, h,hsn,A,fsh,flo,Ta,qa,rain,snow,runo,rsss, & ug,ustar,T_oc,S_oc,H_ML,t,ice_dt,ch,ce,ch_i,ce_i,evap_in,fw,ehf,evap, & - rsf, dhgrowth, dhsngrowth, iflice, hflatow, hfsenow, hflwrdout,lid_clo,subli) + rsf, dhgrowth, dhsngrowth, iflice, hflatow, hfsenow, hflwrdout,lid_clo,geolon, geolat, subli) USE MOD_ICE type(t_ice_thermo), intent(in), target :: ithermp real(kind=WP) h, hsn, A, fsh, flo, Ta, qa, rain, snow, runo, rsss, evap_in, & - ug, ustar, T_oc, S_oc, H_ML, t, ice_dt, ch, ce, ch_i, ce_i, fw, ehf, & + ug, ustar, T_oc, S_oc, H_ML, t, ice_dt, ch, ce, ch_i, ce_i, evap_in, fw, ehf, & dhgrowth, dhsngrowth, ahf, prec, subli, subli_i, rsf, & rhow, show, rhice, shice, sh, thick, thact, lat, & rh, rA, qhst, sn, hsntmp, o2ihf, evap, iflice, hflatow, & - hfsenow, hflwrdout, lid_clo - end subroutine + hfsenow, hflwrdout, lid_clo, geolon, geolat + end subroutine therm_ice end interface -end module +end module ice_therm_interface module ice_budget_interfaces interface @@ -45,22 +45,22 @@ subroutine budget(ithermp, hice, hsn, t, ta, qa, fsh, flo, ug, S_oc, ch_i, ce_i, USE MOD_ICE type(t_ice_thermo), intent(in), target :: ithermp real(kind=WP) hice, hsn, t, ta, qa, fsh, flo, ug, S_oc, ch_i, ce_i, fh, subli - end subroutine + end subroutine budget - subroutine obudget(ithermp, qa, fsh, flo, t, ug, ta, ch, ce, fh, evap, hflatow, hfsenow, hflwrdout) + subroutine obudget(ithermp, qa, fsh, flo, t, ug, ta, ch, ce, geolon, geolat, fh, evap, hflatow, hfsenow, hflwrdout) USE MOD_ICE type(t_ice_thermo), intent(in), target :: ithermp - real(kind=WP) qa, t, ta, fsh, flo, ug, ch, ce, fh, evap, hfsenow, & - hfradow, hflatow, hftotow, hflwrdout - end subroutine + real(kind=WP) qa, t, ta, fsh, flo, ug, ch, ce, geolon, geolat, fh, evap, hfsenow, & + hfsensow, hfradow, hflatow, hftotow, hflwrdout + end subroutine obudget subroutine flooding(ithermp, h, hsn) USE MOD_ICE type(t_ice_thermo), intent(in), target :: ithermp real(kind=WP) h, hsn - end subroutine + end subroutine flooding end interface -end module +end module ice_budget_interfaces ! ! !_______________________________________________________________________________ @@ -94,26 +94,26 @@ subroutine cut_off(ice, partit, mesh) #endif /* (__oifs) */ !___________________________________________________________________________ - ! lower cutoff: a_ice + ! upper cutoff: a_ice !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(n) DO n=1, myDim_nod2D+eDim_nod2D if (a_ice(n) > 1.0_WP) a_ice(n)=1.0_WP - ! upper cutoff: a_ice + ! lower cutoff: a_ice if (a_ice(n) < .1e-8_WP) then - a_ice(n) =0.0_WP + a_ice(n)=0.0_WP +#if defined (__oifs) || defined (__ifsinterface) m_ice(n) =0.0_WP m_snow(n) =0.0_WP -#if defined (__oifs) || defined (__ifsinterface) ice_temp(n)=273.15_WP #endif /* (__oifs) */ end if !___________________________________________________________________________ ! lower cutoff: m_ice if (m_ice(n) < .1e-8_WP) then - m_ice(n) =0.0_WP + m_ice(n)=0.0_WP +#if defined (__oifs) || defined (__ifsinterface) m_snow(n) =0.0_WP a_ice(n) =0.0_WP -#if defined (__oifs) || defined (__ifsinterface) ice_temp(n)=273.15_WP #endif /* (__oifs) */ end if @@ -169,8 +169,9 @@ subroutine thermodynamics(ice, partit, mesh) real(kind=WP) :: h,hsn,A,fsh,flo,Ta,qa,rain,snow,runo,rsss,rsf,evap_in real(kind=WP) :: ug,ustar,T_oc,S_oc,h_ml,t,ch,ce,ch_i,ce_i,fw,ehf,evap real(kind=WP) :: ithdgr, ithdgrsn, iflice, hflatow, hfsenow, hflwrdout, subli - real(kind=WP) :: lid_clo + real(kind=WP) :: lid_clo, o2ihf real(kind=WP) :: lat + real(kind=WP) :: geolon, geolat !_____________________________________________________________________________ ! pointer on necessary derived types @@ -183,6 +184,7 @@ subroutine thermodynamics(ice, partit, mesh) real(kind=WP), dimension(:) , pointer :: thdgr, thdgrsn, thdgr_old, t_skin, ustar_aux real(kind=WP), dimension(:) , pointer :: S_oc_array, T_oc_array, u_w, v_w real(kind=WP), dimension(:) , pointer :: fresh_wa_flux, net_heat_flux + real(kind=WP), external :: TFrez ! Sea water freeze temperature myDim_nod2d => partit%myDim_nod2D eDim_nod2D => partit%eDim_nod2D ulevels_nod2D (1 :myDim_nod2D+eDim_nod2D) => mesh%ulevels_nod2D(:) @@ -246,7 +248,15 @@ subroutine thermodynamics(ice, partit, mesh) fsh = shortwave(i) flo = longwave(i) Ta = Tair(i) - qa = shum(i) + if (l_tdew) then + call + qa = 0. + elseif (l_humi) then + qa = shum(i) + else + stop 'neither specific humidity nor dew-point temperature is defined!' + endif + if (.not. l_snow) then if (Ta>=0.0_WP) then rain=prec_rain(i) @@ -281,12 +291,14 @@ subroutine thermodynamics(ice, partit, mesh) else lid_clo=0.5_WP endif + geolon = geo_coord_nod2D(1, i) + geolat = geo_coord_nod2D(2, i) !_______________________________________________________________________ ! do ice thermodynamics call therm_ice(ice%thermo,h,hsn,A,fsh,flo,Ta,qa,rain,snow,runo,rsss, & ug,ustar,T_oc,S_oc,h_ml,t,ice%ice_dt,ch,ce,ch_i,ce_i,evap_in,fw,ehf,evap, & - rsf, ithdgr, ithdgrsn, iflice, hflatow, hfsenow, hflwrdout,lid_clo,subli) + rsf, ithdgr, ithdgrsn, iflice, hflatow, hfsenow, hflwrdout,lid_clo,geolon, geolat, subli) !_______________________________________________________________________ ! write ice thermodyn. results into arrays @@ -306,6 +318,9 @@ subroutine thermodynamics(ice, partit, mesh) net_heat_flux(i) = ehf !positive down evaporation(i) = evap !negative up ice_sublimation(i)= subli + ice_ocean_hflx(i) = o2ihf + ice_dTfrez(i) = T_oc-TFrez(S_oc) + ice_ustar(i) = ustar thdgr(i) = ithdgr thdgrsn(i) = ithdgrsn @@ -333,7 +348,7 @@ end subroutine thermodynamics subroutine therm_ice(ithermp, h, hsn, A, fsh, flo, Ta, qa, rain, snow, runo, rsss, & ug, ustar, T_oc, S_oc, H_ML, t, ice_dt, ch, ce, ch_i, ce_i, & evap_in, fw, ehf, evap, rsf, dhgrowth, dhsngrowth, iflice, & - hflatow, hfsenow, hflwrdout, lid_clo, subli) + hflatow, hfsenow, hflwrdout, lid_clo, geolon, geolat, subli) ! Ice Thermodynamic growth model ! ! Input parameters: @@ -359,6 +374,8 @@ subroutine therm_ice(ithermp, h, hsn, A, fsh, flo, Ta, qa, rain, snow, runo, rss ! ch_i - transfer coefficient for sensible heat (for ice) ! ce_i - transfer coefficient for evaporation (for ice) ! lid_clo - lid closing parameter + ! geolon - geographical longitude + ! geolat - geographical latitude ! Output parameters: !------------------- ! h - ice mass @@ -367,6 +384,8 @@ subroutine therm_ice(ithermp, h, hsn, A, fsh, flo, Ta, qa, rain, snow, runo, rss ! t - temperature of snow/ice top surface ! fw - freshwater flux due to ice melting [m water/ice_dt] ! ehf - net heat flux at the ocean surface [W/m2] !RTnew + ! subli - sublimatione over ice + ! o2ihf - ocean to ice heat flux [W/m2] USE MOD_ICE use g_forcing_param, only: use_virt_salt @@ -378,14 +397,20 @@ subroutine therm_ice(ithermp, h, hsn, A, fsh, flo, Ta, qa, rain, snow, runo, rss real(kind=WP) h,hsn,A,fsh,flo,Ta,qa,rain,snow,runo,rsss,evap_in real(kind=WP) ug,ustar,T_oc,S_oc,H_ML,t,ice_dt,ch,ce,ch_i,ce_i,fw,ehf real(kind=WP) dhgrowth,dhsngrowth,ahf,prec,subli,subli_i,rsf - real(kind=WP) rhow,show,rhice,shice,sh,thick,thact,lat + real(kind=WP) rhow,show,rhice,shice,sh,snthick,thick,thact,lat real(kind=WP) rh,rA,qhst,sn,hsntmp,o2ihf,evap real(kind=WP) iflice,hflatow,hfsenow,hflwrdout real(kind=WP), external :: TFrez ! Sea water freeze temperature. - real(kind=WP) lid_clo + real(kind=WP) lid_clo, geolon, geolat !___________________________________________________________________________ - real(kind=WP), pointer :: hmin, Sice, Armin, cc, cl, con, consn, rhosno, rhoice, inv_rhowat, inv_rhosno - integer , pointer :: iclasses + logical , pointer :: snowdist, new_iclasses + integer , pointer :: iclasses, open_water_albedo + real(kind=WP), pointer :: hmin, Sice, Armin, cc, cl, con, consn, rhosno, rhoice, inv_rhowat, inv_rhosno, c_melt, h_cutoff + real(kind=WP), pointer, dimension (:) :: hpdf + snowdist => ithermp%snowdist + new_iclasses => ithermp%new_iclasses + iclasses => ithermp%iclasses + open_water_albedo => ithermp%open_water_albedo hmin => ithermp%hmin Armin => ithermp%Armin Sice => ithermp%Sice @@ -398,6 +423,9 @@ subroutine therm_ice(ithermp, h, hsn, A, fsh, flo, Ta, qa, rain, snow, runo, rss rhoice => ithermp%rhoice inv_rhowat => ithermp%inv_rhowat inv_rhosno => ithermp%inv_rhosno + c_melt => ithermp%c_melt + h_cutoff => ithermp%h_cutoff + hpdf => ithermp%hpdf !___________________________________________________________________________ ! Store ice thickness at start of growth routine @@ -408,13 +436,14 @@ subroutine therm_ice(ithermp, h, hsn, A, fsh, flo, Ta, qa, rain, snow, runo, rss ! of conductivities of ice and snow, according to 0-layer approach ! of Semtner (1976). ! thickness at the ice covered part - thick=hsn*(con/consn)/max(A,Armin) ! Effective snow thickness - thick=thick+h/max(A,Armin) ! Effective total snow-ice thickness + snthick=hsn*(con/consn)/max(A,Armin) ! Effective snow thickness + thick=h/max(A,Armin) ! Effective ice thickness + if (snowdist) thick=snthick+thick ! Effective ice and ice thickness ! Growth rate for ice in open ocean rhow=0.0_WP evap=0.0_WP - call obudget(ithermp, qa,fsh,flo,T_oc,ug,ta,ch,ce,rhow,evap,hflatow,hfsenow,hflwrdout) + call obudget(ithermp, qa,fsh,flo,T_oc,ug,ta,ch,ce,geolon, geolat, rhow,evap,hflatow,hfsenow,hflwrdout) hflatow=hflatow*(1.0_WP-A) hfsenow=hfsenow*(1.0_WP-A) hflwrdout=hflwrdout*(1.0_WP-A) @@ -431,13 +460,22 @@ subroutine therm_ice(ithermp, h, hsn, A, fsh, flo, Ta, qa, rain, snow, runo, rss if (thick.gt.hmin) then do k=1,iclasses thact = real((2*k-1),WP)*thick/real(iclasses,WP) ! Thicknesses of actual ice class + if(new_iclasses) thact=h_cutoff/2.*thact ! h_cutoff is variable (originally hcutoff was 2*h => factor 2 singles out) + if(.not. snowdist) thact=thact+snthick ! if snowdist=.true. snow depth is the same on every ice class call budget(ithermp, thact, hsn,t,Ta,qa,fsh,flo,ug,S_oc,ch_i,ce_i,shice,subli_i) !Thick ice K-class growth rate - rhice=rhice+shice ! Add to average heat flux - subli=subli+subli_i + if(new_iclasses) then + rhice=rhice+shice*hpdf(k) + subli=subli+subli_i*hpdf(k) + else + rhice=rhice+shice + subli=subli+subli_i + end if end do - rhice=rhice/real(iclasses,WP) ! Add to average heat flux - subli=subli/real(iclasses,WP) + if(.not. new_iclasses) then + rhice=rhice/real(iclasses,WP) ! Add to average heat flux + subli=subli/real(iclasses,WP) + end if end if ! Convert growth rates [m ice/sec] into growth per time step DT. @@ -546,15 +584,13 @@ subroutine therm_ice(ithermp, h, hsn, A, fsh, flo, Ta, qa, rain, snow, runo, rss rsf= -dhgrowth*rhoice*inv_rhowat*Sice else fw= prec+evap - dhgrowth*rhoice*inv_rhowat*(rsss-Sice)/rsss - dhsngrowth*rhosno*inv_rhowat - ! initialize rsf with zero since it is not used for linfs only in zstar - rsf = 0.0_WP end if ! Changes in compactnesses (equation 16 of Hibler 1979) rh=-min(h,-rh) ! Make sure we do not try to melt more ice than is available rA= rhow - o2ihf*ice_dt/cl !Qiang: it was -(T_oc-TFrez(S_oc))*H_ML*cc/cl, changed in June 2010 !rA= rhow - (T_oc-TFrez(S_oc))*H_ML*cc/cl*(1.0-A) - A=A + 0.5_WP*min(rh,0.0_WP)*A/max(h,hmin) + max(rA,0.0_WP)*(1._WP-A)/lid_clo !/h0 + A=A + c_melt*min(rh,0.0_WP)*A/max(h,hmin) + max(rA,0.0_WP)*(1._WP-A)/lid_clo !/h0 !meaning: melting freezing A=min(A,h*1.e6_WP) ! A -> 0 for h -> 0 @@ -699,7 +735,7 @@ end subroutine budget ! ! !_______________________________________________________________________________ -subroutine obudget (ithermp, qa,fsh,flo,t,ug,ta,ch,ce,fh,evap,hflatow,hfsenow,hflwrdout) +subroutine obudget (ithermp, qa,fsh,flo,t,ug,ta,ch,ce,geolon, geolat,fh,evap,hflatow,hfsenow,hflwrdout) ! Ice growth rate for open ocean [m ice/sec] ! ! INPUT: @@ -711,20 +747,26 @@ subroutine obudget (ithermp, qa,fsh,flo,t,ug,ta,ch,ce,fh,evap,hflatow,hfsenow,hf ! ug - wind speed [m/sec] ! ch - transfer coefficient for sensible heat ! ce - transfer coefficient for evaporation + ! geolon - geographical longitude + ! geolat - geographical latitude ! ! OUTPUT: fh - growth rate ! evap - evaporation use MOD_ICE + use MOD_MESH use o_param, only: WP + use g_clock implicit none type(t_ice_thermo), intent(in), target :: ithermp real(kind=WP) qa,t,ta,fsh,flo,ug,ch,ce,fh,evap real(kind=WP) hfsenow,hfradow,hflatow,hftotow,hflwrdout,b real(kind=WP) q1, q2 ! coefficients for saturated specific humidity - real(kind=WP) c1, c4, c5 + real(kind=WP) c1, c4, c5, coszen, geolon, geolat + real(kind=WP), external :: solar_zenith_cosangle, albw_taylor, albw_briegleb logical :: standard_saturation_shum_formula = .true. integer :: ii !___________________________________________________________________________ + integer, pointer :: open_water_albedo real(kind=WP), pointer :: boltzmann, emiss_wat, inv_rhowat, inv_rhoair, rhoair, & tmelt, cl, clhw, cpair, albw boltzmann => ithermp%boltzmann @@ -737,6 +779,7 @@ subroutine obudget (ithermp, qa,fsh,flo,t,ug,ta,ch,ce,fh,evap,hflatow,hfsenow,hf clhw => ithermp%clhw cpair => ithermp%cpair albw => ithermp%albw + open_water_albedo => ithermp%open_water_albedo !___________________________________________________________________________ c1 = 3.8e-3_WP @@ -744,6 +787,14 @@ subroutine obudget (ithermp, qa,fsh,flo,t,ug,ta,ch,ce,fh,evap,hflatow,hfsenow,hf c5 = 237.3_WP q1 = 640380._WP q2 = -5107.4_WP + if(open_water_albedo > 0)then + coszen=solar_zenith_cosangle(daynew, timenew/3600., geolon, geolat) + if(open_water_albedo > 1)then + albw=albw_briegleb(coszen) + else + albw=albw_taylor(coszen) + endif + endif ! (saturated) surface specific humidity if(standard_saturation_shum_formula) then @@ -815,6 +866,95 @@ function TFrez(S) TFrez= -0.0575_WP*S+1.7105e-3_WP *sqrt(S**3)-2.155e-4_WP *S*S end function TFrez + + +!----------------------------------------------------------------------- +! Copyright by the GOTM-team under the GNU Public License - www.gnu.org +!----------------------------------------------------------------------- +!----------------------------------------------------------------------- +! +! ROUTINE: Calculate the solar zenith angle +! +! INTERFACE: +function solar_zenith_cosangle(yday, hh, rlon, rlat) + ! + ! DESCRIPTION: + ! This subroutine calculates the solar zenith angle as being used both + ! in albedo_water() and short_wave_radiation(). The result is in degrees. + ! + use o_param, only: WP + + IMPLICIT NONE + + integer, intent(in) :: yday + real(kind=WP), intent(in) :: hh + real(kind=WP), intent(in) :: rlon,rlat + ! + ! REVISION HISTORY: + ! Original author(s): Karsten Bolding + ! + real(kind=WP), parameter :: pi=3.14159265358979323846 + real(kind=WP), parameter :: deg2rad=pi/180. + real(kind=WP), parameter :: rad2deg=180./pi + + real(kind=WP) :: yrdays, th0, th02, th03, sundec, thsun, solar_zenith_cosangle + + ! from now on everything in radians + ! rlon = deg2rad*dlon + ! rlat = deg2rad*dlat + + yrdays=365.25_WP + th0 = 2._WP * pi * yday/yrdays + th02 = 2._WP * th0 + th03 = 3._WP * th0 + + ! sun declination : + sundec = 0.006918_WP - 0.399912_WP * cos(th0) + 0.070257_WP * sin(th0) & + - 0.006758_WP * cos(th02) + 0.000907_WP * sin(th02) & + - 0.002697*cos(th03) + 0.001480*sin(th03) + + ! sun hour angle : + thsun = (hh-12._WP) * 15._WP * deg2rad + rlon + + ! cosine of the solar zenith angle : + solar_zenith_cosangle = sin(rlat) * sin(sundec) + cos(rlat) * cos(sundec) * cos(thsun) + if (solar_zenith_cosangle .lt. 0.0_WP) solar_zenith_cosangle = 0.0_WP + return + end function solar_zenith_cosangle + + function albw_taylor(coszen) + ! purpose: zenith depending open water albedo + ! method: taylor et al. (1996) + ! author: frank kauker + ! date: 26. april 2017 + use o_param, only: WP + implicit none + + real(kind=WP), intent(in) :: coszen + real(kind=WP) :: albw_taylor + + albw_taylor = 0.037_WP/(1.1_WP * coszen**1.4_WP + 0.15_WP) + + return + end function albw_taylor + + function albw_briegleb(coszen) + ! purpose: zenith depending open water albedo + ! method: briegleb et al. (1986) + ! author: frank kauker + ! date: 26. april 2017 + + use o_param, only: WP + implicit none + + real(kind=WP), intent(in) :: coszen + real(kind=WP) :: albw_briegleb + + albw_briegleb = 0.026_WP/(1.1_WP * coszen**1.7_WP + 0.065_WP) + 0.15_WP * (coszen-1._WP)**2 * (coszen-0.5_WP) + + return + end function albw_briegleb + ! ! ! !_______________________________________________________________________________ From 7be0379ed7e840f4c6be9a278e6bb87124b140d7 Mon Sep 17 00:00:00 2001 From: Nikolay Koldunov Date: Mon, 3 Mar 2025 10:53:49 +0100 Subject: [PATCH 2/8] running version of Frank's sea ice --- src/MOD_ICE.F90 | 35 +++++++++++++++++++++++++++-------- src/gen_forcing_init.F90 | 5 ++++- src/gen_modules_forcing.F90 | 2 +- src/ice_thermo_oce.F90 | 20 ++++++++++---------- 4 files changed, 42 insertions(+), 20 deletions(-) diff --git a/src/MOD_ICE.F90 b/src/MOD_ICE.F90 index 0375d37d1..de95b48a5 100644 --- a/src/MOD_ICE.F90 +++ b/src/MOD_ICE.F90 @@ -70,7 +70,7 @@ MODULE MOD_ICE ! --- namelist parameter /ice_therm/ real(kind=WP) :: con= 2.1656, consn = 0.31 ! Thermal conductivities: ice & snow; W/m/K real(kind=WP) :: Sice = 4.0 ! Ice salinity 3.2--5.0 ppt. - real(kind=WP) :: h0=1.0 ! Lead closing parameter [m] ! 0.5 + real(kind=WP) :: h0=0.5 ! Lead closing parameter [m] ! 0.5 real(kind=WP) :: emiss_ice=0.97 ! Emissivity of Snow/Ice, real(kind=WP) :: emiss_wat=0.97 ! Emissivity of open water real(kind=WP) :: albsn = 0.81 ! Albedo: frozen snow @@ -78,6 +78,17 @@ MODULE MOD_ICE real(kind=WP) :: albi = 0.70 ! frozen ice real(kind=WP) :: albim = 0.68 ! melting ice real(kind=WP) :: albw = 0.066 ! open water, LY2004 + + + ! --- additional namelist parameters (Frank.Kauker(at)awi.de 2023/04/04) + logical :: snowdist=.true. ! distribution of snow depth according to ice distribution + logical :: new_iclasses=.false. ! ice thickness distribution based on EM observations (Castro-Morales et al., JGR, 2013) + integer :: open_water_albedo=0 ! 0=standard; 1=taylor; 2=briegleb + REAL(kind=WP) :: c_melt=0.5 ! constant in concentration equation for melting conditions + REAL(kind=WP) :: h_cutoff=3.0 ! cutoff thickness of thickness pdf + REAL(kind=WP), DIMENSION(15) :: hpdf = (/ 0.066745491, 0.1462317, 0.17769822, 0.13131106, & + 0.11518432, 0.08514193, 0.06871303, 0.05592151, 0.04428673, 0.03584652, 0.02970195, 0.02469673, & + 0.02001543, 0.01653681, 0.0141026 /) ! pdf of ice thickness based on EM observations contains procedure WRITE_T_ICE_THERMO procedure READ_T_ICE_THERMO @@ -549,11 +560,12 @@ subroutine ice_init(ice, partit, mesh) namelist /ice_dyn/ whichEVP, Pstar, ellipse, c_pressure, delta_min, evp_rheol_steps, & Cd_oce_ice, ice_gamma_fct, ice_diff, theta_io, ice_ave_steps, & alpha_evp, beta_evp, c_aevp - + logical :: snowdist, new_iclasses + integer :: open_water_albedo, iclasses real(kind=WP) :: Sice, h0, emiss_ice, emiss_wat, albsn, albsnm, albi, & - albim, albw, con, consn - namelist /ice_therm/ Sice, h0, emiss_ice, emiss_wat, albsn, albsnm, albi, & - albim, albw, con, consn + albim, albw, con, consn, hmin, armin, c_melt, h_cutoff + namelist /ice_therm/ Sice, iclasses, h0, hmin, armin, emiss_ice, emiss_wat, albsn, albsnm, albi, & + albim, albw, con, consn, snowdist, new_iclasses, open_water_albedo, c_melt, h_cutoff !___________________________________________________________________________ ! pointer on necessary derived types #include "associate_part_def.h" @@ -596,7 +608,10 @@ subroutine ice_init(ice, partit, mesh) ice%thermo%con = con ice%thermo%consn = consn ice%thermo%Sice = Sice + ice%thermo%iclasses = iclasses ice%thermo%h0 = h0 + ice%thermo%hmin = hmin + ice%thermo%armin = armin ice%thermo%emiss_ice= emiss_ice ice%thermo%emiss_wat= emiss_wat ice%thermo%albsn = albsn @@ -604,9 +619,13 @@ subroutine ice_init(ice, partit, mesh) ice%thermo%albi = albi ice%thermo%albim = albim ice%thermo%albw = albw - - ice%thermo%cc=ice%thermo%rhowat*4190.0 ! Volumetr. heat cap. of water [J/m**3/K](cc = rhowat*cp_water) - ice%thermo%cl=ice%thermo%rhoice*3.34e5 ! Volumetr. latent heat of ice fusion [J/m**3](cl=rhoice*Lf) + ice%thermo%snowdist = snowdist + ice%thermo%new_iclasses=new_iclasses + ice%thermo%open_water_albedo=open_water_albedo + ice%thermo%c_melt = c_melt + ice%thermo%h_cutoff = h_cutoff + ice%thermo%cc =ice%thermo%rhowat*4190.0 ! Volumetr. heat cap. of water [J/m**3/K](cc = rhowat*cp_water) + ice%thermo%cl =ice%thermo%rhoice*3.34e5 ! Volumetr. latent heat of ice fusion [J/m**3](cl=rhoice*Lf) !___________________________________________________________________________ ! define local vertice & elem array size diff --git a/src/gen_forcing_init.F90 b/src/gen_forcing_init.F90 index 8d92be366..52dbe8edc 100644 --- a/src/gen_forcing_init.F90 +++ b/src/gen_forcing_init.F90 @@ -137,10 +137,13 @@ subroutine forcing_array_setup(partit, mesh) allocate(Tair(n2), shum(n2)) Tair=0.0_WP shum=0.0_WP - allocate(runoff(n2), evaporation(n2),ice_sublimation(n2)) + allocate(runoff(n2), evaporation(n2),ice_sublimation(n2), ice_ocean_hflx(n2), ice_dTfrez(n2), ice_ustar(n2)) runoff=0.0_WP evaporation = 0.0_WP ice_sublimation = 0.0_WP + ice_ocean_hflx = 0.0_WP + ice_dTfrez = 0.0_WP + ice_ustar = 0.0_WP #if defined (__oasis) || defined (__ifsinterface) allocate(sublimation(n2), evap_no_ifrac(n2)) diff --git a/src/gen_modules_forcing.F90 b/src/gen_modules_forcing.F90 index 2c27b7118..811dc25c7 100755 --- a/src/gen_modules_forcing.F90 +++ b/src/gen_modules_forcing.F90 @@ -70,7 +70,7 @@ module g_forcing_arrays real(kind=WP), allocatable, dimension(:,:) :: Tair_t, shum_t real(kind=WP), allocatable, dimension(:) :: shortwave, longwave real(kind=WP), allocatable, dimension(:) :: prec_rain, prec_snow - real(kind=WP), allocatable, dimension(:) :: runoff, evaporation, ice_sublimation + real(kind=WP), allocatable, dimension(:) :: runoff, evaporation, ice_sublimation, ice_ocean_hflx, ice_dTfrez, ice_ustar real(kind=WP), allocatable, dimension(:) :: cloudiness, press_air !---wiso-code real(kind=WP), allocatable, dimension(:) :: www1,www2,www3,iii1,iii2,iii3 diff --git a/src/ice_thermo_oce.F90 b/src/ice_thermo_oce.F90 index 89d942f76..3eb9e44f3 100755 --- a/src/ice_thermo_oce.F90 +++ b/src/ice_thermo_oce.F90 @@ -29,7 +29,7 @@ subroutine therm_ice(ithermp, h,hsn,A,fsh,flo,Ta,qa,rain,snow,runo,rsss, & rsf, dhgrowth, dhsngrowth, iflice, hflatow, hfsenow, hflwrdout,lid_clo,geolon, geolat, subli) USE MOD_ICE type(t_ice_thermo), intent(in), target :: ithermp - real(kind=WP) h, hsn, A, fsh, flo, Ta, qa, rain, snow, runo, rsss, evap_in, & + real(kind=WP) h, hsn, A, fsh, flo, Ta, qa, rain, snow, runo, rsss, & ug, ustar, T_oc, S_oc, H_ML, t, ice_dt, ch, ce, ch_i, ce_i, evap_in, fw, ehf, & dhgrowth, dhsngrowth, ahf, prec, subli, subli_i, rsf, & rhow, show, rhice, shice, sh, thick, thact, lat, & @@ -37,7 +37,7 @@ subroutine therm_ice(ithermp, h,hsn,A,fsh,flo,Ta,qa,rain,snow,runo,rsss, & hfsenow, hflwrdout, lid_clo, geolon, geolat end subroutine therm_ice end interface -end module ice_therm_interface +end module module ice_budget_interfaces interface @@ -248,14 +248,14 @@ subroutine thermodynamics(ice, partit, mesh) fsh = shortwave(i) flo = longwave(i) Ta = Tair(i) - if (l_tdew) then - call - qa = 0. - elseif (l_humi) then - qa = shum(i) - else - stop 'neither specific humidity nor dew-point temperature is defined!' - endif + ! if (l_tdew) then + ! !call ! I am not sure what is meant here + ! qa = 0. + ! elseif (l_humi) then + ! qa = shum(i) + ! else + ! stop 'neither specific humidity nor dew-point temperature is defined!' + ! endif if (.not. l_snow) then if (Ta>=0.0_WP) then From bd3bdbc232094a9ce0ca869e9896d639b4ac958e Mon Sep 17 00:00:00 2001 From: Nikolay Koldunov Date: Mon, 3 Mar 2025 11:31:38 +0100 Subject: [PATCH 3/8] update namelist --- config/namelist.ice | 8 ++++++++ src/MOD_ICE.F90 | 2 +- src/ice_thermo_oce.F90 | 9 ++++++--- 3 files changed, 15 insertions(+), 4 deletions(-) diff --git a/config/namelist.ice b/config/namelist.ice index bcd86f145..8897a4ebe 100644 --- a/config/namelist.ice +++ b/config/namelist.ice @@ -18,7 +18,10 @@ ice_ave_steps=1 ! ice step=ice_ave_steps*oce_step &ice_therm Sice=4.0 ! Ice salinity 3.2--5.0 ppt. +iclasses=15 ! default = 7; in case of EM distribution ('new_iceclasses=.true.') must be set to 15 h0=.5 ! Lead closing parameter [m] +hmin=0.05 ! default=0.01 +armin=0.15 ! default=0.01 emiss_ice=0.97 ! Emissivity of Snow/Ice, emiss_wat=0.97 ! Emissivity of open water albsn=0.81 ! Albedo: frozen snow @@ -28,4 +31,9 @@ albim=0.68 ! melting ice albw=0.1 ! open water con=2.1656 ! Thermal conductivities: ice; W/m/K consn=0.31 ! snow +snowdist=.true. ! distribution of snow depth according to ice distribution - default: .true. +new_iclasses=.false. ! default=.false.; ice thickness distribution based on EM observations (Castro-Morales et al., JGR, 2013) +open_water_albedo=0 ! 0=default; 1=taylor; 2=briegleb +c_melt=0.5 ! constant in concentration equation for melting conditions - default=0.5 +h_cutoff=3.0 ! only used for new_iclasses=.true. / diff --git a/src/MOD_ICE.F90 b/src/MOD_ICE.F90 index de95b48a5..ec2f786f2 100644 --- a/src/MOD_ICE.F90 +++ b/src/MOD_ICE.F90 @@ -70,7 +70,7 @@ MODULE MOD_ICE ! --- namelist parameter /ice_therm/ real(kind=WP) :: con= 2.1656, consn = 0.31 ! Thermal conductivities: ice & snow; W/m/K real(kind=WP) :: Sice = 4.0 ! Ice salinity 3.2--5.0 ppt. - real(kind=WP) :: h0=0.5 ! Lead closing parameter [m] ! 0.5 + real(kind=WP) :: h0=1.0 ! Lead closing parameter [m] ! 0.5 real(kind=WP) :: emiss_ice=0.97 ! Emissivity of Snow/Ice, real(kind=WP) :: emiss_wat=0.97 ! Emissivity of open water real(kind=WP) :: albsn = 0.81 ! Albedo: frozen snow diff --git a/src/ice_thermo_oce.F90 b/src/ice_thermo_oce.F90 index 3eb9e44f3..b99a23f81 100755 --- a/src/ice_thermo_oce.F90 +++ b/src/ice_thermo_oce.F90 @@ -101,9 +101,10 @@ subroutine cut_off(ice, partit, mesh) ! lower cutoff: a_ice if (a_ice(n) < .1e-8_WP) then a_ice(n)=0.0_WP + m_ice(n) =0.0_WP + m_snow(n) =0.0_WP #if defined (__oifs) || defined (__ifsinterface) - m_ice(n) =0.0_WP - m_snow(n) =0.0_WP + ice_temp(n)=273.15_WP #endif /* (__oifs) */ end if @@ -111,9 +112,9 @@ subroutine cut_off(ice, partit, mesh) ! lower cutoff: m_ice if (m_ice(n) < .1e-8_WP) then m_ice(n)=0.0_WP -#if defined (__oifs) || defined (__ifsinterface) m_snow(n) =0.0_WP a_ice(n) =0.0_WP +#if defined (__oifs) || defined (__ifsinterface) ice_temp(n)=273.15_WP #endif /* (__oifs) */ end if @@ -248,6 +249,7 @@ subroutine thermodynamics(ice, partit, mesh) fsh = shortwave(i) flo = longwave(i) Ta = Tair(i) + qa = shum(i) ! if (l_tdew) then ! !call ! I am not sure what is meant here ! qa = 0. @@ -584,6 +586,7 @@ subroutine therm_ice(ithermp, h, hsn, A, fsh, flo, Ta, qa, rain, snow, runo, rss rsf= -dhgrowth*rhoice*inv_rhowat*Sice else fw= prec+evap - dhgrowth*rhoice*inv_rhowat*(rsss-Sice)/rsss - dhsngrowth*rhosno*inv_rhowat + rsf = 0.0_WP end if ! Changes in compactnesses (equation 16 of Hibler 1979) From dd22bfc37410384c488a36a928caa6859f64253c Mon Sep 17 00:00:00 2001 From: Nikolay Koldunov Date: Mon, 3 Mar 2025 11:37:28 +0100 Subject: [PATCH 4/8] update namelist --- config/namelist.ice | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/config/namelist.ice b/config/namelist.ice index 8897a4ebe..292ebe031 100644 --- a/config/namelist.ice +++ b/config/namelist.ice @@ -20,8 +20,8 @@ ice_ave_steps=1 ! ice step=ice_ave_steps*oce_step Sice=4.0 ! Ice salinity 3.2--5.0 ppt. iclasses=15 ! default = 7; in case of EM distribution ('new_iceclasses=.true.') must be set to 15 h0=.5 ! Lead closing parameter [m] -hmin=0.05 ! default=0.01 -armin=0.15 ! default=0.01 +hmin=0.01 ! default=0.01 +armin=0.01 ! default=0.01 emiss_ice=0.97 ! Emissivity of Snow/Ice, emiss_wat=0.97 ! Emissivity of open water albsn=0.81 ! Albedo: frozen snow From a03ec06ecbdd137763c69715d719f35a7f1963da Mon Sep 17 00:00:00 2001 From: Nikolay Koldunov Date: Mon, 3 Mar 2025 21:37:42 +0100 Subject: [PATCH 5/8] update namelist for iceclasses --- config/namelist.ice | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config/namelist.ice b/config/namelist.ice index 292ebe031..f11dd8b8b 100644 --- a/config/namelist.ice +++ b/config/namelist.ice @@ -18,7 +18,7 @@ ice_ave_steps=1 ! ice step=ice_ave_steps*oce_step &ice_therm Sice=4.0 ! Ice salinity 3.2--5.0 ppt. -iclasses=15 ! default = 7; in case of EM distribution ('new_iceclasses=.true.') must be set to 15 +iclasses=7 ! default = 7; in case of EM distribution ('new_iceclasses=.true.') must be set to 15 h0=.5 ! Lead closing parameter [m] hmin=0.01 ! default=0.01 armin=0.01 ! default=0.01 From eec9052120203d166895f9ec634d6c76f30d1058 Mon Sep 17 00:00:00 2001 From: Nikolay Koldunov Date: Sat, 8 Mar 2025 12:14:45 +0100 Subject: [PATCH 6/8] remove unnecesary variables --- src/gen_forcing_init.F90 | 6 ++---- src/gen_modules_forcing.F90 | 2 +- src/ice_thermo_oce.F90 | 12 +----------- 3 files changed, 4 insertions(+), 16 deletions(-) diff --git a/src/gen_forcing_init.F90 b/src/gen_forcing_init.F90 index 52dbe8edc..73184b362 100644 --- a/src/gen_forcing_init.F90 +++ b/src/gen_forcing_init.F90 @@ -137,13 +137,11 @@ subroutine forcing_array_setup(partit, mesh) allocate(Tair(n2), shum(n2)) Tair=0.0_WP shum=0.0_WP - allocate(runoff(n2), evaporation(n2),ice_sublimation(n2), ice_ocean_hflx(n2), ice_dTfrez(n2), ice_ustar(n2)) + allocate(runoff(n2), evaporation(n2),ice_sublimation(n2)) runoff=0.0_WP evaporation = 0.0_WP ice_sublimation = 0.0_WP - ice_ocean_hflx = 0.0_WP - ice_dTfrez = 0.0_WP - ice_ustar = 0.0_WP + #if defined (__oasis) || defined (__ifsinterface) allocate(sublimation(n2), evap_no_ifrac(n2)) diff --git a/src/gen_modules_forcing.F90 b/src/gen_modules_forcing.F90 index 811dc25c7..2c27b7118 100755 --- a/src/gen_modules_forcing.F90 +++ b/src/gen_modules_forcing.F90 @@ -70,7 +70,7 @@ module g_forcing_arrays real(kind=WP), allocatable, dimension(:,:) :: Tair_t, shum_t real(kind=WP), allocatable, dimension(:) :: shortwave, longwave real(kind=WP), allocatable, dimension(:) :: prec_rain, prec_snow - real(kind=WP), allocatable, dimension(:) :: runoff, evaporation, ice_sublimation, ice_ocean_hflx, ice_dTfrez, ice_ustar + real(kind=WP), allocatable, dimension(:) :: runoff, evaporation, ice_sublimation real(kind=WP), allocatable, dimension(:) :: cloudiness, press_air !---wiso-code real(kind=WP), allocatable, dimension(:) :: www1,www2,www3,iii1,iii2,iii3 diff --git a/src/ice_thermo_oce.F90 b/src/ice_thermo_oce.F90 index b99a23f81..f20b7f590 100755 --- a/src/ice_thermo_oce.F90 +++ b/src/ice_thermo_oce.F90 @@ -250,14 +250,7 @@ subroutine thermodynamics(ice, partit, mesh) flo = longwave(i) Ta = Tair(i) qa = shum(i) - ! if (l_tdew) then - ! !call ! I am not sure what is meant here - ! qa = 0. - ! elseif (l_humi) then - ! qa = shum(i) - ! else - ! stop 'neither specific humidity nor dew-point temperature is defined!' - ! endif + if (.not. l_snow) then if (Ta>=0.0_WP) then @@ -320,9 +313,6 @@ subroutine thermodynamics(ice, partit, mesh) net_heat_flux(i) = ehf !positive down evaporation(i) = evap !negative up ice_sublimation(i)= subli - ice_ocean_hflx(i) = o2ihf - ice_dTfrez(i) = T_oc-TFrez(S_oc) - ice_ustar(i) = ustar thdgr(i) = ithdgr thdgrsn(i) = ithdgrsn From d1fa4f6b88fd69bd72e1660bbd40b7e8437075a3 Mon Sep 17 00:00:00 2001 From: Nikolay Koldunov Date: Sat, 8 Mar 2025 12:32:53 +0100 Subject: [PATCH 7/8] rewrite GPL function to be complient with Apache license --- src/ice_thermo_oce.F90 | 116 +++++++++++++++++++++++------------------ 1 file changed, 65 insertions(+), 51 deletions(-) diff --git a/src/ice_thermo_oce.F90 b/src/ice_thermo_oce.F90 index f20b7f590..bc69dbdb2 100755 --- a/src/ice_thermo_oce.F90 +++ b/src/ice_thermo_oce.F90 @@ -755,7 +755,7 @@ subroutine obudget (ithermp, qa,fsh,flo,t,ug,ta,ch,ce,geolon, geolat,fh,evap,hfl real(kind=WP) hfsenow,hfradow,hflatow,hftotow,hflwrdout,b real(kind=WP) q1, q2 ! coefficients for saturated specific humidity real(kind=WP) c1, c4, c5, coszen, geolon, geolat - real(kind=WP), external :: solar_zenith_cosangle, albw_taylor, albw_briegleb + real(kind=WP), external :: compute_solar_zenith_angle, albw_taylor, albw_briegleb logical :: standard_saturation_shum_formula = .true. integer :: ii !___________________________________________________________________________ @@ -781,7 +781,7 @@ subroutine obudget (ithermp, qa,fsh,flo,t,ug,ta,ch,ce,geolon, geolat,fh,evap,hfl q1 = 640380._WP q2 = -5107.4_WP if(open_water_albedo > 0)then - coszen=solar_zenith_cosangle(daynew, timenew/3600., geolon, geolat) + coszen=compute_solar_zenith_angle(daynew, timenew/3600., geolon, geolat) if(open_water_albedo > 1)then albw=albw_briegleb(coszen) else @@ -861,59 +861,73 @@ function TFrez(S) end function TFrez -!----------------------------------------------------------------------- -! Copyright by the GOTM-team under the GNU Public License - www.gnu.org -!----------------------------------------------------------------------- -!----------------------------------------------------------------------- -! -! ROUTINE: Calculate the solar zenith angle -! -! INTERFACE: -function solar_zenith_cosangle(yday, hh, rlon, rlat) +function compute_solar_zenith_angle(day_of_year, hour_utc, longitude, latitude) result(cos_zenith) + !----------------------------------------------------------------------- + ! Purpose: Compute the cosine of the solar zenith angle given the day of the year, + ! time in UTC, and geographic coordinates. ! - ! DESCRIPTION: - ! This subroutine calculates the solar zenith angle as being used both - ! in albedo_water() and short_wave_radiation(). The result is in degrees. + ! Inputs: + ! - day_of_year: Integer, day of the year (1 to 365 or 366) + ! - hour_utc: Real, time in UTC (decimal hours) + ! - longitude: Real, longitude in radians + ! - latitude: Real, latitude in radians ! - use o_param, only: WP - - IMPLICIT NONE - - integer, intent(in) :: yday - real(kind=WP), intent(in) :: hh - real(kind=WP), intent(in) :: rlon,rlat - ! - ! REVISION HISTORY: - ! Original author(s): Karsten Bolding + ! Output: + ! - cos_zenith: Real, cosine of the solar zenith angle ! - real(kind=WP), parameter :: pi=3.14159265358979323846 - real(kind=WP), parameter :: deg2rad=pi/180. - real(kind=WP), parameter :: rad2deg=180./pi - - real(kind=WP) :: yrdays, th0, th02, th03, sundec, thsun, solar_zenith_cosangle - - ! from now on everything in radians - ! rlon = deg2rad*dlon - ! rlat = deg2rad*dlat - - yrdays=365.25_WP - th0 = 2._WP * pi * yday/yrdays - th02 = 2._WP * th0 - th03 = 3._WP * th0 - - ! sun declination : - sundec = 0.006918_WP - 0.399912_WP * cos(th0) + 0.070257_WP * sin(th0) & - - 0.006758_WP * cos(th02) + 0.000907_WP * sin(th02) & - - 0.002697*cos(th03) + 0.001480*sin(th03) - - ! sun hour angle : - thsun = (hh-12._WP) * 15._WP * deg2rad + rlon - - ! cosine of the solar zenith angle : - solar_zenith_cosangle = sin(rlat) * sin(sundec) + cos(rlat) * cos(sundec) * cos(thsun) - if (solar_zenith_cosangle .lt. 0.0_WP) solar_zenith_cosangle = 0.0_WP + ! Notes: + ! - Based on standard solar position calculations + ! - Uses an approximation for solar declination + ! + ! Written for Apache 2.0 licensed projects. + !----------------------------------------------------------------------- + + use o_param, only: WP ! Ensure precision consistency + + implicit none + + ! Input variables + integer, intent(in) :: day_of_year + real(kind=WP), intent(in) :: hour_utc + real(kind=WP), intent(in) :: longitude, latitude + + ! Output variable + real(kind=WP) :: cos_zenith + + ! Constants + real(kind=WP), parameter :: PI = 3.141592653589793_WP + real(kind=WP), parameter :: DEG_TO_RAD = PI / 180.0_WP + real(kind=WP), parameter :: RAD_TO_DEG = 180.0_WP / PI + real(kind=WP), parameter :: DAYS_PER_YEAR = 365.25_WP + + ! Computed values + real(kind=WP) :: solar_declination, hour_angle + real(kind=WP) :: solar_fraction + + ! Compute solar position parameters + solar_fraction = 2.0_WP * PI * day_of_year / DAYS_PER_YEAR + + ! Approximate solar declination using Fourier terms + solar_declination = 0.006918_WP - 0.399912_WP * cos(solar_fraction) & + + 0.070257_WP * sin(solar_fraction) & + - 0.006758_WP * cos(2.0_WP * solar_fraction) & + + 0.000907_WP * sin(2.0_WP * solar_fraction) & + - 0.002697_WP * cos(3.0_WP * solar_fraction) & + + 0.001480_WP * sin(3.0_WP * solar_fraction) + + ! Compute hour angle (relative to noon, converted to radians) + hour_angle = (hour_utc - 12.0_WP) * 15.0_WP * DEG_TO_RAD + longitude + + ! Compute cosine of the solar zenith angle + cos_zenith = sin(latitude) * sin(solar_declination) & + + cos(latitude) * cos(solar_declination) * cos(hour_angle) + + ! Ensure cosine is non-negative (no below-horizon values) + if (cos_zenith < 0.0_WP) cos_zenith = 0.0_WP + return - end function solar_zenith_cosangle +end function compute_solar_zenith_angle + function albw_taylor(coszen) ! purpose: zenith depending open water albedo From 83b673f1f5787e6d0bdcecc5c36267d246389834 Mon Sep 17 00:00:00 2001 From: Nikolay Koldunov Date: Sun, 9 Mar 2025 12:09:20 +0100 Subject: [PATCH 8/8] split h0 to NH and SH, also add h0_s to namelist and include namelist.ice_ERA5 to config --- config/namelist.ice | 3 ++- config/namelist.ice_ERA5 | 40 ++++++++++++++++++++++++++++++++++++++++ src/MOD_ICE.F90 | 8 +++++--- src/ice_thermo_oce.F90 | 15 ++++++++------- 4 files changed, 55 insertions(+), 11 deletions(-) create mode 100644 config/namelist.ice_ERA5 diff --git a/config/namelist.ice b/config/namelist.ice index f11dd8b8b..bf7ea69bd 100644 --- a/config/namelist.ice +++ b/config/namelist.ice @@ -19,7 +19,8 @@ ice_ave_steps=1 ! ice step=ice_ave_steps*oce_step &ice_therm Sice=4.0 ! Ice salinity 3.2--5.0 ppt. iclasses=7 ! default = 7; in case of EM distribution ('new_iceclasses=.true.') must be set to 15 -h0=.5 ! Lead closing parameter [m] +h0=0.5 ! Lead closing parameter for Nothern Hemisphere [m], default 0.5 +h0_s=0.5 ! Lead closing parameter [m] for Southern Hemisphere, default 0.5 hmin=0.01 ! default=0.01 armin=0.01 ! default=0.01 emiss_ice=0.97 ! Emissivity of Snow/Ice, diff --git a/config/namelist.ice_ERA5 b/config/namelist.ice_ERA5 new file mode 100644 index 000000000..9a57e0032 --- /dev/null +++ b/config/namelist.ice_ERA5 @@ -0,0 +1,40 @@ +! Ice namelist +&ice_dyn +whichEVP=0 ! 0=standard; 1=mEVP; 2=aEVP +Pstar=25000 ! standard: 30000 [N/m^2] +ellipse=2.0 +c_pressure=8. ! ice concentration parameter used in ice strength computation - standard: 20 +delta_min=1.0e-11 ! [s^(-1)] +evp_rheol_steps=120 ! number of EVP subcycles +alpha_evp=500 ! constant that control numerical stability of mEVP. Adjust with resolution. +beta_evp=500 ! constant that control numerical stability of mEVP. Adjust with resolution. +c_aevp=0.15 ! a tuning constant in aEVP. Adjust with resolution. +Cd_oce_ice=0.0085 ! drag coef. oce - ice - standard: 0.0055 +ice_gamma_fct=0.5 ! smoothing parameter +ice_diff=0.0 ! diffusion to stabilize +theta_io=0.0 ! rotation angle +ice_ave_steps=1 ! ice step=ice_ave_steps*oce_step +/ + +&ice_therm +Sice=4.0 ! Ice salinity 3.2--5.0 ppt. +iclasses=15 ! default = 7; in case of EM distribution ('new_iceclasses=.true.') must be set to 15 +h0=0.8 ! Lead closing parameter for Nothern Hemisphere [m], default 0.5 +h0_s=0.5 ! Lead closing parameter [m] for Southern Hemisphere, default 0.5 +hmin=0.05 ! default=0.01 +armin=0.15 ! default=0.01 +emiss_ice=0.97 ! Emissivity of Snow/Ice, +emiss_wat=0.97 ! Emissivity of open water +albsn=0.81 ! Albedo: frozen snow (default: AOMIP albedos from Perovic (1986)) +albsnm=0.77 ! melting snow +albi=0.7 ! frozen ice +albim=0.68 ! melting ice +albw=0.1 ! open water +con=2.1656 ! Thermal conductivities: ice; W/m/K +consn=0.31 ! snow +snowdist=.false. ! distribution of snow depth according to ice distribution - default: .true. +new_iclasses=.true. ! default=.false.; ice thickness distribution based on EM observations (Castro-Morales et al., JGR, 2013) +open_water_albedo=2 ! 0=default; 1=taylor; 2=briegleb +c_melt=0.3 ! constant in concentration equation for melting conditions - default=0.5 +h_cutoff=3.0 ! only used for new_iclasses=.true. +/ diff --git a/src/MOD_ICE.F90 b/src/MOD_ICE.F90 index ec2f786f2..b1d175e3b 100644 --- a/src/MOD_ICE.F90 +++ b/src/MOD_ICE.F90 @@ -70,7 +70,8 @@ MODULE MOD_ICE ! --- namelist parameter /ice_therm/ real(kind=WP) :: con= 2.1656, consn = 0.31 ! Thermal conductivities: ice & snow; W/m/K real(kind=WP) :: Sice = 4.0 ! Ice salinity 3.2--5.0 ppt. - real(kind=WP) :: h0=1.0 ! Lead closing parameter [m] ! 0.5 + real(kind=WP) :: h0=0.5 ! Lead closing parameter [m] for Nothern Hemisphere! 0.5 + real(kind=WP) :: h0_s=0.5 ! Lead closing parameter [m] for Southern Hemisphere! 0.5 real(kind=WP) :: emiss_ice=0.97 ! Emissivity of Snow/Ice, real(kind=WP) :: emiss_wat=0.97 ! Emissivity of open water real(kind=WP) :: albsn = 0.81 ! Albedo: frozen snow @@ -562,9 +563,9 @@ subroutine ice_init(ice, partit, mesh) alpha_evp, beta_evp, c_aevp logical :: snowdist, new_iclasses integer :: open_water_albedo, iclasses - real(kind=WP) :: Sice, h0, emiss_ice, emiss_wat, albsn, albsnm, albi, & + real(kind=WP) :: Sice, h0, h0_s, emiss_ice, emiss_wat, albsn, albsnm, albi, & albim, albw, con, consn, hmin, armin, c_melt, h_cutoff - namelist /ice_therm/ Sice, iclasses, h0, hmin, armin, emiss_ice, emiss_wat, albsn, albsnm, albi, & + namelist /ice_therm/ Sice, iclasses, h0, h0_s, hmin, armin, emiss_ice, emiss_wat, albsn, albsnm, albi, & albim, albw, con, consn, snowdist, new_iclasses, open_water_albedo, c_melt, h_cutoff !___________________________________________________________________________ ! pointer on necessary derived types @@ -610,6 +611,7 @@ subroutine ice_init(ice, partit, mesh) ice%thermo%Sice = Sice ice%thermo%iclasses = iclasses ice%thermo%h0 = h0 + ice%thermo%h0_s = h0_s ice%thermo%hmin = hmin ice%thermo%armin = armin ice%thermo%emiss_ice= emiss_ice diff --git a/src/ice_thermo_oce.F90 b/src/ice_thermo_oce.F90 index bc69dbdb2..9299dc48a 100755 --- a/src/ice_thermo_oce.F90 +++ b/src/ice_thermo_oce.F90 @@ -220,7 +220,7 @@ subroutine thermodynamics(ice, partit, mesh) ! Friction velocity !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i, j, elem, h, hsn, A, fsh, flo, Ta, qa, rain, snow, runo, rsss, rsf, evap_in, ug, ustar, T_oc, S_oc, & !$OMP h_ml, t, ch, ce, ch_i, ce_i, fw, ehf, evap, ithdgr, ithdgrsn, iflice, hflatow, hfsenow, hflwrdout, & -!$OMP subli, lid_clo, lat) +!$OMP subli, lid_clo, lat, geolon, geolat, o2ihf) !$OMP DO do i=1, myDim_nod2D ustar=0.0_WP @@ -280,15 +280,16 @@ subroutine thermodynamics(ice, partit, mesh) h_ml = 2.5_WP ! 10.0 or 30. used previously fw = 0.0_WP ehf = 0.0_WP - lid_Clo=ice%thermo%h0 - if (geo_coord_nod2D(2, i)>0) then !TODO 2 separate pars for each hemisphere - lid_clo=0.5_WP - else - lid_clo=0.5_WP - endif geolon = geo_coord_nod2D(1, i) geolat = geo_coord_nod2D(2, i) + if (geolat>0) then !TODO 2 separate pars for each hemisphere + lid_clo=ice%thermo%h0 + else + lid_clo=ice%thermo%h0_s + endif + + !_______________________________________________________________________ ! do ice thermodynamics call therm_ice(ice%thermo,h,hsn,A,fsh,flo,Ta,qa,rain,snow,runo,rsss, &