Skip to content

Commit

Permalink
Merge branch 'user/mjh/SIS2_nudging' into dev/master
Browse files Browse the repository at this point in the history
  • Loading branch information
@adcroft
adcroft committed Oct 9, 2015
2 parents 0a93815 + aad727f commit 5b77bfc
Show file tree
Hide file tree
Showing 2 changed files with 103 additions and 5 deletions.
95 changes: 91 additions & 4 deletions ice_model.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,8 +71,9 @@ module ice_model_mod
use MOM_time_manager, only : operator(>), operator(*), operator(/), operator(/=)
use astronomy_mod, only: astronomy_init, astronomy_end
use astronomy_mod, only: universal_time, orbital_time, diurnal_solar, daily_mean_solar
use coupler_types_mod,only: coupler_3d_bc_type
use constants_mod, only: hlv, hlf, T_0degC=>Tfreeze, grav, STEFAN
use coupler_types_mod,only: coupler_3d_bc_type
use constants_mod, only: hlv, hlf, T_0degC=>Tfreeze, grav, STEFAN
use data_override_mod, only : data_override
use ocean_albedo_mod, only: compute_ocean_albedo ! ice sets ocean surface
use ocean_rough_mod, only: compute_ocean_roughness ! properties over water

Expand All @@ -92,7 +93,7 @@ module ice_model_mod
use SIS_tracer_registry, only : register_SIS_tracer

use ice_thm_mod, only: slab_ice_optics, ice_thm_param, ice5lay_temp, ice5lay_resize
use ice_thm_mod, only: MU_TS, TFI, CI, e_to_melt, get_thermo_coefs
use ice_thm_mod, only: MU_TS, TFI, CI, e_to_melt, get_thermo_coefs
use SIS2_ice_thm, only: ice_temp_SIS2, ice_resize_SIS2, SIS2_ice_thm_init, SIS2_ice_thm_end
use SIS2_ice_thm, only: ice_optics_SIS2, get_SIS2_thermo_coefs, enthalpy_liquid_freeze
use SIS2_ice_thm, only: enthalpy_from_TS, enth_from_TS, Temp_from_En_S, Temp_from_Enth_S
Expand Down Expand Up @@ -721,6 +722,7 @@ subroutine set_ice_surface_state(Ice, IST, t_surf_ice_bot, u_surf_ice_bot, v_sur

real, dimension(G%isc:G%iec,G%jsc:G%jec) :: m_ice_tot
real, dimension(G%isc:G%iec,G%jsc:G%jec) :: h_ice_input
real, dimension(G%isc:G%iec,G%jsc:G%jec) :: icec, icec_obs
real, dimension(SZI_(G),SZJ_(G)) :: u_nonsym, v_nonsym
real, dimension(G%NkIce) :: sw_abs_lay
real :: u, v
Expand All @@ -731,7 +733,9 @@ subroutine set_ice_surface_state(Ice, IST, t_surf_ice_bot, u_surf_ice_bot, v_sur
logical :: sent
real :: H_to_m_ice ! The specific volumes of ice and snow times the
real :: H_to_m_snow ! conversion factor from thickness units, in m H-1.
real, parameter :: LI = hlf
real :: cp_inv, drho_dT,drho_dS
real, parameter :: LI = hlf

isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; ncat = G%CatIce
i_off = LBOUND(Ice%t_surf,1) - G%isc ; j_off = LBOUND(Ice%t_surf,2) - G%jsc
I_Nk = 1.0 / G%NkIce ; kg_H_Nk = G%H_to_kg_m2 * I_Nk
Expand Down Expand Up @@ -784,6 +788,26 @@ subroutine set_ice_surface_state(Ice, IST, t_surf_ice_bot, u_surf_ice_bot, v_sur
call chksum(v_surf_ice_bot(isc:iec,jsc:jec), "Start IB2IT v_surf_ice_bot")
endif

if (IST%nudge_sea_ice) then
cp_inv=1.0/CI
IST%frazil_nudge(isc:iec,jsc:jec)=0.0
IST%melt_nudge(isc:iec,jsc:jec)=0.0
call data_override('ICE','icec',icec_obs,Ice%Time)
icec = sum(IST%part_size(isc:iec,jsc:jec,2:),dim=3)
do j = jsc,jec
do i = isc,iec
if (icec(i,j) < icec_obs(i,j)) then
IST%frazil_nudge(i,j) = hlf*(1.0-exp(-(icec_obs(i,j)-icec(i,j))**2.0))*IST%nudge_sea_ice_coeff ! J/m2
if (IST%t_ocn(i,j) > TFI .and. IST%s_surf(i,j) > 1.0 ) then
call calculate_density_derivs_wright(IST%t_ocn(i,j),IST%s_surf(i,j),0.0,drho_dT,drho_dS)
IST%melt_nudge(i,j) = -IST%frazil_nudge(i,j)*drho_dT/drho_dS*cp_inv/IST%s_surf(i,j)
endif
IST%frazil(i,j) = IST%frazil(i,j) + IST%frazil_nudge(i,j)
endif
enddo
enddo
endif

! Transfer the ocean state for extra tracer fluxes.
do n=1,OIB%fields%num_bcs ; do m=1,OIB%fields%bc(n)%num_fields
Ice%ocean_fields%bc(n)%field(m)%values(:,:,1) = OIB%fields%bc(n)%field(m)%values
Expand Down Expand Up @@ -1645,6 +1669,9 @@ subroutine update_ice_model_slow(Ice, IST, G, runoff, calving, &
if (IST%id_frazil>0) &
call post_data(IST%id_frazil, IST%frazil*Idt_slow, IST%diag, mask=G%Lmask2dT)

if (IST%nudge_sea_ice .and. IST%id_fwnudge>0) &
call post_data(IST%id_fwnudge, IST%melt_nudge, IST%diag, mask=Ice%mask )

call avg_top_quantities(Ice, IST, G) ! average fluxes from update_ice_model_fast

!$OMP parallel do default(none) shared(isc,iec,jsc,jec,IST)
Expand Down Expand Up @@ -2631,6 +2658,12 @@ subroutine SIS1_5L_thermodynamics(Ice, IST, G) !, runoff, calving, &
IST%lprec_ocn_top(i,j) = IST%lprec_ocn_top(i,j) + &
(h2o_to_ocn*IST%part_size(i,j,k)) / dt_slow


if (IST%nudge_sea_ice) then
IST%lprec_top(i,j,:) = IST%lprec_top(i,j,:) + IST%melt_nudge(i,j)*IST%part_size(i,j,k)/dt_slow
IST%lprec_ocn_top(i,j) = IST%lprec_ocn_top(i,j) + IST%melt_nudge(i,j)*Ice%part_size(i,j,k)/dt_slow

This comment has been minimized.

Sign in to view

Copy link
@nikizadehgfdl

nikizadehgfdl Dec 31, 2015

Contributor

Just double checking, is "Ice%" correct on line 2664? If so, "Ice" has to be added to line 2585 (OMP shared directive ) for sis2 to compile with -openmp.
endif

IST%enth_snow(i,j,k,1) = enth_from_TS(T_col(0), 0.0, IST%ITV)
do m=1,NkIce ; IST%enth_ice(i,j,k,m) = enth_from_TS(T_col(m), S_col(m), IST%ITV) ; enddo

Expand Down Expand Up @@ -3677,6 +3710,12 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow )
call get_param(param_file, mod, "MAX_ICE_THICK_LIMIT", IST%max_ice_limit, &
"The maximum permitted sea ice thickness when \n"//&
"APPLY_ICE_LIMIT is true.", units="m", default=4.0)


call get_param(param_file, mod, "NUDGE_SEA_ICE", IST%nudge_sea_ice, &
"If true, constrain the sea ice concentrations using observations.", &
default=.false.)

call get_param(param_file, mod, "APPLY_SLP_TO_OCEAN", IST%slp2ocean, &
"If true, apply the atmospheric sea level pressure to \n"//&
"the ocean.", default=.false.)
Expand Down Expand Up @@ -3768,6 +3807,19 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow )
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec
i_off = LBOUND(Ice%t_surf,1) - G%isc ; j_off = LBOUND(Ice%t_surf,2) - G%jsc

if (IST%nudge_sea_ice) then
allocate(IST%frazil_nudge(isc:iec,jsc:jec)); IST%frazil_nudge(:,:)=0.0
allocate(IST%melt_nudge(isc:iec,jsc:jec)); IST%melt_nudge(:,:)=0.0


call get_param(param_file, mod, "NUDGE_SEA_ICE_COEFF", IST%nudge_sea_ice_coeff, &
"dimensional coefficient controls how strongly sea ice \n"//&
"is constrained to observations. A suggested value is 1.e2.", &
units = "kg m-2", default=0.0)

endif


IST%coszen(:,:) = cos(3.14*67.0/180.0) ! NP summer solstice.

do j=jsc,jec ; do i=isc,iec ; i2 = i+i_off ; j2 = j+j_off
Expand Down Expand Up @@ -4135,4 +4187,39 @@ subroutine ice_aging(G, mi, age, mi_old, dt)

end subroutine ice_aging

!> Calculates the derivatives of seawater density with respect to potential
!! temperature and salinity
subroutine calculate_density_derivs_wright(T, S, pressure, drho_dT, drho_dS)
real, intent(in) :: T !< Potential temperature relative to the surface in C
real, intent(in) :: S !< Salinity in PSU
real, intent(in) :: pressure !< Pressure in Pa
real, intent(out) :: drho_dT !< Partial derivative of density with potential
!! tempetature, in kg m-3 K-1
real, intent(out) :: drho_dS !< Partial derivative of density with salinity,
!! in kg m-3 psu-1
! Local variables
real :: al0, p0, lambda, I_denom2
integer :: j
! Following are the values for the reduced range formula.
real, parameter :: a0 = 7.057924e-4, a1 = 3.480336e-7, a2 = -1.112733e-7
real, parameter :: b0 = 5.790749e8, b1 = 3.516535e6, b2 = -4.002714e4
real, parameter :: b3 = 2.084372e2, b4 = 5.944068e5, b5 = -9.643486e3
real, parameter :: c0 = 1.704853e5, c1 = 7.904722e2, c2 = -7.984422
real, parameter :: c3 = 5.140652e-2, c4 = -2.302158e2, c5 = -3.079464

al0 = a0 + a1*T + a2*S
p0 = b0 + b4*S + T * (b1 + T*(b2 + b3*T) + b5*S)
lambda = c0 +c4*S + T * (c1 + T*(c2 + c3*T) + c5*S)

I_denom2 = 1.0 / (lambda + al0*(pressure + p0))
I_denom2 = I_denom2 *I_denom2
drho_dT = I_denom2 * &
(lambda* (b1 + T*(2.0*b2 + 3.0*b3*T) + b5*S) - &
(pressure+p0) * ( (pressure+p0)*a1 + &
(c1 + T*(c2*2.0 + c3*3.0*T) + c5*S) ))
drho_dS = I_denom2 * (lambda* (b4 + b5*T) - &
(pressure+p0) * ( (pressure+p0)*a2 + (c4 + c5*T) ))

end subroutine calculate_density_derivs_wright

end module ice_model_mod
13 changes: 12 additions & 1 deletion ice_type.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -186,6 +186,10 @@ module ice_type_mod

real, pointer, dimension(:,:) :: frazil =>NULL()
real, pointer, dimension(:,:) :: frazil_input =>NULL()

real, pointer, dimension(:,:) :: frazil_nudge =>NULL()
real, pointer, dimension(:,:) :: melt_nudge =>NULL()

real, pointer, dimension(:,:) :: bheat =>NULL()
real, pointer, dimension(:,:) :: mi =>NULL() ! The total ice+snow mass, in kg m-2.
logical :: slab_ice ! If true, do the old style GFDL slab ice.
Expand Down Expand Up @@ -259,6 +263,11 @@ module ice_type_mod
! that needs to be done.
logical :: first_time = .true. ! If true, this is the first call to
! update_ice_model_slow_up
logical :: nudge_sea_ice = .false. ! If true, nudge sea ice concentrations towards observations.
real :: nudge_sea_ice_coeff = 0.0 ! Dimensional coefficient controls how strongly sea ice
! is constrained to observations. Units are kg m-2. A suggested value
! is 1.e2

integer :: num_tr_fluxes = -1 ! The number of tracer flux fields
integer, allocatable, dimension(:,:) :: tr_flux_index
real, allocatable, dimension(:,:,:,:) :: tr_flux_top
Expand All @@ -273,7 +282,7 @@ module ice_type_mod
integer :: id_lh=-1, id_sw=-1, id_lw=-1, id_snofl=-1, id_rain=-1, id_runoff=-1
integer :: id_calving=-1, id_runoff_hflx=-1, id_calving_hflx=-1, id_evap=-1
integer :: id_saltf=-1, id_tmelt=-1, id_bmelt=-1, id_bheat=-1, id_e2m=-1
integer :: id_rdgr=-1,id_rdgf=-1,id_rdgo=-1,id_rdgv=-1,id_age=-1
integer :: id_rdgr=-1, id_rdgf=-1, id_rdgo=-1, id_rdgv=-1, id_age=-1, id_fwnudge=-1
integer :: id_frazil=-1, id_alb=-1, id_xprt=-1, id_lsrc=-1, id_lsnk=-1, id_bsnk=-1
integer :: id_strna=-1, id_fax=-1, id_fay=-1, id_swdn=-1, id_lwdn=-1, id_sn2ic=-1
integer :: id_slp=-1, id_ext=-1, id_sst=-1, id_sss=-1, id_ssh=-1, id_uo=-1, id_vo=-1
Expand Down Expand Up @@ -1044,6 +1053,8 @@ subroutine ice_diagnostics_init(Ice, IST, G, diag, Time)
'rate of rain fall', 'kg/(m^2*s)', missing_value=missing)
IST%id_runoff = register_SIS_diag_field('ice_model','RUNOFF' ,diag%axesT1, Time, &
'liquid runoff', 'kg/(m^2*s)', missing_value=missing)
IST%id_fwnudge = register_SIS_diag_field('ice_model','FW_NUDGE' ,diag%axesT1, Time, &
'nudging freshwater flux', 'kg/(m^2*s)', missing_value=missing)
IST%id_calving = register_SIS_diag_field('ice_model','CALVING',diag%axesT1, Time, &
'frozen runoff', 'kg/(m^2*s)', missing_value=missing)
IST%id_runoff_hflx = register_SIS_diag_field('ice_model','RUNOFF_HFLX' ,diag%axesT1, Time, &
Expand Down

0 comments on commit 5b77bfc

Please sign in to comment.