16 implicit none ;
private
18 #include <MOM_memory.h>
20 public user_initialize_shelf_mass, user_update_shelf_mass
21 public user_init_ice_thickness
33 real :: flat_shelf_width
34 real :: shelf_slope_scale
35 real :: pos_shelf_edge_0
37 logical :: first_call = .true.
43 subroutine user_initialize_shelf_mass(mass_shelf, area_shelf_h, h_shelf, hmask, G, US, CS, param_file, new_sim)
46 real,
dimension(SZDI_(G),SZDJ_(G)), &
47 intent(out) :: mass_shelf
49 real,
dimension(SZDI_(G),SZDJ_(G)), &
50 intent(out) :: h_shelf
51 real,
dimension(SZDI_(G),SZDJ_(G)), &
52 intent(out) :: area_shelf_h
53 real,
dimension(SZDI_(G),SZDJ_(G)), &
59 logical,
intent(in) :: new_sim
65 real :: flat_shelf_width
67 character(len=40) :: mdl =
"USER_initialize_shelf_mass"
73 if (.not.
associated(cs))
allocate(cs)
76 if (cs%first_call)
call write_user_log(param_file)
77 cs%first_call = .false.
78 call get_param(param_file, mdl,
"RHO_0", cs%Rho_ocean, &
79 "The mean ocean density used with BOUSSINESQ true to "//&
80 "calculate accelerations and the mass for conservation "//&
81 "properties, or with BOUSSINSEQ false to convert some "//&
82 "parameters from vertical units of m to kg m-2.", &
83 units=
"kg m-3", default=1035.0, scale=us%kg_m3_to_R)
84 call get_param(param_file, mdl,
"SHELF_MAX_DRAFT", cs%max_draft, &
85 units=
"m", default=1.0, scale=us%m_to_Z)
86 call get_param(param_file, mdl,
"SHELF_MIN_DRAFT", cs%min_draft, &
87 units=
"m", default=1.0, scale=us%m_to_Z)
88 call get_param(param_file, mdl,
"FLAT_SHELF_WIDTH", cs%flat_shelf_width, &
89 units=
"axis_units", default=0.0)
90 call get_param(param_file, mdl,
"SHELF_SLOPE_SCALE", cs%shelf_slope_scale, &
91 units=
"axis_units", default=0.0)
92 call get_param(param_file, mdl,
"SHELF_EDGE_POS_0", cs%pos_shelf_edge_0, &
93 units=
"axis_units", default=0.0)
94 call get_param(param_file, mdl,
"SHELF_SPEED", cs%shelf_speed, &
95 units=
"axis_units day-1", default=0.0)
97 call user_update_shelf_mass(mass_shelf, area_shelf_h, h_shelf, hmask, g, cs, set_time(0,0), new_sim)
99 end subroutine user_initialize_shelf_mass
102 subroutine user_init_ice_thickness(h_shelf, area_shelf_h, hmask, G, US, param_file)
104 real,
dimension(SZDI_(G),SZDJ_(G)), &
105 intent(out) :: h_shelf
106 real,
dimension(SZDI_(G),SZDJ_(G)), &
107 intent(out) :: area_shelf_h
108 real,
dimension(SZDI_(G),SZDJ_(G)), &
116 real,
dimension(SZI_(G),SZJ_(G)) :: mass_shelf
119 call user_initialize_shelf_mass(mass_shelf, area_shelf_h, h_shelf, hmask, g, us, cs, param_file, .true.)
121 end subroutine user_init_ice_thickness
124 subroutine user_update_shelf_mass(mass_shelf, area_shelf_h, h_shelf, hmask, G, CS, Time, new_sim)
126 real,
dimension(SZDI_(G),SZDJ_(G)), &
127 intent(inout) :: mass_shelf
129 real,
dimension(SZDI_(G),SZDJ_(G)), &
130 intent(inout) :: area_shelf_h
131 real,
dimension(SZDI_(G),SZDJ_(G)), &
132 intent(inout) :: h_shelf
133 real,
dimension(SZDI_(G),SZDJ_(G)), &
134 intent(inout) :: hmask
137 type(time_type),
intent(in) :: time
138 logical,
intent(in) :: new_sim
141 real :: c1, edge_pos, slope_pos
144 edge_pos = cs%pos_shelf_edge_0 + cs%shelf_speed*(time_type_to_real(time) / 86400.0)
146 slope_pos = edge_pos - cs%flat_shelf_width
147 c1 = 0.0 ;
if (cs%shelf_slope_scale > 0.0) c1 = 1.0 / cs%shelf_slope_scale
152 if (((j+g%jdg_offset) <= g%domain%njglobal+g%domain%njhalo) .AND. &
153 ((j+g%jdg_offset) >= g%domain%njhalo+1))
then
160 if ((j >= g%jsc) .and. (j <= g%jec))
then
162 if (new_sim)
then ;
if (g%geoLonCu(i-1,j) >= edge_pos)
then
164 mass_shelf(i,j) = 0.0
165 area_shelf_h(i,j) = 0.0
169 if (g%geoLonCu(i,j) > edge_pos)
then
170 area_shelf_h(i,j) = g%areaT(i,j) * (edge_pos - g%geoLonCu(i-1,j)) / &
171 (g%geoLonCu(i,j) - g%geoLonCu(i-1,j))
174 area_shelf_h(i,j) = g%areaT(i,j)
178 if (g%geoLonT(i,j) > slope_pos)
then
179 h_shelf(i,j) = cs%min_draft
180 mass_shelf(i,j) = cs%Rho_ocean * cs%min_draft
182 mass_shelf(i,j) = cs%Rho_ocean * (cs%min_draft + &
183 (cs%max_draft - cs%min_draft) * &
184 min(1.0, (c1*(slope_pos - g%geoLonT(i,j)))**2) )
185 h_shelf(i,j) = (cs%min_draft + &
186 (cs%max_draft - cs%min_draft) * &
187 min(1.0, (c1*(slope_pos - g%geoLonT(i,j)))**2) )
190 endif ;
endif ;
endif
192 if ((i+g%idg_offset) == g%domain%nihalo+1)
then
196 enddo ;
endif ;
enddo
198 end subroutine user_update_shelf_mass
201 subroutine write_user_log(param_file)
204 character(len=128) :: version =
'$Id: user_shelf_init.F90,v 1.1.2.7 2012/06/19 22:15:52 Robert.Hallberg Exp $'
205 character(len=128) :: tagname =
'$Name: MOM_ogrp $'
206 character(len=40) :: mdl =
"user_shelf_init"
208 call log_version(param_file, mdl, version, tagname)
210 end subroutine write_user_log