20 implicit none ;
private
22 #include <MOM_memory.h>
24 public iceberg_forces, iceberg_fluxes, marine_ice_init
29 real :: berg_area_threshold
32 real :: latent_heat_fusion
33 real :: density_iceberg
35 type(time_type),
pointer :: time
44 subroutine iceberg_forces(G, forces, use_ice_shelf, sfc_state, &
48 type(
surface),
intent(inout) :: sfc_state
50 logical,
intent(in) :: use_ice_shelf
51 real,
intent(in) :: time_step
55 integer :: i, j, is, ie, js, je
56 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
62 if (.not.
associated(cs))
return
64 if (.not.(
associated(forces%area_berg) .and.
associated(forces%mass_berg) ) )
return
66 if (.not.(
associated(forces%frac_shelf_u) .and.
associated(forces%frac_shelf_v) .and. &
67 associated(forces%rigidity_ice_u) .and.
associated(forces%rigidity_ice_v)) )
return
70 if (.not. use_ice_shelf)
then
71 forces%frac_shelf_u(:,:) = 0.0 ; forces%frac_shelf_v(:,:) = 0.0
73 if (.not. forces%accumulate_rigidity)
then
74 forces%rigidity_ice_u(:,:) = 0.0 ; forces%rigidity_ice_v(:,:) = 0.0
77 call pass_var(forces%area_berg, g%domain, to_all+omit_corners, halo=1, complete=.false.)
78 call pass_var(forces%mass_berg, g%domain, to_all+omit_corners, halo=1, complete=.true.)
79 kv_rho_ice = cs%kv_iceberg / cs%density_iceberg
80 do j=js,je ;
do i=is-1,ie
81 if ((g%areaT(i,j) + g%areaT(i+1,j) > 0.0)) &
82 forces%frac_shelf_u(i,j) = forces%frac_shelf_u(i,j) + &
83 (((forces%area_berg(i,j)*g%areaT(i,j)) + &
84 (forces%area_berg(i+1,j)*g%areaT(i+1,j))) / &
85 (g%areaT(i,j) + g%areaT(i+1,j)) )
86 forces%rigidity_ice_u(i,j) = forces%rigidity_ice_u(i,j) + kv_rho_ice * &
87 min(forces%mass_berg(i,j), forces%mass_berg(i+1,j))
89 do j=js-1,je ;
do i=is,ie
90 if ((g%areaT(i,j) + g%areaT(i,j+1) > 0.0)) &
91 forces%frac_shelf_v(i,j) = forces%frac_shelf_v(i,j) + &
92 (((forces%area_berg(i,j)*g%areaT(i,j)) + &
93 (forces%area_berg(i,j+1)*g%areaT(i,j+1))) / &
94 (g%areaT(i,j) + g%areaT(i,j+1)) )
95 forces%rigidity_ice_v(i,j) = forces%rigidity_ice_v(i,j) + kv_rho_ice * &
96 min(forces%mass_berg(i,j), forces%mass_berg(i,j+1))
99 call pass_vector(forces%frac_shelf_u, forces%frac_shelf_v, g%domain, to_all, cgrid_ne)
101 end subroutine iceberg_forces
105 subroutine iceberg_fluxes(G, fluxes, use_ice_shelf, sfc_state, &
108 type(
forcing),
intent(inout) :: fluxes
110 type(
surface),
intent(inout) :: sfc_state
112 logical,
intent(in) :: use_ice_shelf
113 real,
intent(in) :: time_step
118 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
119 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
120 isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
126 if (.not.
associated(cs))
return
127 if (.not.(
associated(fluxes%area_berg) .and.
associated(fluxes%ustar_berg) .and. &
128 associated(fluxes%mass_berg) ) )
return
129 if (.not.(
associated(fluxes%frac_shelf_h) .and.
associated(fluxes%ustar_shelf)) )
return
132 if (.not.(
associated(fluxes%area_berg) .and.
associated(fluxes%ustar_berg) .and. &
133 associated(fluxes%mass_berg) ) )
return
134 if (.not. use_ice_shelf)
then
135 fluxes%frac_shelf_h(:,:) = 0.
136 fluxes%ustar_shelf(:,:) = 0.
138 do j=jsd,jed ;
do i=isd,ied ;
if (g%areaT(i,j) > 0.0)
then
139 fluxes%frac_shelf_h(i,j) = fluxes%frac_shelf_h(i,j) + fluxes%area_berg(i,j)
140 fluxes%ustar_shelf(i,j) = fluxes%ustar_shelf(i,j) + fluxes%ustar_berg(i,j)
141 endif ;
enddo ;
enddo
144 if (cs%berg_area_threshold >= 0.)
then
145 i_dt_lhf = 1.0 / (time_step * cs%latent_heat_fusion)
146 do j=jsd,jed ;
do i=isd,ied
147 if (fluxes%frac_shelf_h(i,j) > cs%berg_area_threshold)
then
150 if (
associated(fluxes%sw)) fluxes%sw(i,j) = 0.0
151 if (
associated(fluxes%lw)) fluxes%lw(i,j) = 0.0
152 if (
associated(fluxes%latent)) fluxes%latent(i,j) = 0.0
153 if (
associated(fluxes%evap)) fluxes%evap(i,j) = 0.0
159 if (
associated(sfc_state%frazil))
then
160 fraz = sfc_state%frazil(i,j) * i_dt_lhf
161 if (
associated(fluxes%evap)) fluxes%evap(i,j) = fluxes%evap(i,j) - fraz
163 sfc_state%frazil(i,j) = 0.0
167 if (
associated(fluxes%sens)) fluxes%sens(i,j) = 0.0
168 if (
associated(fluxes%salt_flux)) fluxes%salt_flux(i,j) = 0.0
169 if (
associated(fluxes%lprec)) fluxes%lprec(i,j) = 0.0
174 end subroutine iceberg_fluxes
177 subroutine marine_ice_init(Time, G, param_file, diag, CS)
178 type(time_type),
target,
intent(in) :: time
181 type(
diag_ctrl),
target,
intent(inout) :: diag
184 #include "version_variable.h"
185 character(len=40) :: mdl =
"MOM_marine_ice"
187 if (
associated(cs))
then
188 call mom_error(warning,
"marine_ice_init called with an associated control structure.")
190 else ;
allocate(cs) ;
endif
195 call get_param(param_file, mdl,
"KV_ICEBERG", cs%kv_iceberg, &
196 "The viscosity of the icebergs", units=
"m2 s-1",default=1.0e10)
197 call get_param(param_file, mdl,
"DENSITY_ICEBERGS", cs%density_iceberg, &
198 "A typical density of icebergs.", units=
"kg m-3", default=917.0)
199 call get_param(param_file, mdl,
"LATENT_HEAT_FUSION", cs%latent_heat_fusion, &
200 "The latent heat of fusion.", units=
"J/kg", default=hlf)
201 call get_param(param_file, mdl,
"BERG_AREA_THRESHOLD", cs%berg_area_threshold, &
202 "Fraction of grid cell which iceberg must occupy, so that fluxes "//&
203 "below berg are set to zero. Not applied for negative "//&
204 "values.", units=
"non-dim", default=-1.0)
206 end subroutine marine_ice_init