7 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
8 use mom_cpu_clock,
only : clock_module_driver, clock_module, clock_routine
25 use time_interp_external_mod,
only : init_external_field, time_interp_external
26 use time_interp_external_mod,
only : time_interp_external_init
28 implicit none ;
private
30 #include <MOM_memory.h>
43 logical :: do_rivermix = .false.
45 real :: rivermix_depth = 0.0
46 logical :: reclaim_frazil
49 logical :: pressure_dependent_frazil
53 logical :: ignore_fluxes_over_land
57 logical :: use_river_heat_content
60 logical :: use_calving_heat_content
67 logical :: chl_from_file
69 type(time_type),
pointer :: time => null()
73 integer :: id_createdh = -1
74 integer :: id_brine_lay = -1
75 integer :: id_pensw_diag = -1
76 integer :: id_penswflux_diag = -1
77 integer :: id_nonpensw_diag = -1
78 integer :: id_chl = -1
81 real,
allocatable,
dimension(:,:) :: createdh
83 real,
allocatable,
dimension(:,:,:) :: pensw_diag
85 real,
allocatable,
dimension(:,:,:) :: penswflux_diag
87 real,
allocatable,
dimension(:,:) :: nonpensw_diag
102 subroutine make_frazil(h, tv, G, GV, CS, p_surf, halo)
105 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
111 real,
dimension(SZI_(G),SZJ_(G)), &
112 optional,
intent(in) :: p_surf
113 integer,
optional,
intent(in) :: halo
116 real,
dimension(SZI_(G)) :: &
120 real,
dimension(SZI_(G),SZK_(G)) :: &
125 integer :: i, j, k, is, ie, js, je, nz
127 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
128 if (
present(halo))
then
129 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
134 if (.not.cs%pressure_dependent_frazil)
then
135 do k=1,nz ;
do i=is,ie ; pressure(i,k) = 0.0 ;
enddo ;
enddo
142 if (
PRESENT(p_surf))
then ;
do i=is,ie
146 do i=is,ie ; fraz_col(i) = 0.0 ;
enddo
148 if (cs%pressure_dependent_frazil)
then
150 pressure(i,1) = ps(i) + (0.5*gv%H_to_Pa)*h(i,j,1)
152 do k=2,nz ;
do i=is,ie
153 pressure(i,k) = pressure(i,k-1) + &
154 (0.5*gv%H_to_Pa) * (h(i,j,k) + h(i,j,k-1))
158 if (cs%reclaim_frazil)
then
160 do i=is,ie ;
if (tv%frazil(i,j) > 0.0)
then
161 if (.not.t_fr_set)
then
163 1, ie-i+1, tv%eqn_of_state)
167 if (tv%T(i,j,1) > t_freeze(i))
then
170 hc = (tv%C_p*gv%H_to_kg_m2) * h(i,j,1)
171 if (tv%frazil(i,j) - hc * (tv%T(i,j,1) - t_freeze(i)) <= 0.0)
then
172 tv%T(i,j,1) = tv%T(i,j,1) - tv%frazil(i,j)/hc
175 tv%frazil(i,j) = tv%frazil(i,j) - hc * (tv%T(i,j,1) - t_freeze(i))
176 tv%T(i,j,1) = t_freeze(i)
185 if ((g%mask2dT(i,j) > 0.0) .and. &
186 ((tv%T(i,j,k) < 0.0) .or. (fraz_col(i) > 0.0)))
then
187 if (.not.t_fr_set)
then
189 1, ie-i+1, tv%eqn_of_state)
193 hc = (tv%C_p*gv%H_to_kg_m2) * h(i,j,k)
194 if (h(i,j,k) <= 10.0*gv%Angstrom_H)
then
196 if (tv%T(i,j,k) < t_freeze(i))
then
197 fraz_col(i) = fraz_col(i) + hc * (t_freeze(i) - tv%T(i,j,k))
198 tv%T(i,j,k) = t_freeze(i)
201 if (fraz_col(i) + hc * (t_freeze(i) - tv%T(i,j,k)) <= 0.0)
then
202 tv%T(i,j,k) = tv%T(i,j,k) - fraz_col(i)/hc
205 fraz_col(i) = fraz_col(i) + hc * (t_freeze(i) - tv%T(i,j,k))
206 tv%T(i,j,k) = t_freeze(i)
213 tv%frazil(i,j) = tv%frazil(i,j) + fraz_col(i)
225 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
231 real,
intent(in) :: dt
234 real,
dimension(SZI_(G)) :: &
237 real,
dimension(SZI_(G),SZK_(G)) :: &
239 real,
dimension(SZI_(G),SZK_(G)+1) :: &
249 real,
dimension(:,:,:),
pointer :: t=>null(), s=>null()
250 real,
dimension(:,:,:),
pointer :: kd_t=>null(), kd_s=>null()
251 integer :: i, j, k, is, ie, js, je, nz
253 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
254 h_neglect = gv%H_subroundoff
256 if (.not.
associated(tv%T))
call mom_error(fatal, &
257 "differential_diffuse_T_S: Called with an unassociated tv%T")
258 if (.not.
associated(tv%S))
call mom_error(fatal, &
259 "differential_diffuse_T_S: Called with an unassociated tv%S")
260 if (.not.
associated(visc%Kd_extra_T))
call mom_error(fatal, &
261 "differential_diffuse_T_S: Called with an unassociated visc%Kd_extra_T")
262 if (.not.
associated(visc%Kd_extra_S))
call mom_error(fatal, &
263 "differential_diffuse_T_S: Called with an unassociated visc%Kd_extra_S")
265 t => tv%T ; s => tv%S
266 kd_t => visc%Kd_extra_T ; kd_s => visc%Kd_extra_S
272 i_h_int = 1.0 / (0.5 * (h(i,j,1) + h(i,j,2)) + h_neglect)
273 mix_t(i,2) = ((dt * kd_t(i,j,2)) * gv%Z_to_H**2) * i_h_int
274 mix_s(i,2) = ((dt * kd_s(i,j,2)) * gv%Z_to_H**2) * i_h_int
276 h_tr = h(i,j,1) + h_neglect
277 b1_t(i) = 1.0 / (h_tr + mix_t(i,2))
278 b1_s(i) = 1.0 / (h_tr + mix_s(i,2))
279 d1_t(i) = h_tr * b1_t(i)
280 d1_s(i) = h_tr * b1_s(i)
281 t(i,j,1) = (b1_t(i)*h_tr)*t(i,j,1)
282 s(i,j,1) = (b1_s(i)*h_tr)*s(i,j,1)
284 do k=2,nz-1 ;
do i=is,ie
286 i_h_int = 1.0 / (0.5 * (h(i,j,k) + h(i,j,k+1)) + h_neglect)
287 mix_t(i,k+1) = ((dt * kd_t(i,j,k+1)) * gv%Z_to_H**2) * i_h_int
288 mix_s(i,k+1) = ((dt * kd_s(i,j,k+1)) * gv%Z_to_H**2) * i_h_int
290 c1_t(i,k) = mix_t(i,k) * b1_t(i)
291 c1_s(i,k) = mix_s(i,k) * b1_s(i)
293 h_tr = h(i,j,k) + h_neglect
294 b_denom_t = h_tr + d1_t(i)*mix_t(i,k)
295 b_denom_s = h_tr + d1_s(i)*mix_s(i,k)
296 b1_t(i) = 1.0 / (b_denom_t + mix_t(i,k+1))
297 b1_s(i) = 1.0 / (b_denom_s + mix_s(i,k+1))
298 d1_t(i) = b_denom_t * b1_t(i)
299 d1_s(i) = b_denom_s * b1_s(i)
301 t(i,j,k) = b1_t(i) * (h_tr*t(i,j,k) + mix_t(i,k)*t(i,j,k-1))
302 s(i,j,k) = b1_s(i) * (h_tr*s(i,j,k) + mix_s(i,k)*s(i,j,k-1))
305 c1_t(i,nz) = mix_t(i,nz) * b1_t(i)
306 c1_s(i,nz) = mix_s(i,nz) * b1_s(i)
308 h_tr = h(i,j,nz) + h_neglect
309 b1_t(i) = 1.0 / (h_tr + d1_t(i)*mix_t(i,nz))
310 b1_s(i) = 1.0 / (h_tr + d1_s(i)*mix_s(i,nz))
312 t(i,j,nz) = b1_t(i) * (h_tr*t(i,j,nz) + mix_t(i,nz)*t(i,j,nz-1))
313 s(i,j,nz) = b1_s(i) * (h_tr*s(i,j,nz) + mix_s(i,nz)*s(i,j,nz-1))
315 do k=nz-1,1,-1 ;
do i=is,ie
316 t(i,j,k) = t(i,j,k) + c1_t(i,k+1)*t(i,j,k+1)
317 s(i,j,k) = s(i,j,k) + c1_s(i,k+1)*s(i,j,k+1)
328 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
334 integer,
optional,
intent(in) :: halo
337 real :: salt_add_col(szi_(g),szj_(g))
340 integer :: i, j, k, is, ie, js, je, nz
342 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
343 if (
present(halo))
then
344 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
349 s_min = tv%min_salinity
351 salt_add_col(:,:) = 0.0
355 do k=nz,1,-1 ;
do i=is,ie
356 if ( (g%mask2dT(i,j) > 0.0) .and. &
357 ((tv%S(i,j,k) < s_min) .or. (salt_add_col(i,j) > 0.0)) )
then
358 mc = gv%H_to_RZ * h(i,j,k)
359 if (h(i,j,k) <= 10.0*gv%Angstrom_H)
then
361 if (tv%S(i,j,k) < s_min)
then
362 salt_add_col(i,j) = salt_add_col(i,j) + mc * (s_min - tv%S(i,j,k))
365 elseif (salt_add_col(i,j) + mc * (s_min - tv%S(i,j,k)) <= 0.0)
then
366 tv%S(i,j,k) = tv%S(i,j,k) - salt_add_col(i,j) / mc
367 salt_add_col(i,j) = 0.0
369 salt_add_col(i,j) = salt_add_col(i,j) + mc * (s_min - tv%S(i,j,k))
375 tv%salt_deficit(i,j) = tv%salt_deficit(i,j) + salt_add_col(i,j)
385 subroutine insert_brine(h, tv, G, GV, US, fluxes, nkmb, CS, dt, id_brine_lay)
388 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
393 type(
forcing),
intent(in) :: fluxes
394 integer,
intent(in) :: nkmb
397 real,
intent(in) :: dt
398 integer,
intent(in) :: id_brine_lay
402 real :: salt(szi_(g))
403 real :: dzbr(szi_(g))
404 real :: inject_layer(szi_(g),szj_(g))
406 real :: p_ref_cv(szi_(g))
407 real :: t(szi_(g),szk_(g))
408 real :: s(szi_(g),szk_(g))
409 real :: h_2d(szi_(g),szk_(g))
410 real :: rcv(szi_(g),szk_(g))
411 real :: s_new,r_new,t0,scale, cdz
412 integer :: i, j, k, is, ie, js, je, nz, ks
415 real,
parameter :: s_max = 45.0
417 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
419 if (.not.
associated(fluxes%salt_flux))
return
425 p_ref_cv(:) = tv%P_ref
426 brine_dz = 1.0*gv%m_to_H
428 inject_layer(:,:) = nz
432 salt(:)=0.0 ; dzbr(:)=0.0
434 do i=is,ie ;
if (g%mask2dT(i,j) > 0.)
then
435 salt(i) = dt * (1000. * fluxes%salt_flux(i,j))
440 t(i,k) = tv%T(i,j,k) ; s(i,k) = tv%S(i,j,k)
442 h_2d(i,k) = max(h(i,j,k), gv%Angstrom_H)
446 ie-is+1, tv%eqn_of_state)
453 do k=nkmb+1,nz-1 ;
do i=is,ie
454 if ((g%mask2dT(i,j) > 0.0) .and. dzbr(i) < brine_dz .and. salt(i) > 0.)
then
455 s_new = s(i,k) + salt(i) / (gv%H_to_RZ * h_2d(i,k))
458 if (r_new < 0.5*(rcv(i,k)+rcv(i,k+1)) .and. s_new<s_max)
then
459 dzbr(i) = dzbr(i)+h_2d(i,k)
460 inject_layer(i,j) = min(inject_layer(i,j),real(k))
466 do k=nkmb,gv%nkml+1,-1 ;
do i=is,ie
467 if ((g%mask2dT(i,j) > 0.0) .and. dzbr(i) < brine_dz .and. salt(i) > 0.)
then
468 dzbr(i) = dzbr(i) + h_2d(i,k)
469 inject_layer(i,j) = min(inject_layer(i,j), real(k))
475 do k=1,gv%nkml ;
do i=is,ie
476 if ((g%mask2dT(i,j) > 0.0) .and. dzbr(i) < brine_dz .and. salt(i) > 0.)
then
477 dzbr(i) = dzbr(i) + h_2d(i,k)
478 inject_layer(i,j) = min(inject_layer(i,j), real(k))
484 if ((g%mask2dT(i,j) > 0.0) .and. salt(i) > 0.)
then
486 ks = inject_layer(i,j)
489 scale = h_2d(i,k) / dzbr(i)
490 cdz = cdz + h_2d(i,k)
493 if (cdz > brine_dz)
exit
494 tv%S(i,j,k) = tv%S(i,j,k) + scale*salt(i) / (gv%H_to_RZ * h_2d(i,k))
501 if (cs%id_brine_lay > 0)
call post_data(cs%id_brine_lay, inject_layer, cs%diag)
507 subroutine tridiagts(G, GV, is, ie, js, je, hold, ea, eb, T, S)
510 integer,
intent(in) :: is
511 integer,
intent(in) :: ie
512 integer,
intent(in) :: js
513 integer,
intent(in) :: je
514 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: hold
516 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: ea
518 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: eb
520 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: t
521 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: s
524 real :: b1(szib_(g)), d1(szib_(g))
525 real :: c1(szib_(g),szk_(g))
526 real :: h_tr, b_denom_1
532 h_tr = hold(i,j,1) + gv%H_subroundoff
533 b1(i) = 1.0 / (h_tr + eb(i,j,1))
535 t(i,j,1) = (b1(i)*h_tr)*t(i,j,1)
536 s(i,j,1) = (b1(i)*h_tr)*s(i,j,1)
538 do k=2,g%ke ;
do i=is,ie
539 c1(i,k) = eb(i,j,k-1) * b1(i)
540 h_tr = hold(i,j,k) + gv%H_subroundoff
541 b_denom_1 = h_tr + d1(i)*ea(i,j,k)
542 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
543 d1(i) = b_denom_1 * b1(i)
544 t(i,j,k) = b1(i) * (h_tr*t(i,j,k) + ea(i,j,k)*t(i,j,k-1))
545 s(i,j,k) = b1(i) * (h_tr*s(i,j,k) + ea(i,j,k)*s(i,j,k-1))
547 do k=g%ke-1,1,-1 ;
do i=is,ie
548 t(i,j,k) = t(i,j,k) + c1(i,k+1)*t(i,j,k+1)
549 s(i,j,k) = s(i,j,k) + c1(i,k+1)*s(i,j,k+1)
556 subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb)
560 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
562 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
564 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
566 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
568 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
570 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
571 optional,
intent(in) :: ea
574 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
575 optional,
intent(in) :: eb
583 real :: b1(szi_(g)), d1(szi_(g)), c1(szi_(g),szk_(g))
584 real :: a_n(szi_(g)), a_s(szi_(g))
585 real :: a_e(szi_(g)), a_w(szi_(g))
587 real :: sum_area, idenom
588 logical :: mix_vertically
589 integer :: i, j, k, is, ie, js, je, nz
590 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
592 h_neglect = gv%H_subroundoff
594 mix_vertically =
present(ea)
595 if (
present(ea) .neqv.
present(eb))
call mom_error(fatal, &
596 "find_uv_at_h: Either both ea and eb or neither one must be present "// &
597 "in call to find_uv_at_h.")
603 sum_area = g%areaCu(i-1,j) + g%areaCu(i,j)
604 if (sum_area>0.0)
then
605 idenom = sqrt(0.5*g%IareaT(i,j) / sum_area)
606 a_w(i) = g%areaCu(i-1,j) * idenom
607 a_e(i) = g%areaCu(i,j) * idenom
609 a_w(i) = 0.0 ; a_e(i) = 0.0
612 sum_area = g%areaCv(i,j-1) + g%areaCv(i,j)
613 if (sum_area>0.0)
then
614 idenom = sqrt(0.5*g%IareaT(i,j) / sum_area)
615 a_s(i) = g%areaCv(i,j-1) * idenom
616 a_n(i) = g%areaCv(i,j) * idenom
618 a_s(i) = 0.0 ; a_n(i) = 0.0
622 if (mix_vertically)
then
624 b_denom_1 = h(i,j,1) + h_neglect
625 b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
626 d1(i) = b_denom_1 * b1(i)
627 u_h(i,j,1) = (h(i,j,1)*b1(i)) * (a_e(i)*u(i,j,1) + a_w(i)*u(i-1,j,1))
628 v_h(i,j,1) = (h(i,j,1)*b1(i)) * (a_n(i)*v(i,j,1) + a_s(i)*v(i,j-1,1))
630 do k=2,nz ;
do i=is,ie
631 c1(i,k) = eb(i,j,k-1) * b1(i)
632 b_denom_1 = h(i,j,k) + d1(i)*ea(i,j,k) + h_neglect
633 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
634 d1(i) = b_denom_1 * b1(i)
635 u_h(i,j,k) = (h(i,j,k) * (a_e(i)*u(i,j,k) + a_w(i)*u(i-1,j,k)) + &
636 ea(i,j,k)*u_h(i,j,k-1))*b1(i)
637 v_h(i,j,k) = (h(i,j,k) * (a_n(i)*v(i,j,k) + a_s(i)*v(i,j-1,k)) + &
638 ea(i,j,k)*v_h(i,j,k-1))*b1(i)
640 do k=nz-1,1,-1 ;
do i=is,ie
641 u_h(i,j,k) = u_h(i,j,k) + c1(i,k+1)*u_h(i,j,k+1)
642 v_h(i,j,k) = v_h(i,j,k) + c1(i,k+1)*v_h(i,j,k+1)
645 do k=1,nz ;
do i=is,ie
646 u_h(i,j,k) = a_e(i)*u(i,j,k) + a_w(i)*u(i-1,j,k)
647 v_h(i,j,k) = a_n(i)*v(i,j,k) + a_s(i)*v(i,j-1,k)
656 subroutine set_pen_shortwave(optics, fluxes, G, GV, CS, opacity_CSp, tracer_flow_CSp)
659 type(
forcing),
intent(inout) :: fluxes
668 real,
dimension(SZI_(G),SZJ_(G)) :: chl_2d
669 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: chl_3d
670 character(len=128) :: mesg
671 integer :: i, j, k, is, ie, js, je
672 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
674 if (.not.
associated(optics))
return
676 if (cs%var_pen_sw)
then
677 if (cs%chl_from_file)
then
680 call time_interp_external(cs%sbc_chl, cs%Time, chl_2d)
681 do j=js,je ;
do i=is,ie
682 if ((g%mask2dT(i,j) > 0.5) .and. (chl_2d(i,j) < 0.0))
then
683 write(mesg,
'(" Time_interp negative chl of ",(1pe12.4)," at i,j = ",&
684 & 2(i3), "lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
685 chl_2d(i,j), i, j, g%geoLonT(i,j), g%geoLatT(i,j)
686 call mom_error(fatal,
"MOM_diabatic_aux set_pen_shortwave: "//trim(mesg))
690 if (cs%id_chl > 0)
call post_data(cs%id_chl, chl_2d, cs%diag)
692 call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
693 fluxes%sw_nir_dir, fluxes%sw_nir_dif, g, gv, opacity_csp, chl_2d=chl_2d)
695 if (.not.
associated(tracer_flow_csp))
call mom_error(fatal, &
696 "The tracer flow control structure must be associated when the model sets "//&
697 "the chlorophyll internally in set_pen_shortwave.")
700 if (cs%id_chl > 0)
call post_data(cs%id_chl, chl_3d(:,:,1), cs%diag)
702 call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
703 fluxes%sw_nir_dir, fluxes%sw_nir_dif, g, gv, opacity_csp, chl_3d=chl_3d)
706 call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
707 fluxes%sw_nir_dir, fluxes%sw_nir_dif, g, gv, opacity_csp)
716 id_N2subML, id_MLDsq, dz_subML)
720 integer,
intent(in) :: id_mld
721 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
725 real,
intent(in) :: densitydiff
727 integer,
optional,
intent(in) :: id_n2subml
728 integer,
optional,
intent(in) :: id_mldsq
729 real,
optional,
intent(in) :: dz_subml
733 real,
dimension(SZI_(G)) :: deltarhoatkm1, deltarhoatk
734 real,
dimension(SZI_(G)) :: pref_mld, pref_n2
735 real,
dimension(SZI_(G)) :: h_subml, dh_n2
736 real,
dimension(SZI_(G)) :: t_subml, t_deeper
737 real,
dimension(SZI_(G)) :: s_subml, s_deeper
738 real,
dimension(SZI_(G)) :: rho_subml, rho_deeper
739 real,
dimension(SZI_(G)) :: dk, dkm1
740 real,
dimension(SZI_(G)) :: rhosurf
741 real,
dimension(SZI_(G), SZJ_(G)) :: mld
742 real,
dimension(SZI_(G), SZJ_(G)) :: submln2
743 real,
dimension(SZI_(G), SZJ_(G)) :: mld2
744 logical,
dimension(SZI_(G)) :: n2_region_set
750 integer :: i, j, is, ie, js, je, k, nz, id_n2, id_sq
752 id_n2 = -1 ;
if (
PRESENT(id_n2subml)) id_n2 = id_n2subml
754 id_sq = -1 ;
if (
PRESENT(id_mldsq)) id_sq = id_mldsq
756 ge_rho0 = us%L_to_Z**2*gv%g_Earth / (gv%Rho0)
757 dh_subml = 50.*gv%m_to_H ;
if (
present(dz_subml)) dh_subml = gv%Z_to_H*dz_subml
759 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
763 do i=is,ie ; dk(i) = 0.5 * h(i,j,1) * gv%H_to_Z ;
enddo
764 call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pref_mld, rhosurf, is, ie-is+1, &
765 tv%eqn_of_state, scale=us%kg_m3_to_R)
771 h_subml(i) = h(i,j,1) ; dh_n2(i) = 0.0
772 t_subml(i) = 0.0 ; s_subml(i) = 0.0 ; t_deeper(i) = 0.0 ; s_deeper(i) = 0.0
773 n2_region_set(i) = (g%mask2dT(i,j)<0.5)
779 dk(i) = dk(i) + 0.5 * ( h(i,j,k) + h(i,j,k-1) ) * gv%H_to_Z
786 if (mld(i,j)==0.0)
then
787 h_subml(i) = h_subml(i) + h(i,j,k)
788 elseif (.not.n2_region_set(i))
then
789 if (dh_n2(i)==0.0)
then
790 t_subml(i) = tv%T(i,j,k) ; s_subml(i) = tv%S(i,j,k)
791 h_subml(i) = h_subml(i) + 0.5 * h(i,j,k)
792 dh_n2(i) = 0.5 * h(i,j,k)
793 elseif (dh_n2(i) + h(i,j,k) < dh_subml)
then
794 dh_n2(i) = dh_n2(i) + h(i,j,k)
796 t_deeper(i) = tv%T(i,j,k) ; s_deeper(i) = tv%S(i,j,k)
797 dh_n2(i) = dh_n2(i) + 0.5 * h(i,j,k)
798 n2_region_set(i) = .true.
805 do i=is,ie ; deltarhoatkm1(i) = deltarhoatk(i) ;
enddo
806 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pref_mld, deltarhoatk, is, ie-is+1, &
807 tv%eqn_of_state, scale=us%kg_m3_to_R)
809 deltarhoatk(i) = deltarhoatk(i) - rhosurf(i)
810 ddrho = deltarhoatk(i) - deltarhoatkm1(i)
811 if ((mld(i,j)==0.) .and. (ddrho>0.) .and. &
812 (deltarhoatkm1(i)<densitydiff) .and. (deltarhoatk(i)>=densitydiff))
then
813 afac = ( densitydiff - deltarhoatkm1(i) ) / ddrho
814 mld(i,j) = dk(i) * afac + dkm1(i) * (1. - afac)
816 if (id_sq > 0) mld2(i,j) = mld(i,j)**2
820 if ((mld(i,j)==0.) .and. (deltarhoatk(i)<densitydiff)) mld(i,j) = dk(i)
824 do i=is,ie ; pref_n2(i) = gv%H_to_Pa * (h_subml(i) + 0.5*dh_n2(i)) ;
enddo
831 tv%eqn_of_state, scale=us%kg_m3_to_R)
833 tv%eqn_of_state, scale=us%kg_m3_to_R)
834 do i=is,ie ;
if ((g%mask2dT(i,j)>0.5) .and. n2_region_set(i))
then
835 submln2(i,j) = ge_rho0 * (rho_deeper(i) - rho_subml(i)) / (gv%H_to_z * dh_n2(i))
840 if (id_mld > 0)
call post_data(id_mld, mld, diagptr)
841 if (id_n2 > 0)
call post_data(id_n2, submln2, diagptr)
842 if (id_sq > 0)
call post_data(id_sq, mld2, diagptr)
849 subroutine applyboundaryfluxesinout(CS, G, GV, US, dt, fluxes, optics, nsw, h, tv, &
850 aggregate_FW_forcing, evap_CFL_limit, &
851 minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, &
857 real,
intent(in) :: dt
858 type(
forcing),
intent(inout) :: fluxes
860 integer,
intent(in) :: nsw
862 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
866 logical,
intent(in) :: aggregate_fw_forcing
867 real,
intent(in) :: evap_cfl_limit
869 real,
intent(in) :: minimum_forcing_depth
871 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
872 optional,
intent(out) :: ctke
874 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
875 optional,
intent(out) :: dsv_dt
877 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
878 optional,
intent(out) :: dsv_ds
880 real,
dimension(SZI_(G),SZJ_(G)), &
881 optional,
intent(out) :: skinbuoyflux
884 integer,
parameter :: maxgroundings = 5
885 integer :: numberofgroundings, iground(maxgroundings), jground(maxgroundings)
886 real :: h_limit_fluxes
887 real :: iforcingdepthscale
889 real :: dthickness, dtemp, dsalt
890 real :: fractionofforcing, hold, ithickness
891 real :: rivermixconst
893 real,
dimension(SZI_(G)) :: &
913 real,
dimension(SZI_(G), SZK_(G)) :: &
919 real,
dimension(SZI_(G)) :: &
922 real,
dimension(max(nsw,1),SZI_(G)) :: &
927 real,
dimension(max(nsw,1),SZI_(G),SZK_(G)) :: &
930 real,
dimension(maxGroundings) :: hgrounding
931 real :: temp_in, salin_in
936 logical :: calculate_energetics
937 logical :: calculate_buoyancy
938 integer :: i, j, is, ie, js, je, k, nz, n, nb
939 integer :: start, npts
940 character(len=45) :: mesg
942 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
945 if (.not.
associated(fluxes%sw))
return
950 calculate_energetics = (
present(ctke) .and.
present(dsv_dt) .and.
present(dsv_ds))
951 calculate_buoyancy =
present(skinbuoyflux)
952 if (calculate_buoyancy) skinbuoyflux(:,:) = 0.0
953 g_hconv2 = (us%L_to_Z**2*gv%g_Earth * gv%H_to_RZ) * gv%H_to_RZ
955 if (
present(ctke)) ctke(:,:,:) = 0.0
956 if (calculate_buoyancy)
then
957 surfpressure(:) = 0.0
958 gorho = us%L_to_Z**2*gv%g_Earth / gv%Rho0
959 start = 1 + g%isc - g%isd
960 npts = 1 + g%iec - g%isc
968 h_limit_fluxes = max(gv%Angstrom_H, 1.e-30*gv%m_to_H)
971 if (cs%id_createdH>0) cs%createdH(:,:) = 0.
972 numberofgroundings = 0
993 do k=1,nz ;
do i=is,ie
995 t2d(i,k) = tv%T(i,j,k)
997 if (nsw>0)
call extract_optics_slice(optics, j, g, gv, opacity=opacityband, opacity_scale=(1.0/gv%m_to_H))
999 if (calculate_energetics)
then
1003 do i=is,ie ; pres(i) = 0.0 ;
enddo
1006 d_pres(i) = gv%H_to_Pa * h2d(i,k)
1007 p_lay(i) = pres(i) + 0.5*d_pres(i)
1008 pres(i) = pres(i) + d_pres(i)
1011 dsv_dt(:,j,k), dsv_ds(:,j,k), is, ie-is+1, tv%eqn_of_state, scale=us%R_to_kg_m3)
1012 do i=is,ie ; dsv_dt_2d(i,k) = dsv_dt(i,j,k) ;
enddo
1014 pen_tke_2d(:,:) = 0.0
1060 if (calculate_buoyancy)
then
1062 h_limit_fluxes, cs%use_river_heat_content, cs%use_calving_heat_content, &
1063 h2d, t2d, netmassinout, netmassout, netheat, netsalt, &
1064 pen_sw_bnd, tv, aggregate_fw_forcing, nonpensw=nonpensw, &
1065 net_heat_rate=netheat_rate, net_salt_rate=netsalt_rate, &
1066 netmassinout_rate=netmassinout_rate, pen_sw_bnd_rate=pen_sw_bnd_rate)
1069 h_limit_fluxes, cs%use_river_heat_content, cs%use_calving_heat_content, &
1070 h2d, t2d, netmassinout, netmassout, netheat, netsalt, &
1071 pen_sw_bnd, tv, aggregate_fw_forcing, nonpensw=nonpensw)
1076 if (aggregate_fw_forcing)
then
1077 netmassout(i) = netmassinout(i)
1080 netmassin(i) = netmassinout(i) - netmassout(i)
1082 if (g%mask2dT(i,j)>0.0)
then
1083 fluxes%netMassOut(i,j) = netmassout(i)
1084 fluxes%netMassIn(i,j) = netmassin(i)
1086 fluxes%netMassOut(i,j) = 0.0
1087 fluxes%netMassIn(i,j) = 0.0
1097 if (g%mask2dT(i,j)>0.)
then
1103 dthickness = netmassin(i)
1109 netmassin(i) = netmassin(i) - dthickness
1113 dtemp = dtemp + dthickness*temp_in
1116 if (
associated(fluxes%heat_content_massin)) &
1117 fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + &
1118 t2d(i,k) * max(0.,dthickness) * gv%H_to_RZ * fluxes%C_p * idt
1119 if (
associated(fluxes%heat_content_massout)) &
1120 fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
1121 t2d(i,k) * min(0.,dthickness) * gv%H_to_RZ * fluxes%C_p * idt
1122 if (
associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
1123 t2d(i,k) * dthickness * gv%H_to_RZ
1126 if (calculate_energetics .and.
associated(fluxes%lrunoff) .and. cs%do_rivermix)
then
1137 if (gv%Boussinesq)
then
1138 rivermixconst = -0.5*(cs%rivermix_depth*dt) * ( us%L_to_Z**2*gv%g_Earth ) * gv%Rho0
1140 rivermixconst = -0.5*(cs%rivermix_depth*dt) * gv%Rho0 * ( us%L_to_Z**2*gv%g_Earth )
1142 ctke(i,j,k) = ctke(i,j,k) + max(0.0, rivermixconst*dsv_ds(i,j,1) * &
1143 (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) * tv%S(i,j,1))
1148 h2d(i,k) = h2d(i,k) + dthickness
1149 if (h2d(i,k) > 0.0)
then
1150 if (calculate_energetics .and. (dthickness > 0.))
then
1153 ctke(i,j,k) = ctke(i,j,k) + 0.5*g_hconv2*(hold*dthickness) * &
1154 ((t2d(i,k) - temp_in) * dsv_dt(i,j,k) + (tv%S(i,j,k) - salin_in) * dsv_ds(i,j,k))
1156 ithickness = 1.0/h2d(i,k)
1158 if (dthickness /= 0. .or. dtemp /= 0.) t2d(i,k) = (hold*t2d(i,k) + dtemp)*ithickness
1159 if (dthickness /= 0. .or. dsalt /= 0.) tv%S(i,j,k) = (hold*tv%S(i,j,k) + dsalt)*ithickness
1171 iforcingdepthscale = 1. / max(gv%H_subroundoff, minimum_forcing_depth - netmassout(i) )
1173 fractionofforcing = min(1.0, h2d(i,k)*iforcingdepthscale)
1178 if (-fractionofforcing*netmassout(i) > evap_cfl_limit*h2d(i,k))
then
1179 fractionofforcing = -evap_cfl_limit*h2d(i,k)/netmassout(i)
1184 dthickness = max( fractionofforcing*netmassout(i), -h2d(i,k) )
1185 dtemp = fractionofforcing*netheat(i)
1187 dsalt = max( fractionofforcing*netsalt(i), -0.9999*h2d(i,k)*tv%S(i,j,k))
1191 netmassout(i) = netmassout(i) - dthickness
1192 netheat(i) = netheat(i) - dtemp
1193 netsalt(i) = netsalt(i) - dsalt
1196 dtemp = dtemp + dthickness*t2d(i,k)
1199 if (
associated(fluxes%heat_content_massin)) &
1200 fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + &
1201 t2d(i,k) * max(0.,dthickness) * gv%H_to_RZ * fluxes%C_p * idt
1202 if (
associated(fluxes%heat_content_massout)) &
1203 fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
1204 t2d(i,k) * min(0.,dthickness) * gv%H_to_RZ * fluxes%C_p * idt
1205 if (
associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
1206 t2d(i,k) * dthickness * gv%H_to_RZ
1210 h2d(i,k) = h2d(i,k) + dthickness
1212 if (h2d(i,k) > 0.)
then
1213 if (calculate_energetics)
then
1219 ctke(i,j,k) = ctke(i,j,k) - (0.5*h2d(i,k)*g_hconv2) * &
1220 ((dtemp - dthickness*t2d(i,k)) * dsv_dt(i,j,k) + &
1221 (dsalt - dthickness*tv%S(i,j,k)) * dsv_ds(i,j,k))
1223 ithickness = 1.0/h2d(i,k)
1224 t2d(i,k) = (hold*t2d(i,k) + dtemp)*ithickness
1225 tv%S(i,j,k) = (hold*tv%S(i,j,k) + dsalt)*ithickness
1226 elseif (h2d(i,k) < 0.0)
then
1228 write(0,*)
'applyBoundaryFluxesInOut(): lon,lat=',g%geoLonT(i,j),g%geoLatT(i,j)
1229 write(0,*)
'applyBoundaryFluxesInOut(): netT,netS,netH=',netheat(i),netsalt(i),netmassinout(i)
1230 write(0,*)
'applyBoundaryFluxesInOut(): dT,dS,dH=',dtemp,dsalt,dthickness
1231 write(0,*)
'applyBoundaryFluxesInOut(): h(n),h(n+1),k=',hold,h2d(i,k),k
1232 call mom_error(fatal,
"MOM_diabatic_driver.F90, applyBoundaryFluxesInOut(): "//&
1233 "Complete mass loss in column!")
1239 elseif ((abs(netheat(i))+abs(netsalt(i))+abs(netmassin(i))+abs(netmassout(i)))>0.)
then
1241 if (.not. cs%ignore_fluxes_over_land)
then
1243 write(0,*)
'applyBoundaryFluxesInOut(): lon,lat=',g%geoLonT(i,j),g%geoLatT(i,j)
1244 write(0,*)
'applyBoundaryFluxesInOut(): netHeat,netSalt,netMassIn,netMassOut=',&
1245 netheat(i),netsalt(i),netmassin(i),netmassout(i)
1247 call mom_error(fatal,
"MOM_diabatic_driver.F90, applyBoundaryFluxesInOut(): "//&
1248 "Mass loss over land?")
1254 if (netmassin(i)+netmassout(i) /= 0.0)
then
1256 numberofgroundings = numberofgroundings +1
1257 if (numberofgroundings<=maxgroundings)
then
1258 iground(numberofgroundings) = i
1259 jground(numberofgroundings) = j
1260 hgrounding(numberofgroundings) = netmassin(i)+netmassout(i)
1263 if (cs%id_createdH>0) cs%createdH(i,j) = cs%createdH(i,j) - (netmassin(i)+netmassout(i))/dt
1274 if (cs%id_penSW_diag > 0 .or. cs%id_penSWflux_diag > 0)
then
1275 do k=1,nz ;
do i=is,ie
1276 cs%penSW_diag(i,j,k) = t2d(i,k)
1277 cs%penSWflux_diag(i,j,k) = 0.0
1280 cs%penSWflux_diag(i,j,k) = 0.0
1284 if (calculate_energetics)
then
1285 call absorbremainingsw(g, gv, us, h2d, opacityband, nsw, optics, j, dt, h_limit_fluxes, &
1286 .false., .true., t2d, pen_sw_bnd, tke=pen_tke_2d, dsv_dt=dsv_dt_2d)
1288 do k=1,nz ;
do i=is,ie
1289 ctke(i,j,k) = ctke(i,j,k) + pen_tke_2d(i,k)
1292 call absorbremainingsw(g, gv, us, h2d, opacityband, nsw, optics, j, dt, h_limit_fluxes, &
1293 .false., .true., t2d, pen_sw_bnd)
1299 do k=1,nz ;
do i=is,ie
1301 tv%T(i,j,k) = t2d(i,k)
1306 if (cs%id_penSW_diag > 0 .or. cs%id_penSWflux_diag > 0)
then
1309 do k=1,nz ;
do i=is,ie
1310 cs%penSW_diag(i,j,k) = (t2d(i,k)-cs%penSW_diag(i,j,k))*h(i,j,k) * us%s_to_T*idt * tv%C_p * gv%H_to_kg_m2
1319 if (cs%id_penSWflux_diag > 0)
then
1320 do k=nz,1,-1 ;
do i=is,ie
1321 cs%penSWflux_diag(i,j,k) = cs%penSW_diag(i,j,k) + cs%penSWflux_diag(i,j,k+1)
1328 if (cs%id_nonpenSW_diag > 0)
then
1330 cs%nonpenSW_diag(i,j) = nonpensw(i)
1339 if (calculate_buoyancy)
then
1342 netpen_rate(:) = 0.0
1349 do i=is,ie ;
do nb=1,nsw ; netpen_rate(i) = netpen_rate(i) + pen_sw_bnd_rate(nb,i) ;
enddo ;
enddo
1353 drhodt, drhods, start, npts, tv%eqn_of_state, scale=us%kg_m3_to_R)
1360 skinbuoyflux(i,j) = - gorho * gv%H_to_Z * &
1361 (drhods(i) * (netsalt_rate(i) - tv%S(i,j,1)*netmassinout_rate(i)) + &
1362 drhodt(i) * ( netheat_rate(i) + netpen_rate(i)) )
1369 if (cs%id_createdH > 0)
call post_data(cs%id_createdH , cs%createdH , cs%diag)
1370 if (cs%id_penSW_diag > 0)
call post_data(cs%id_penSW_diag , cs%penSW_diag , cs%diag)
1371 if (cs%id_penSWflux_diag > 0)
call post_data(cs%id_penSWflux_diag, cs%penSWflux_diag, cs%diag)
1372 if (cs%id_nonpenSW_diag > 0)
call post_data(cs%id_nonpenSW_diag , cs%nonpenSW_diag , cs%diag)
1375 if (numberofgroundings>0 .and. .not. cs%ignore_fluxes_over_land)
then
1376 do i = 1, min(numberofgroundings, maxgroundings)
1378 write(mesg(1:45),
'(3es15.3)') g%geoLonT( iground(i), jground(i) ), &
1379 g%geoLatT( iground(i), jground(i)) , hgrounding(i)
1380 call mom_error(warning,
"MOM_diabatic_driver.F90, applyBoundaryFluxesInOut(): "//&
1381 "Mass created. x,y,dh= "//trim(mesg), all_print=.true.)
1384 if (numberofgroundings - maxgroundings > 0)
then
1385 write(mesg,
'(i4)') numberofgroundings - maxgroundings
1386 call mom_error(warning,
"MOM_diabatic_driver:F90, applyBoundaryFluxesInOut(): "//&
1387 trim(mesg) //
" groundings remaining")
1394 subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgorithm, use_ePBL)
1395 type(time_type),
target,
intent(in) :: time
1400 type(
diag_ctrl),
target,
intent(inout) :: diag
1403 logical,
intent(in) :: usealealgorithm
1405 logical,
intent(in) :: use_epbl
1410 #include "version_variable.h"
1411 character(len=40) :: mdl =
"MOM_diabatic_aux"
1412 character(len=48) :: thickness_units
1413 character(len=200) :: inputdir
1414 character(len=240) :: chl_filename
1415 character(len=128) :: chl_file
1417 character(len=32) :: chl_varname
1418 logical :: use_temperature
1419 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz, nbands
1420 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
1421 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1423 if (
associated(cs))
then
1424 call mom_error(warning,
"diabatic_aux_init called with an "// &
1425 "associated control structure.")
1436 "The following parameters are used for auxiliary diabatic processes.")
1438 call get_param(param_file, mdl,
"ENABLE_THERMODYNAMICS", use_temperature, &
1439 "If true, temperature and salinity are used as state "//&
1440 "variables.", default=.true.)
1442 call get_param(param_file, mdl,
"RECLAIM_FRAZIL", cs%reclaim_frazil, &
1443 "If true, try to use any frazil heat deficit to cool any "//&
1444 "overlying layers down to the freezing point, thereby "//&
1445 "avoiding the creation of thin ice when the SST is above "//&
1446 "the freezing point.", default=.true.)
1447 call get_param(param_file, mdl,
"PRESSURE_DEPENDENT_FRAZIL", &
1448 cs%pressure_dependent_frazil, &
1449 "If true, use a pressure dependent freezing temperature "//&
1450 "when making frazil. The default is false, which will be "//&
1451 "faster but is inappropriate with ice-shelf cavities.", &
1455 call get_param(param_file, mdl,
"IGNORE_FLUXES_OVER_LAND", cs%ignore_fluxes_over_land,&
1456 "If true, the model does not check if fluxes are being applied "//&
1457 "over land points. This is needed when the ocean is coupled "//&
1458 "with ice shelves and sea ice, since the sea ice mask needs to "//&
1459 "be different than the ocean mask to avoid sea ice formation "//&
1460 "under ice shelves. This flag only works when use_ePBL = True.", default=.false.)
1461 call get_param(param_file, mdl,
"DO_RIVERMIX", cs%do_rivermix, &
1462 "If true, apply additional mixing wherever there is "//&
1463 "runoff, so that it is mixed down to RIVERMIX_DEPTH "//&
1464 "if the ocean is that deep.", default=.false.)
1465 if (cs%do_rivermix) &
1466 call get_param(param_file, mdl,
"RIVERMIX_DEPTH", cs%rivermix_depth, &
1467 "The depth to which rivers are mixed if DO_RIVERMIX is "//&
1468 "defined.", units=
"m", default=0.0, scale=us%m_to_Z)
1470 cs%do_rivermix = .false. ; cs%rivermix_depth = 0.0 ; cs%ignore_fluxes_over_land = .false.
1473 if (gv%nkml == 0)
then
1474 call get_param(param_file, mdl,
"USE_RIVER_HEAT_CONTENT", cs%use_river_heat_content, &
1475 "If true, use the fluxes%runoff_Hflx field to set the "//&
1476 "heat carried by runoff, instead of using SST*CP*liq_runoff.", &
1478 call get_param(param_file, mdl,
"USE_CALVING_HEAT_CONTENT", cs%use_calving_heat_content, &
1479 "If true, use the fluxes%calving_Hflx field to set the "//&
1480 "heat carried by runoff, instead of using SST*CP*froz_runoff.", &
1483 cs%use_river_heat_content = .false.
1484 cs%use_calving_heat_content = .false.
1487 if (usealealgorithm)
then
1488 cs%id_createdH = register_diag_field(
'ocean_model',
"created_H",diag%axesT1, &
1489 time,
"The volume flux added to stop the ocean from drying out and becoming negative in depth", &
1490 "m s-1", conversion=gv%H_to_m*us%s_to_T)
1491 if (cs%id_createdH>0)
allocate(cs%createdH(isd:ied,jsd:jed))
1494 cs%id_penSW_diag = register_diag_field(
'ocean_model',
'rsdoabsorb', &
1495 diag%axesTL, time,
'Convergence of Penetrative Shortwave Flux in Sea Water Layer',&
1496 'W m-2', standard_name=
'net_rate_of_absorption_of_shortwave_energy_in_ocean_layer',v_extensive=.true.)
1500 cs%id_penSWflux_diag = register_diag_field(
'ocean_model',
'rsdo', &
1501 diag%axesTi, time,
'Downwelling Shortwave Flux in Sea Water at Grid Cell Upper Interface',&
1502 'W m-2', standard_name=
'downwelling_shortwave_flux_in_sea_water')
1505 if (cs%id_penSW_diag>0 .or. cs%id_penSWflux_diag>0)
then
1506 allocate(cs%penSW_diag(isd:ied,jsd:jed,nz))
1507 cs%penSW_diag(:,:,:) = 0.0
1508 allocate(cs%penSWflux_diag(isd:ied,jsd:jed,nz+1))
1509 cs%penSWflux_diag(:,:,:) = 0.0
1513 cs%id_nonpenSW_diag = register_diag_field(
'ocean_model',
'nonpenSW', &
1514 diag%axesT1, time, &
1515 'Non-downwelling SW radiation (i.e., SW absorbed in ocean surface with LW,SENS,LAT)',&
1516 'W m-2', standard_name=
'nondownwelling_shortwave_flux_in_sea_water')
1517 if (cs%id_nonpenSW_diag > 0)
then
1518 allocate(cs%nonpenSW_diag(isd:ied,jsd:jed))
1519 cs%nonpenSW_diag(:,:) = 0.0
1523 if (use_temperature)
then
1524 call get_param(param_file, mdl,
"VAR_PEN_SW", cs%var_pen_sw, &
1525 "If true, use one of the CHL_A schemes specified by "//&
1526 "OPACITY_SCHEME to determine the e-folding depth of "//&
1527 "incoming short wave radiation.", default=.false.)
1528 if (cs%var_pen_sw)
then
1530 call get_param(param_file, mdl,
"CHL_FROM_FILE", cs%chl_from_file, &
1531 "If true, chl_a is read from a file.", default=.true.)
1532 if (cs%chl_from_file)
then
1533 call time_interp_external_init()
1535 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
1536 call get_param(param_file, mdl,
"CHL_FILE", chl_file, &
1537 "CHL_FILE is the file containing chl_a concentrations in "//&
1538 "the variable CHL_A. It is used when VAR_PEN_SW and "//&
1539 "CHL_FROM_FILE are true.", fail_if_missing=.true.)
1540 chl_filename = trim(slasher(inputdir))//trim(chl_file)
1541 call log_param(param_file, mdl,
"INPUTDIR/CHL_FILE", chl_filename)
1542 call get_param(param_file, mdl,
"CHL_VARNAME", chl_varname, &
1543 "Name of CHL_A variable in CHL_FILE.", default=
'CHL_A')
1544 cs%sbc_chl = init_external_field(chl_filename, trim(chl_varname), domain=g%Domain%mpp_domain)
1547 cs%id_chl = register_diag_field(
'ocean_model',
'Chl_opac', diag%axesT1, time, &
1548 'Surface chlorophyll A concentration used to find opacity',
'mg m-3')
1552 id_clock_uv_at_h = cpu_clock_id(
'(Ocean find_uv_at_h)', grain=clock_routine)
1563 if (.not.
associated(cs))
return
1565 if (cs%id_createdH >0)
deallocate(cs%createdH)
1566 if (cs%id_penSW_diag >0)
deallocate(cs%penSW_diag)
1567 if (cs%id_penSWflux_diag >0)
deallocate(cs%penSWflux_diag)
1568 if (cs%id_nonpenSW_diag >0)
deallocate(cs%nonpenSW_diag)
1570 if (
associated(cs))
deallocate(cs)