diff --git a/CMakeLists.txt b/CMakeLists.txt index cd76251..67ce38f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -33,6 +33,7 @@ endif(STIM_EMBED_VERSION) option(STIM_LEBEDEV "Lebedev ICE model" ON) option(STIM_MYLAKE "MyLake ICE model" ON) option(STIM_WINTON "Winton ICE model" ON) +option(STIM_FLATO "Flato ICE model" ON) include(src/CMakeLists.txt) include(tests/CMakeLists.txt) diff --git a/src/drivers/getm/getm_stim.F90 b/src/drivers/getm/getm_stim.F90 index 278658c..9f8785c 100644 --- a/src/drivers/getm/getm_stim.F90 +++ b/src/drivers/getm/getm_stim.F90 @@ -214,7 +214,9 @@ end subroutine post_init_ice #ifdef _TESTING_ subroutine do_ice(dz,dt,Tw,S,Ta) #else - subroutine do_ice(dz,dt,Tw,S,Ta,precip,Qsw,Qfluxes) + subroutine do_ice(dz,dt,Tw,S,Ta,precip,Qsw,Qfluxes,julianday,secondsofday,longitude, & + latitude,I_0,airt,airp,hum,u10,v10,cloud,rho,rho_0,longwave_radiation, & + hum_method,fluxes_method,albedo,heat) #endif ! ! !DESCRIPTION: @@ -224,9 +226,13 @@ subroutine do_ice(dz,dt,Tw,S,Ta,precip,Qsw,Qfluxes) ! ! !INPUT PARAMETERS: #ifdef _TESTING_ - REALTYPE, intent(in) :: dz(:,:),dt,Ta(:,:),S(:,:) + REALTYPE, intent(in) :: dz(:,:),dt,S(:,:),Ta(:,:) #else REALTYPE, intent(in) :: dz(:,:),dt,Ta(:,:),S(:,:),precip(:,:),Qsw(:,:) + REALTYPE, intent(in) :: longitude,latitude,I_0,airt,airp,hum,u10,v10,cloud,rho,rho_0,albedo,heat + integer, intent(in) :: julianday,secondsofday + integer, intent(in) :: longwave_radiation,hum_method,fluxes_method + #endif ! ! !INPUT/OUTPUT PARAMETERS: diff --git a/src/drivers/gotm/gotm_stim.F90 b/src/drivers/gotm/gotm_stim.F90 index ac99460..63ac95d 100644 --- a/src/drivers/gotm/gotm_stim.F90 +++ b/src/drivers/gotm/gotm_stim.F90 @@ -72,6 +72,8 @@ subroutine init_stim_yaml() ! See log for the ice module ! ! !LOCAL VARIABLES: + integer :: k,rc + class (type_gotm_settings), pointer :: branch !EOP !----------------------------------------------------------------------- @@ -84,16 +86,111 @@ subroutine init_stim_yaml() (/option(0, 'none', 'none'), & option(1, 'Lebedev (1938)', 'Lebedev'), & option(2, 'MyLake', 'MyLake'), & - option(3, 'Winton', 'Winton')/)) + option(3, 'Winton', 'Winton'), & + option(4,'Flato', 'Flato')/)) + call branch%get(Hice, 'H', 'initial ice thickness', 'm',default=0._rk) call branch%get(ocean_ice_flux, 'ocean_ice_flux', & 'ocean->ice heat flux','W/m^2',default=0._rk, display=display_hidden) + + !flato + call branch%get(nilay, 'nilay', 'number of ice layers', '', default=0) + call branch%get(sfall_method, 'sfall_method', 'define how snow fall is determined ','', default=0) + call branch%get(const_sfall,'const_sfall ', 'constant snow fall rate', 'm d^-1', default=0._rk) + call branch%get(dfact,'dfact', 'drift factor', '', default=0._rk) + call branch%get(depmix ,'depmix', 'prescribed mixed layer depth', '', default=0._rk) + call branch%get(sice_method, 'sice_method','define how sea-ice salinity is to be calculated','', default=0) + call branch%get(const_Sice, 'const_Sice', 'prescribed sea ice salinity', 'ppt', default=0._rk) + call branch%get(snow_dist,'snow_dist ','logical switch between uniform and Weibull-distributed snow', default=.false.) + call branch%get(distr_type,'distr_type', 'integer to chose the type of distribution', '', default=-1) + call branch%get(meltpond,'meltpond', 'If true meltponds are included If false only bare ice is included', default=.false.) + call branch%get(Ameltmax,'Ameltmax ', 'Maximum meltpond area fraction allowed', '', default=0._rk) + call branch%get(drainrate,'drainrate', 'Melt pond drainage rate in ','m/d', default=0._rk) + call branch%get(hh0, 'hh0', 'initial thickness for S calculation','', default=0._rk) + call branch%get(ice_hi_i ,'ice_hi_i ', 'initial ice thickness', '', default=0._rk) + call branch%get(ice_hs_i ,'ice_hs_i ',' initial snow thickness', '', default=0._rk) + call branch%get(albice_method,'albice_method ', 'define how ice albedo is determined ', '', default=0) + call branch%get(albice_f, 'albice_f', 'freezing ice albedo','', default=0._rk) + call branch%get(albmelt,'albmelt', 'melt pond albedo','', default=0._rk) + call branch%get(albsnow_f,'albsnow_f', 'freezing snow albedo','', default=0._rk) + call branch%get(albice_m ,'albice_m ', 'melting ice albedo','', default=0._rk) + call branch%get(albsnow_m, 'albsnow_m', 'melting snow albedo','', default=0._rk) + call branch%get(transsf,'transsf', 'freezing snow transmission coefficient','', default=0._rk) + call branch%get(transsm ,'transsm ', 'melting snow transmission coefficient','', default=0._rk) + call branch%get(transif, 'transif ', 'freezing ice transmission coefficient','', default=0._rk) + call branch%get(transim, 'transim ', 'melting ice transmission coefficient','', default=0._rk) + call branch%get(transm,'transm', 'melt pond transmision coefficient','', default=0._rk) + call branch%get(swkappasm,'swkappasm', 'melting snow extinction coefficient','', default=0._rk) + call branch%get(swkappasf,'swkappasf ', 'freezing snow extinction coefficient','', default=0._rk) + call branch%get(swkappaim,'swkappaim ', 'melting ice extinction coefficient','', default=0._rk) + call branch%get(swkappaif, 'swkappaif', 'freezing ice extinction coefficient','', default=0._rk) + LEVEL2 'done.' allocate(Tice(2)) + +#ifdef STIM_FLATO + allocate(ice_uvic_Tice(nilay+1),stat=rc) + if (rc /= 0) STOP 'init_ice: Error allocating (ice_uvic_Tice)' + ice_uvic_Tice=0 + do k=1,nilay+1 + ice_uvic_Tice(k)=245.+(Tfreezi-245.)*float(k-1)/float(nilay) + enddo + allocate(ice_uvic_Cond(nilay),stat=rc) + if (rc /= 0) STOP 'init_ice: Error allocating (ice_uvic_Cond)' + ice_uvic_Cond =0 + allocate(ice_uvic_rhoCp(nilay),stat=rc) + if (rc /= 0) STOP 'init_ice: Error allocating (ice_uvic_rhoCp)' + ice_uvic_rhoCp =0 + allocate(ice_uvic_Sint(nilay+1),stat=rc) + if (rc /= 0) STOP 'init_ice: Error allocating (ice_uvic_Sint)' + ice_uvic_Sint =0 + allocate(ice_uvic_dzi(nilay),stat=rc) + if (rc /= 0) STOP 'init_ice: Error allocating (ice_uvic_dzi)' + ice_uvic_dzi =0 + allocate(ice_uvic_zi(nilay+1),stat=rc) + if (rc /= 0) STOP 'init_ice: Error allocating (ice_uvic_zi)' + ice_uvic_zi =0 + allocate(ice_uvic_Told(nilay+1),stat=rc) + if (rc /= 0) STOP 'init_ice: Error allocating (ice_uvic_Told)' + ice_uvic_Told =0 + allocate(ice_uvic_Pari(nilay+1),stat=rc) + if (rc /= 0) STOP 'init_ice: Error allocating (ice_uvic_Pari)' + ice_uvic_Pari =0 + allocate(ice_uvic_dum(nilay+1),stat=rc) + if (rc /= 0) STOP 'init_ice: Error allocating (ice_uvic_dum)' + allocate(ice_uvic_dzice(nilay),stat=rc) + if (rc /= 0) STOP 'init_ice: Error allocating (ice_uvic_dzi)' + do k=1,nilay + ice_uvic_dzice(k)=float(k) + enddo + allocate(ice_uvic_zice(nilay+1),stat=rc) + if (rc /= 0) STOP 'init_ice: Error allocating (ice_uvic_zice)' + do k=1,nilay+1 + ice_uvic_zice(k)=float(k) + enddo + + ice_uvic_dum =0 + hsnow=0;hice=0;ice_uvic_hm=0 + ice_uvic_ts=273.16D+00;ice_uvic_tb=273.16D+00;ice_uvic_Fh=0 + ice_uvic_swr_0=0;ice_uvic_precip_i=0;ice_uvic_sfall_i=0 + ice_uvic_parb=0;ice_uvic_parui=0; + ice_uvic_Ff=0;ice_uvic_Fs=0 + ice_uvic_Sicebulk=6.0D+00 + ice_uvic_topmelt=0;ice_uvic_botmelt=0 + ice_uvic_termelt=0;ice_uvic_topgrowth=0 + ice_uvic_botgrowth=0;ice_uvic_Hmix=0 + ice_uvic_Aice=0;ice_uvic_Asnow=0 + ice_uvic_Amelt=0 + +#endif return end subroutine init_stim_yaml !EOC + + + + !----------------------------------------------------------------------- !BOP ! @@ -116,6 +213,7 @@ subroutine post_init_stim(Ta,S) !----------------------------------------------------------------------- !BOC LEVEL1 'post_init_stim' + print *, 'ice_model: ', ice_model if(Hice .gt. _ZERO_ .and. ice_model /= 0) then ice_cover=2 @@ -146,6 +244,19 @@ subroutine post_init_stim(Ta,S) LEVEL0 ".... in STIM" stop 'post_init_stim(): init_stim_winton()' #endif +#endif +#ifdef STIM_FLATO + case(4) +#if 1 +!KB allocate(Tice(2)) + call init_stim_flato() +#else + LEVEL0 "Flato model is compiled - but execution is disabled" + LEVEL0 "change line 138 in gotm_stim.F90 - then recompile - " + LEVEL0 "then do some work to make the Flato ice model work ...." + LEVEL0 ".... in STIM" + stop 'post_init_stim(): init_stim_flato()' +#endif #endif case default stop 'invalid ice model' @@ -162,17 +273,21 @@ end subroutine post_init_stim !BOP ! ! !IROUTINE: do the ice calculations -! + ! !INTERFACE: - subroutine do_stim(dz,dt,Tw,S,Ta,precip,Qsw,Qfluxes) -! + subroutine do_stim(dz,dt,Tw,S,Ta,precip,Qsw,Qfluxes,Qfluxes_uvic,julianday,secondsofday, & + I_0,airt,rho,rho_0,albedo,heat) + ! !DESCRIPTION: ! ! !USES: IMPLICIT NONE ! ! !INPUT PARAMETERS: - REALTYPE, intent(in) :: dz,dt,Ta,S,precip,Qsw + + REALTYPE, intent(in) :: dz,dt,S,Qsw,airt,rho,rho_0 + REALTYPE, intent(inout) :: Ta,precip,I_0,albedo,heat + integer, intent(in) :: julianday,secondsofday ! ! !INPUT/OUTPUT PARAMETERS: REALTYPE, intent(inout) :: Tw @@ -183,6 +298,13 @@ subroutine Qfluxes(T,qh,qe,qb) REALTYPE, intent(out) :: qh,qe,qb end subroutine end interface +!jpnote +interface + subroutine Qfluxes_uvic(T,qh,qe,qb) + REALTYPE, intent(in) :: T + REALTYPE, intent(out) :: qh,qe,qb + end subroutine +end interface ! ! !REVISION HISTORY: ! Original author(s): Karsten Bolding @@ -191,6 +313,7 @@ subroutine Qfluxes(T,qh,qe,qb) ! ! !LOCAL VARIABLES: REALTYPE :: Tf + !EOP !----------------------------------------------------------------------- !BOC @@ -224,6 +347,39 @@ subroutine Qfluxes(T,qh,qe,qb) else call do_stim_winton(ice_cover,dz,dt,Tw,S,Ta,precip,Qsw,Qfluxes) end if +#endif +#ifdef STIM_FLATO + case(4) + + if (S .lt. 0.01) then + LEVEL0 'The FLato ice model is developed for oceanic conditions.' + LEVEL0 'Very low salinity is not supported - and the principle' + LEVEL0 'advantage of the model (brine contribution to latent' + LEVEL0 'heat calculation) is not met.' + LEVEL0 'Please select another ice model.' + stop 'do_stim()' + else + call do_ice_uvic(dt,dz,julianday,secondsofday, & + I_0,airt,precip,Tw,S,rho,rho_0, & + ice_hi,ice_hs,ice_uvic_hm,ice_uvic_Tice, & + ice_uvic_Cond,ice_uvic_rhoCp, & + ice_uvic_Sint,ice_uvic_dzi,ice_uvic_zi, & + ice_uvic_Pari,ice_uvic_Told,albedo,heat,& + ice_uvic_Fh,ice_uvic_Ff,ice_uvic_Fs,& + ice_uvic_Sicebulk, ice_uvic_topmelt, & + ice_uvic_botmelt,ice_uvic_termelt, & + ice_uvic_topgrowth,ice_uvic_botgrowth,& + ice_uvic_Hmix,ice_uvic_Aice,ice_uvic_Asnow,& + ice_uvic_Amelt,ice_uvic_swr_0,ice_uvic_precip_i,ice_uvic_sfall_i,Qfluxes_uvic) + + ice_uvic_ts=ice_uvic_Tice(1) + ice_uvic_tb=ice_uvic_Tice(nilay) + ice_uvic_parb=ice_uvic_Pari(nilay) + ice_uvic_parui=ice_uvic_Pari(nilay+1) + + end if + + #endif case default stop 'invalid ice model' diff --git a/src/models/CMakeLists.txt b/src/models/CMakeLists.txt index bb5ac1a..cfbc8d2 100644 --- a/src/models/CMakeLists.txt +++ b/src/models/CMakeLists.txt @@ -9,6 +9,7 @@ set_property(TARGET stim PROPERTY Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINAR include(${CMAKE_CURRENT_LIST_DIR}/lebedev/CMakeLists.txt) include(${CMAKE_CURRENT_LIST_DIR}/mylake/CMakeLists.txt) include(${CMAKE_CURRENT_LIST_DIR}/winton/CMakeLists.txt) +include(${CMAKE_CURRENT_LIST_DIR}/flato/CMakeLists.txt) target_include_directories(stim PUBLIC $ $ diff --git a/src/models/flato/CMakeLists.txt b/src/models/flato/CMakeLists.txt new file mode 100644 index 0000000..a745d47 --- /dev/null +++ b/src/models/flato/CMakeLists.txt @@ -0,0 +1,7 @@ +if(STIM_FLATO) + target_sources(stim + PRIVATE + ${CMAKE_CURRENT_LIST_DIR}/stim_flato.F90 + ) + target_compile_definitions(stim PUBLIC STIM_FLATO) +endif(STIM_FLATO) diff --git a/src/models/flato/README.md b/src/models/flato/README.md new file mode 100644 index 0000000..6fe87e3 --- /dev/null +++ b/src/models/flato/README.md @@ -0,0 +1,6 @@ +# Flato ice model + +This is an ice model based on the 1-D thermodynamic sea ice model by Greg Flato (Flato and Brown, 1996) and modified to fit the GOTM structure, as well as containing the snow distribution and meltponds subgrid parametrization (Abraham et al. 2014). + +It performs a surface energy budget calculation to get net flux at ice, snow or open water surface, solves the 1-D heat conduction problem using an implicit finite difference scheme with the ice/snow slab discretized into an arbitrary number of thickness layers (presently configured for only one snow layer). A Newton-Raphson iterative scheme is used to solve for the surface temperature which appears non-linearly in both surface and conductive fluxes. Surface melt is calculated to use up the net incoming energy flux at the surface with the surface temperature fixed at the melting point. Additional melt may occur to restore temperatures to the freezing point where it is exceeded. Under cooling conditions, the upper surface has a flux boundary condition with surface growth occurring only if snow load causes surface flooding. Growth at the ice underside is calculated so as to balance the flux at the bottom where the temperature remains at the freezing point of sea water. If ice melts completely, incoming heat is stored in the mixed layer which must cool to the freezing point before new ice can form. Initial growth of very thin ice is calculated assuming a linear temperature profile. + diff --git a/src/models/flato/stim_flato.F90 b/src/models/flato/stim_flato.F90 new file mode 100644 index 0000000..7e79d67 --- /dev/null +++ b/src/models/flato/stim_flato.F90 @@ -0,0 +1,2736 @@ +!----------------------------------------------------------------------- +!BOP +! +! MODULE: stim_flato --- Flato&Brown thermodynamic ice model +! \label{sec:stim_flato} +! +! INTERFACE +module stim_flato +! +! DESCRIPTION: +!! This module is based on the 1-D thermodynamic sea ice model +!! (ONE_D_THERMOv1.0) by Greg Flato (Flato and Brown, 1996. JGR, 101(C10)) +!! and modified to fit the GOTM structure. +!! It performs a surface energy budget calculation to get net flux at ice, +!! snow or open water surface, solves the 1-D heat conduction problem using +!! an implicit finite difference scheme with the ice/snow slab discretized +!! into an arbitrary number of thickness layers (presently configured for only +!! one snow layer). A Newton-Raphson iterative scheme is used to solve for the +!! surface temperature which appears non-linearly in both surface and +!! conductive fluxes. Surface melt is calculated to use up the net incoming +!! energy flux at the surface with the surface temperature fixed at the +!! melting point. Additional melt may occur to restore temperatures to the +!! freezing point where it is exceeded. Under cooling conditions, the upper +!! surface has a flux boundary condition with surface growth occuring only +!! if snowload causes surface flooding. Growth at the ice underside is +!! calculated so as to balance the flux at the bottom where the temperature +!! remains at the freezing point of sea water. If ice melts completely, +!! incoming heat is stored in the mixed layer {\sl fixed MLD in original code} +!! which must cool to the freezing point before new ice can form. Initial +!! growth of very thin ice is calculated assuming a linear temperature profile. +! +!! Further documentation is provided by comments in the code and headers +!! in each subroutine. +!!----------------------------------------------------------------------- +!! original code written by Gregory M. Flato, Canadian Centre for Climate +!! Modelling and Analysis, Environment Canada +!! +!! ver. 1.05 - G.Flato - 29/Oct/99 - last change before transfer to +!! implement module in GOTM +!! ver. 2.0 - N.Steiner - Nov/08 - conversion to subroutine as +!! part of GOTM +!! ver. 3.0 - N. Steiner - Jul/14 - recodng for gotm git structure +!! +!! ver. 3.1 - N. Steiner - Nov/14 - include saltflux from +!! Vancoppenolle et al 2009 - in prog. +!! ver. 3.2 - N. Steiner - Nov/14 - correct heat transfer for open water +!! ( remove double calc inside and outside the ice module, +!! recorrect I_0 for ice case +!! ver. 3.3 - N. Steiner - Dec/14 - Implement Abraham et al. 2014 +!! radiative transfers +!! - include functions for extinction +!! and transmissivity +!! +!! +! !USES: + use stim_variables + + IMPLICIT NONE +! +! default: all is private. + private +! +! !PUBLIC MEMBER FUNCTIONS: + public :: init_stim_flato + public :: do_ice_uvic +! +! +! LOCAL VARIABLES: + + real(rk), public :: rhosnow + !! snow density (kg m-3) + +!----fixed size variables +! + real(rk), dimension(2), public :: Iceflux + !! a 2-element array of time-step averaged boundary fluxes + !! returned from heat conduction scheme, defined as positive + !! downward in accordance with vertical coordinate sign + !! convention (W m-2) + real(rk), dimension(2) :: bctype + !! 2 element array defining upper and lower boundary + !! condition type. bctype(1) refers to upper boundary, + !! bctype(2) refers to lower boundary. Negative value + !! indicates temperature boundary condition, positive + !! value indicates flux boundary condition. + real(rk), dimension(2) :: bcs + !! 2 element array containing boundary condition values. + !! bcs(1) refers to upper boundary, bcs(2) refers to lower + !! boundary. If temperature boundary condition, units: (K), + !! if flux boundary condition, units: (W m-2) + real(rk) :: dti + !! timestep in the ice model (usually set to ocean time step) + +! +!----fluxes + real(rk), public :: qb + !! long wave back radiation (in-out) (W m-2) + real(rk) :: qe + !! latent heat flux into ice (W m-2) + real(rk) :: qh + !! sensible heat flux into ice (W m-2) + real(rk) :: tx,ty + !! surface stress components in x and y direction (Pa)!check + real(rk), public :: PenSW + !! short wave radiation that penetrates the surface (W m-2) + real(rk) :: fluxt + !! net flux at surface of ice/snow slab calculated from + !! surface energy budget and used as upper boundary condition + !! in heat conduction solution. NOTE: this does not include + !! 'PenSW' which is taken as a distributed source within the + !! slab (W m-2) + real(rk), public :: simass + !! ice mass per unit area (kg m-2) + real(rk), public :: snmass + !! snow mass per unit area (kg m-2) + real(rk) :: simasso + !! ice mass per unit area at previous timestep(kg m-2) + real(rk) :: snmasso + !! ice mass per unit area at previous timestep (kg m-2) !snow mass? + real(rk) :: Ts + !! upper surface temperature (K) + real(rk) :: Tsav + !! average snow layer temperature (K) + + +! Ice salinity + logical, public :: ice_salt=.false. + !! logical variable to turn on/off the salt profile scheme + +! Atmospheric forcing + real(rk), public :: sfall + !! snow fall rate (m s-1) + real(rk) :: airtk + !! surface air temperature (K) + real(rk), dimension(5,3) :: C + !! coefficients for + real(rk), dimension(5,3) :: R + !! coefficients for RHS vector + real(rk) :: dto + !! time step from GOTM (dt) + integer, public :: nslay + !! nslay number of snow layers + +!Snow_dist: Snowdistribution variables: + real(rk), public :: Asnow + !! Area which is covered with snow + real(rk), public :: Aice + !! Area which is covered with ice + real(rk), public :: Amelt + !! Area which is covered with melt pond + real(rk) :: hsmax = hsmin + !! Maximal height of snow for the calculations of Weibull-distributed snow + real(rk) :: albice + !! Albedo of ice + real(rk):: albsnow + !! Albedo of snow + real(rk), public :: meltmass + !! melt pond mass per unit area (kg m-2) + real(rk) :: meltmasso + !! melt pond mass per unit area at previous timestep(kg m-2) + real(rk), parameter :: pi=4.D+00*atan(1.D+00) + !! Pi + !NSnote read from gotm? + +! !REVISION HISTORY: +! Original author: Michael Winton +! Author(s): Adolf Stips, Jesper Larsen and Karsten Bolding +! +!----------------------------------------------------------------------- + + contains + +!----------------------------------------------------------------------- +!BOP +! +!IROUTINE: Initialise the uvic--ice model \label{sec:init-uvic-ice} +! !INTERFACE: +!KB subroutine init_stim_flato(ice_cover,dz,dt,Tw,S,Ta,precip) + subroutine init_stim_flato() +!! Initialise the uvic-ice model +! +! !DESCRIPTION: +! +! !USES: + IMPLICIT NONE +! +! !LOCAL VARIABLES: +! + integer :: k,rc +! !LOCAL PARAMETERS: + +!EOP +!----------------------------------------------------------------------- +!BOC +! +! +!------------------------------------------------------------------------------------- +! copy paste from init_ice_uvic FROM ice_uvic.F90 +!------------------------------------------------------------------------------------- + +! !DESCRIPTION: !copy paste from ice_uvic.f90 without the namelist initialization +! +! !USES: + ! IMPLICIT NONE +! +! !INPUT PARAMETERS: + +! !REVISION HISTORY: +!! Original author(s): Karsten Bolding, Nadja Steiner based on original model from Greg Flato +! + + ! check if the snow distribution is set + if(snow_dist .and. distr_type .eq. -1) then + print*,'No specific snow distribution was chosen' + stop + end if + +! The snow fall rate + select case (sfall_method) + case (1) + ! LEVEL2 'Using snow fall rate= ', const_sfall + sfall=const_sfall/86400.D+00 !convert to [m/s] + case (2) + sfall = 0.00D+00 ! initialize snowfall from precipitation to 0.0 + case default + end select + + snmass=ice_hs_i*rhoscold + simass=ice_hi_i*rhoice + meltmass=0.0 + drainrate=drainrate*rhowaterfresh/86400.D+00 + !conversion from [m/d] to [kg/m^2/s] + +! The sea ice salinity calculation + select case (sice_method) + case (1) +! LEVEL2 'Using constant bulk salinity',const_Sice + case (2) +! LEVEL2 'Using simple ice salinity calculation (Vancoppenolle 2009)' + ice_salt=.true. + case default + end select +! check if the snow distribution is set snow_dist + if(snow_dist .and. distr_type .eq. -1) then + !LEVEL2 'No specific snow distribution was chosen' + stop + end if +!return + + + return + + end subroutine init_stim_flato +!EOC + +!----------------------------------------------------------------------- + +!BOP +! +! ROUTINE: Calculate ice thermodynamics \label{sec:do_ice_uvic} +! +! !INTERFACE: + subroutine do_ice_uvic(dto,h,julianday,secondsofday, & + I_0,airt,precip,TSS,SSS,rhowater,rho_0, & + ice_hi,ice_hs,ice_hm,Tice,Cond,rhoCp,Sint,dzi,zi, & + Pari,Told,alb,heat,Fh,Ff,Fs,Sice_bulk,TopMelt,BotMelt,& + TerMelt,TopGrowth,BotGrowth,Hmix,Aice_i,Asnow_i,Amelt_i,& + swr_0,precip_i,sfall_i,Qfluxes_uvic) +! DESCRIPTION: +!! Update the sea-ice conditions: if no ice, call open_water and determine +!! new ice growth, if ice exists calculate growth or melt. The model has +!! nilay ice layers and one snow layer. +!! This subroutine updates the sea ice prognostic variables. The updated +!! variables are ice_hi, ice_hs,ice_hm, Tice(nilay+1), Cond(nilay),rhoCp(nilay), +!! Sint(nilay+1), dzi(nilay), zi(nilay+1), Told(nilay+1),ts,alb, heat +!! melt and growth terms per timestep are transferred out for potential use in +!! ice algae model: TopMelt,BotMelt,TerMelt,TopGrowth,BotGrowth, note these are +!! not accumulated over the run as in the original code from GF +! REVISION HISTORY: +!! Original author(s): Nadja Steiner, modified from ice_thermo by Greg Flato +!! based on GOTM template + +! !USES: + implicit none +! +! !INPUT PARAMETERS: + real(rk), intent(in) :: dto + !! ocean timestep (sec) + real(rk), intent(in) :: h + !! sea surface layer thickness + integer, intent(in) :: julianday + !! this julian day + integer, intent(in) :: secondsofday + !! seconds for this day + real(rk), intent(inout) :: I_0 + !! shortwave radiation at sea surface + real(rk), intent(in) :: airt + !! 2m temperature + real(rk), intent(inout) :: precip + !! freshwater precipitation (m/s) + ! real(rk), intent(in) :: cloud ! cloud cover + real(rk), intent(inout) :: TSS + !! sea surface temperature + real(rk), intent(in) :: SSS + !! sea surface salinity + real(rk), intent(in) :: rhowater + !! sea surface layer density + real(rk), intent(in) :: rho_0 + !! reference density + + ! !INPUT/OUTPUT PARAMETERS: + real(rk), intent(inout) :: ice_hi + !! ice thickness (m) + real(rk), intent(inout) :: ice_hs + !! snow thickness (m) + real(rk), intent(inout) :: ice_hm + !! meltpond thickness (m) + real(rk), intent(inout) :: Tice(nilay+1) + !! ice layer temperature Tice(nilay +1)(deg-C) + real(rk), intent(inout) :: Cond(nilay) + !! thermal conductivities defined at the centre of each layer Cond(nilay)(W m-1 K-1) + real(rk), intent(inout) :: rhoCp(nilay) + !! volumetric heat capacities defined at the centre of each layer rhoCp(nilay)(J m-3 K-1) + real(rk), intent(inout) :: Sint(nilay+1) + !! internal heat source due to penetrating short wave radiation Sint(nilay)(W m-3) + real(rk), intent(inout) :: dzi(nilay) + !! layer thicknesses dzi(nilay)(m) + real(rk), intent(inout) :: zi(nlmax) + !! layer interface depths zi(nilay+1)(m) + real(rk), intent(inout) :: Told(nilay+1) + !! ice temperature two time steps previous to calculation + !! of outgoing long-wave flux in SEBUDGET Told (nilay+1) (K) + real(rk), intent(inout) :: alb + !! surface albedo - water, ice/snow-total + real(rk), intent(inout) :: heat + !! surface heat flux + +!NSnote, check - maybe adjust units... + real(rk), intent(inout) :: Sice_bulk + !! bulk ice salinity (ppt) + real(rk), intent(inout) :: Fh + !! interface heat flux (W/m2) + real(rk), intent(out) :: Ff + !! interface freshwater flux (m s-1) + real(rk), intent(out) :: Fs + !! interface salt flux - (ppt m s-1) + real(rk), intent(out) :: Pari(nilay+1) + !! photosynthatically available radiation in ice (W m-2) + +! !OUTPUT PARAMETERS: + real(rk), intent(out) :: TopMelt + !! top melting - ice mass melted at the surface (snow+ice) at time step(m) + real(rk), intent(out) :: Botmelt + !! bottom melting - ice mass melted at the ice bottom at time step(m) + real(rk), intent(out) :: TerMelt + !! internal melting - ice mass melted in the ice interior at time step (m) + real(rk), intent(out) :: TopGrowth + !! top growth ice mass growth at slab surface due to snow submersion (m) + real(rk), intent(out) :: BotGrowth + !! bottom growth - ice mass growth at the ice bottom at time step (m) + real(rk), intent(out) :: Hmix + !! transferred energy - check (m) +! Hmix - mixed layer heat storage (J m-2) =======> accounts only for +! the SWR which crosses the ice slab and reach the water. keep it for now + real(rk), intent(out) :: Aice_i,Asnow_i,Amelt_i + !! ice area fraction which is : open ice, snow and + !! meltpond, respectively + real(rk), intent(out) :: swr_0,precip_i,sfall_i !H! + !! incidental SWR,snowfall + + interface + subroutine Qfluxes_uvic(T,qh,qe,qb) + integer, parameter :: rk = kind(1.d0) + real(rk), intent(in) :: T + real(rk), intent(out) :: qh,qe,qb + end subroutine + end interface + +! LOCAL VARIABLES + real(rk) :: dmsi ! dmsi - new ice formation at open water [kg m-2] + integer :: yy,mm,dd,j,k + real(rk) :: p_tmp ! p_tmp ! for rain to snow conversion +!KB + real(rk) :: evap ! evaporation of ice (m/s) + real(rk) :: ohflux !heat flux from ocean into ice underside (W m-2) +!NSnote, not sure what evap is for +! +!BOC +! LEVEL0 'do_ice_uvic' +! Calculate seawater freezing temperature +! Tfreezi = (-0.0575D00*SSS)+kelvin +!KB STDERR T,S +! +! set ice timestep to ocean timestep + + dti=dto +! initialize new ice formation in open water, melt/growth for this timestep + dmsi=0 + TopMelt = 0 + Botmelt = 0 + TerMelt = 0 + TopGrowth = 0 + BotGrowth = 0 + Ff= 0 + Fs= 0 +! degC to K unit conversion + +! jpnote: add option for if kelvin is given instead of celcius + if (airt .lt. 100.) then + airtk = airt + kelvin ! airt is given in celcius, need to convert to kelvin + else + airtk = airt !airt is given in kelvin + end if + + !airtk=airt+kelvin + + if(sfall_method .eq. 2) then +!get snowfall from precipitation + if((airtk-Tmelts).lt.(-5.0)) then + p_tmp = 1.0 + elseif((airtk-Tmelts).ge.(-5.0).and.(airtk-Tmelts).le.(5.0)) then + p_tmp = 1.0-(airtk-Tmelts+5.0)*0.1 + else + p_tmp = 0.0 + endif + ! - add factor for drifting snow... + p_tmp=p_tmp*dfact + sfall=precip*p_tmp*rhowaterfresh/rhoscold ! m/s + endif +! set old snow and ice mass and meltpond mass + +!nsnote, are all the mass variables known from the last timestep? if not, set here...??? + snmasso=snmass + simasso=simass + meltmasso=meltmass + +! add snowfall to existing snow if ice is present (no accumulation over open water!) + +! and if air temp. is below freezing. + if (simass.gt.rhoice*hsmin .and. airtk.lt. kelvin) then + snmass=snmass+rhoscold*sfall*dti +! Remove snow falling on ice from precip + precip_i=precip !H! Save precip before set to zero. + precip=precip-sfall*rhoscold/rhowaterfresh + precip =max(0.0D0,precip) + !NSnote, make sure no double counting of precip + ! if no ice then precip (or if precip not given, water equivalent for + ! snow fall) should go in the ocean, if ice then sfall, but no + ! additional precip in the ocean !!! + endif + +! calculate snow density - depends on mean snow layer temperature + + Tsav=(Tice(1)+Tice(2))/2.D+00 + if(Tsav.lt.Tmelti) then + rhosnow=rhoscold + else + rhosnow=rhoswarm + endif + + + swr_0 = I_0/(1-alb) !H! save incidental swr. +! if ice is present, calculate growth or melt, if not calculate potential new growth +! + if(simass.eq.0.D+00) then + Hmix=Hmix+(heat+I_0)*dti + if(depmix.eq.0.0) depmix=h + call open_water(nilay,I_0,Sice_bulk,Hmix,Tice,depmix,TSS,Fh,heat,precip,precip_i) + else +! reset I_0 so it can be recalculated for ice + I_0 = I_0/(1-alb) ! Here alb is water albedo from GOTM + +! do surface energy budget and heat conduction calculations via N-R +! iteration on surface temperature +! This will update the ice temperature at each layer + + call nr_iterate(nilay,I_0,Told,Tice,Pari,& + Sice_bulk,ice_hi,ice_hs,dzi,Cond,rhoCp,zi, & + Sint,evap,alb,Qfluxes_uvic) + +! construct depth array by adding up dzi's + zi(1)= 0 + do k=2,nilay+1 + zi(k)=zi(k-1)+dzi(k-1) + enddo + +! add short wave radiation that penetrates snow/ice slab to mixed layer +!nsnote, check if ok that if snow_dist can be removed and calculation done with PAR for all cases +!if (snow_dist) then + Hmix=Hmix+Pari(size(Pari))*dti +!else +! Hmix=Hmix+PenSW*exp(-swkappa*zi(nilay+1))*dti +!endif + +!add mixed layer heat storage to oceanic heat flux - + + ohflux=Hmix/dti + Hmix=0.D+00 + +!...Calculate sea-ice turbulent heat flux and input it to the ohflux + Fh=ohflux !NSnote check + if(ice_hi .gt. 0) then +! heat=-Fh !NSnote, don't think so, if, ohflux is used for ice growth +! heat=0.D00 + end if ! NSnote: What if no ice?, all melted..what heat then? +! I_0=PenSW + I_0=Pari(nilay+1) !NScheck + +!...........calculate growth or melt amount at top and bottom surface, +!...........converting snow to ice if ice surface is submerged + +! initialize top, interior and bottom melt and growth +! + + + + call growthtb(rhowater,nilay,rhoCp,dzi,Tice,TopGrowth,TerMelt,TopMelt,& + BotGrowth,BotMelt,Hmix,ohflux) +!NSnote, Hmix contains extra heat left from melting, +!this should be transferred back into ocean to heat the water + heat=Hmix/dti + + if (ice_salt) & + call saltice_prof_simple(dti,nilay,SSS,simasso,snmasso,meltmasso, & + rhoice,rhosnow,rhowater,rhowaterfresh,zi,Tice,Sice_bulk,Ff,Fs, & + TopGrowth,BotGrowth,TopMelt,BotMelt,TerMelt) + + + endif + +!...Update ice and snow thickness from mass and density so accurate for end of timestep +! + ice_hs=snmass/rhosnow + ice_hi=simass/rhoice + ice_hm=meltmass/rhowaterfresh + +! Change units from m (per timestep) to m/s for potential transfer to FABM +! (ice algae, since FABM is not aware of the timestep) + + TopMelt=TopMelt/dti + BotMelt=BotMelt/dti + TerMelt=TerMelt/dti + TopGrowth=TopGrowth/dti + BotGrowth=BotGrowth/dti + +! set area fractions of ice so they can be written out + Aice_i=Aice + Asnow_i=Asnow + Amelt_i=Amelt + sfall_i=sfall + + return + +!#endif + + end subroutine do_ice_uvic +!EOC +!----------------------------------------------------------------------- +!----------------------------------------------------------------------- + ! put any internal subroutines below here +!----------------------------------------------------------------------- +!----------------------------------------------------------------------- +!BOP +! !IROUTINE: +! +! !INTERFACE: + subroutine therm1d(nilay,Sice_bulk,ice_hi,ice_hs,dzi,Cond,rhoCp,zi,Sint,Pari,Tice,I_0) + +! !DESCRIPTION: +! +!! Subroutine to calculate heat conduction through a 1-dimensional slab +!! of sea ice of a given thickness for specified boundary conditions. +!! Parameterization of conductivity and heat capacity follows closely +!! that of Ebert and Curry (J. Geophys. Res.,98:10085-10109, 1993) and +!! uses a generalized Crank-Nicholson or implicit algorithm to solve the +!! 1-dimensional heat conduction problem. + +! !USES: + + IMPLICIT NONE + +! !INPUT PARAMETERS: + + ! !INPUT PARAMETERS: + integer, intent(in) :: nilay + real(rk), intent(in) :: I_0 !incoming swr for snow_dist + +! !OUTPUT PARAMETERS: + real(rk),intent(out) :: ice_hi,ice_hs + real(rk),intent(out) :: dzi(nilay),zi(nilay+1) + real(rk),intent(out) :: Sint(nilay+1),Cond(nilay),rhoCp(nilay) + real(rk),intent(out) :: Pari(nilay+1) + +! !INOUT PARAMETERS: + real(rk),intent(inout) :: Tice(nilay+1) + real(rk),intent(inout) :: Sice_bulk + +! +! !REVISION HISTORY: +!! adjusted to GOTM F90 Nadja Steiner Dec 2008 +!! Original author(s): Gregory M. Flato - Canadian Centre for +!! Climate Modelling and Analysis + +! +!BOC +! !LOCAL VARIABLES: +! bctype - 2 element array defining upper and lower boundary +! condition type. bctype(1) refers to upper boundary, +! bctype(2) refers to lower boundary. Negative value +! indicates temperature boundary condition, positive +! value indicates flux boundary condition. + real(rk) :: bctype(2) +! bcs - 2 element array containing boundary condition values. +! bcs(1) refers to upper boundary, bcs(2) refers to lower +! boundary. If temperature boundary condition, units: (K), +! if flux boundary condition, units: (W m-2) +! NOTE: flux is positive downward due to sign convention. + real(rk) :: bcs(2) +! hlay - ice or ice/snow thickness partitioned by ice layers (m) + real(rk) :: hlay +! Ticeav - average ice temperature (interim) + real(rk) :: Ticeav +! nslay - number of snow layers +! integer :: nslay +! Sice_bulk - ice salinity used in parameterization of heat capacity +! and conductivity (ppt) + real(rk) :: snlaythick,Tbot,thicklay,Ttop +! fluxtsplit - averaged surface flux over all 3 surface conditions + real(rk) :: fluxtsplit + integer :: k +! hzero - minimum ice thickness to avoid dividing by zero elsewhere in code + real(rk), parameter :: hzero=1.D-03 + +!EOP + !----------------------------------------------------------------------- + + + + ! LEVEL1'therm1d' + +! +!...Calculate ice and snow thickness from mass and density +!...(minimum thickness to avoid dividing by zero elsewhere in code) +! + ice_hs=max(hzero*0.1D+00,snmass/rhosnow) + ice_hi=max(hzero,simass/rhoice) + +!...snow_dist: max. snow height within a year to get the ratio of current snow +! depth to maximum snow depth for the Weibull-distributed snow +! if(snow_dist .and. snmass .eq. 0.D+00) then + if(snmass .eq. 0.D+00) then + hsmax=hsmin +! else if (snow_dist .and. snmass .gt. 0.D+00) then + else if (snmass .gt. 0.D+00) then + hsmax=max(ice_hs,hsmax) + end if +!...end snow_dist +! +!...Calculate array of layer thicknesses from ice and snow thickness +!...and calculate conductivity and heat capacity arrays. Parameterizations +!...based on Ebert and Curry (1993) +! +!...Total slab thickness +! + hlay=ice_hi+ice_hs +! +!...Check if sufficient snow to warrant snow layers +! + if(ice_hs.gt.hsmin) then +! +!.......estimate layer thickness (i.e. if evenly-spaced) + thicklay=hlay/float(nilay) +! +!.......number of snow layers +! nslay=max(1,nint(ice_hs/thicklay)) +! nslay=min(nslay,nilay-1) + +!---Temporary fix--- + +! --- ONLY ONE SNOW LAYER USED +! --- Do not change this until a conservative method for reapportioning +! --- temperature when number of snow layers changes is developed +! --- G. Flato 14/Sept/94 + + nslay=1 + +!---Temporary fix--- +! +!.......thickness of snow layers + snlaythick=ice_hs/float(nslay) +! +!.......loop over snow layers + do k=1,nslay + dzi(k)=snlaythick +!.........conductivity of snow, including parameterization of water vapor +!.........diffusion + Ticeav=(Tice(k)+Tice(k+1))/2.D+00 + Cond(k)=(2.845D-06*rhosnow**2)+(2.7D-04*2.D+00**((Ticeav-233.D+00)/5.D+00)) +!....snow_dist + if (snow_dist .and. distr_type .eq. 0 .and. ice_hs .ge. hsmax) then + Cond(k)=Cond(k)*pi/4.D+00 + end if + +!end snow_dist +!.........heat capacity of snow + rhoCp(k)=rhosnow*(92.88D+00+7.364D+00*Ticeav) + enddo +!.......reset 'hlay' to be thickness spanned by ice layers + hlay=ice_hi + else +!....if not, number of snow layers equals zero + nslay=0 + endif +! +!...Now do ice layers +! + do k=nslay+1,nilay +!.......ice layer thickness + dzi(k)=hlay/float(nilay-nslay) +!.......salinity of ice + if (ice_salt) then + Sice_bulk=Sice_bulk + else + if(ice_hi.ge.0.57D+00) then + Sice_bulk=3.2 + else + Sice_bulk=(14.24D+00-19.39D+00*ice_hi) + endif + endif + Ticeav=min(Tmelti-0.1D+00,(Tice(k)+Tice(k+1))/2.D+00) +!.......conductivity of ice + Cond(k)=Condfi+0.1172D+00*Sice_bulk/(Ticeav-Tmelti) +!.......minimum value to avoid zero or negative conductivities + Cond(k)=max(Cond(k),Condfi/5.D+00) +!.......heat capacity of ice + rhoCp(k)=rhoCpfi+1.715D+07*Sice_bulk/(Ticeav-Tmelti)**2 +!.......limit value to 100 times fresh water value to avoid +!.......extreme values when temperature nears melting + rhoCp(k)=min(rhoCpfi*100.D+00,rhoCp(k)) + enddo +! +!...Calculate internal sources due to penetrating shortwave radiation +!...Note, internal sources specified at layer interfaces +! + zi(1)=0.D+00 + zi(2)=zi(1)+dzi(1)/2.D+00 + zi(3)=zi(2)+dzi(1)/2.D+00+dzi(2)/2.D+00 + Sint(1)=0.D+00 !internal sources must be zero at + Sint(nilay+1)=0.D+00 !top and bottom surfaces + Pari(1)=PenSW !Pari at top of snow set to PenSW +!....if no ice, put all short wave radiation into surface flux +! NSnote: not sure why this is here, since this routine should not +! be called if no ice + + if(hlay.lt.hzero) then + fluxt=fluxt+PenSW + else +!....snow_dist + if (snow_dist .and. distr_type .eq. 0 .and. ice_hs .ge. hsmax) then +! absorption at surface, use PenSW from albedo_ice since Anow=1.0 + fluxt=fluxt+PenSW*swkappas(Tsav)*(1.D+00-swkappas(Tsav)*ice_hs *erfc(swkappas(Tsav)& + *ice_hs/sqrt(pi))*exp((swkappas(Tsav)*ice_hs)**2.D+00/pi)) + Pari(2)=PenSW*(1-swkappas(Tsav)*ice_hs*erfc(swkappas(Tsav)*ice_hs/sqrt(pi))& + *exp((swkappas(Tsav)*ice_hs)**2.D+00/pi)) + else if(snow_dist .and. distr_type .eq. 0 .and. ice_hs .gt. 0.D+00 & + .and. ice_hs .lt. hsmax .and. Asnow.gt.0.D+00) then +! absortion at surface redone here to accurately account for different surface types + fluxt=fluxt+Asnow*I_0*(1.D+00-albsnow)*transs(Tsav)*swkappas(Tsav) & + *(exp(-(inverfc(ice_hs/hsmax))**2.D+00-swkappas(Tsav)*2.D+00 & + *hsmax/sqrt(pi)*inverfc(ice_hs/hsmax))& + -swkappas(Tsav)*hsmax*erfc(swkappas(Tsav)*hsmax/sqrt(pi) & + +inverfc(ice_hs/hsmax))*exp((swkappas(Tsav)*hsmax)**2.D+00/pi)) & + +Aice*I_0*(1.D+00-albice)*transi(Tsav) + + Pari(2)=Asnow*I_0*(1.D+00-albsnow)*transs(Tsav) & + *(exp(-(inverfc(ice_hs/hsmax))**2.D+00-swkappas(Tsav)*2.D+00 & + *hsmax/sqrt(pi)*inverfc(ice_hs/hsmax))& + -swkappas(Tsav)*hsmax*erfc(swkappas(Tsav)*hsmax/sqrt(pi) & + +inverfc(ice_hs/hsmax))*exp((swkappas(Tsav)*hsmax)**2.D+00/pi)) & + +Aice*I_0*(1.D+00-albice)*transi(Tsav) + else if(snow_dist .and. distr_type .eq. 1 .and. ice_hs .ge. hsmax) then +! absorption at surface, use PenSW from albedo_ice since Anow=1.0 + fluxt=fluxt+PenSW*swkappas(Tsav)*4.D+00/ice_hs**2.D+00*1.D+00 & + /(2.D+00/ice_hs+swkappas(Tsav))**2 + Pari(2)=PenSW*4.D+00/ice_hs**2.D+00*1.D+00/(2.D+00/ice_hs+swkappas(Tsav))**2 + else if(snow_dist .and. distr_type .eq. 1 .and. ice_hs .gt. 0.D+00 .and. ice_hs .lt. hsmax) then +! absortion at surface redone, here to accurately account for different surface types + fluxt=fluxt+Asnow*I_0*(1.D+00-albsnow)*transs(Tsav)*swkappas(Tsav) & + *4.D+00/hsmax**2.D+00*1.D+00/(2.D+00/hsmax+swkappas(Tsav)) & + *exp(hsmax/2.D+00*(W(-2.D+00*ice_hs/hsmax*exp(-2.D+00))+2.D+00) & + *(2.D+00/hsmax+swkappas(Tsav)))*(hsmax/2.D+00*(-W(-2.D+00*ice_hs/hsmax & + *exp(-2.D+00))-2.D+00)+1.D+00/(2.D+00/hsmax+swkappas(Tsav))) & + + Aice*I_0*(1.D+00-albice)*transi(Tsav) + Pari(2)=Asnow*I_0*(1.D+00-albsnow)*transs(Tsav) & + *4.D+00/hsmax**2.D+00*1.D+00/(2.D+00/hsmax+swkappas(Tsav)) & + *exp(hsmax/2.D+00*(W(-2.D+00*ice_hs/hsmax*exp(-2.D+00))+2.D+00) & + *(2.D+00/hsmax+swkappas(Tsav)))*(hsmax/2.D+00*(-W(-2.D+00*ice_hs/hsmax & + *exp(-2.D+00))-2.D+00)+1.D+00/(2.D+00/hsmax+swkappas(Tsav))) & + + Aice*I_0*(1.D+00-albice)*transi(Tsav) + else if (.not. snow_dist .and. ice_hs .gt. 1.D-02) then + fluxt=fluxt+PenSW*swkappas(Tsav)*(exp(-swkappas(Tsav)*zi(1)) & + -exp(-swkappas(Tsav)*zi(2))) + Pari(2)=PenSW*exp(-swkappas(Tsav)*ice_hs) + else + fluxt=fluxt+PenSW*swkappai(Tsav)*(exp(-swkappai(Tsav)*zi(1)) & + -exp(-swkappai(Tsav)*zi(2))) + Pari(2)=PenSW*exp(-swkappai(Tsav)*zi(2)) + end if + + if (meltpond) then + !note for ice_hs ge hsmax Asnow is 1.0 and Amelt is 0.0, so no doublecounting + fluxt=fluxt+Amelt*I_0*(1.D+00-albmelt)*transm + Pari(2)=Pari(2)+Amelt*I_0*(1.D+00-albmelt)*transm + end if + +!......put short wave radiation absorbed in upper half of top layer +!......into surface warming + Sint(2)=Pari(1)*swkappai(Tice(2))*(exp(-swkappai(Tice(2))*zi(2))- & + exp(-swkappai(Tice(2))*zi(3)))/(0.5D+00*(dzi(1)+dzi(2))) + +!......do ice layers + do k=3,nilay-1 + zi(k+1)=zi(k)+dzi(k-1)/2.D+00+dzi(k)/2.D+00 + + Pari(k)=Pari(k-1)*exp(-swkappai(Tice(k))*dzi(k-1)) + Sint(k)=Pari(1)*swkappai(Tice(k))*(exp(-swkappai(Tice(k))*zi(k))- & + exp(-swkappai(Tice(k))*zi(k+1)))/(0.5D+00*(dzi(k-1)+dzi(k))) + enddo + + zi(nilay+1)=zi(nilay)+dzi(nilay-1)/2.D+00+dzi(nilay) +! zi(nilay+1)=zi(nilay)+dzi(nilay-1)/2.D+00+dzi(nilay)/2.D+00 !NScheck +!......all of short wave radiation absorbed in bottom-most layer +!......applied at interface between bottom two layers for simplicity + Pari(nilay)=Pari(nilay-1)*exp(-swkappai(Tice(nilay))*dzi(nilay-1)) + Sint(nilay)=Pari(1)*swkappai(Tice(k))*(exp(-swkappai(Tice(k))*zi(nilay)) & + -exp(-swkappai(Tice(k))*zi(nilay+1)))/(0.5D+00*dzi(nilay-1)+dzi(nilay)) + Pari(nilay+1)=Pari(nilay)*exp(-swkappai(Tice(nilay+1))*dzi(nilay)) !NScheck + endif +!...end snow_dist + + +! +!...Set up boundary conditions for 1-D heat conduction equation +! +!....if upper surface is melting, use temp. boundary condition + if(fluxt.ge.0.D+00.and.Tice(1).ge.Tmelts) then + bctype(1)=-1.D+00 !indicates temp. boundary condition + bcs(1)=Tmelts !Note, under melting conditions, surface +! temperature is taken to be melting point +! of snow (fresh water) + else +! +!....if upper surface is cooling, use flux boundary condition + bctype(1)=1.D+00 !indicates flux boundary condition + bcs(1)=fluxt !Note, under cooling conditions, fluxt is +! the outgoing flux calculated in SEBUDGET + endif +! +!....bottom boundary condition always freezing temp. of sea water + bctype(2)=-1.D+00 !indicates lower boundary has temp specified + bcs(2)=Tfreezi + +! +!...Solve 1-D heat conduction equation to get new temperature profile +!...and conductive flux at top and bottom +! + if(hlay.ge.hlaymin) then + call cndiffus(bctype,bcs,nilay,dzi,rhoCp,Cond,Sint,Tice) + +!...Check if layer is thinner than HLAYMIN but greater than 0. If so, +!...assume linear temperature profile and calculate conduction from +!...constant conductivity rather than solving transient diffusion +!...equation. If ice and snow have vanished, calculate ocean surface +!...temperature assuming specified mixed-layer depth. +! + elseif(hlay.lt.hlaymin .and. hlay.gt.0.D+00)then +!.......thin layer of ice: assume linear temp. profile and calculate +!.......new surface temperature + Tbot=Tfreezi + Ttop=(hlay/Condfi)*fluxt+Tbot + if(Ttop.gt.Tmelts) then + Ttop=Tmelts + do k=1,nilay + Tice(k)=Ttop + enddo + Tice(nilay+1)=Tbot + Iceflux(1)=0.D+00 + Iceflux(2)=0.D+00 + else + do k=1,nilay+1 + Tice(k)=Ttop+(Tbot-Ttop)*float(k-1)/float(nilay) + enddo + Iceflux(1)=fluxt + Iceflux(2)=fluxt + endif + endif +! +! +! LEVEL1'end therm1d' + + return + + end subroutine therm1d + +!EOC----------------------------------------------------------------------- + +!BOP +! !IROUTINE: +! +! !INTERFACE: + subroutine cndiffus(bctype,bcs,nilay,dzi,rhoCp,Cond,Sint,Tice) + +! !DESCRIPTION: +! +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! +!! Subroutine to solve 1-dimensional heat diffusion equation using a +!! generalized Crank-Nicholson or implicit algorithm which allows +!! arbitrary spacing of layers and spatially varying diffusion +!! coefficients. Parameter 'theta' determines whether standard +!! Crank-Nicholson or implicit algorithm is used. +!! +!! +!! sketch of grid layout: +!! +!! dzi(i-1) dzi(i) +!! ----+---------+--------------+----- +!! Tice(i-1) Tice(i) Tice(i+1) +!! Cond(i-1) Cond(i) +!! rhoCp(i-1) rhoCp(i) +!! rCpav(i-1) rCpav(i) rCpav(i+1) +!! +!! NOTE: minimum number of layers is 2 (this is required so that flux +!! boundary conditions can be specified at both boundaries). +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! +! +! !USES: + IMPLICIT NONE + +! !INPUT PARAMETERS: + real(rk), intent(in) :: bctype(2),bcs(2) + integer, intent(in) :: nilay + real(rk), intent(in) :: dzi(nilay),rhoCp(nilay),Cond(nilay) + real(rk), intent(in) :: Sint(nilay+1) +! !INOUT PARAMETERS: + real(rk), intent(inout) :: Tice(nilay+1) + +! !REVISION HISTORY: +! adjusted to GOTM F90 Nadja Steiner Dec 2008 +! Original author(s): Gregory M. Flato - Canadian Centre for +! Climate Modelling and Analysis +!BOC +! !LOCAL VARIABLES: +! coefficients for tridiagonal matrix + real(rk) :: C(nlmax+1,3) +! coefficients for RHS vector + real(rk) :: R(nlmax+1) +! ratio of conductivity to heat capacity + real(rk) :: rCpav(nlmax+1) +! array to store previous temperature + real(rk) :: To(nlmax+1) + integer :: j,l + +!EOP +!----------------------------------------------------------------------- + +! LEVEL1'cndiffus' +! +!...check if nilay is greater than maximum allowed +! + if(nilay.gt.nlmax) then + print*,'!!!! FATAL ERROR IN CNDIFFUS - NILAY > NLMAX' + stop + endif +! +!...check if nilay is less than minimum allowed + if(nilay.lt.2) then + print*,'!!!! FATAL ERROR IN CNDIFFUS - NILAY < 2' + stop + endif +! +!...calculate average volumetric heat capacity at temperature grid points +!...from layer values +! + do j=1,nilay-1 + rCpav(j+1)=(rhoCp(j)*dzi(j)+rhoCp(j+1)*dzi(j+1))/(dzi(j)+dzi(j+1)) + enddo + rCpav(1)=rhoCp(1) + rCpav(nilay+1)=rhoCp(nilay) +! +!...set up tridiagonal matrix coefficients for Crank-Nicholson scheme +! + do l=2,nilay + C(l,1)=-theta*(2.D+00*Cond(l-1)*dti/dzi(l-1))/(dzi(l)+dzi(l-1))/ & + rCpav(l) + C(l,2)=1.D+00+(theta*2.D+00*Cond(l-1)*dti/dzi(l-1)+theta*2.D+00* & + Cond(l)*dti/dzi(l))/(dzi(l)+dzi(l-1))/rCpav(l) + C(l,3)=-theta*(2.D+00*Cond(l)*dti/dzi(l))/(dzi(l)+dzi(l-1))/ & + rCpav(l) + enddo +! +!...set up right hand side vector +! + do l=2,nilay + R(l)=(1.D+00-theta)*(((2.D+00*Cond(l)*dti/dzi(l))/(dzi(l)+ & + dzi(l-1)))*(Tice(l+1)-Tice(l))-((2.D+00*Cond(l-1)*dti/dzi(l-1))/ & + (dzi(l)+dzi(l-1)))*(Tice(l)-Tice(l-1)))/rCpav(l)+Tice(l) +!......add internal heat sources (if any) + R(l)=R(l)+Sint(l)*dti/rCpav(l) + enddo +! +!...add temperature boundary conditions (if any) to RHS +! +!.....upper boundary + if(bctype(1).lt.0.D+00) then + C(1,1)=0.D+00 + C(1,2)=1.D+00 + C(1,3)=0.D+00 + R(1)=bcs(1) + endif +!.....lower boundary + if(bctype(2).lt.0.D+00) then + C(nilay+1,1)=0.D+00 + C(nilay+1,2)=1.D+00 + C(nilay+1,3)=0.D+00 + R(nilay+1)=bcs(2) + endif +! +!...add flux boundary conditions (if any) to RHS +! +!.....upper boundary + if(bctype(1).gt.0.D+00) then + C(1,1)=0.D+00 + C(1,2)=1.D+00+(2.D+00*dti*Cond(1)*theta/dzi(1)**2)/rCpav(1) + C(1,3)=-(2.D+00*dti*Cond(1)*theta/dzi(1)**2)/rCpav(1) + R(1)=Tice(1)+(Tice(2)-Tice(1))*(2.D+00*dti*Cond(1)*(1.D+00-theta)/ & + dzi(1)**2)/rCpav(1)+2.D+00*bcs(1)*dti/rCpav(1)/dzi(1) + endif +!.....lower boundary + if(bctype(2).gt.0.D+00) then + C(nilay+1,1)=-(2.D+00*dti*Cond(nilay)*theta/dzi(nilay)**2)/ & + rCpav(nilay+1) + C(nilay+1,2)=1.D+00+(2.D+00*dti*Cond(nilay)*theta/dzi(nilay)**2)/ & + rCpav(nilay+1) + C(nilay+1,3)=0.D+00 + R(nilay+1)=Tice(nilay+1)-((Tice(nilay+1)-Tice(nilay))*2.D+00*dti*Cond(nilay)* & + (1.D+00-theta)/dzi(nilay)**2)/rCpav(nilay+1)-(2.D+00* & + bcs(2)*dti/rCpav(nilay+1)/dzi(nilay)) + endif +! +!...save temperatures from previous time step for use in flux +!...calculation later +! + do l=1,nilay+1 + To(l)=Tice(l) + enddo +! +!...solve system of equations using efficient tridiagonal +!...matrix solver. Subroutine returns new values of temperature +!...in array 'Tice' - overwriting previous values. +! + call trisol(C,R,nilay+1,Tice) +! +!...calculate fluxes at boundaries where temperature boundary +!...conditions were specified. Flux returned in array 'flux' - +!...overwriting previous value. +! +!.....upper boundary + Iceflux(1)=-Cond(1)*((1.D+00-theta)*(To(2)-To(1))+theta* & + (Tice(2)-Tice(1)))/dzi(1) +!.....lower boundary + Iceflux(2)=-Cond(nilay)*((1.D+00-theta)*(To(nilay+1)-To(nilay)) & + +theta*(Tice(nilay+1)-Tice(nilay)))/dzi(nilay) +! + +! +! LEVEL1'end cndiffus' + + return + + end subroutine cndiffus +!EOC +!----------------------------------------------------------------------- + +!BOP +! +! !IROUTINE: +! +! !INTERFACE: + subroutine trisol(C,R,n,Tice) +! +! !DESCRIPTION: +! +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! +!! Subroutine to solve tri-diagonal matrix. Based on algorithm +!! described in Press et al. ("Numerical Recipies", Cambridge +!! University Press, 1986, pp. 40-41). + +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! +! +! !USES: + IMPLICIT NONE + +! !INPUT PARAMETERS: + real(rk), intent(in) :: C(nlmax+1,3) + real(rk), intent(in) :: R(nlmax+1) +! n - size of system of equations + integer, intent(in) :: n +! !INOUT PARAMETERS: + real(rk), intent(inout) :: Tice(n) + + +! !REVISION HISTORY: +! adjusted to GOTM F90 Nadja Steiner Dec 2008 +! Original author(s): Gregory M. Flato - Canadian Centre for +! Climate Modelling and Analysis +! !LOCAL VARIABLES: +! coefficients for + integer,parameter :: nmax=100 + real(rk),parameter :: eps=1.D-12 + real(rk) :: D1 + real(rk) :: temp(nmax) + integer :: j + +!EOP +!----------------------------------------------------------------------- + +! LEVEL1'trisol' + +!...dimension arrays +! +! dimension Tice(nilay+1),C(nmax,3),R(nmax),temp(nmax) +! +!...check that 'n' is not greater than nmax +! + if(n.gt.nmax) then + print*,'!!!! FATAL ERROR IN TRISOL - N > NMAX' + stop + endif +! +!...check that first diagonal element is not too small +! + if(C(1,2).lt.eps) then + print*,'!!!! FATAL ERROR IN TRISOL ' + print*,' - FIRST DIAGONAL ELEMENT < EPS' + stop + endif +! +!..."decomposition" phase of solution +! + D1=C(1,2) + Tice(1)=R(1)/D1 + do j=2,n + temp(j)=C(j-1,3)/D1 + D1=C(j,2)-C(j,1)*temp(j) + if(D1.lt.eps) then + print*,'!!!! FATAL ERROR IN TRISOL ' + print*,' - DIVISOR, D1 < EPS' + print*,' - Try making hlaymin larger' + stop + endif + Tice(j)=(R(j)-C(j,1)*Tice(j-1))/D1 + enddo +! +!..."back-substitution" phase of solution +! + do j=n-1,1,-1 + Tice(j)=Tice(j)-temp(j+1)*Tice(j+1) + enddo +! +! +! LEVEL1'end trisol' + + return + + + end subroutine trisol +!EOC!----------------------------------------------------------------------- + +!BOP +! +! !IROUTINE: +! +! !INTERFACE: + subroutine growthtb(rhowater,nilay,rhoCp,dzi,Tice,TopGrowth,TerMelt, & + TopMelt,BotGrowth,BotMelt,Hmix,ohflux) +! +! !DESCRIPTION: +! +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! +!! Subroutine to calculate growth or melt at the top and bottom of a +!! 1-dimensional slab of snow/sea ice given the calculated surface and +!! bottom fluxes. Also checks for submergence of ice surface and forms +!! 'snow ice' if neccessary. +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! +! +! !USES: + IMPLICIT NONE +! +!INPUT PARAMETERS: + real(rk), intent(in) :: rhowater + integer, intent(in) :: nilay + real(rk), intent(in) :: rhoCp(nilay),dzi(nilay) + real(rk), intent(in) :: ohflux +!OUTPUT PARAMETERS: + +!INOUT PARAMETERS: + real(rk), intent(inout) :: BotGrowth,BotMelt + real(rk), intent(inout) :: TopGrowth,TerMelt,TopMelt,Hmix + real(rk), intent(inout) :: Tice(nilay+1) +! +! !REVISION HISTORY: +! adjusted to GOTM F90 Nadja Steiner Dec 2008 +! Original author(s): Gregory M. Flato - Canadian Centre for +! Climate Modelling and Analysis +!BOC +! !LOCAL VARIABLES: +! coefficients for + real(rk) :: dhib,dhst,dhs,dhi,Ts +! coefficients for + real(rk) :: Tl +! snsub - submerged snow mass per unit area (kg m-2) + real(rk) :: snsub +! porewat - mass per unit area or pore water frozen into submerged +! snow (kg m-2) + real(rk) :: porewat + integer :: k + +!EOP +!----------------------------------------------------------------------- + +! LEVEL1'growthtb' +! +!...Calculate surface melt, if any +! + call surfmelt(Tice(1),TopMelt) + + ! + !...Calculate change in thickness due to bottom growth or melt + ! + + dhib=-(Iceflux(2)+ohflux)*dti/(rhoice*Hfi) + simass=simass+dhib*rhoice + if(dhib.gt.0.D+00) BotGrowth=BotGrowth+dhib + if(dhib.lt.0.D+00) BotMelt=BotMelt-dhib + + ! + !...Check if sufficient snow to depress ice surface below sea level, if + !...so, convert some snow to ice ("flooding"). + !...Add latent heat released by sea water + !...freezing in snow pores to upper ice layer. (Note, if this heat raises + !...temperature above freezing, some ice will subsequently be melted by + !...the calculations which follow + ! + + + snsub=snmass-(rhowater-rhoice)*max(0.D+00,simass)/rhoice + if(snsub.gt.0.D+00) then + !......limit submerged snow to that available + if(snmass-snsub.lt.0.D+00) then + snsub=snmass + snmass=0.D+00 + else + snmass=snmass-snsub + endif + + !......ice mass increase enhanced due to sea water filling pores in snow + !......and freezing + simass=simass+snsub*(1.D+00+(rhoice-rhosnow)/rhoice) + TopGrowth=TopGrowth+snsub*(1.D+00+(rhoice-rhosnow)/ & + rhoice)/rhoice + !......add latent heat released by freezing sea water to ice (this may + !......subsequently cause some melt if temp. is already near melting) + porewat=snsub*(rhoice-rhosnow)/rhoice + Tice(nslay+1)=Tice(nslay+1)+porewat*Hfi/(0.5D+00*(rhoCp(nslay)* & + dzi(nslay)+rhoCp(nslay+1)*dzi(nslay+1))) + + endif + ! + !...Calculate melt if temperature has risen above melting point). Note: + !...don't need to do bottom temperatures because it is always at the + !...freezing temp. of sea water. However, have to do surface temperature + !...in case conductive flux causes melting even when surface flux is + !...outward. + ! + !....if there is snow cover, check for snow melt + if(nslay.gt.0) then + !......first do surface layer + Ts=Tice(1) !nsnote this is not in CA code? + if(Ts.gt.Tmelts) then + !...snow_dist + if(snow_dist) then + dhst=Asnow*(Ts-Tmelts)*rhoCp(1)*dzi(1)*0.5D+00/(Hfw*rhosnow) + dhi=(Aice+Amelt)*(Ts-Tmelts)*rhoCp(1)*dzi(2)*0.5D+00/(Hfi*rhoice) + simass=simass-dhi*rhoice + Topmelt=Topmelt+dhi + else + dhst=(Ts-Tmelts)*rhoCp(1)*dzi(1)*0.5D+00/(Hfw*rhosnow) + end if + !...end snow_dist + Ts=Tmelts + !........reduce snow amount accordingly and replace heat that went unused + snmass=snmass-dhst*rhosnow + if(snmass.lt.0.D+00) then + Ts=Ts-snmass*Hfw/(rhoCp(1)*dzi(1)*0.5D+00) + snmass=0.D+00 + nslay=0 + endif + Tice(1)=Ts + endif + ! + ! + !......now check internal points + if(nslay.gt.1) then + do k=2,nslay + Tl=Tice(k) + if(Tl.gt.Tmelts) then + dhs=(Tl-Tmelts)*0.5D+00*(rhoCp(k-1)*dzi(k-1)+rhoCp(k) & + *dzi(k))/(Hfw*rhosnow) + Tl=Tmelts + !..........reduce snow thickness and replace heat that went unused + snmass=snmass-dhs*rhosnow + TerMelt=TerMelt+dhs + if(snmass.lt.0.D+00) then + TerMelt=TerMelt+snmass/rhosnow + Tl=Tl-snmass*Hfw/(0.5D+00*(rhoCp(k-1)*dzi(k-1) & + +rhoCp(k)*dzi(k))) + snmass=0.D+00 + endif + Tice(k)=Tl + endif + enddo + endif + endif + ! + !....do top level if no snow (NOTE: surface temp. taken as Tmelts) + + if(nslay.eq.0) then + Tl=Tice(1) + if(Tl.gt.Tmelts) then + dhi=(Tl-Tmelts)*rhoCp(1)*dzi(1)*0.5D+00/(Hfi*rhoice) + Tl=Tmelts + !........first use heat to melt remaining snow, if any + snmass=snmass-dhi*rhoice + if(snmass.lt.0.D+00) then + dhi=-snmass/rhoice + snmass=0.D+00 + else + dhi=0.D+00 + endif + !........reduce ice mass accordingly and replace heat that went unused + simass=simass-dhi*rhoice + TopMelt=TopMelt+dhi + Tice(1)=Tl + endif + endif + ! + !....now check internal points (NOTE: internal points melt at Tmelti) + do k=max(2,nslay+1),nilay + Tl=Tice(k) + if(Tl.gt.Tmelti) then + dhi=(Tl-Tmelti)*0.5D+00*(rhoCp(k-1)*dzi(k-1)+rhoCp(k) & + *dzi(k))/(Hfi*rhoice) + Tl=Tmelti + !........reduce ice thickness accordingly and replace heat that went unused + simass=simass-dhi*rhoice + TerMelt=TerMelt+dhi + Tice(k)=Tl + endif + enddo + ! + !...Check if negative ice has been created. If so, put heat into + !...mixed layer (Note, negative ice is used to account for additional + !...heat which may remain after all ice is melted) + ! + if(simass.le.0.D+00) then + !......add negative ice to mixed layer heat !NSnote, this is the dtemp comp from Winton? => should go into heat, is this done??? + Hmix=Hmix-min(0.D+00,simass)*Hfi + simass=0.D+00 + endif + ! + ! LEVEL1'end growthtb' + + return + + end subroutine growthtb + +!EOC!----------------------------------------------------------------------- + +!BOP +! +! !IROUTINE: +! +! !INTERFACE: + subroutine sebudget(TTss,ice_hi,ice_hs,evap,Qfluxes_uvic) +! +! !DESCRIPTION: +! +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! +!! Subroutine to calculate surface energy budget. Includes calculation of +!! surface fluxes from meteorological forcing data. +! NSnote: replaced original call fluxform,which calculated humidity, +! sensible, latent and incoming longwave heat flux with gotm call +! sequence: humidity, back_radiation, airsea_fluxes - this also meant +! calling airt, rather than airtk in sebudget and hence open_water +! and nr_iterate and the transfer of extra variables, e.g. methods + +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! +! +! !USES: + IMPLICIT NONE +! +! !INPUT PARAMETERS: + + real(rk), intent(in) :: TTss +! !INOUT PARAMETERS: + real(rk), intent(inout) :: evap ! freshwater evaporation (m/s) + real(rk), intent(inout) :: ice_hi,ice_hs + + interface + subroutine Qfluxes_uvic(T,qh,qe,qb) + integer, parameter :: rk = kind(1.d0) + real(rk), intent(in) :: T + real(rk), intent(out) :: qh,qe,qb + end subroutine + end interface + +! !REVISION HISTORY: +! adjusted to GOTM F90 Nadja Steiner Dec 2008 +! Original author(s): Gregory M. Flato - Canadian Centre for +! Climate Modelling and Analysis +!BOC +! !LOCAL VARIABLES: + + +!EOP +!----------------------------------------------------------------------- +! LEVEL1'sebudget' + + +!...Calculate ice and snow thickness from mass and density + ice_hs=snmass/rhosnow + ice_hi=simass/rhoice +! +!...Calculate surface fluxes from meteorological data using bulk formulae +! Recalculate surface fluxes for sea ice conditions +! NSnote TTss is in Kelvin!!! + + call Qfluxes_uvic(TTss,qh,qe,qb) + +! evap is set to zero for now....as for ice_winton, not sure what it is... + evap = 0 + +!...Calculate energy flux into surface +! +!.... long wave flux + fluxt=qb +! +!....sensible heat flux + fluxt=fluxt+qh +! +!....latent heat flux + fluxt=fluxt+qe +! +! +! LEVEL1'end sebudget' + + return + + end subroutine sebudget + +!EOC!----------------------------------------------------------------------- + +!BOP +! +! !IROUTINE: +! +! !INTERFACE: + subroutine surfmelt(Ts,TopMelt) + +! !DESCRIPTION: +! +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! +!! Subroutine to calculate melt from flux divergence at surface. Available +!! energy is first used to melt snow, and when snow has disappeared, +!! remaining energy is used to melt ice. +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! +! +! !USES: + IMPLICIT NONE +! +! !INPUT PARAMETERS: + real(rk), intent(in) :: Ts + real(rk) :: hmeltinput + +! !OUTPUT PARAMETERS: + +! !INOUT PARAMETERS: + real(rk), intent(inout) :: TopMelt + +! !REVISION HISTORY: +! adjusted to GOTM F90 Nadja Steiner Dec 2008 +! Original author(s): Gregory M. Flato - Canadian Centre for +! Climate Modelling and Analysis +!BOC +! !LOCAL VARIABLES: +! coefficients for + real(rk) :: Fdiv,dmsw,dmi,Fdivice,Fdivsnow + +!EOP +!----------------------------------------------------------------------- + +! LEVEL1'surfmelt' + +! +!...Calculate flux divergence at surface. Note: Iceflux(1) is the downward +!...positive heat flux defined at the centre of the uppermost layer +! + Fdiv=Iceflux(1)-fluxt + ! + !...If no melt is implied, return; otherwise, calculate surface melt + ! still apply constant drainrate for meltponds + ! + if(Fdiv.ge.0.D+00 .or. fluxt.lt.0.D+00 .or. Ts.lt.Tmelts) then + if (meltpond) then + meltmass=meltmass-drainrate*dti*Amelt + if (meltmass .lt. 0.D+00) then + meltmass= 0.D+00 + Amelt=0.D+00 + Aice=1-Asnow + end if + end if + return + else + ! + !......use heat to melt available snow, then ice + dmsw=-Fdiv*dti/Hfw + fluxt=0.D+00 + snmass=snmass-Asnow*dmsw + simass=simass-(Aice+Amelt)*dmsw*Hfw/Hfi + !nsnote again double snow check! rm one? + if (meltpond .and. snmass .gt. 0.D+00 .and. snmass/rhosnow .gt. hsmin) then + meltmass=meltmass+dmsw*rhosnow/rhowaterfresh-drainrate*dti*Amelt + + if (meltmass .lt. 0.D+00) then + meltmass=0.D+00 + Amelt=0.D+00 + Aice=1-Asnow + end if + else if (meltpond .and. snmass .lt. 0.D+00) then + meltmass=meltmass+dmsw*rhoice/rhowaterfresh-drainrate*dti*Amelt + if (meltmass .lt. 0.D+00) then + meltmass=0.D+00 + Amelt=0.D+00 + Aice=1-Asnow + end if + end if + + if(snmass.lt.0.D+00) then + ! + !........if all of snow melted, use remaining heat to melt ice + if (snow_dist) then + hsmax=hsmin + end if + dmi=-snmass*Hfw/Hfi + snmass=0.D+00 + simass=simass-dmi + TopMelt=TopMelt+dmi/rhoice + if(simass.lt.0.D+00) then + ! + !..........if all ice is melted, restore remaining surface flux to + !..........be added to mixed layer heat storage + TopMelt=TopMelt+simass/rhoice + fluxt=-simass*Hfi/dti + simass=0.D+00 + endif + endif + endif + + ! + ! LEVEL1'end surfmelt' + + return + + + end subroutine surfmelt + +!----------------------------------------------------------------------- + + + subroutine nr_iterate(nilay,I_0,Told,Tice,Pari,& + Sice_bulk,ice_hi,ice_hs,dzi,Cond,rhoCp,zi, & + Sint,evap,alb,Qfluxes_uvic) + +! !DESCRIPTION: +! +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! +!! Subroutine to perform Newton-Raphson iteration to solve for a +!! surface temperature consistent with the surface energy budget +!! including conductive heat flux. If no solution is found in +!! initial Newton-Raphson iteration, a bisection algorithm is +!! used. NOTE: if uppermost layer is greater than 0.1m thick, +!! surface temperature is defined as the temperature at the surface +!! Tice(1); however, if uppermost layer is less than 0.1m thick, surface +!! temperature is defined as the average temperature of the top layer. +!! This prevents rapid fluctuations in surface temperature and so +!! allows a longer time step to be taken. +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! + + +! !USES: + IMPLICIT NONE + +! !INPUT PARAMETERS: + integer, intent(in) :: nilay +! !OUTPUT PARAMETERS: + real(rk), intent(out) :: Sice_bulk + real(rk), intent(out) :: Told(nilay+1) + real(rk), intent(out) :: dzi(nilay),zi(nilay+1) + real(rk), intent(out) :: Sint(nilay+1),Cond(nilay),rhoCp(nilay) + real(rk), intent(out) :: alb + real(rk), intent(out) :: Pari(nilay+1) +! !INOUT PARAMETERS: + real(rk), intent(inout) :: Tice(nilay+1) + real(rk), intent(inout) :: ice_hi,ice_hs + real(rk), intent(inout) :: I_0 + real(rk), intent(inout) :: evap! freshwater evaporation (m/s) + +! +! !REVISION HISTORY: +! adjusted to GOTM F90 Nadja Steiner Dec 2008 +! Original author(s): Gregory M. Flato - Canadian Centre for +! Climate Modelling and Analysis + interface + subroutine Qfluxes_uvic(T,qh,qe,qb) + integer, parameter :: rk = kind(1.d0) + real(rk), intent(in) :: T + real(rk), intent(out) :: qh,qe,qb + end subroutine + end interface + + ! +!BOC +! !LOCAL VARIABLES: +! coefficients for + real(rk) :: fTs,dTemp,Tsp,fTsdT1,fTsdT2 + real(rk) :: fprime,Tsnew,Error,Ts1,fTs1 + real(rk) :: dTs,Tsu,fTsu,Ts2,Tsm,fTsm,fTsl, Tsl,Ts + integer :: nrit,ksearch,kb_uvic,l + real(rk),parameter :: toler=1.D-02 +!EOP +!----------------------------------------------------------------------- + + + +! LEVEL1'nr_iterate' +! +! +!...Save initial temperature profile +! + do l=1,nilay+1 + Told(l)=Tice(l) + enddo +! +!...First guess at surface temperature +! + Ts=Tice(1) + +! +!----------------------------------------------------------------------------- +!--- Start initial Newton-Raphson iteration - do a maximum of 5 iterations --- +!----------------------------------------------------------------------------- +! + do nrit=1,5 +! +!......calculate surface energy budget terms + call sebudget(Ts,ice_hi,ice_hs,evap,Qfluxes_uvic) + call albedo_ice_uvic(fluxt,I_0,PenSW,alb,ice_hs,ice_hi,Ts) +! +! +!......restore intial temperature profile as therm1d takes a forward time step + do l=1,nilay+1 + Tice(l)=Told(l) + enddo + +!......solve unsteady heat conduction equation to get new temperature profile +!......and bottom flux after one time step + call therm1d(nilay,Sice_bulk,ice_hi,ice_hs,dzi,Cond,& + rhoCp,zi,Sint,Pari,Tice,I_0) +! +!......now find mismatch in assumed and computed surface temperature + fTs=Ts-Tice(1) +! +!......stop iterating if mismatch is less than tolerance + if(abs(fTs).lt.toler) then +! print*,'Newton-Raphson scheme met tolerance after', & +! nrit-1,'iterations' + go to 987 + endif +! +!......re-do above calculations with slightly increased and decreased +!......initial temperature to calculate first derivative by centred +!......finite difference + dTemp=0.1 + Tsp=Ts+dTemp + + call sebudget(Tsp,ice_hi,ice_hs,evap,Qfluxes_uvic) + call albedo_ice_uvic(fluxt,I_0,PenSW,alb,ice_hs,ice_hi,Tsp) +! + do l=1,nilay+1 + Tice(l)=Told(l) + enddo + call therm1d(nilay,Sice_bulk,ice_hi,ice_hs,dzi, & + Cond,rhoCp,zi,Sint,Pari,Tice,I_0) + fTsdT1=Tsp-Tice(1) + Tsp=Ts-dTemp + + call sebudget(Tsp,ice_hi,ice_hs,evap,Qfluxes_uvic) + call albedo_ice_uvic(fluxt,I_0,PenSW,alb,ice_hs,ice_hi,Tsp) +! + do l=1,nilay+1 + Tice(l)=Told(l) + enddo + call therm1d(nilay,Sice_bulk,ice_hi,ice_hs,dzi, & + Cond,rhoCp,zi,Sint,Pari,Tice,I_0) + fTsdT2=Tsp-Tice(1) +! +!......if derivative vanishes, take solution obtained at this point +!......and don't iterate any further + if(abs(fTsdT1-fTsdT2).lt.1.D-02) then + print*,'fTsdT1~=fTsdT2 after ',nrit-1,' iterations' + go to 987 + endif +! +!......otherwise, compute derivative and hence next guess from +!......Newton-Raphson formula + fprime=(fTsdT1-fTsdT2)/(2.D+00*dTemp) + Tsnew=Ts-fTs/fprime +! +!......check error and stop iterating if tolerance is met + Error=Tsnew-Ts + if(abs(Tsnew-Ts).lt.toler) then +! print*,'Newton-Raphson scheme met tolerance after', & +! nrit,'iterations' + go to 987 + endif +! +!......if another iteration is required, update surface temperature +!......and try again + Ts=Tsnew + enddo +! +!...If iteration has not reached a solution at this point, (which is +!...rare) we have to resort to a cruder, brute force approach. So, ... +! + print*,'***Newton-Raphson scheme failed: Error =',Error + print*,'***Starting method of Bisection' +! +!---------------------------------------------------------------------- +!--- If first try at Newton-Raphson iteration fails, try method of --- +!--- bisection which is less efficient, but more robust --- +!---------------------------------------------------------------------- +! +!...To start, find where Ts-Tice(1) changes sign +! +!....first, do initial guess again + Ts1=Told(1) + + call sebudget(Ts1,ice_hi,ice_hs,evap,Qfluxes_uvic) + call albedo_ice_uvic(fluxt,I_0,PenSW,alb,ice_hs,ice_hi,Ts1) +! + do l=1,nilay+1 + Tice(l)=Told(l) + enddo + call therm1d(nilay,Sice_bulk,ice_hi,ice_hs,dzi, & + Cond,rhoCp,zi,Sint,Pari,Tice,I_0) + fTs1=Ts1-Tice(1) +! +!....increase and decrease Ts in 1 degree intervals until zero-crossing +!....is bracketed + dTs=1. + do ksearch=1,10 + Tsu=Ts1+dTs*float(ksearch-1) + + call sebudget(Tsu,ice_hi,ice_hs,evap,Qfluxes_uvic) + call albedo_ice_uvic(fluxt,I_0,PenSW,alb,ice_hs,ice_hi,Tsu) +! + do l=1,nilay+1 + Tice(l)=Told(l) + enddo + call therm1d(nilay,Sice_bulk,ice_hi,ice_hs,dzi, & + Cond,rhoCp,zi,Sint,Pari,Tice,I_0) + fTsu=Tsu-Tice(1) +!......check if crossing is between Ts1 and Tsu>Ts1 + if(fTsu.gt.0.D+00.and.fTs1.le.0.D+00 & + .or. & + fTsu.lt.0.D+00.and.fTs1.ge.0.D+00) then + Ts2=Tsu + go to 887 + endif +! + Tsl=Ts1-dTs*float(ksearch-1) + + call sebudget(Tsl,ice_hi,ice_hs,evap,Qfluxes_uvic) + call albedo_ice_uvic(fluxt,I_0,PenSW,alb,ice_hs,ice_hi,Tsl) +! + do l=1,nilay+1 + Tice(l)=Told(l) + enddo + call therm1d(nilay,Sice_bulk,ice_hi,ice_hs,dzi, & + Cond,rhoCp,zi,Sint,Pari,Tice,I_0) + fTsl=Tsl-Tice(1) +!......check if crossing is between Ts1 and Tsl> linear profile + if(Sice_bulk .lt. S1) alpha=1.D+00 +!...for a high bulk salinity (Sice_bulk > S2 = 4.5 ppt) >> const profile + if(Sice_bulk .gt. S2) alpha=0.D+00 +!...for intermediate bulk salinity >> intermediate values between previous cases + if(Sice_bulk .gt. S1 .and. Sice_bulk .lt. S2) then + alpha=(Sice_bulk-S2)/(S1-S2) + end if + + do k=1,nilay+1 + Sice(k) = alpha*Szero(k)+(1.D+00-alpha)*Sice_bulk + end do + + +! Ice-ocean salt flux (Tartinville et al. 2001) +! Terms representing salt rejection of the sea ice and +! salt input to the ocean. +! +! Freshwater flux + +! Flux due to snow melt + if (meltpond.and.Amelt.ne.0.0) then + rates=meltmass-meltmasso + else + rates=snmass-snmasso + rates=rates*(rhosnow/rhof) + endif +! to be summed to ppt,evap and river inflow in GOTM + Ff = (-1.D+00)* min(0.D+00,rates)/dti ![kg m-2 s-1] !sign positive if any freshwater goes into the ocean + Ff = Ff/rhow ![m s-1] + rates=0 + +! Flux due to brine drainage +! + Fb=simass*(x3+x4)/dti + +! Basal accretion +! + xx2=-rhoice*(SSS-Sice(nilay+1))*max(0.D+00,higb)/dti !Vancoppenole multiplies it by the ice fractionation at the mth categorie, +!check how it would work properly using Hibler scheme (check 'g' in the black book) + +! Snow ice formation +! + xx3=-rhoice*(SSS-Ssi)*(max(0.D+00,higs))/dti ![ppt kg m-2 s-1] + +! Melt of saline ice +! + xx4=rhoice*(SSS-Sice_bulk)*(hims+himb+himi)/dti + + Feq=xx2+xx3+xx4 + + Fs = (Fb + Feq)/rhow ! [ppt m s-1] - sign with respect to ice ! neg => S transfer from ice to water + ! pos => S transfer from water to ice ! need to invert sign in S equation + + + return + end subroutine saltice_prof_simple + +!----------------------------------------------------------------------- +!EOP + +!BOP +! +! !IROUTINE: Calculation of ice/snow albedo +! +! !INTERFACE: + subroutine albedo_ice_uvic(fluxt,I_0,PenSW,alb,ice_hs,ice_hi,TTss) +! +! !DESCRIPTION: +! calculation of ice/snow albedo taken out of original sebudget routine. +! Adjust incoming shortwave with ice albedo + +! !USES: + + IMPLICIT NONE +! +! !INPUT PARAMETERS: + real(rk), intent(in) :: I_0,ice_hs,ice_hi,TTss +! !OUTPUT PARAMETERS: + real(rk), intent(out) :: alb,PenSW +! !INPUT/OUTPUT PARAMETERS: + real(rk), intent(inout) :: fluxt + +! +! !REVISION HISTORY: +! Original author(s): N. Steiner +! +! See log for the icemodel module +! !LOCAL VARIABLES: +! snow_dist + real(rk) :: I_0_tmp +! +!snow_dist +! +!EOP +!----------------------------------------------------------------------------- +!BOP +!----------------------------------------------------------------------------- +! LEVEL1 'albedo_ice' + +!...snow_dist + if (ice_hs .ge. hsmax.and.ice_hs.ne.0.D+00) then + !nsnote isn't hsmax always gt 0.0 ? + Asnow=1.D+00 + Aice=0.D+00 + Amelt=0.D+00 + meltmass=0.D+00 + else if (ice_hs .lt. hsmin .and. ice_hi .ge. hsmin) then + if (meltpond .and. meltmass .gt. 0.D+00) then + Amelt=Ameltmax + Aice=1.D+00-Amelt + else + Amelt=0.D+00 + Aice=1.D+00 + end if + Asnow=0.D+00 + else if (ice_hi .lt. hsmin) then + Aice=0.D+00 + Asnow=0.D+00 + Amelt=1.D+00 + meltmass=0.D+00 + else + if(snow_dist .and. distr_type .eq. 0) then + Asnow=exp(-1.D+00*(inverfc(ice_hs/hsmax)**2)) + else if(snow_dist .and. distr_type .eq. 1) then + Asnow=exp(W(-2.D+00*ice_hs/hsmax*exp(-2.D+00))+2.D+00) & + *(-W(-2.D+00*ice_hs/hsmax*exp(-2.D+00))-2.D+00+1.D+00) + else!Nfix + Asnow=1.D+00 !Nfix + end if + if (meltpond .and. meltmass .gt. 0.D+00) then + if(snow_dist) then + Amelt=min(Ameltmax,1.D+00-Asnow) + Aice=1.D+00-Asnow-Amelt + else + Amelt=Ameltmax !uniform snow with meltponds + Asnow=1.0D+00-Amelt !uniform snow + Aice=0.0D+00 !Nfix + end if + else + Aice=1.D+00-Asnow + Amelt=0.D+00 + end if + endif + !H!begin + ! ice and snow albedo + if (albice_method.eq.1) then + albice=albice_f + albsnow=albsnow_f + else if (albice_method.eq.2) then + albice=(0.44D+00*ice_hi**0.28D+00)+0.08D+00 + albice=max(albice,albmelt) + albsnow=albsnow_f + else if (albice_method.eq.3) then + if (airtk.lt.273.05) then + albice=max(albmelt,(0.44D+00*ice_hi**0.28D+00)+0.08D+00) + else + albice=min(albice_m,(0.075D+00*ice_hi**2.D+00)+albmelt) + end if + if (airtk.lt.273.15) then + albsnow=albsnow_f + else + albsnow=albsnow_m + end if + else if (albice_method.eq.4) then + albice=((albice_m+albice_f)+(albice_m-albice_f)*tanh(TTss-273.15D+00))/2.D+00 + albsnow=((albsnow_m+albsnow_f)+(albsnow_m-albsnow_f)*tanh(TTss-273.15D+00))/2.D+00 + end if + !H!end + !...Calculate penetrating SW flux + ! + ! + + + PenSW=Asnow*I_0*(1.D+00-albsnow)*transs(TTss)+ & + Aice*I_0*(1.D+00-albice)*transi(TTss)+ & + Amelt*I_0*(1.D+00-albmelt)*transm + + + ! Calculate short wave flux absorbed at the surface + + + I_0_tmp=Asnow*I_0*(1.D+00-albsnow)*(1.D+00-transs(TTss))+ & + Aice*I_0*(1.D+00-albice)*(1.D+00-transi(TTss))+ & + Amelt*I_0*(1.D+00-albmelt)*(1.D+00-transm) + + + !....Add short wave flux absorbed at surface to total flux + fluxt=fluxt+I_0_tmp + + alb=Asnow*albsnow+Aice*albice+Amelt*albmelt +! +!...end snow_dist +! + + return + + end subroutine albedo_ice_uvic + + +!----------------------------------------------------------------------- + +! subroutine ice_interpol(var_grid2,var_grid1,nilay_int,nint_total) ! jpnote: not used + +!----------------------------------------------------------------------- + +!BOP +! +! !IROUTINE: Cleaning up the mean flow variables +! +! !INTERFACE: + subroutine clean_ice_uvic() +! +! !DESCRIPTION: +! De-allocates all memory allocated via init\_icemodel() +! +! !USES: +! use ncdfout, only: close_ncdf + + IMPLICIT NONE +! +! !INPUT PARAMETERS: +! +! !REVISION HISTORY: +! Original author(s): +! +! See log for the icemodel module +! +!EOP +!----------------------------------------------------------------------- +!BOC + !LEVEL1 'clean_ice_uvic' + + !LEVEL3 'closing ice.nc file...' +! call close_ncdf() + + return + + end subroutine clean_ice_uvic +!EOC +!---------------------------------------------------------------------- +! end internal subroutines +!---------------------------------------------------------------------- +! +! +!----------------------------------------------------------------------- +! begin internal functions +!----------------------------------------------------------------------- +!BOF +! !IROUTINE: inverse complementary error function +! +! !INTERFACE: + real function erfc(x) + +!DESCRIPTION: Note: for FORTRAN2008 and higher ERFC is an +! intrinsic Fortran function! + +!USES: + implicit none + +!INPUT PARAMETERS: + real(rk), intent(in) :: x + +!LOCAL VARIABLES: + real(rk) :: t,z +! !REVISION HISTORY: +! Original author(s): +! +! See log for the icemodel module +! +!----------------------------------------------------------------------- +!BOC + z = abs(x) + t = 1.0 / ( 1.0 + 0.5 * z ) + + erfc = t * exp( -z * z - 1.26551223 + t * & + ( 1.00002368 + t * ( 0.37409196 + t * & + ( 0.09678418 + t * (-0.18628806 + t * & + ( 0.27886807 + t * (-1.13520398 + t * & + ( 1.48851587 + t * (-0.82215223 + t * 0.17087277 ))))))))) + + if ( x.lt.0.0 ) erfc = 2.0 - erfc + + return +end function erfc + +!EOF +!---------------------------------------------------------------------- + +!BOF +! !IROUTINE: inverse complementary error function +! +! !INTERFACE: + real function inverfc(y) + +!DESCRIPTION: + +!USES: + implicit none + +!INPUT PARAMETERS: + real(rk), intent(in) :: y + +!LOCAL VARIABLES: + real(rk), parameter :: qa = 9.16461398268964d-01, & + & qb = 2.31729200323405d-01, & + & qc = 4.88826640273108d-01, & + & qd = 1.24610454613712d-01, & + & q0 = 4.99999303439796d-01, & + & q1 = 1.16065025341614d-01, & + & q2 = 1.50689047360223d-01, & + & q3 = 2.69999308670029d-01, & + & q4 = -7.28846765585675d-02, & + + & pa = 3.97886080735226000d+00, & + & pb = 1.20782237635245222d-01, & + & p0 = 2.44044510593190935d-01, & + & p1 = 4.34397492331430115d-01, & + & p2 = 6.86265948274097816d-01, & + & p3 = 9.56464974744799006d-01, & + & p4 = 1.16374581931560831d+00, & + & p5 = 1.21448730779995237d+00, & + & p6 = 1.05375024970847138d+00, & + & p7 = 7.13657635868730364d-01, & + & p8 = 3.16847638520135944d-01, & + & p9 = 1.47297938331485121d-02, & + & p10 = -1.05872177941595488d-01, & + & p11 = -7.43424357241784861d-02, & + + & p12 = 2.20995927012179067d-03, & + & p13 = 3.46494207789099922d-02, & + & p14 = 1.42961988697898018d-02, & + & p15 = -1.18598117047771104d-02, & + & p16 = -1.12749169332504870d-02, & + & p17 = 3.39721910367775861d-03, & + & p18 = 6.85649426074558612d-03, & + & p19 = -7.71708358954120939d-04, & + & p20 = -3.51287146129100025d-03, & + & p21 = 1.05739299623423047d-04, & + & p22 = 1.12648096188977922d-03 + + + real(rk) :: z, w, u, s, t +! !REVISION HISTORY: +! Original author(s): +! +! See log for the icemodel module +! +!----------------------------------------------------------------------- +!BOC + z = y + if (y .gt. 1d0) then + z = 2d0 - y + end if + w = qa - log(z) + u = sqrt(w) + s = (qc + log(u)) / w + t = 1d0 / (u + qb) + inverfc = u * (1d0 - s * (0.5d0 + s * qd)) - & + ((((q4 * t + q3) * t + q2) * t + q1) * t + q0) * t + t = pa / (pa + inverfc) + u = t - 0.5d0 + s = (((((((((p22*u+p21)*u+p20)*u+p19)*u+p18)*u+p17)*u+p16)*u+ & + p15)*u+p14)*u+p13)*u+p12 + s = ((((((((((((s*u+p11)*u+p10)*u+p9)*u+p8)*u+p7)*u+ & + p6)*u+p5)*u+p4)*u+p3)*u+p2)*u+p1)*u+p0)*t- & + z*exp(inverfc*inverfc-pb) + inverfc = inverfc + s * (1d0 + inverfc * s) + if (y .gt. 1d0) then + inverfc = -inverfc + end if + + return + end function inverfc + +!EOF +!----------------------------------------------------------------------- + +!BOF +! !IROUTINE: +! +! !INTERFACE: + real function W(y) +! +!DESCRIPTION: W??? +! +!USES: + implicit none +!INPUT PARAMETERS: + real(rk), intent(in) :: y +!LOCAL VARIABLES: + real(rk), parameter :: M1=0.3361D+00, M2=-0.0042D+00 + real(rk), parameter :: M3=-0.0201D+00 + real(rk) :: sigma +! !REVISION HISTORY: +! Original author(s): +! +! See log for the icemodel module +! +!----------------------------------------------------------------------- +!BOC + sigma=-1.D+00-log(-y) + + W=-1.D+00-sigma-2.D+00/M1*(1.D+00-1.D+00/(1.D+00+ & + ((M1*sqrt(sigma/2.D+00))/(1.D+00+M2*sigma*exp(M3*sqrt(sigma)))))) + + end function W +!EOF +!---------------------------------------------------------------------- + +!BOF +! !IROUTINE: snow extinction +! +! !INTERFACE: + + real function swkappas(Tin) +!DESCRIPTION: + +!USES: + IMPLICIT NONE +!INPUT PARAMETERS: + real(rk), intent(in) :: Tin +! !REVISION HISTORY: +! Original author(s): +! +! See log for the icemodel module +! +!----------------------------------------------------------------------- +!BOC + swkappas=((swkappasm+swkappasf)+(swkappasm-swkappasf)* & + tanh(Tin-273.15D+00))/2.D+00 + +!test +! swkappas=1.5 + + end function swkappas +!EOF +!---------------------------------------------------------------------- +!BOF +! !IROUTINE: +! +! !INTERFACE: + real function swkappai(Tin) +!DESCRIPTION: + +!USES: + implicit none +!INPUT PARAMETERS: + real(rk), intent(in) :: Tin +! !REVISION HISTORY: +! Original author(s): +! +! See log for the icemodel module +! +!----------------------------------------------------------------------- +!BOC + swkappai=((swkappaim+swkappaif)+(swkappaim-swkappaif)* & + tanh(Tin-273.15D+00))/2.D+00 +!test +! swkappai=1.5 + + end function swkappai +!EOF +!---------------------------------------------------------------------- + +!BOF +! !IROUTINE: snow transmissivity +! +! !INTERFACE: + + real function transs(Tin) +!DESCRIPTION: + +!USES: + implicit none +!INPUT PARAMETERS: + real(rk), intent(in) :: Tin +! !REVISION HISTORY: +! Original author(s): +! +! See log for the icemodel module +! +!----------------------------------------------------------------------- +!BOC + transs=((transsm+transsf)+(transsm-transsf)* & + tanh(Tin-273.15D+00))/2.D+00 + !test +! transs=0.05 + + end function transs +!EOF +!---------------------------------------------------------------------- +!BOF +! !IROUTINE: ice transmissivity +! +! !INTERFACE: + + real function transi(Tin) +!DESCRIPTION: + +!USES: + +!INPUT PARAMETERS: + real(rk), intent(in) :: Tin + +!LOCAL VARIABLES: +! !REVISION HISTORY: +! Original author(s): +! +! See log for the icemodel module +! +!----------------------------------------------------------------------- +!BOC + +!H!begin +! transi=0.5 + transi=((transim+transif)+(transim-transif)*tanh(Tin-273.15D+00))/2.D+00 +!H!end + + end function transi +!EOF + + +end module stim_flato + +!----------------------------------------------------------------------- +! Copyright by the GETM-team under the GNU Public License - www.gnu.org +!----------------------------------------------------------------------- diff --git a/src/models/ice_flato.F90 b/src/models/ice_flato.F90 deleted file mode 100644 index e69de29..0000000 diff --git a/src/models/stim_models.F90 b/src/models/stim_models.F90 index dde531c..c575a1f 100644 --- a/src/models/stim_models.F90 +++ b/src/models/stim_models.F90 @@ -18,6 +18,9 @@ module stim_models #endif #ifdef STIM_WINTON use stim_winton, only: init_stim_winton, do_stim_winton +#endif +#ifdef STIM_FLATO + use stim_flato, only: init_stim_flato, do_ice_uvic #endif IMPLICIT NONE ! diff --git a/src/models/stim_variables.F90 b/src/models/stim_variables.F90 index dd317b9..2b6b3d1 100644 --- a/src/models/stim_variables.F90 +++ b/src/models/stim_variables.F90 @@ -31,7 +31,7 @@ module stim_variables real(rk), target, public :: Hfrazil = 0._rk ! Total frazil ice thickness real(rk), target, public :: Hsnow = 0._rk ! Total snow thickness real(rk), target, public :: Hice = 0._rk ! Total ice thickness - real(rk), target, public :: dHis = 0._rk ! surface ice growth + real(rk), target, public :: dHis = 0._rk ! surface ice growth real(rk), target, public :: dHib = 0._rk ! bottom ice growth real(rk), target, public :: albedo_ice = 0._rk real(rk), target, public :: attenuation_ice = 0._rk @@ -40,6 +40,208 @@ module stim_variables ! Lebedev real(rk), target, public :: fdd + + ! Flato +!---------------------------------------------------------------------------------------------------------------------------------- +! Flato PUBLIC DATA MEMBERS from uvic_ice.F90 +!---------------------------------------------------------------------------------------------------------------------------------- + +! parameters + public :: hlaymin,rhoice,Tfreezi,rCpmix,Hfi + +! !DEFINED PARAMETERS: +! hsmin - minimum snow thickness required to have a separate snow +! layer. If 'ice_hs' is less than this snow is ignored in +! heat conduction scheme (m) + real(rk), parameter :: hsmin = 0.001D+00 +! theta - a parameter between 0.5 and 1. which determines how +! 'implicit' the scheme is. theta=0.5 is the conventional +! Crank-Nicholson algorithm and theta=1. is the fully +! implicit scheme (the method of choice in this model because +! C-R scheme allows spurious temperature oscillations when +! layer thickness becomes small - an instability in the +! numerical scheme which is avoided by using theta=1). + real(rk), parameter :: theta = 1.D+00 +! sigma - Stefan-Boltzmann constant (W m-2 K-4) + real(rk), parameter :: sigma = 5.67D-08 +! epsilon - emissivity of ice (dimensionless) + real(rk), parameter :: epsilon = 0.99D+00 +! PenFrac - fraction of incoming short wave radiation that penetrates the surface + real(rk), parameter :: PenFrac = 0.17D+00 +! hlaymin - thickness below which a linear temperature profile is assumed (m) + real(rk), parameter :: hlaymin = 0.2D+00 +! rhoscold - specified 'cold' snow density (kg m-3) + real(rk), parameter :: rhoscold = 330.D+00 +! rhoswarm - specified 'warm' snow density (kg m-3) + real(rk), parameter :: rhoswarm = 450.D+00 +! rhowaterfresh - fresh water density (kg m-3) + real(rk), parameter :: rhowaterfresh = 1000.D+00 +! rhoice - ice density (kg m-3) + real(rk), parameter :: rhoice = 913.D+00 +! kelvin - zero deg Celsius (K) + real(rk), parameter :: kelvin = 273.15D+00 +! Tmelts - melting temperature of snow (fresh water) (K) + real(rk), parameter :: Tmelts = 273.15D+00 +! Tmelti - melting temperature of sea ice (K) + real(rk), parameter :: Tmelti=273.05D+00 +! Condfi - conductivity of pure ice (W m-1 K-1) + real(rk), parameter :: Condfi = 2.034D+00 +! rhoCpfi - heat capacity of pure ice (J m-3 K-1) + real(rk), parameter :: rhoCpfi = 1.883D+06 +! rCpmix - volumetric heat capacity of sea water (J m-3 K-1) + real(rk), parameter :: rCpmix = 1025.D+00*3.99D+03 +! Hfi - latent heat of fusion of sea ice (J kg-1) + real(rk), parameter :: Hfi = 2.93D+05 +! Hfw - latent heat of fusion of fresh water (J kg-1) + real(rk), parameter :: Hfw = 3.335D+05 +! swkappa - bulk short-wave extinction coefficient (m-1) + real(rk), parameter :: swkappa = 1.5 +! Tfreezi - freezing temperature of sea water (K) + real(rk), parameter :: Tfreezi = 271.2D+00 + +! !Maximum snow and ice layers + integer, public, parameter :: nlmax=99 + +!----------------------------------------------------------- +! YAML VARIABLES from old code nml +!----------------------------------------------------------- + +! !initializing yaml variables to reasonable defaults +! nilay - number of snow & ice layers + integer, public :: nilay = 0 +! sfall_method - define how snow fall is determined +! 1:constant snow fall +! 2:calculate snowfall from precipitation + integer, public :: sfall_method = 1 +! const_sfall - constant snow fall rate (m d^-1) + real(rk), public :: const_sfall = 0.0017D+00 +! dfact - drift factor allowing a factor to increase snow fall from +! precipitation via drifing snow (only for sfall_method=2) + real(rk), public :: dfact = 1.0D+00 +! depmix - prescribed mixed layer depth for open_water +! calculation. If set to zero top layer thickness is used. + real(rk), public :: depmix = 10.0D+00 +! sice_method - define how sea-ice salinity is to be calculated +! 1: constant ice salinity +! 2: Simple ice salinity profile (Vancoppenolle et al. 2009?) + integer, public :: sice_method = 1 +! snow_dist - logical switch between uniform and Weibull-distributed snow +! if true a snow distribution function will be applied +! if false a uniform snow thickness will be applied + logical, public :: snow_dist = .false. +! const_Sice - prescribed sea ice salinity (ppt) + real(rk), public :: const_Sice = 6.0D+00 +! distr_type - integer to chose the type of distribution + integer, public :: distr_type = -1 +! meltpond - logical to switch on meltpond + logical, public :: meltpond = .false. +! Ameltmax - maximum meltpond area + real(rk), public :: Ameltmax = 0.25D+00 +! drainrate - fixed meltpond drainrate - (m/d) converted to (kg/m^2/s) + real(rk), public :: drainrate = 0.0175D+00 +! hh0 - initial thickness for S calculation + real(rk), public :: hh0 = 0.15D+00 +! ice_hi_i - initial ice thickness (m) + real(rk), public :: ice_hi_i = 0.2D+00 +! ice_hs_i - initial snow thickness (m) + real(rk), public :: ice_hs_i = 0.01D+00 +! albice_method - 1:constant albedo (albice_f and albsnow_f) +! 2:albice dependent on ice thickness - attempt to +! capture the evolution of melt ponds (Flato&Brown 1996) +! albsnow is set to albsnow_f (constant) +! 3:albice and albsnow based on eq12&13 of Flato&Brown1996 +! 4:albice and albsnow dependent on temperature (same as transs/transi formulation) + integer, public :: albice_method = 1 +! albice_f - freezing ice albedo + real(rk), public :: albice_f = 0.55D+00 +! albmelt - melt pond albedo + real(rk), public :: albmelt = 0.2D+00 +! albsnow_f - freezing snow albedo + real(rk), public :: albsnow_f = 0.75D+00 +! albice_m - melting ice albedo + real(rk), public :: albice_m = 0.75D+00 +! albsnow_m - melting snow albedo + real(rk), public :: albsnow_m = 0.55D+00 +! transsf - freezing snow transmission coefficient + real(rk), public :: transsf = 0.05D+00 +! transsm - melting snow transmission coefficient + real(rk), public :: transsm = 0.08D+00 +! transif - freezing ice transmission coefficient + real(rk), public :: transif = 0.5D+00 +! transim - melting ice transmission coefficient + real(rk), public :: transim = 0.5D+00 +! transm - melt pond transmision coefficient + real(rk), public :: transm = 0.8D+00 +! swkappasm - melting snow extinction coefficient + real(rk), public :: swkappasm = 7.5D+00 +! swkappasf - freezing snow extinction coefficient + real(rk), public :: swkappasf = 14.0D+00 +! swkappaim - melting ice extinction coefficient + real(rk), public :: swkappaim = 0.8D+00 +! swkappaif - freezing ice extinction coefficient + real(rk), public :: swkappaif = 1.2D+00 + + +!------------------------------------------------------------------------------------------------------------- +! Flato PUBLIC DATA MEMBERS from ice.F90 +!------------------------------------------------------------------------------------------------------------- + +! University of Victoria ice model +! ts,tb are the upper surface temperature of the ice/snow and tb is the lowest layer temperature (not the bottom temperature, which is set to Tfreezi) all in (K) + real(rk), public, target :: ice_hs,ice_hi + real(rk), public, target :: ice_uvic_ts + real(rk), public, target :: ice_uvic_tb + real(rk), public, target :: ice_uvic_swr_0 + real(rk), public, target :: ice_uvic_precip_i + real(rk), public, target :: ice_uvic_sfall_i + real(rk), public, target :: ice_uvic_parb + real(rk), public, target :: ice_uvic_parui + real(rk), public, target :: ice_uvic_topmelt ! top melting - ice mass melted at the surface (snow+ice) at time step(m) + real(rk), public, target :: ice_uvic_botmelt ! bottom melting - ice mass melted at the ice bottom at time step(m) + real(rk), public, target :: ice_uvic_topgrowth ! top growth ice mass growth at slab surface due to snow submersion (m) + real(rk), public, target :: ice_uvic_botgrowth ! bottom growth - ice mass growth at the ice bottom at time step (m) + real(rk), public, target :: ice_uvic_termelt ! internal melting - ice mass melted in the ice interior at time step (m) + + real(rk), public, target :: ice_uvic_Fh = 0._rk ! interface heat flux (W/m2) + real(rk), public, target :: ice_uvic_Ff = 0._rk ! interface freshwater flux (m s-1) + real(rk), public, target :: ice_uvic_Fs = 0._rk ! interface salt flux - (ppt m s-1) + + + + real(rk), public, target :: ice_uvic_Sicebulk + real(rk), public, target :: ice_uvic_Hmix ! transferred energy - check (m) + ! Hmix - mixed layer heat storage (J m-2) =======> accounts only for + ! the SWR which crosses the ice slab and reach the water. keep it for now + real(rk), public, target :: ice_uvic_Aice ! ice area fraction which is : open + real(rk), public, target :: ice_uvic_Asnow ! ice area fraction which is : snow + real(rk), public, target :: ice_uvic_Amelt ! ice area fraction which is : meltpond + real(rk), public , target :: ice_uvic_hm +! University of Victoria ice model: allocatable variables +! ice_uvic_Tice - temperature array (dimensioned nlay+1), Tice(1) is the +! upper surface temperature, Tice(nlay+1) is the bottom +! temperature, remaining values at layer interfaces (K) + + + real(rk), public, dimension(:), allocatable, target :: ice_uvic_Tice ! ice layer temperature Tice(nilay +1)(deg-C) + + real(rk), public, dimension(:), allocatable, target :: ice_uvic_Cond ! thermal conductivities defined at the + ! centre of each layer Cond(nilay)(W m-1 K-1) + real(rk), public, dimension(:), allocatable, target :: ice_uvic_rhoCp ! volumetric heat capacities defined at + ! the centre of each layer rhoCp(nilay)(J m-3 K-1) + real(rk), public, dimension(:), allocatable, target :: ice_uvic_Sint ! internal heat source due to penetrating + ! short wave radiation Sint(nilay)(W m-3) + real(rk), public, dimension(:), allocatable, target :: ice_uvic_dzi ! layer thicknesses dzi(nilay)(m) + real(rk), public, dimension(:), allocatable, target :: ice_uvic_zi ! layer interface depths zi(nilay+1)(m) + real(rk), public, dimension(:), allocatable, target :: ice_uvic_Told ! ice temperature two time steps + ! previous to calculation of outgoing + ! long-wave flux in SEBUDGET Told (nilay+1) (K) + real(rk), public, dimension(:), allocatable, target :: ice_uvic_Pari !photosynthatically available radiation in ice (W m-2) + + real(rk), public, dimension(:), allocatable :: ice_uvic_dum + real(rk), public, dimension(:), allocatable :: ice_uvic_dzice + real(rk), public, dimension(:), allocatable :: ice_uvic_zice + +!------------------------------------------------------------------------------------------------------------- ! ! !REVISION HISTORY: ! Original author(s): Karsten Bolding diff --git a/tests/test_stim.F90 b/tests/test_stim.F90 index c7468eb..74e44fa 100644 --- a/tests/test_stim.F90 +++ b/tests/test_stim.F90 @@ -34,6 +34,11 @@ program test_stim #else write(*,*) 'Winton: off' #endif +#ifdef STIM_FLATO + write(*,*) 'Flato: on' +#else + write(*,*) 'Flato: off' +#endif end program test_stim !EOC