21 implicit none ;
private
28 logical :: use_temperature
29 logical :: restorebuoy
35 real,
dimension(:,:),
pointer :: &
36 t_restore(:,:) => null(), &
37 s_restore(:,:) => null(), &
40 real,
dimension(:,:),
pointer :: heat(:,:) => null()
42 character(len=200) :: inputdir
43 character(len=200) :: salinityrestore_file
44 character(len=200) :: sstrestore_file
45 character(len=200) :: solar_file
46 character(len=200) :: heating_file
47 character(len=200) :: pme_file
61 type(
forcing),
intent(inout) :: fluxes
62 type(time_type),
intent(in) :: day
63 real,
intent(in) :: dt
82 real :: density_restore
85 real :: buoy_rest_const
88 integer :: i, j, is, ie, js, je
89 integer :: isd, ied, jsd, jed
91 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
92 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
99 if (cs%use_temperature)
then
100 call safe_alloc_ptr(fluxes%evap, isd, ied, jsd, jed)
101 call safe_alloc_ptr(fluxes%lprec, isd, ied, jsd, jed)
102 call safe_alloc_ptr(fluxes%fprec, isd, ied, jsd, jed)
103 call safe_alloc_ptr(fluxes%lrunoff, isd, ied, jsd, jed)
104 call safe_alloc_ptr(fluxes%frunoff, isd, ied, jsd, jed)
105 call safe_alloc_ptr(fluxes%vprec, isd, ied, jsd, jed)
107 call safe_alloc_ptr(fluxes%sw, isd, ied, jsd, jed)
108 call safe_alloc_ptr(fluxes%lw, isd, ied, jsd, jed)
109 call safe_alloc_ptr(fluxes%latent, isd, ied, jsd, jed)
110 call safe_alloc_ptr(fluxes%sens, isd, ied, jsd, jed)
111 call safe_alloc_ptr(fluxes%heat_content_lprec, isd, ied, jsd, jed)
113 call safe_alloc_ptr(fluxes%buoy, isd, ied, jsd, jed)
119 call safe_alloc_ptr(cs%T_Restore, isd, ied, jsd, jed)
120 call safe_alloc_ptr(cs%S_Restore, isd, ied, jsd, jed)
121 call safe_alloc_ptr(cs%Heat, isd, ied, jsd, jed)
122 call safe_alloc_ptr(cs%PmE, isd, ied, jsd, jed)
123 call safe_alloc_ptr(cs%Solar, isd, ied, jsd, jed)
125 call mom_read_data(trim(cs%inputdir)//trim(cs%SSTrestore_file),
"SST", &
126 cs%T_Restore(:,:), g%Domain)
127 call mom_read_data(trim(cs%inputdir)//trim(cs%salinityrestore_file),
"SAL", &
128 cs%S_Restore(:,:), g%Domain)
129 call mom_read_data(trim(cs%inputdir)//trim(cs%heating_file),
"Heat", &
130 cs%Heat(:,:), g%Domain)
131 call mom_read_data(trim(cs%inputdir)//trim(cs%PmE_file),
"PmE", &
132 cs%PmE(:,:), g%Domain)
133 call mom_read_data(trim(cs%inputdir)//trim(cs%Solar_file),
"NET_SOL", &
134 cs%Solar(:,:), g%Domain)
138 if ( cs%use_temperature )
then
141 do j=js,je ;
do i=is,ie
144 fluxes%evap(i,j) = -0.0 * g%mask2dT(i,j)
145 fluxes%lprec(i,j) = us%m_to_Z*us%T_to_s * cs%PmE(i,j) * cs%Rho0 * g%mask2dT(i,j)
148 fluxes%vprec(i,j) = 0.0
151 fluxes%lw(i,j) = 0.0 * g%mask2dT(i,j)
152 fluxes%latent(i,j) = 0.0 * g%mask2dT(i,j)
153 fluxes%sens(i,j) = cs%Heat(i,j) * g%mask2dT(i,j)
154 fluxes%sw(i,j) = cs%Solar(i,j) * g%mask2dT(i,j)
157 do j=js,je ;
do i=is,ie
160 fluxes%buoy(i,j) = 0.0 * g%mask2dT(i,j)
164 if (cs%restorebuoy)
then
165 if (cs%use_temperature)
then
166 call safe_alloc_ptr(fluxes%heat_added, isd, ied, jsd, jed)
172 rhoxcp = us%R_to_kg_m3*cs%Rho0 * fluxes%C_p
173 do j=js,je ;
do i=is,ie
176 if (g%mask2dT(i,j) > 0)
then
177 fluxes%heat_added(i,j) = g%mask2dT(i,j) * &
178 ((cs%T_Restore(i,j) - sfc_state%SST(i,j)) * rhoxcp * us%Z_to_m*us%s_to_T*cs%Flux_const)
179 fluxes%vprec(i,j) = - (cs%Rho0 * cs%Flux_const) * &
180 (cs%S_Restore(i,j) - sfc_state%SSS(i,j)) / &
181 (0.5*(sfc_state%SSS(i,j) + cs%S_Restore(i,j)))
183 fluxes%heat_added(i,j) = 0.0
184 fluxes%vprec(i,j) = 0.0
190 call mom_error(fatal,
"MESO_buoyancy_surface_forcing: " // &
191 "Buoyancy restoring used without modification." )
194 buoy_rest_const = -1.0 * (cs%G_Earth * cs%Flux_const) / cs%Rho0
195 do j=js,je ;
do i=is,ie
198 density_restore = 1030.0
200 fluxes%buoy(i,j) = g%mask2dT(i,j) * buoy_rest_const * &
201 us%kg_m3_to_R * (density_restore - sfc_state%sfc_density(i,j))
211 type(time_type),
intent(in) :: time
215 type(
diag_ctrl),
target,
intent(inout) :: diag
220 #include "version_variable.h"
221 character(len=40) :: mdl =
"MESO_surface_forcing"
223 if (
associated(cs))
then
224 call mom_error(warning,
"MESO_surface_forcing_init called with an associated "// &
225 "control structure.")
233 call get_param(param_file, mdl,
"ENABLE_THERMODYNAMICS", cs%use_temperature, &
234 "If true, Temperature and salinity are used as state "//&
235 "variables.", default=.true.)
237 call get_param(param_file, mdl,
"G_EARTH", cs%G_Earth, &
238 "The gravitational acceleration of the Earth.", &
239 units=
"m s-2", default = 9.80, scale=us%m_to_L**2*us%Z_to_m*us%T_to_s**2)
240 call get_param(param_file, mdl,
"RHO_0", cs%Rho0, &
241 "The mean ocean density used with BOUSSINESQ true to "//&
242 "calculate accelerations and the mass for conservation "//&
243 "properties, or with BOUSSINSEQ false to convert some "//&
244 "parameters from vertical units of m to kg m-2.", &
245 units=
"kg m-3", default=1035.0, scale=us%kg_m3_to_R)
246 call get_param(param_file, mdl,
"GUST_CONST", cs%gust_const, &
247 "The background gustiness in the winds.", units=
"Pa", &
250 call get_param(param_file, mdl,
"RESTOREBUOY", cs%restorebuoy, &
251 "If true, the buoyancy fluxes drive the model back "//&
252 "toward some specified surface state with a rate "//&
253 "given by FLUXCONST.", default= .false.)
255 if (cs%restorebuoy)
then
256 call get_param(param_file, mdl,
"FLUXCONST", cs%Flux_const, &
257 "The constant that relates the restoring surface fluxes "//&
258 "to the relative surface anomalies (akin to a piston "//&
259 "velocity). Note the non-MKS units.", &
260 units=
"m day-1", scale=us%m_to_Z/(86400.0*us%s_to_T), &
261 fail_if_missing=.true.)
263 call get_param(param_file, mdl,
"SSTRESTORE_FILE", cs%SSTrestore_file, &
264 "The file with the SST toward which to restore in "//&
265 "variable TEMP.", fail_if_missing=.true.)
266 call get_param(param_file, mdl,
"SALINITYRESTORE_FILE", cs%salinityrestore_file, &
267 "The file with the surface salinity toward which to "//&
268 "restore in variable SALT.", fail_if_missing=.true.)
269 call get_param(param_file, mdl,
"SENSIBLEHEAT_FILE", cs%heating_file, &
270 "The file with the non-shortwave heat flux in "//&
271 "variable Heat.", fail_if_missing=.true.)
272 call get_param(param_file, mdl,
"PRECIP_FILE", cs%PmE_file, &
273 "The file with the net precipiation minus evaporation "//&
274 "in variable PmE.", fail_if_missing=.true.)
275 call get_param(param_file, mdl,
"SHORTWAVE_FILE", cs%Solar_file, &
276 "The file with the shortwave heat flux in "//&
277 "variable NET_SOL.", fail_if_missing=.true.)
278 call get_param(param_file, mdl,
"INPUTDIR", cs%inputdir, default=
".")
279 cs%inputdir = slasher(cs%inputdir)