21 implicit none ;
private
23 #include <MOM_memory.h>
30 real :: berg_area_threshold
33 real :: latent_heat_fusion
34 real :: density_iceberg
36 type(time_type),
pointer :: time
49 type(
surface),
intent(inout) :: sfc_state
51 logical,
intent(in) :: use_ice_shelf
52 real,
intent(in) :: time_step
56 integer :: i, j, is, ie, js, je
57 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
63 if (.not.
associated(cs))
return
65 if (.not.(
associated(forces%area_berg) .and.
associated(forces%mass_berg) ) )
return
67 if (.not.(
associated(forces%frac_shelf_u) .and.
associated(forces%frac_shelf_v) .and. &
68 associated(forces%rigidity_ice_u) .and.
associated(forces%rigidity_ice_v)) )
return
71 if (.not. use_ice_shelf)
then
72 forces%frac_shelf_u(:,:) = 0.0 ; forces%frac_shelf_v(:,:) = 0.0
74 if (.not. forces%accumulate_rigidity)
then
75 forces%rigidity_ice_u(:,:) = 0.0 ; forces%rigidity_ice_v(:,:) = 0.0
78 call pass_var(forces%area_berg, g%domain,
to_all+omit_corners, halo=1, complete=.false.)
79 call pass_var(forces%mass_berg, g%domain,
to_all+omit_corners, halo=1, complete=.true.)
80 kv_rho_ice = cs%kv_iceberg / cs%density_iceberg
81 do j=js,je ;
do i=is-1,ie
82 if ((g%areaT(i,j) + g%areaT(i+1,j) > 0.0)) &
83 forces%frac_shelf_u(i,j) = forces%frac_shelf_u(i,j) + &
84 (((forces%area_berg(i,j)*g%US%L_to_m**2*g%areaT(i,j)) + &
85 (forces%area_berg(i+1,j)*g%US%L_to_m**2*g%areaT(i+1,j))) / &
86 (g%US%L_to_m**2*g%areaT(i,j) + g%US%L_to_m**2*g%areaT(i+1,j)) )
87 forces%rigidity_ice_u(i,j) = forces%rigidity_ice_u(i,j) + kv_rho_ice * &
88 min(forces%mass_berg(i,j), forces%mass_berg(i+1,j))
90 do j=js-1,je ;
do i=is,ie
91 if ((g%areaT(i,j) + g%areaT(i,j+1) > 0.0)) &
92 forces%frac_shelf_v(i,j) = forces%frac_shelf_v(i,j) + &
93 (((forces%area_berg(i,j)*g%US%L_to_m**2*g%areaT(i,j)) + &
94 (forces%area_berg(i,j+1)*g%US%L_to_m**2*g%areaT(i,j+1))) / &
95 (g%US%L_to_m**2*g%areaT(i,j) + g%US%L_to_m**2*g%areaT(i,j+1)) )
96 forces%rigidity_ice_v(i,j) = forces%rigidity_ice_v(i,j) + kv_rho_ice * &
97 min(forces%mass_berg(i,j), forces%mass_berg(i,j+1))
100 call pass_vector(forces%frac_shelf_u, forces%frac_shelf_v, g%domain,
to_all, cgrid_ne)
106 subroutine iceberg_fluxes(G, US, fluxes, use_ice_shelf, sfc_state, &
110 type(
forcing),
intent(inout) :: fluxes
112 type(
surface),
intent(inout) :: sfc_state
114 logical,
intent(in) :: use_ice_shelf
115 real,
intent(in) :: time_step
120 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
121 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
122 isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
128 if (.not.
associated(cs))
return
129 if (.not.(
associated(fluxes%area_berg) .and.
associated(fluxes%ustar_berg) .and. &
130 associated(fluxes%mass_berg) ) )
return
131 if (.not.(
associated(fluxes%frac_shelf_h) .and.
associated(fluxes%ustar_shelf)) )
return
134 if (.not.(
associated(fluxes%area_berg) .and.
associated(fluxes%ustar_berg) .and. &
135 associated(fluxes%mass_berg) ) )
return
136 if (.not. use_ice_shelf)
then
137 fluxes%frac_shelf_h(:,:) = 0.
138 fluxes%ustar_shelf(:,:) = 0.
140 do j=jsd,jed ;
do i=isd,ied ;
if (g%areaT(i,j) > 0.0)
then
141 fluxes%frac_shelf_h(i,j) = fluxes%frac_shelf_h(i,j) + fluxes%area_berg(i,j)
142 fluxes%ustar_shelf(i,j) = fluxes%ustar_shelf(i,j) + fluxes%ustar_berg(i,j)
143 endif ;
enddo ;
enddo
146 if (cs%berg_area_threshold >= 0.)
then
147 i_dt_lhf = 1.0 / (us%s_to_T*time_step * cs%latent_heat_fusion)
148 do j=jsd,jed ;
do i=isd,ied
149 if (fluxes%frac_shelf_h(i,j) > cs%berg_area_threshold)
then
152 if (
associated(fluxes%sw)) fluxes%sw(i,j) = 0.0
153 if (
associated(fluxes%lw)) fluxes%lw(i,j) = 0.0
154 if (
associated(fluxes%latent)) fluxes%latent(i,j) = 0.0
155 if (
associated(fluxes%evap)) fluxes%evap(i,j) = 0.0
161 if (
associated(sfc_state%frazil))
then
162 fraz = us%kg_m3_to_R*us%m_to_Z*sfc_state%frazil(i,j) * i_dt_lhf
163 if (
associated(fluxes%evap)) &
164 fluxes%evap(i,j) = fluxes%evap(i,j) - fraz
166 sfc_state%frazil(i,j) = 0.0
170 if (
associated(fluxes%sens)) fluxes%sens(i,j) = 0.0
171 if (
associated(fluxes%salt_flux)) fluxes%salt_flux(i,j) = 0.0
172 if (
associated(fluxes%lprec)) fluxes%lprec(i,j) = 0.0
181 type(time_type),
target,
intent(in) :: time
184 type(
diag_ctrl),
target,
intent(inout) :: diag
187 #include "version_variable.h"
188 character(len=40) :: mdl =
"MOM_marine_ice"
190 if (
associated(cs))
then
191 call mom_error(warning,
"marine_ice_init called with an associated control structure.")
193 else ;
allocate(cs) ;
endif
198 call get_param(param_file, mdl,
"KV_ICEBERG", cs%kv_iceberg, &
199 "The viscosity of the icebergs", units=
"m2 s-1",default=1.0e10)
200 call get_param(param_file, mdl,
"DENSITY_ICEBERGS", cs%density_iceberg, &
201 "A typical density of icebergs.", units=
"kg m-3", default=917.0)
202 call get_param(param_file, mdl,
"LATENT_HEAT_FUSION", cs%latent_heat_fusion, &
203 "The latent heat of fusion.", units=
"J/kg", default=hlf)
204 call get_param(param_file, mdl,
"BERG_AREA_THRESHOLD", cs%berg_area_threshold, &
205 "Fraction of grid cell which iceberg must occupy, so that fluxes "//&
206 "below berg are set to zero. Not applied for negative "//&
207 "values.", units=
"non-dim", default=-1.0)