7 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
12 use mom_domains,
only : global_field_sum, bitwise_exact_sum
14 use mom_domains,
only : to_north, to_east, omit_corners
33 use coupler_types_mod,
only : coupler_2d_bc_type, coupler_type_write_chksums
34 use coupler_types_mod,
only : coupler_type_initialized, coupler_type_spawn
35 use coupler_types_mod,
only : coupler_type_copy_data
36 use data_override_mod,
only : data_override_init, data_override
37 use fms_mod,
only : stdout
38 use mpp_mod,
only : mpp_chksum
39 use time_interp_external_mod,
only : init_external_field, time_interp_external
40 use time_interp_external_mod,
only : time_interp_external_init
42 implicit none ;
private
44 #include <MOM_memory.h>
59 integer :: wind_stagger
63 logical :: use_temperature
64 real :: wind_stress_multiplier
67 real :: area_surf = -1.0
68 real :: latent_heat_fusion
69 real :: latent_heat_vapor
77 logical :: use_limited_p_ssh
82 logical :: read_gust_2d
84 real,
pointer,
dimension(:,:) :: &
85 tke_tidal => null(), &
94 logical :: read_tideamp
95 logical :: rigid_sea_ice
100 real :: density_sea_ice
103 real :: rigid_sea_ice_mass
106 logical :: allow_flux_adjustments
108 logical :: salt_restore_as_sflux
109 logical :: adjust_net_srestore_to_zero
110 logical :: adjust_net_srestore_by_scaling
111 logical :: adjust_net_fresh_water_to_zero
112 logical :: use_net_fw_adjustment_sign_bug
113 logical :: adjust_net_fresh_water_by_scaling
114 logical :: mask_srestore_under_ice
116 real :: ice_salt_concentration
117 logical :: mask_srestore_marginal_seas
118 real :: max_delta_srestore
119 real :: max_delta_trestore
120 real,
pointer,
dimension(:,:) :: basin_mask => null()
122 character(len=200) :: inputdir
123 character(len=200) :: salt_restore_file
124 character(len=30) :: salt_restore_var_name
125 logical :: mask_srestore
129 real,
pointer,
dimension(:,:) :: srestore_mask => null()
130 character(len=200) :: temp_restore_file
131 character(len=30) :: temp_restore_var_name
132 logical :: mask_trestore
136 real,
pointer,
dimension(:,:) :: trestore_mask => null()
137 integer :: id_srestore = -1
138 integer :: id_trestore = -1
148 real,
pointer,
dimension(:,:) :: rofl_flux =>null()
149 real,
pointer,
dimension(:,:) :: rofi_flux =>null()
150 real,
pointer,
dimension(:,:) :: u_flux =>null()
151 real,
pointer,
dimension(:,:) :: v_flux =>null()
152 real,
pointer,
dimension(:,:) :: t_flux =>null()
153 real,
pointer,
dimension(:,:) :: q_flux =>null()
154 real,
pointer,
dimension(:,:) :: salt_flux =>null()
155 real,
pointer,
dimension(:,:) :: seaice_melt_heat =>null()
156 real,
pointer,
dimension(:,:) :: seaice_melt =>null()
157 real,
pointer,
dimension(:,:) :: lw_flux =>null()
158 real,
pointer,
dimension(:,:) :: sw_flux_vis_dir =>null()
159 real,
pointer,
dimension(:,:) :: sw_flux_vis_dif =>null()
160 real,
pointer,
dimension(:,:) :: sw_flux_nir_dir =>null()
161 real,
pointer,
dimension(:,:) :: sw_flux_nir_dif =>null()
162 real,
pointer,
dimension(:,:) :: lprec =>null()
163 real,
pointer,
dimension(:,:) :: fprec =>null()
164 real,
pointer,
dimension(:,:) :: runoff =>null()
165 real,
pointer,
dimension(:,:) :: calving =>null()
166 real,
pointer,
dimension(:,:) :: ustar_berg =>null()
167 real,
pointer,
dimension(:,:) :: area_berg =>null()
168 real,
pointer,
dimension(:,:) :: mass_berg =>null()
169 real,
pointer,
dimension(:,:) :: runoff_hflx =>null()
170 real,
pointer,
dimension(:,:) :: calving_hflx =>null()
171 real,
pointer,
dimension(:,:) :: p =>null()
173 real,
pointer,
dimension(:,:) :: mi =>null()
174 real,
pointer,
dimension(:,:) :: ice_rigidity =>null()
179 type(coupler_2d_bc_type) :: fluxes
181 integer :: wind_stagger = -999
196 sfc_state, restore_salt, restore_temp)
199 target,
intent(in) :: iob
202 type(
forcing),
intent(inout) :: fluxes
205 integer,
dimension(4),
intent(in) :: index_bounds
206 type(time_type),
intent(in) :: time
208 real,
intent(in) :: valid_time
214 type(
surface),
intent(in) :: sfc_state
216 logical,
optional,
intent(in) :: restore_salt
217 logical,
optional,
intent(in) :: restore_temp
220 real,
dimension(SZI_(G),SZJ_(G)) :: &
221 data_restore, & !< The surface value toward which to restore [g/kg or degC]
222 sst_anom, & !< Instantaneous sea surface temperature anomalies from a target value [deg C]
223 sss_anom, & !< Instantaneous sea surface salinity anomalies from a target value [g/kg]
224 sss_mean, & !< A (mean?) salinity about which to normalize local salinity
234 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq, i0, j0
235 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, isr, ier, jsr, jer
236 integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd
238 logical :: restore_salinity
240 logical :: restore_sst
245 real :: kg_m2_s_conversion
248 real :: sign_for_net_fw_bug
252 isc_bnd = index_bounds(1) ; iec_bnd = index_bounds(2)
253 jsc_bnd = index_bounds(3) ; jec_bnd = index_bounds(4)
254 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
255 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
256 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
257 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
258 isr = is-isd+1 ; ier = ie-isd+1 ; jsr = js-jsd+1 ; jer = je-jsd+1
260 kg_m2_s_conversion = us%kg_m3_to_R*us%m_to_Z*us%T_to_s
262 open_ocn_mask(:,:) = 1.0
264 fluxes%vPrecGlobalAdj = 0.0
265 fluxes%vPrecGlobalScl = 0.0
266 fluxes%saltFluxGlobalAdj = 0.0
267 fluxes%saltFluxGlobalScl = 0.0
268 fluxes%netFWGlobalAdj = 0.0
269 fluxes%netFWGlobalScl = 0.0
271 restore_salinity = .false.
272 if (
present(restore_salt)) restore_salinity = restore_salt
273 restore_sst = .false.
274 if (
present(restore_temp)) restore_sst = restore_temp
278 if (fluxes%dt_buoy_accum < 0)
then
279 call allocate_forcing_type(g, fluxes, water=.true., heat=.true., &
280 ustar=.true., press=.true.)
282 call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed)
283 call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed)
284 call safe_alloc_ptr(fluxes%sw_nir_dir,isd,ied,jsd,jed)
285 call safe_alloc_ptr(fluxes%sw_nir_dif,isd,ied,jsd,jed)
287 call safe_alloc_ptr(fluxes%p_surf,isd,ied,jsd,jed)
288 call safe_alloc_ptr(fluxes%p_surf_full,isd,ied,jsd,jed)
289 if (cs%use_limited_P_SSH)
then
290 fluxes%p_surf_SSH => fluxes%p_surf
292 fluxes%p_surf_SSH => fluxes%p_surf_full
295 call safe_alloc_ptr(fluxes%salt_flux,isd,ied,jsd,jed)
296 call safe_alloc_ptr(fluxes%salt_flux_in,isd,ied,jsd,jed)
297 call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed)
299 call safe_alloc_ptr(fluxes%TKE_tidal,isd,ied,jsd,jed)
300 call safe_alloc_ptr(fluxes%ustar_tidal,isd,ied,jsd,jed)
302 if (cs%allow_flux_adjustments)
then
303 call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed)
304 call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed)
307 do j=js-2,je+2 ;
do i=is-2,ie+2
308 fluxes%TKE_tidal(i,j) = cs%TKE_tidal(i,j)
309 fluxes%ustar_tidal(i,j) = cs%ustar_tidal(i,j)
312 if (restore_temp)
call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed)
317 if (((
associated(iob%ustar_berg) .and. (.not.
associated(fluxes%ustar_berg))) &
318 .or. (
associated(iob%area_berg) .and. (.not.
associated(fluxes%area_berg)))) &
319 .or. (
associated(iob%mass_berg) .and. (.not.
associated(fluxes%mass_berg)))) &
320 call allocate_forcing_type(g, fluxes, iceberg=.true.)
322 if ((.not.coupler_type_initialized(fluxes%tr_fluxes)) .and. &
323 coupler_type_initialized(iob%fluxes)) &
324 call coupler_type_spawn(iob%fluxes, fluxes%tr_fluxes, &
325 (/is,is,ie,ie/), (/js,js,je,je/))
331 if (cs%area_surf < 0.0)
then
332 do j=js,je ;
do i=is,ie
333 work_sum(i,j) = us%L_to_m**2*g%areaT(i,j) * g%mask2dT(i,j)
340 fluxes%fluxes_used = .false.
341 fluxes%dt_buoy_accum = us%s_to_T*valid_time
343 if (cs%allow_flux_adjustments)
then
344 fluxes%heat_added(:,:)=0.0
345 fluxes%salt_flux_added(:,:)=0.0
348 do j=js,je ;
do i=is,ie
349 fluxes%salt_flux(i,j) = 0.0
350 fluxes%vprec(i,j) = 0.0
354 if (restore_salinity)
then
355 call time_interp_external(cs%id_srestore,time,data_restore)
357 open_ocn_mask(:,:) = 1.0
358 if (cs%mask_srestore_under_ice)
then
359 do j=js,je ;
do i=is,ie
360 if (sfc_state%SST(i,j) <= -0.0539*sfc_state%SSS(i,j)) open_ocn_mask(i,j)=0.0
363 if (cs%salt_restore_as_sflux)
then
364 do j=js,je ;
do i=is,ie
365 delta_sss = data_restore(i,j)- sfc_state%SSS(i,j)
366 delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),cs%max_delta_srestore)
367 fluxes%salt_flux(i,j) = 1.e-3*g%mask2dT(i,j) * (cs%Rho0*us%m_to_Z*us%T_to_s*cs%Flux_const)* &
368 (cs%basin_mask(i,j)*open_ocn_mask(i,j)*cs%srestore_mask(i,j)) *delta_sss
370 if (cs%adjust_net_srestore_to_zero)
then
371 if (cs%adjust_net_srestore_by_scaling)
then
373 unit_scale=us%R_to_kg_m3*us%Z_to_m*us%s_to_T)
374 fluxes%saltFluxGlobalAdj = 0.
376 work_sum(is:ie,js:je) = us%L_to_m**2*us%R_to_kg_m3*us%Z_to_m*us%s_to_T * &
377 g%areaT(is:ie,js:je)*fluxes%salt_flux(is:ie,js:je)
378 fluxes%saltFluxGlobalAdj =
reproducing_sum(work_sum(:,:), isr,ier, jsr,jer)/cs%area_surf
379 fluxes%salt_flux(is:ie,js:je) = fluxes%salt_flux(is:ie,js:je) - kg_m2_s_conversion * fluxes%saltFluxGlobalAdj
382 fluxes%salt_flux_added(is:ie,js:je) = fluxes%salt_flux(is:ie,js:je)
384 do j=js,je ;
do i=is,ie
385 if (g%mask2dT(i,j) > 0.5)
then
386 delta_sss = sfc_state%SSS(i,j) - data_restore(i,j)
387 delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),cs%max_delta_srestore)
388 fluxes%vprec(i,j) = (cs%basin_mask(i,j)*open_ocn_mask(i,j)*cs%srestore_mask(i,j))* &
389 (us%m_to_Z*us%T_to_s * cs%Rho0*cs%Flux_const) * &
390 delta_sss / (0.5*(sfc_state%SSS(i,j) + data_restore(i,j)))
393 if (cs%adjust_net_srestore_to_zero)
then
394 if (cs%adjust_net_srestore_by_scaling)
then
396 unit_scale=us%R_to_kg_m3*us%Z_to_m*us%s_to_T)
397 fluxes%vPrecGlobalAdj = 0.
399 work_sum(is:ie,js:je) = us%L_to_m**2*g%areaT(is:ie,js:je) * &
400 us%R_to_kg_m3*us%Z_to_m*us%s_to_T*fluxes%vprec(is:ie,js:je)
401 fluxes%vPrecGlobalAdj =
reproducing_sum(work_sum(:,:), isr, ier, jsr, jer) / cs%area_surf
402 do j=js,je ;
do i=is,ie
403 fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - kg_m2_s_conversion*fluxes%vPrecGlobalAdj ) * g%mask2dT(i,j)
411 if (restore_sst)
then
412 call time_interp_external(cs%id_trestore,time,data_restore)
413 do j=js,je ;
do i=is,ie
414 delta_sst = data_restore(i,j)- sfc_state%SST(i,j)
415 delta_sst = sign(1.0,delta_sst)*min(abs(delta_sst),cs%max_delta_trestore)
416 fluxes%heat_added(i,j) = g%mask2dT(i,j) * cs%trestore_mask(i,j) * &
417 (us%R_to_kg_m3*cs%Rho0*fluxes%C_p) * delta_sst * cs%Flux_const
423 do j=js,je ;
do i=is,ie
425 if (
associated(iob%lprec)) &
426 fluxes%lprec(i,j) = kg_m2_s_conversion * iob%lprec(i-i0,j-j0) * g%mask2dT(i,j)
429 if (
associated(iob%fprec)) &
430 fluxes%fprec(i,j) = kg_m2_s_conversion * iob%fprec(i-i0,j-j0) * g%mask2dT(i,j)
433 if (
associated(iob%q_flux)) &
434 fluxes%evap(i,j) = kg_m2_s_conversion * iob%q_flux(i-i0,j-j0) * g%mask2dT(i,j)
437 if (
associated(iob%rofl_flux))
then
438 fluxes%lrunoff(i,j) = kg_m2_s_conversion * iob%rofl_flux(i-i0,j-j0) * g%mask2dT(i,j)
439 else if (
associated(iob%runoff))
then
440 fluxes%lrunoff(i,j) = kg_m2_s_conversion * iob%runoff(i-i0,j-j0) * g%mask2dT(i,j)
444 if (
associated(iob%rofi_flux))
then
445 fluxes%frunoff(i,j) = kg_m2_s_conversion * iob%rofi_flux(i-i0,j-j0) * g%mask2dT(i,j)
446 else if (
associated(iob%calving))
then
447 fluxes%frunoff(i,j) = kg_m2_s_conversion * iob%calving(i-i0,j-j0) * g%mask2dT(i,j)
450 if (
associated(iob%ustar_berg)) &
451 fluxes%ustar_berg(i,j) = us%m_to_Z * iob%ustar_berg(i-i0,j-j0) * g%mask2dT(i,j)
453 if (
associated(iob%area_berg)) &
454 fluxes%area_berg(i,j) = iob%area_berg(i-i0,j-j0) * g%mask2dT(i,j)
456 if (
associated(iob%mass_berg)) &
457 fluxes%mass_berg(i,j) = iob%mass_berg(i-i0,j-j0) * g%mask2dT(i,j)
461 if (
associated(fluxes%heat_content_lrunoff)) &
462 fluxes%heat_content_lrunoff(i,j) = 0.0 * g%mask2dT(i,j)
463 if (
associated(fluxes%heat_content_frunoff)) &
464 fluxes%heat_content_frunoff(i,j) = 0.0 * g%mask2dT(i,j)
466 if (
associated(iob%calving_hflx)) &
467 fluxes%heat_content_frunoff(i,j) = kg_m2_s_conversion * iob%calving_hflx(i-i0,j-j0) * g%mask2dT(i,j)
470 if (
associated(iob%lw_flux)) &
471 fluxes%LW(i,j) = iob%lw_flux(i-i0,j-j0) * g%mask2dT(i,j)
474 if (
associated(iob%t_flux)) &
475 fluxes%sens(i,j) = iob%t_flux(i-i0,j-j0) * g%mask2dT(i,j)
478 if (
associated(iob%seaice_melt_heat)) &
479 fluxes%seaice_melt_heat(i,j) = g%mask2dT(i,j) * iob%seaice_melt_heat(i-i0,j-j0)
482 if (
associated(iob%seaice_melt)) &
483 fluxes%seaice_melt(i,j) = g%mask2dT(i,j) * kg_m2_s_conversion * iob%seaice_melt(i-i0,j-j0)
486 fluxes%latent(i,j) = 0.0
488 if (
associated(iob%fprec))
then
489 fluxes%latent(i,j) = fluxes%latent(i,j) + iob%fprec(i-i0,j-j0)*cs%latent_heat_fusion
490 fluxes%latent_fprec_diag(i,j) = g%mask2dT(i,j) * iob%fprec(i-i0,j-j0)*cs%latent_heat_fusion
493 if (
associated(fluxes%frunoff))
then
494 fluxes%latent(i,j) = fluxes%latent(i,j) + iob%rofi_flux(i-i0,j-j0)*cs%latent_heat_fusion
495 fluxes%latent_frunoff_diag(i,j) = g%mask2dT(i,j) * iob%rofi_flux(i-i0,j-j0)*cs%latent_heat_fusion
498 if (
associated(iob%q_flux))
then
499 fluxes%latent(i,j) = fluxes%latent(i,j) + iob%q_flux(i-i0,j-j0)*cs%latent_heat_vapor
500 fluxes%latent_evap_diag(i,j) = g%mask2dT(i,j) * iob%q_flux(i-i0,j-j0)*cs%latent_heat_vapor
502 fluxes%latent(i,j) = g%mask2dT(i,j) * fluxes%latent(i,j)
504 if (
associated(iob%sw_flux_vis_dir)) &
505 fluxes%sw_vis_dir(i,j) = g%mask2dT(i,j) * iob%sw_flux_vis_dir(i-i0,j-j0)
507 if (
associated(iob%sw_flux_vis_dif)) &
508 fluxes%sw_vis_dif(i,j) = g%mask2dT(i,j) * iob%sw_flux_vis_dif(i-i0,j-j0)
510 if (
associated(iob%sw_flux_nir_dir)) &
511 fluxes%sw_nir_dir(i,j) = g%mask2dT(i,j) * iob%sw_flux_nir_dir(i-i0,j-j0)
513 if (
associated(iob%sw_flux_nir_dif)) &
514 fluxes%sw_nir_dif(i,j) = g%mask2dT(i,j) * iob%sw_flux_nir_dif(i-i0,j-j0)
516 fluxes%sw(i,j) = fluxes%sw_vis_dir(i,j) + fluxes%sw_vis_dif(i,j) + &
517 fluxes%sw_nir_dir(i,j) + fluxes%sw_nir_dif(i,j)
522 if (
associated(iob%p))
then
523 if (cs%max_p_surf >= 0.0)
then
524 do j=js,je ;
do i=is,ie
525 fluxes%p_surf_full(i,j) = g%mask2dT(i,j) * iob%p(i-i0,j-j0)
526 fluxes%p_surf(i,j) = min(fluxes%p_surf_full(i,j),cs%max_p_surf)
529 do j=js,je ;
do i=is,ie
530 fluxes%p_surf_full(i,j) = g%mask2dT(i,j) * iob%p(i-i0,j-j0)
531 fluxes%p_surf(i,j) = fluxes%p_surf_full(i,j)
534 fluxes%accumulate_p_surf = .true.
537 if (
associated(iob%salt_flux))
then
538 do j=js,je ;
do i=is,ie
539 fluxes%salt_flux(i,j) = g%mask2dT(i,j)*(fluxes%salt_flux(i,j) + kg_m2_s_conversion*iob%salt_flux(i-i0,j-j0))
540 fluxes%salt_flux_in(i,j) = g%mask2dT(i,j)*( kg_m2_s_conversion*iob%salt_flux(i-i0,j-j0) )
545 if (cs%adjust_net_fresh_water_to_zero)
then
546 sign_for_net_fw_bug = 1.
547 if (cs%use_net_FW_adjustment_sign_bug) sign_for_net_fw_bug = -1.
548 do j=js,je ;
do i=is,ie
549 net_fw(i,j) = us%R_to_kg_m3*us%Z_to_m*us%s_to_T * &
550 (((fluxes%lprec(i,j) + fluxes%fprec(i,j) + fluxes%seaice_melt(i,j)) + &
551 (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j))) + &
552 (fluxes%evap(i,j) + fluxes%vprec(i,j)) ) * us%L_to_m**2*g%areaT(i,j)
554 net_fw2(i,j) = net_fw(i,j) / (us%L_to_m**2*g%areaT(i,j))
557 if (cs%adjust_net_fresh_water_by_scaling)
then
559 do j=js,je ;
do i=is,ie
560 fluxes%vprec(i,j) = fluxes%vprec(i,j) + kg_m2_s_conversion * &
561 (net_fw2(i,j) - net_fw(i,j)/(us%L_to_m**2*g%areaT(i,j))) * g%mask2dT(i,j)
564 fluxes%netFWGlobalAdj =
reproducing_sum(net_fw(:,:), isr, ier, jsr, jer) / cs%area_surf
565 do j=js,je ;
do i=is,ie
566 fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - kg_m2_s_conversion * fluxes%netFWGlobalAdj ) * g%mask2dT(i,j)
571 if (coupler_type_initialized(fluxes%tr_fluxes) .and. &
572 coupler_type_initialized(iob%fluxes)) &
573 call coupler_type_copy_data(iob%fluxes, fluxes%tr_fluxes)
575 if (cs%allow_flux_adjustments)
then
581 call user_alter_forcing(sfc_state, fluxes, time, g, cs%urf_CS)
592 target,
intent(in) :: iob
595 integer,
dimension(4),
intent(in) :: index_bounds
596 type(time_type),
intent(in) :: time
604 real,
dimension(SZIB_(G),SZJB_(G)) :: &
605 taux_at_q, & !< Zonal wind stresses at q points [R Z L T-2 ~> Pa]
608 real,
dimension(SZI_(G),SZJ_(G)) :: &
609 rigidity_at_h, & !< Ice rigidity at tracer points [m3 s-1]
610 taux_at_h, & !< Zonal wind stresses at h points [R Z L T-2 ~> Pa]
617 real :: pa_conversion
618 real :: stress_conversion
624 integer :: wind_stagger
625 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq, i0, j0
626 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, isr, ier, jsr, jer
627 integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd
631 isc_bnd = index_bounds(1) ; iec_bnd = index_bounds(2)
632 jsc_bnd = index_bounds(3) ; jec_bnd = index_bounds(4)
633 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
634 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
635 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
636 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
637 isr = is-isd+1 ; ier = ie-isd+1 ; jsr = js-jsd+1 ; jer = je-jsd+1
641 irho0 = us%L_to_Z / cs%Rho0
642 pa_conversion = us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z
643 stress_conversion = pa_conversion * cs%wind_stress_multiplier
647 if (.not.forces%initialized)
then
651 call safe_alloc_ptr(forces%p_surf,isd,ied,jsd,jed)
652 call safe_alloc_ptr(forces%p_surf_full,isd,ied,jsd,jed)
654 if (cs%rigid_sea_ice)
then
655 call safe_alloc_ptr(forces%rigidity_ice_u,isdb,iedb,jsd,jed)
656 call safe_alloc_ptr(forces%rigidity_ice_v,isd,ied,jsdb,jedb)
659 forces%initialized = .true.
662 if ( (
associated(iob%area_berg) .and. (.not.
associated(forces%area_berg))) .or. &
663 (
associated(iob%mass_berg) .and. (.not.
associated(forces%mass_berg))) ) &
665 if (
associated(iob%ice_rigidity))
then
666 rigidity_at_h(:,:) = 0.0
667 call safe_alloc_ptr(forces%rigidity_ice_u,isdb,iedb,jsd,jed)
668 call safe_alloc_ptr(forces%rigidity_ice_v,isd,ied,jsdb,jedb)
671 forces%accumulate_rigidity = .true.
672 if (
associated(forces%rigidity_ice_u)) forces%rigidity_ice_u(:,:) = 0.0
673 if (
associated(forces%rigidity_ice_v)) forces%rigidity_ice_v(:,:) = 0.0
676 if (cs%use_limited_P_SSH)
then
677 forces%p_surf_SSH => forces%p_surf
679 forces%p_surf_SSH => forces%p_surf_full
681 if (
associated(iob%p))
then
682 if (cs%max_p_surf >= 0.0)
then
683 do j=js,je ;
do i=is,ie
684 forces%p_surf_full(i,j) = g%mask2dT(i,j) * iob%p(i-i0,j-j0)
685 forces%p_surf(i,j) = min(forces%p_surf_full(i,j),cs%max_p_surf)
688 do j=js,je ;
do i=is,ie
689 forces%p_surf_full(i,j) = g%mask2dT(i,j) * iob%p(i-i0,j-j0)
690 forces%p_surf(i,j) = forces%p_surf_full(i,j)
694 do j=js,je ;
do i=is,ie
695 forces%p_surf_full(i,j) = 0.0
696 forces%p_surf(i,j) = 0.0
699 forces%accumulate_p_surf = .true.
704 if (wind_stagger == bgrid_ne)
then
706 taux_at_q(:,:) = 0.0 ; tauy_at_q(:,:) = 0.0
708 if (wind_stagger == agrid)
then
710 taux_at_h(:,:) = 0.0 ; tauy_at_h(:,:) = 0.0
714 do j=js,je ;
do i=is,ie
715 if (
associated(iob%area_berg)) &
716 forces%area_berg(i,j) = iob%area_berg(i-i0,j-j0) * g%mask2dT(i,j)
718 if (
associated(iob%mass_berg)) &
719 forces%mass_berg(i,j) = iob%mass_berg(i-i0,j-j0) * g%mask2dT(i,j)
721 if (
associated(iob%ice_rigidity)) &
722 rigidity_at_h(i,j) = iob%ice_rigidity(i-i0,j-j0) * g%mask2dT(i,j)
724 if (wind_stagger == bgrid_ne)
then
725 if (
associated(iob%u_flux)) taux_at_q(i,j) = iob%u_flux(i-i0,j-j0) * stress_conversion
726 if (
associated(iob%v_flux)) tauy_at_q(i,j) = iob%v_flux(i-i0,j-j0) * stress_conversion
727 elseif (wind_stagger == agrid)
then
728 if (
associated(iob%u_flux)) taux_at_h(i,j) = iob%u_flux(i-i0,j-j0) * stress_conversion
729 if (
associated(iob%v_flux)) tauy_at_h(i,j) = iob%v_flux(i-i0,j-j0) * stress_conversion
731 if (
associated(iob%u_flux)) forces%taux(i,j) = iob%u_flux(i-i0,j-j0) * stress_conversion
732 if (
associated(iob%v_flux)) forces%tauy(i,j) = iob%v_flux(i-i0,j-j0) * stress_conversion
738 if (wind_stagger == bgrid_ne)
then
741 call pass_vector(taux_at_q, tauy_at_q, g%Domain, stagger=bgrid_ne, halo=1)
743 do j=js,je ;
do i=isq,ieq
744 forces%taux(i,j) = 0.0
745 if ((g%mask2dBu(i,j) + g%mask2dBu(i,j-1)) > 0) &
746 forces%taux(i,j) = (g%mask2dBu(i,j)*taux_at_q(i,j) + &
747 g%mask2dBu(i,j-1)*taux_at_q(i,j-1)) / &
748 (g%mask2dBu(i,j) + g%mask2dBu(i,j-1))
751 do j=jsq,jeq ;
do i=is,ie
752 forces%tauy(i,j) = 0.0
753 if ((g%mask2dBu(i,j) + g%mask2dBu(i-1,j)) > 0) &
754 forces%tauy(i,j) = (g%mask2dBu(i,j)*tauy_at_q(i,j) + &
755 g%mask2dBu(i-1,j)*tauy_at_q(i-1,j)) / &
756 (g%mask2dBu(i,j) + g%mask2dBu(i-1,j))
763 do j=js,je ;
do i=is,ie
764 tau_mag = 0.0 ; gustiness = cs%gust_const
765 if (((g%mask2dBu(i,j) + g%mask2dBu(i-1,j-1)) + &
766 (g%mask2dBu(i,j-1) + g%mask2dBu(i-1,j))) > 0)
then
767 tau_mag = sqrt(((g%mask2dBu(i,j)*(taux_at_q(i,j)**2 + tauy_at_q(i,j)**2) + &
768 g%mask2dBu(i-1,j-1)*(taux_at_q(i-1,j-1)**2 + tauy_at_q(i-1,j-1)**2)) + &
769 (g%mask2dBu(i,j-1)*(taux_at_q(i,j-1)**2 + tauy_at_q(i,j-1)**2) + &
770 g%mask2dBu(i-1,j)*(taux_at_q(i-1,j)**2 + tauy_at_q(i-1,j)**2)) ) / &
771 ((g%mask2dBu(i,j) + g%mask2dBu(i-1,j-1)) + (g%mask2dBu(i,j-1) + g%mask2dBu(i-1,j))) )
772 if (cs%read_gust_2d) gustiness = cs%gust(i,j)
774 forces%ustar(i,j) = sqrt(gustiness*irho0 + irho0*tau_mag)
777 elseif (wind_stagger == agrid)
then
778 call pass_vector(taux_at_h, tauy_at_h, g%Domain, to_all+omit_corners, stagger=agrid, halo=1)
780 do j=js,je ;
do i=isq,ieq
781 forces%taux(i,j) = 0.0
782 if ((g%mask2dT(i,j) + g%mask2dT(i+1,j)) > 0) &
783 forces%taux(i,j) = (g%mask2dT(i,j)*taux_at_h(i,j) + &
784 g%mask2dT(i+1,j)*taux_at_h(i+1,j)) / &
785 (g%mask2dT(i,j) + g%mask2dT(i+1,j))
788 do j=jsq,jeq ;
do i=is,ie
789 forces%tauy(i,j) = 0.0
790 if ((g%mask2dT(i,j) + g%mask2dT(i,j+1)) > 0) &
791 forces%tauy(i,j) = (g%mask2dT(i,j)*tauy_at_h(i,j) + &
792 g%mask2dT(i,j+1)*tauy_at_h(i,j+1)) / &
793 (g%mask2dT(i,j) + g%mask2dT(i,j+1))
796 do j=js,je ;
do i=is,ie
797 gustiness = cs%gust_const
798 if (cs%read_gust_2d .and. (g%mask2dT(i,j) > 0)) gustiness = cs%gust(i,j)
799 forces%ustar(i,j) = sqrt(gustiness*irho0 + irho0 * g%mask2dT(i,j) * &
800 sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2))
806 call pass_vector(forces%taux, forces%tauy, g%Domain, halo=1)
808 do j=js,je ;
do i=is,ie
810 if ((g%mask2dCu(i-1,j) + g%mask2dCu(i,j)) > 0) &
811 taux2 = (g%mask2dCu(i-1,j)*forces%taux(i-1,j)**2 + &
812 g%mask2dCu(i,j)*forces%taux(i,j)**2) / (g%mask2dCu(i-1,j) + g%mask2dCu(i,j))
815 if ((g%mask2dCv(i,j-1) + g%mask2dCv(i,j)) > 0) &
816 tauy2 = (g%mask2dCv(i,j-1)*forces%tauy(i,j-1)**2 + &
817 g%mask2dCv(i,j)*forces%tauy(i,j)**2) / (g%mask2dCv(i,j-1) + g%mask2dCv(i,j))
819 if (cs%read_gust_2d)
then
820 forces%ustar(i,j) = sqrt(cs%gust(i,j)*irho0 + irho0*sqrt(taux2 + tauy2))
822 forces%ustar(i,j) = sqrt(cs%gust_const*irho0 + irho0*sqrt(taux2 + tauy2))
829 if (
associated(iob%ice_rigidity))
then
830 call pass_var(rigidity_at_h, g%Domain, halo=1)
831 do i=is-1,ie ;
do j=js,je
832 forces%rigidity_ice_u(i,j) = forces%rigidity_ice_u(i,j) + &
833 min(rigidity_at_h(i,j), rigidity_at_h(i+1,j))
835 do i=is,ie ;
do j=js-1,je
836 forces%rigidity_ice_v(i,j) = forces%rigidity_ice_v(i,j) + &
837 min(rigidity_at_h(i,j), rigidity_at_h(i,j+1))
841 if (cs%rigid_sea_ice)
then
842 call pass_var(forces%p_surf_full, g%Domain, halo=1)
843 i_gearth = 1.0 / cs%G_Earth
844 kv_rho_ice = (cs%kv_sea_ice / cs%density_sea_ice)
845 do i=is-1,ie ;
do j=js,je
846 mass_ice = min(forces%p_surf_full(i,j), forces%p_surf_full(i+1,j)) * i_gearth
848 if (mass_ice > cs%rigid_sea_ice_mass)
then
849 mass_eff = (mass_ice - cs%rigid_sea_ice_mass) **2 / &
850 (mass_ice + cs%rigid_sea_ice_mass)
852 forces%rigidity_ice_u(i,j) = forces%rigidity_ice_u(i,j) + kv_rho_ice * mass_eff
854 do i=is,ie ;
do j=js-1,je
855 mass_ice = min(forces%p_surf_full(i,j), forces%p_surf_full(i,j+1)) * i_gearth
857 if (mass_ice > cs%rigid_sea_ice_mass)
then
858 mass_eff = (mass_ice - cs%rigid_sea_ice_mass) **2 / &
859 (mass_ice + cs%rigid_sea_ice_mass)
861 forces%rigidity_ice_v(i,j) = forces%rigidity_ice_v(i,j) + kv_rho_ice * mass_eff
865 if (cs%allow_flux_adjustments)
then
886 type(time_type),
intent(in) :: time
887 type(
forcing),
intent(inout) :: fluxes
890 real,
dimension(SZI_(G),SZJ_(G)) :: temp_at_h
892 integer :: isc, iec, jsc, jec, i, j
893 logical :: overrode_h
895 isc = g%isc; iec = g%iec ; jsc = g%jsc; jec = g%jec
898 call data_override(
'OCN',
'hflx_adj', temp_at_h(isc:iec,jsc:jec), time, override=overrode_h)
900 if (overrode_h)
then ;
do j=jsc,jec ;
do i=isc,iec
901 fluxes%heat_added(i,j) = fluxes%heat_added(i,j) + temp_at_h(i,j)* g%mask2dT(i,j)
902 enddo ;
enddo ;
endif
906 call data_override(
'OCN',
'sflx_adj', temp_at_h(isc:iec,jsc:jec), time, override=overrode_h)
908 if (overrode_h)
then ;
do j=jsc,jec ;
do i=isc,iec
909 fluxes%salt_flux_added(i,j) = fluxes%salt_flux_added(i,j) + &
910 us%kg_m3_to_R*us%m_to_Z*us%T_to_s * temp_at_h(i,j)* g%mask2dT(i,j)
911 enddo ;
enddo ;
endif
915 call data_override(
'OCN',
'prcme_adj', temp_at_h(isc:iec,jsc:jec), time, override=overrode_h)
917 if (overrode_h)
then ;
do j=jsc,jec ;
do i=isc,iec
918 fluxes%vprec(i,j) = fluxes%vprec(i,j) + us%kg_m3_to_R*us%m_to_Z*us%T_to_s * temp_at_h(i,j)* g%mask2dT(i,j)
919 enddo ;
enddo ;
endif
933 type(time_type),
intent(in) :: time
937 real,
dimension(SZI_(G),SZJ_(G)) :: tempx_at_h
938 real,
dimension(SZI_(G),SZJ_(G)) :: tempy_at_h
940 integer :: isc, iec, jsc, jec, i, j
941 real :: dlondx, dlondy, rdlon, cosa, sina, zonal_tau, merid_tau
942 real :: pa_conversion
943 logical :: overrode_x, overrode_y
945 isc = g%isc; iec = g%iec ; jsc = g%jsc; jec = g%jec
946 pa_conversion = us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z
948 tempx_at_h(:,:) = 0.0 ; tempy_at_h(:,:) = 0.0
950 overrode_x = .false. ; overrode_y = .false.
951 call data_override(
'OCN',
'taux_adj', tempx_at_h(isc:iec,jsc:jec), time, override=overrode_x)
952 call data_override(
'OCN',
'tauy_adj', tempy_at_h(isc:iec,jsc:jec), time, override=overrode_y)
954 if (overrode_x .or. overrode_y)
then
955 if (.not. (overrode_x .and. overrode_y))
call mom_error(fatal,
"apply_flux_adjustments: "//&
956 "Both taux_adj and tauy_adj must be specified, or neither, in data_table")
959 call pass_vector(tempx_at_h, tempy_at_h, g%Domain, to_all, agrid, halo=1)
960 do j=jsc-1,jec+1 ;
do i=isc-1,iec+1
961 dlondx = g%geoLonCu(i,j) - g%geoLonCu(i-1,j)
962 dlondy = g%geoLonCv(i,j) - g%geoLonCv(i,j-1)
963 rdlon = sqrt( dlondx * dlondx + dlondy * dlondy )
964 if (rdlon > 0.) rdlon = 1. / rdlon
965 cosa = dlondx * rdlon
966 sina = dlondy * rdlon
967 zonal_tau = pa_conversion * tempx_at_h(i,j)
968 merid_tau = pa_conversion * tempy_at_h(i,j)
969 tempx_at_h(i,j) = cosa * zonal_tau - sina * merid_tau
970 tempy_at_h(i,j) = sina * zonal_tau + cosa * merid_tau
974 do j=jsc,jec ;
do i=isc-1,iec
975 forces%taux(i,j) = forces%taux(i,j) + 0.5 * ( tempx_at_h(i,j) + tempx_at_h(i+1,j) )
978 do j=jsc-1,jec ;
do i=isc,iec
979 forces%tauy(i,j) = forces%tauy(i,j) + 0.5 * ( tempy_at_h(i,j) + tempy_at_h(i,j+1) )
991 type(time_type),
intent(in) :: time
992 character(len=*),
intent(in) :: directory
994 logical,
optional,
intent(in) :: time_stamped
996 character(len=*),
optional,
intent(in) :: filename_suffix
999 if (.not.
associated(cs))
return
1000 if (.not.
associated(cs%restart_CSp))
return
1001 call save_restart(directory, time, g, cs%restart_CSp, time_stamped)
1007 type(time_type),
intent(in) :: time
1011 type(
diag_ctrl),
target,
intent(inout) :: diag
1015 logical,
optional,
intent(in) :: restore_salt
1017 logical,
optional,
intent(in) :: restore_temp
1023 logical :: new_sim, iceberg_flux_diags
1024 type(time_type) :: time_frc
1025 character(len=200) :: tideamp_file, gust_file, salt_file, temp_file
1027 #include "version_variable.h"
1028 character(len=40) :: mdl =
"MOM_surface_forcing_mct"
1029 character(len=48) :: stagger
1030 character(len=48) :: flnam
1031 character(len=240) :: basin_file
1032 integer :: i, j, isd, ied, jsd, jed
1034 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1036 if (
associated(cs))
then
1037 call mom_error(warning,
"surface_forcing_init called with an associated "// &
1038 "control structure.")
1043 id_clock_forcing=cpu_clock_id(
'Ocean surface forcing', grain=clock_subcomponent)
1048 call write_version_number(version)
1052 call get_param(param_file, mdl,
"INPUTDIR", cs%inputdir, &
1053 "The directory in which all input files are found.", &
1055 cs%inputdir = slasher(cs%inputdir)
1056 call get_param(param_file, mdl,
"ENABLE_THERMODYNAMICS", cs%use_temperature, &
1057 "If true, Temperature and salinity are used as state "//&
1058 "variables.", default=.true.)
1059 call get_param(param_file, mdl,
"RHO_0", cs%Rho0, &
1060 "The mean ocean density used with BOUSSINESQ true to "//&
1061 "calculate accelerations and the mass for conservation "//&
1062 "properties, or with BOUSSINSEQ false to convert some "//&
1063 "parameters from vertical units of m to kg m-2.", &
1064 units=
"kg m-3", default=1035.0, scale=us%kg_m3_to_R)
1065 call get_param(param_file, mdl,
"LATENT_HEAT_FUSION", cs%latent_heat_fusion, &
1066 "The latent heat of fusion.", units=
"J/kg", default=hlf)
1067 call get_param(param_file, mdl,
"LATENT_HEAT_VAPORIZATION", cs%latent_heat_vapor, &
1068 "The latent heat of fusion.", units=
"J/kg", default=hlv)
1069 call get_param(param_file, mdl,
"MAX_P_SURF", cs%max_p_surf, &
1070 "The maximum surface pressure that can be exerted by the "//&
1071 "atmosphere and floating sea-ice or ice shelves. This is "//&
1072 "needed because the FMS coupling structure does not "//&
1073 "limit the water that can be frozen out of the ocean and "//&
1074 "the ice-ocean heat fluxes are treated explicitly. No "//&
1075 "limit is applied if a negative value is used.", units=
"Pa", &
1077 call get_param(param_file, mdl,
"ADJUST_NET_SRESTORE_TO_ZERO", &
1078 cs%adjust_net_srestore_to_zero, &
1079 "If true, adjusts the salinity restoring seen to zero "//&
1080 "whether restoring is via a salt flux or virtual precip.",&
1081 default=restore_salt)
1082 call get_param(param_file, mdl,
"ADJUST_NET_SRESTORE_BY_SCALING", &
1083 cs%adjust_net_srestore_by_scaling, &
1084 "If true, adjustments to salt restoring to achieve zero net are "//&
1085 "made by scaling values without moving the zero contour.",&
1087 call get_param(param_file, mdl,
"ADJUST_NET_FRESH_WATER_TO_ZERO", &
1088 cs%adjust_net_fresh_water_to_zero, &
1089 "If true, adjusts the net fresh-water forcing seen "//&
1090 "by the ocean (including restoring) to zero.", default=.false.)
1091 if (cs%adjust_net_fresh_water_to_zero) &
1092 call get_param(param_file, mdl,
"USE_NET_FW_ADJUSTMENT_SIGN_BUG", &
1093 cs%use_net_FW_adjustment_sign_bug, &
1094 "If true, use the wrong sign for the adjustment to "//&
1095 "the net fresh-water.", default=.false.)
1096 call get_param(param_file, mdl,
"ADJUST_NET_FRESH_WATER_BY_SCALING", &
1097 cs%adjust_net_fresh_water_by_scaling, &
1098 "If true, adjustments to net fresh water to achieve zero net are "//&
1099 "made by scaling values without moving the zero contour.",&
1101 call get_param(param_file, mdl,
"ICE_SALT_CONCENTRATION", &
1102 cs%ice_salt_concentration, &
1103 "The assumed sea-ice salinity needed to reverse engineer the "//&
1104 "melt flux (or ice-ocean fresh-water flux).", &
1105 units=
"kg/kg", default=0.005)
1106 call get_param(param_file, mdl,
"USE_LIMITED_PATM_SSH", cs%use_limited_P_SSH, &
1107 "If true, return the sea surface height with the "//&
1108 "correction for the atmospheric (and sea-ice) pressure "//&
1109 "limited by max_p_surf instead of the full atmospheric "//&
1110 "pressure.", default=.true.)
1112 call get_param(param_file, mdl,
"WIND_STAGGER", stagger, &
1113 "A case-insensitive character string to indicate the "//&
1114 "staggering of the input wind stress field. Valid "//&
1115 "values are 'A', 'B', or 'C'.", default=
"C")
1116 if (
uppercase(stagger(1:1)) ==
'A')
then ; cs%wind_stagger = agrid
1117 elseif (
uppercase(stagger(1:1)) ==
'B')
then ; cs%wind_stagger = bgrid_ne
1118 elseif (
uppercase(stagger(1:1)) ==
'C')
then ; cs%wind_stagger = cgrid_ne
1119 else ;
call mom_error(fatal,
"surface_forcing_init: WIND_STAGGER = "// &
1120 trim(stagger)//
" is invalid.") ;
endif
1121 call get_param(param_file, mdl,
"WIND_STRESS_MULTIPLIER", cs%wind_stress_multiplier, &
1122 "A factor multiplying the wind-stress given to the ocean by the "//&
1123 "coupler. This is used for testing and should be =1.0 for any "//&
1124 "production runs.", default=1.0)
1126 if (restore_salt)
then
1127 call get_param(param_file, mdl,
"FLUXCONST", cs%Flux_const, &
1128 "The constant that relates the restoring surface fluxes "//&
1129 "to the relative surface anomalies (akin to a piston "//&
1130 "velocity). Note the non-MKS units.", units=
"m day-1", &
1131 fail_if_missing=.true.)
1132 call get_param(param_file, mdl,
"SALT_RESTORE_FILE", cs%salt_restore_file, &
1133 "A file in which to find the surface salinity to use for restoring.", &
1134 default=
"salt_restore.nc")
1135 call get_param(param_file, mdl,
"SALT_RESTORE_VARIABLE", cs%salt_restore_var_name, &
1136 "The name of the surface salinity variable to read from "//&
1137 "SALT_RESTORE_FILE for restoring salinity.", &
1140 cs%Flux_const = cs%Flux_const / 86400.0
1142 call get_param(param_file, mdl,
"SRESTORE_AS_SFLUX", cs%salt_restore_as_sflux, &
1143 "If true, the restoring of salinity is applied as a salt "//&
1144 "flux instead of as a freshwater flux.", default=.false.)
1145 call get_param(param_file, mdl,
"MAX_DELTA_SRESTORE", cs%max_delta_srestore, &
1146 "The maximum salinity difference used in restoring terms.", &
1147 units=
"PSU or g kg-1", default=999.0)
1148 call get_param(param_file, mdl,
"MASK_SRESTORE_UNDER_ICE", &
1149 cs%mask_srestore_under_ice, &
1150 "If true, disables SSS restoring under sea-ice based on a frazil "//&
1151 "criteria (SST<=Tf). Only used when RESTORE_SALINITY is True.", &
1153 call get_param(param_file, mdl,
"MASK_SRESTORE_MARGINAL_SEAS", &
1154 cs%mask_srestore_marginal_seas, &
1155 "If true, disable SSS restoring in marginal seas. Only used when "//&
1156 "RESTORE_SALINITY is True.", default=.false.)
1157 call get_param(param_file, mdl,
"BASIN_FILE", basin_file, &
1158 "A file in which to find the basin masks, in variable 'basin'.", &
1160 basin_file = trim(cs%inputdir) // trim(basin_file)
1161 call safe_alloc_ptr(cs%basin_mask,isd,ied,jsd,jed) ; cs%basin_mask(:,:) = 1.0
1162 if (cs%mask_srestore_marginal_seas)
then
1163 call mom_read_data(basin_file,
'basin',cs%basin_mask,g%domain, timelevel=1)
1164 do j=jsd,jed ;
do i=isd,ied
1165 if (cs%basin_mask(i,j) >= 6.0)
then ; cs%basin_mask(i,j) = 0.0
1166 else ; cs%basin_mask(i,j) = 1.0 ;
endif
1169 call get_param(param_file, mdl,
"MASK_SRESTORE", cs%mask_srestore, &
1170 "If true, read a file (salt_restore_mask) containing "//&
1171 "a mask for SSS restoring.", default=.false.)
1174 if (restore_temp)
then
1175 call get_param(param_file, mdl,
"FLUXCONST", cs%Flux_const, &
1176 "The constant that relates the restoring surface fluxes "//&
1177 "to the relative surface anomalies (akin to a piston "//&
1178 "velocity). Note the non-MKS units.", units=
"m day-1", &
1179 fail_if_missing=.true.)
1180 call get_param(param_file, mdl,
"SST_RESTORE_FILE", cs%temp_restore_file, &
1181 "A file in which to find the surface temperature to use for restoring.", &
1182 default=
"temp_restore.nc")
1183 call get_param(param_file, mdl,
"SST_RESTORE_VARIABLE", cs%temp_restore_var_name, &
1184 "The name of the surface temperature variable to read from "//&
1185 "SST_RESTORE_FILE for restoring sst.", &
1188 cs%Flux_const = cs%Flux_const / 86400.0
1190 call get_param(param_file, mdl,
"MAX_DELTA_TRESTORE", cs%max_delta_trestore, &
1191 "The maximum sst difference used in restoring terms.", &
1192 units=
"degC ", default=999.0)
1194 call get_param(param_file, mdl,
"MASK_TRESTORE", cs%mask_trestore, &
1195 "If true, read a file (temp_restore_mask) containing "//&
1196 "a mask for SST restoring.", default=.false.)
1203 call get_param(param_file, mdl,
"CD_TIDES", cs%cd_tides, &
1204 "The drag coefficient that applies to the tides.", &
1205 units=
"nondim", default=1.0e-4)
1206 call get_param(param_file, mdl,
"READ_TIDEAMP", cs%read_TIDEAMP, &
1207 "If true, read a file (given by TIDEAMP_FILE) containing "//&
1208 "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
1209 if (cs%read_TIDEAMP)
then
1210 call get_param(param_file, mdl,
"TIDEAMP_FILE", tideamp_file, &
1211 "The path to the file containing the spatially varying "//&
1212 "tidal amplitudes with INT_TIDE_DISSIPATION.", &
1213 default=
"tideamp.nc")
1216 call get_param(param_file, mdl,
"UTIDE", cs%utide, &
1217 "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
1218 units=
"m s-1", default=0.0, scale=us%m_to_Z*us%T_to_s)
1221 call safe_alloc_ptr(cs%TKE_tidal,isd,ied,jsd,jed)
1222 call safe_alloc_ptr(cs%ustar_tidal,isd,ied,jsd,jed)
1224 if (cs%read_TIDEAMP)
then
1225 tideamp_file = trim(cs%inputdir) // trim(tideamp_file)
1226 call mom_read_data(tideamp_file,
'tideamp',cs%TKE_tidal,g%domain,timelevel=1, scale=us%m_to_Z*us%T_to_s)
1227 do j=jsd, jed;
do i=isd, ied
1228 utide = cs%TKE_tidal(i,j)
1229 cs%TKE_tidal(i,j) = g%mask2dT(i,j)*cs%Rho0*cs%cd_tides*(utide*utide*utide)
1230 cs%ustar_tidal(i,j) = sqrt(cs%cd_tides)*utide
1233 do j=jsd,jed;
do i=isd,ied
1235 cs%TKE_tidal(i,j) = cs%Rho0*cs%cd_tides*(utide*utide*utide)
1236 cs%ustar_tidal(i,j) = sqrt(cs%cd_tides)*utide
1240 call time_interp_external_init
1245 call get_param(param_file, mdl,
"READ_GUST_2D", cs%read_gust_2d, &
1246 "If true, use a 2-dimensional gustiness supplied from "//&
1247 "an input file", default=.false.)
1248 call get_param(param_file, mdl,
"GUST_CONST", cs%gust_const, &
1249 "The background gustiness in the winds.", units=
"Pa", &
1251 if (cs%read_gust_2d)
then
1252 call get_param(param_file, mdl,
"GUST_2D_FILE", gust_file, &
1253 "The file in which the wind gustiness is found in "//&
1254 "variable gustiness.")
1256 call safe_alloc_ptr(cs%gust,isd,ied,jsd,jed)
1257 gust_file = trim(cs%inputdir) // trim(gust_file)
1258 call mom_read_data(gust_file,
'gustiness',cs%gust,g%domain, timelevel=1, &
1259 scale=us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z)
1263 call get_param(param_file, mdl,
"USE_RIGID_SEA_ICE", cs%rigid_sea_ice, &
1264 "If true, sea-ice is rigid enough to exert a "//&
1265 "nonhydrostatic pressure that resist vertical motion.", &
1267 if (cs%rigid_sea_ice)
then
1268 call get_param(param_file, mdl,
"G_EARTH", cs%g_Earth, &
1269 "The gravitational acceleration of the Earth.", &
1270 units=
"m s-2", default = 9.80)
1271 call get_param(param_file, mdl,
"SEA_ICE_MEAN_DENSITY", cs%density_sea_ice, &
1272 "A typical density of sea ice, used with the kinematic "//&
1273 "viscosity, when USE_RIGID_SEA_ICE is true.", units=
"kg m-3", &
1275 call get_param(param_file, mdl,
"SEA_ICE_VISCOSITY", cs%Kv_sea_ice, &
1276 "The kinematic viscosity of sufficiently thick sea ice "//&
1277 "for use in calculating the rigidity of sea ice.", &
1278 units=
"m2 s-1", default=1.0e9)
1279 call get_param(param_file, mdl,
"SEA_ICE_RIGID_MASS", cs%rigid_sea_ice_mass, &
1280 "The mass of sea-ice per unit area at which the sea-ice "//&
1281 "starts to exhibit rigidity", units=
"kg m-2", default=1000.0)
1284 call get_param(param_file, mdl,
"ALLOW_ICEBERG_FLUX_DIAGNOSTICS", iceberg_flux_diags, &
1285 "If true, makes available diagnostics of fluxes from icebergs "//&
1286 "as seen by MOM6.", default=.false.)
1287 call register_forcing_type_diags(time, diag, us, cs%use_temperature, cs%handles, &
1288 use_berg_fluxes=iceberg_flux_diags)
1290 call get_param(param_file, mdl,
"ALLOW_FLUX_ADJUSTMENTS", cs%allow_flux_adjustments, &
1291 "If true, allows flux adjustments to specified via the "//&
1292 "data_table using the component name 'OCN'.", default=.false.)
1293 if (cs%allow_flux_adjustments)
then
1294 call data_override_init(ocean_domain_in=g%Domain%mpp_domain)
1297 if (
present(restore_salt))
then ;
if (restore_salt)
then
1298 salt_file = trim(cs%inputdir) // trim(cs%salt_restore_file)
1299 cs%id_srestore = init_external_field(salt_file, cs%salt_restore_var_name, domain=g%Domain%mpp_domain)
1300 call safe_alloc_ptr(cs%srestore_mask,isd,ied,jsd,jed); cs%srestore_mask(:,:) = 1.0
1301 if (cs%mask_srestore)
then
1302 flnam = trim(cs%inputdir) //
'salt_restore_mask.nc'
1303 call mom_read_data(flnam,
'mask', cs%srestore_mask, g%domain, timelevel=1)
1307 if (
present(restore_temp))
then ;
if (restore_temp)
then
1308 temp_file = trim(cs%inputdir) // trim(cs%temp_restore_file)
1309 cs%id_trestore = init_external_field(temp_file, cs%temp_restore_var_name, domain=g%Domain%mpp_domain)
1310 call safe_alloc_ptr(cs%trestore_mask,isd,ied,jsd,jed); cs%trestore_mask(:,:) = 1.0
1311 if (cs%mask_trestore)
then
1312 flnam = trim(cs%inputdir) //
'temp_restore_mask.nc'
1313 call mom_read_data(flnam,
'mask', cs%trestore_mask, g%domain, timelevel=1)
1318 call restart_init(param_file, cs%restart_CSp,
"MOM_forcing.res")
1321 if (
associated(cs%restart_CSp))
then
1325 if ((dirs%input_filename(1:1) ==
'n') .and. &
1326 (len_trim(dirs%input_filename) == 1)) new_sim = .true.
1327 if (.not.new_sim)
then
1328 call restore_state(dirs%input_filename, dirs%restart_input_dir, time_frc, &
1333 call user_revise_forcing_init(param_file, cs%urf_CS)
1343 type(
forcing),
optional,
intent(inout) :: fluxes
1347 if (
present(fluxes))
call deallocate_forcing_type(fluxes)
1349 if (
associated(cs))
deallocate(cs)
1357 character(len=*),
intent(in) :: id
1358 integer,
intent(in) :: timestep
1364 integer :: n,m, outunit
1368 write(outunit,*)
"BEGIN CHECKSUM(ice_ocean_boundary_type):: ", id, timestep
1369 write(outunit,100)
'iobt%u_flux ' , mpp_chksum( iobt%u_flux )
1370 write(outunit,100)
'iobt%v_flux ' , mpp_chksum( iobt%v_flux )
1371 write(outunit,100)
'iobt%t_flux ' , mpp_chksum( iobt%t_flux )
1372 write(outunit,100)
'iobt%q_flux ' , mpp_chksum( iobt%q_flux )
1373 write(outunit,100)
'iobt%salt_flux ' , mpp_chksum( iobt%salt_flux )
1374 write(outunit,100)
'iobt%seaice_melt_heat' , mpp_chksum( iobt%seaice_melt_heat)
1375 write(outunit,100)
'iobt%seaice_melt ' , mpp_chksum( iobt%seaice_melt )
1376 write(outunit,100)
'iobt%lw_flux ' , mpp_chksum( iobt%lw_flux )
1377 write(outunit,100)
'iobt%sw_flux_vis_dir' , mpp_chksum( iobt%sw_flux_vis_dir)
1378 write(outunit,100)
'iobt%sw_flux_vis_dif' , mpp_chksum( iobt%sw_flux_vis_dif)
1379 write(outunit,100)
'iobt%sw_flux_nir_dir' , mpp_chksum( iobt%sw_flux_nir_dir)
1380 write(outunit,100)
'iobt%sw_flux_nir_dif' , mpp_chksum( iobt%sw_flux_nir_dif)
1381 write(outunit,100)
'iobt%lprec ' , mpp_chksum( iobt%lprec )
1382 write(outunit,100)
'iobt%fprec ' , mpp_chksum( iobt%fprec )
1383 write(outunit,100)
'iobt%runoff ' , mpp_chksum( iobt%runoff )
1384 write(outunit,100)
'iobt%calving ' , mpp_chksum( iobt%calving )
1385 write(outunit,100)
'iobt%p ' , mpp_chksum( iobt%p )
1386 if (
associated(iobt%ustar_berg)) &
1387 write(outunit,100)
'iobt%ustar_berg ' , mpp_chksum( iobt%ustar_berg )
1388 if (
associated(iobt%area_berg)) &
1389 write(outunit,100)
'iobt%area_berg ' , mpp_chksum( iobt%area_berg )
1390 if (
associated(iobt%mass_berg)) &
1391 write(outunit,100)
'iobt%mass_berg ' , mpp_chksum( iobt%mass_berg )
1392 100
FORMAT(
" CHECKSUM::",a20,
" = ",z20)
1394 call coupler_type_write_chksums(iobt%fluxes, outunit,
'iobt%')