16 implicit none ;
private
18 #include <MOM_memory.h>
28 real,
pointer,
dimension(:,:,:,:) :: opacity_band => null()
31 real,
pointer,
dimension(:,:,:) :: sw_pen_band => null()
35 real,
pointer,
dimension(:) :: &
36 min_wavelength_band => null(), &
37 max_wavelength_band => null()
39 real :: pensw_flux_absorb
41 real :: pensw_absorb_invlen
43 logical :: answers_2018
53 integer :: opacity_scheme
58 real :: pen_sw_scale_2nd
60 real :: sw_1st_exp_ratio
65 real :: opacity_land_value
71 integer :: id_sw_pen = -1, id_sw_vis_pen = -1
72 integer,
pointer :: id_opacity(:) => null()
91 subroutine set_opacity(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_nir_dif, &
92 G, GV, CS, chl_2d, chl_3d)
95 real,
dimension(:,:),
pointer :: sw_total
96 real,
dimension(:,:),
pointer :: sw_vis_dir
97 real,
dimension(:,:),
pointer :: sw_vis_dif
98 real,
dimension(:,:),
pointer :: sw_nir_dir
99 real,
dimension(:,:),
pointer :: sw_nir_dif
104 real,
dimension(SZI_(G),SZJ_(G)), &
105 optional,
intent(in) :: chl_2d
106 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
107 optional,
intent(in) :: chl_3d
110 integer :: i, j, k, n, is, ie, js, je, nz
111 real :: inv_sw_pen_scale
114 logical :: call_for_surface
115 real :: tmp(szi_(g),szj_(g),szk_(gv))
116 real :: chl(szi_(g),szj_(g),szk_(gv))
117 real :: pen_sw_tot(szi_(g),szj_(g))
119 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
121 if (.not.
associated(cs))
call mom_error(fatal,
"set_opacity: "// &
122 "Module must be initialized via opacity_init before it is used.")
124 if (
present(chl_2d) .or.
present(chl_3d))
then
126 call opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_nir_dif, &
127 g, gv, cs, chl_2d, chl_3d)
129 if (optics%nbands <= 1)
then ; inv_nbands = 1.0
130 else ; inv_nbands = 1.0 / real(optics%nbands) ;
endif
133 inv_sw_pen_scale = 1.0 / max(cs%pen_sw_scale, 0.1*gv%Angstrom_m, &
134 gv%H_to_m*gv%H_subroundoff)
137 do k=1,nz ;
do j=js,je ;
do i=is,ie
138 optics%opacity_band(1,i,j,k) = inv_sw_pen_scale
139 optics%opacity_band(2,i,j,k) = 1.0 / max(cs%pen_sw_scale_2nd, &
140 0.1*gv%Angstrom_m,gv%H_to_m*gv%H_subroundoff)
141 enddo ;
enddo ;
enddo
142 if (.not.
associated(sw_total) .or. (cs%pen_SW_scale <= 0.0))
then
144 do j=js,je ;
do i=is,ie ;
do n=1,optics%nbands
145 optics%sw_pen_band(n,i,j) = 0.0
146 enddo ;
enddo ;
enddo
149 do j=js,je ;
do i=is,ie
150 optics%sw_pen_band(1,i,j) = (cs%SW_1st_EXP_RATIO) * sw_total(i,j)
151 optics%sw_pen_band(2,i,j) = (1.-cs%SW_1st_EXP_RATIO) * sw_total(i,j)
155 do k=1,nz ;
do j=js,je ;
do i=is,ie ;
do n=1,optics%nbands
156 optics%opacity_band(n,i,j,k) = inv_sw_pen_scale
157 enddo ;
enddo ;
enddo ;
enddo
158 if (.not.
associated(sw_total) .or. (cs%pen_SW_scale <= 0.0))
then
160 do j=js,je ;
do i=is,ie ;
do n=1,optics%nbands
161 optics%sw_pen_band(n,i,j) = 0.0
162 enddo ;
enddo ;
enddo
165 do j=js,je ;
do i=is,ie ;
do n=1,optics%nbands
166 optics%sw_pen_band(n,i,j) = cs%pen_SW_frac * inv_nbands * sw_total(i,j)
167 enddo ;
enddo ;
enddo
173 if (cs%id_sw_pen > 0)
then
175 do j=js,je ;
do i=is,ie
176 pen_sw_tot(i,j) = 0.0
178 pen_sw_tot(i,j) = pen_sw_tot(i,j) + optics%sw_pen_band(n,i,j)
181 call post_data(cs%id_sw_pen, pen_sw_tot, cs%diag)
183 if (cs%id_sw_vis_pen > 0)
then
186 do j=js,je ;
do i=is,ie
187 pen_sw_tot(i,j) = 0.0
188 do n=1,min(optics%nbands,2)
189 pen_sw_tot(i,j) = pen_sw_tot(i,j) + optics%sw_pen_band(n,i,j)
194 do j=js,je ;
do i=is,ie
195 pen_sw_tot(i,j) = 0.0
197 pen_sw_tot(i,j) = pen_sw_tot(i,j) + optics%sw_pen_band(n,i,j)
201 call post_data(cs%id_sw_vis_pen, pen_sw_tot, cs%diag)
203 do n=1,optics%nbands ;
if (cs%id_opacity(n) > 0)
then
205 do k=1,nz ;
do j=js,je ;
do i=is,ie
210 enddo ;
enddo ;
enddo
211 call post_data(cs%id_opacity(n), tmp, cs%diag)
220 subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_nir_dif, &
221 G, GV, CS, chl_2d, chl_3d)
224 real,
dimension(:,:),
pointer :: sw_total
225 real,
dimension(:,:),
pointer :: sw_vis_dir
226 real,
dimension(:,:),
pointer :: sw_vis_dif
227 real,
dimension(:,:),
pointer :: sw_nir_dir
228 real,
dimension(:,:),
pointer :: sw_nir_dif
232 real,
dimension(SZI_(G),SZJ_(G)), &
233 optional,
intent(in) :: chl_2d
234 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
235 optional,
intent(in) :: chl_3d
237 real :: chl_data(SZI_(G),SZJ_(G))
240 real :: Inv_nbands_nir
248 type(time_type) :: day
249 character(len=128) :: mesg
250 integer :: i, j, k, n, is, ie, js, je, nz, nbands
251 logical :: multiband_vis_input, multiband_nir_input
253 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
267 nbands = optics%nbands
269 if (nbands <= 1)
then ; inv_nbands = 1.0
270 else ; inv_nbands = 1.0 / real(nbands) ;
endif
272 if (nbands <= 2)
then ; inv_nbands_nir = 0.0
273 else ; inv_nbands_nir = 1.0 / real(nbands - 2.0) ;
endif
275 multiband_vis_input = (
associated(sw_vis_dir) .and. &
276 associated(sw_vis_dif))
277 multiband_nir_input = (
associated(sw_nir_dir) .and. &
278 associated(sw_nir_dif))
281 if (
present(chl_3d))
then
282 do j=js,je ;
do i=is,ie ; chl_data(i,j) = chl_3d(i,j,1) ;
enddo ;
enddo
283 do k=1,nz;
do j=js,je ;
do i=is,ie
284 if ((g%mask2dT(i,j) > 0.5) .and. (chl_3d(i,j,k) < 0.0))
then
285 write(mesg,
'(" Negative chl_3d of ",(1pe12.4)," found at i,j,k = ", &
286 & 3(1x,i3), " lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
287 chl_3d(i,j,k), i, j, k, g%geoLonT(i,j), g%geoLatT(i,j)
288 call mom_error(fatal,
"MOM_opacity opacity_from_chl: "//trim(mesg))
290 enddo ;
enddo ;
enddo
291 elseif (
present(chl_2d))
then
292 do j=js,je ;
do i=is,ie ; chl_data(i,j) = chl_2d(i,j) ;
enddo ;
enddo
293 do j=js,je ;
do i=is,ie
294 if ((g%mask2dT(i,j) > 0.5) .and. (chl_2d(i,j) < 0.0))
then
295 write(mesg,
'(" Negative chl_2d of ",(1pe12.4)," at i,j = ", &
296 & 2(i3), "lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
297 chl_data(i,j), i, j, g%geoLonT(i,j), g%geoLatT(i,j)
298 call mom_error(fatal,
"MOM_opacity opacity_from_chl: "//trim(mesg))
302 call mom_error(fatal,
"Either chl_2d or chl_3d must be preesnt in a call to opacity_form_chl.")
305 select case (cs%opacity_scheme)
308 do j=js,je ;
do i=is,ie
309 sw_vis_tot = 0.0 ; sw_nir_tot = 0.0
310 if (g%mask2dT(i,j) > 0.5)
then
311 if (multiband_vis_input)
then
312 sw_vis_tot = sw_vis_dir(i,j) + sw_vis_dif(i,j)
314 sw_vis_tot = 0.42 * sw_total(i,j)
316 if (multiband_nir_input)
then
317 sw_nir_tot = sw_nir_dir(i,j) + sw_nir_dif(i,j)
319 sw_nir_tot = sw_total(i,j) - sw_vis_tot
324 optics%sw_pen_band(1,i,j) = cs%blue_frac*sw_vis_tot
327 optics%sw_pen_band(2,i,j) = (1.0-cs%blue_frac)*sw_vis_tot
330 optics%sw_pen_band(n,i,j) = inv_nbands_nir * sw_nir_tot
335 do j=js,je ;
do i=is,ie
337 if (g%mask2dT(i,j) > 0.5)
then ;
if (multiband_vis_input)
then
339 (sw_vis_dir(i,j) + sw_vis_dif(i,j))
346 optics%sw_pen_band(n,i,j) = inv_nbands*sw_pen_tot
350 call mom_error(fatal,
"opacity_from_chl: CS%opacity_scheme is not valid.")
355 if (
present(chl_3d))
then
356 do j=js,je ;
do i=is,ie ; chl_data(i,j) = chl_3d(i,j,k) ;
enddo ;
enddo
359 select case (cs%opacity_scheme)
361 do j=js,je ;
do i=is,ie
362 if (g%mask2dT(i,j) <= 0.5)
then
364 optics%opacity_band(n,i,j,k) = cs%opacity_land_value
368 optics%opacity_band(1,i,j,k) = 0.0232 + 0.074*chl_data(i,j)**0.674
370 optics%opacity_band(2,i,j,k) = 0.225 + 0.037*chl_data(i,j)**0.629
372 do n=3,nbands ; optics%opacity_band(n,i,j,k) = 2.86 ;
enddo
376 do j=js,je ;
do i=is,ie
377 optics%opacity_band(1,i,j,k) = cs%opacity_land_value
378 if (g%mask2dT(i,j) > 0.5) &
382 optics%opacity_band(n,i,j,k) = optics%opacity_band(1,i,j,k)
387 call mom_error(fatal,
"opacity_from_chl: CS%opacity_scheme is not valid.")
397 real,
intent(in) :: chl_data
405 real,
dimension(6),
parameter :: &
406 z2_coef=(/7.925, -6.644, 3.662, -1.815, -0.218, 0.502/)
409 chl = log10(min(max(chl_data,0.02),60.0)) ; chl2 = chl*chl
410 opacity_morel = 1.0 / ( (z2_coef(1) + z2_coef(2)*chl) + chl2 * &
411 ((z2_coef(3) + chl*z2_coef(4)) + chl2*(z2_coef(5) + chl*z2_coef(6))) )
417 real,
intent(in) :: chl_data
426 real,
dimension(6),
parameter :: &
427 v1_coef=(/0.321, 0.008, 0.132, 0.038, -0.017, -0.007/)
429 chl = log10(min(max(chl_data,0.02),60.0)) ; chl2 = chl*chl
431 ((v1_coef(3) + chl*v1_coef(4)) + chl2*(v1_coef(5) + chl*v1_coef(6))) )
437 real,
intent(in) :: chl_data
446 subroutine extract_optics_slice(optics, j, G, GV, opacity, opacity_scale, penSW_top, penSW_scale)
449 integer,
intent(in) :: j
452 real,
dimension(max(optics%nbands,1),SZI_(G),SZK_(GV)), &
453 optional,
intent(out) :: opacity
454 real,
optional,
intent(in) :: opacity_scale
455 real,
dimension(max(optics%nbands,1),SZI_(G)), &
456 optional,
intent(out) :: pensw_top
459 real,
optional,
intent(in) :: pensw_scale
462 real :: scale_opacity, scale_pensw
463 integer :: i, is, ie, k, nz, n
464 is = g%isc ; ie = g%iec ; nz = g%ke
466 scale_opacity = 1.0 ;
if (
present(opacity_scale)) scale_opacity = opacity_scale
467 scale_pensw = 1.0 ;
if (
present(pensw_scale)) scale_pensw = pensw_scale
469 if (
present(opacity))
then ;
do k=1,nz ;
do i=is,ie
471 opacity(n,i,k) = scale_opacity * optics%opacity_band(n,i,j,k)
473 enddo ;
enddo ;
endif
475 if (
present(pensw_top))
then ;
do k=1,nz ;
do i=is,ie
477 pensw_top(n,i) = scale_pensw * optics%SW_pen_band(n,i,j)
479 enddo ;
enddo ;
endif
487 integer,
optional,
intent(out) :: nbands
489 if (
present(nbands)) nbands = optics%nbands
509 subroutine absorbremainingsw(G, GV, US, h, opacity_band, nsw, optics, j, dt, H_limit_fluxes, &
510 adjustAbsorptionProfile, absorbAllSW, T, Pen_SW_bnd, &
511 eps, ksort, htot, Ttot, TKE, dSV_dT)
516 integer,
intent(in) :: nsw
518 real,
dimension(SZI_(G),SZK_(GV)),
intent(in) :: h
519 real,
dimension(max(1,nsw),SZI_(G),SZK_(GV)),
intent(in) :: opacity_band
524 integer,
intent(in) :: j
525 real,
intent(in) :: dt
526 real,
intent(in) :: h_limit_fluxes
532 logical,
intent(in) :: adjustabsorptionprofile
538 logical,
intent(in) :: absorballsw
544 real,
dimension(SZI_(G),SZK_(GV)),
intent(inout) :: t
546 real,
dimension(max(1,nsw),SZI_(G)),
intent(inout) :: pen_sw_bnd
551 real,
dimension(SZI_(G),SZK_(GV)),
optional,
intent(in) :: eps
554 integer,
dimension(SZI_(G),SZK_(GV)),
optional,
intent(in) :: ksort
555 real,
dimension(SZI_(G)),
optional,
intent(in) :: htot
556 real,
dimension(SZI_(G)),
optional,
intent(inout) :: ttot
558 real,
dimension(SZI_(G),SZK_(GV)),
optional,
intent(in) :: dsv_dt
560 real,
dimension(SZI_(G),SZK_(GV)),
optional,
intent(inout) :: tke
564 real,
dimension(SZI_(G),SZK_(GV)) :: &
572 real,
dimension(SZI_(G)) :: &
603 logical :: sw_remains
608 integer :: is, ie, nz, i, k, ks, n
611 min_sw_heat = optics%PenSW_flux_absorb * dt
612 i_habs = optics%PenSW_absorb_Invlen
614 h_min_heat = 2.0*gv%Angstrom_H + gv%H_subroundoff
615 is = g%isc ; ie = g%iec ; nz = g%ke
616 c1_6 = 1.0 / 6.0 ; c1_60 = 1.0 / 60.0
618 tke_calc = (
present(tke) .and.
present(dsv_dt))
620 if (optics%answers_2018)
then
621 g_hconv2 = (us%L_to_Z**2*gv%g_Earth * gv%H_to_RZ) * gv%H_to_RZ
623 g_hconv2 = us%L_to_Z**2*gv%g_Earth * gv%H_to_RZ**2
627 if (
present(htot))
then ;
do i=is,ie ; h_heat(i) = htot(i) ;
enddo ;
endif
631 do ks=1,nz ;
do i=is,ie
633 if (
present(ksort))
then
634 if (ksort(i,ks) <= 0) cycle
637 epsilon = 0.0 ;
if (
present(eps)) epsilon = eps(i,k)
639 t_chg_above(i,k) = 0.0
641 if (h(i,k) > 1.5*epsilon)
then
642 do n=1,nsw ;
if (pen_sw_bnd(n,i) > 0.0)
then
644 opt_depth = h(i,k) * opacity_band(n,i,k)
645 exp_od = exp(-opt_depth)
650 if (optics%answers_2018)
then
651 if (nsw*pen_sw_bnd(n,i)*sw_trans < min_sw_heat*min(1.0, i_habs*h(i,k)) ) sw_trans = 0.0
652 elseif ((nsw*pen_sw_bnd(n,i)*sw_trans < min_sw_heat) .and. (h(i,k) > h_min_heat))
then
653 if (nsw*pen_sw_bnd(n,i) <= min_sw_heat * (i_habs*(h(i,k) - h_min_heat)))
then
656 sw_trans = min(sw_trans, &
657 1.0 - (min_sw_heat*(i_habs*(h(i,k) - h_min_heat))) / (nsw*pen_sw_bnd(n,i)))
661 heat_bnd = pen_sw_bnd(n,i) * (1.0 - sw_trans)
662 if (adjustabsorptionprofile .and. (h_heat(i) > 0.0))
then
673 if (opt_depth > 1e-5)
then
674 swa = ((opt_depth + (opt_depth + 2.0)*exp_od) - 2.0) / &
675 ((opt_depth + opacity_band(n,i,k) * h_heat(i)) * &
680 swa = h(i,k) * (opt_depth * (1.0 - opt_depth)) / &
681 ((h_heat(i) + h(i,k)) * (6.0 - 3.0*opt_depth))
684 if (swa*(h_heat(i) + h(i,k)) > h_heat(i))
then
685 coswa_frac = (swa*(h_heat(i) + h(i,k)) - h_heat(i) ) / &
686 (swa*(h_heat(i) + h(i,k)))
687 swa = h_heat(i) / (h_heat(i) + h(i,k))
690 t_chg_above(i,k) = t_chg_above(i,k) + (swa * heat_bnd) / h_heat(i)
691 t(i,k) = t(i,k) + ((1.0 - swa) * heat_bnd) / h(i,k)
694 t(i,k) = t(i,k) + pen_sw_bnd(n,i) * (1.0 - sw_trans) / h(i,k)
698 if (opt_depth > 1e-2)
then
699 tke(i,k) = tke(i,k) - coswa_frac*heat_bnd*dsv_dt(i,k)* &
700 (0.5*h(i,k)*g_hconv2) * &
701 (opt_depth*(1.0+exp_od) - 2.0*(1.0-exp_od)) / (opt_depth*(1.0-exp_od))
705 tke(i,k) = tke(i,k) - coswa_frac*heat_bnd*dsv_dt(i,k)* &
706 (0.5*h(i,k)*g_hconv2) * &
707 (c1_6*opt_depth * (1.0 - c1_60*opt_depth**2))
711 pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
717 if (h(i,k) >= 2.0*h_min_heat)
then
718 h_heat(i) = h_heat(i) + h(i,k)
719 elseif (h(i,k) > h_min_heat)
then
720 h_heat(i) = h_heat(i) + (2.0*h(i,k) - 2.0*h_min_heat)
728 do i=is,ie ; t_chg(i) = 0.0 ;
enddo
730 if (absorballsw)
then
735 pen_sw_rem(i) = pen_sw_bnd(1,i)
736 do n=2,nsw ; pen_sw_rem(i) = pen_sw_rem(i) + pen_sw_bnd(n,i) ;
enddo
738 do i=is,ie ;
if (pen_sw_rem(i) > 0.0) sw_remains = .true. ;
enddo
740 ih_limit = 1.0 / h_limit_fluxes
741 do i=is,ie ;
if ((pen_sw_rem(i) > 0.0) .and. (h_heat(i) > 0.0))
then
742 if (h_heat(i)*ih_limit >= 1.0)
then
743 t_chg(i) = pen_sw_rem(i) / h_heat(i) ; unabsorbed = 0.0
745 t_chg(i) = pen_sw_rem(i) * ih_limit
746 unabsorbed = 1.0 - h_heat(i)*ih_limit
748 do n=1,nsw ; pen_sw_bnd(n,i) = unabsorbed * pen_sw_bnd(n,i) ;
enddo
752 if (absorballsw .or. adjustabsorptionprofile)
then
753 do ks=nz,1,-1 ;
do i=is,ie
755 if (
present(ksort))
then
756 if (ksort(i,ks) <= 0) cycle
760 if (t_chg(i) > 0.0)
then
762 if (h(i,k) >= 2.0*h_min_heat)
then ; t(i,k) = t(i,k) + t_chg(i)
763 elseif (h(i,k) > h_min_heat)
then
764 t(i,k) = t(i,k) + t_chg(i) * (2.0 - 2.0*h_min_heat/h(i,k))
768 t_chg(i) = t_chg(i) + t_chg_above(i,k)
770 if (
present(htot) .and.
present(ttot))
then
771 do i=is,ie ; ttot(i) = ttot(i) + t_chg(i) * htot(i) ;
enddo
782 H_limit_fluxes, absorbAllSW, iPen_SW_bnd, netPen)
786 real,
dimension(SZI_(G),SZK_(GV)), &
788 integer,
intent(in) :: nsw
792 integer,
intent(in) :: j
793 real,
intent(in) :: dt
794 real,
intent(in) :: h_limit_fluxes
797 logical,
intent(in) :: absorballsw
799 real,
dimension(max(nsw,1),SZI_(G)),
intent(in) :: ipen_sw_bnd
803 real,
dimension(SZI_(G),SZK_(GV)+1), &
804 intent(inout) :: netpen
808 real :: h_heat(szi_(g))
810 real :: pen_sw_rem(szi_(g))
815 real,
dimension(max(nsw,1),SZI_(G)) :: pen_sw_bnd
830 logical :: sw_remains
833 integer :: is, ie, nz, i, k, ks, n
836 min_sw_heat = optics%PenSW_flux_absorb*dt
837 i_habs = 1e3*gv%H_to_m
839 h_min_heat = 2.0*gv%Angstrom_H + gv%H_subroundoff
840 is = g%isc ; ie = g%iec ; nz = g%ke
842 pen_sw_bnd(:,:) = ipen_sw_bnd(:,:)
843 do i=is,ie ; h_heat(i) = 0.0 ;
enddo
847 netpen(i,1) = netpen(i,1) + pen_sw_bnd(n,i)
858 if (h(i,k) > 0.0)
then
859 do n=1,nsw ;
if (pen_sw_bnd(n,i) > 0.0)
then
861 opt_depth = h(i,k)*gv%H_to_m * optics%opacity_band(n,i,j,k)
862 exp_od = exp(-opt_depth)
867 if (optics%answers_2018)
then
868 if (nsw*pen_sw_bnd(n,i)*sw_trans < min_sw_heat*min(1.0, i_habs*h(i,k)) ) sw_trans = 0.0
869 elseif ((nsw*pen_sw_bnd(n,i)*sw_trans < min_sw_heat) .and. (h(i,k) > h_min_heat))
then
870 if (nsw*pen_sw_bnd(n,i) <= min_sw_heat * (i_habs*(h(i,k) - h_min_heat)))
then
873 sw_trans = min(sw_trans, &
874 1.0 - (min_sw_heat*(i_habs*(h(i,k) - h_min_heat))) / (nsw*pen_sw_bnd(n,i)))
878 pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
879 netpen(i,k+1) = netpen(i,k+1) + pen_sw_bnd(n,i)
885 if (h(i,k) >= 2.0*h_min_heat)
then
886 h_heat(i) = h_heat(i) + h(i,k)
887 elseif (h(i,k) > h_min_heat)
then
888 h_heat(i) = h_heat(i) + (2.0*h(i,k) - 2.0*h_min_heat)
893 if (absorballsw)
then
899 pen_sw_rem(i) = pen_sw_bnd(1,i)
900 do n=2,nsw ; pen_sw_rem(i) = pen_sw_rem(i) + pen_sw_bnd(n,i) ;
enddo
902 do i=is,ie ;
if (pen_sw_rem(i) > 0.0) sw_remains = .true. ;
enddo
904 ih_limit = 1.0 / h_limit_fluxes
905 do i=is,ie ;
if ((pen_sw_rem(i) > 0.0) .and. (h_heat(i) > 0.0))
then
906 if (h_heat(i)*ih_limit < 1.0)
then
907 unabsorbed = 1.0 - h_heat(i)*ih_limit
911 do n=1,nsw ; pen_sw_bnd(n,i) = unabsorbed * pen_sw_bnd(n,i) ;
enddo
921 subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics)
922 type(time_type),
target,
intent(in) :: time
928 type(
diag_ctrl),
target,
intent(inout) :: diag
935 character(len=200) :: tmpstr
936 character(len=40) :: mdl =
"MOM_opacity"
937 character(len=40) :: bandnum, shortname
938 character(len=200) :: longname
939 character(len=40) :: scheme_string
941 # include "version_variable.h"
942 real :: pensw_absorb_minthick
944 real :: pensw_minthick_dflt
945 logical :: default_2018_answers
946 logical :: use_scheme
947 integer :: isd, ied, jsd, jed, nz, n
948 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
950 if (
associated(cs))
then
951 call mom_error(warning,
"opacity_init called with an associated"// &
952 "associated control structure.")
954 else ;
allocate(cs) ;
endif
962 call get_param(param_file, mdl,
"VAR_PEN_SW", cs%var_pen_sw, &
963 "If true, use one of the CHL_A schemes specified by "//&
964 "OPACITY_SCHEME to determine the e-folding depth of "//&
965 "incoming short wave radiation.", default=.false.)
967 cs%opacity_scheme =
no_scheme ; scheme_string =
''
968 if (cs%var_pen_sw)
then
969 call get_param(param_file, mdl,
"OPACITY_SCHEME", tmpstr, &
970 "This character string specifies how chlorophyll "//&
971 "concentrations are translated into opacities. Currently "//&
972 "valid options include:\n"//&
973 " \t\t MANIZZA_05 - Use Manizza et al., GRL, 2005. \n"//&
974 " \t\t MOREL_88 - Use Morel, JGR, 1988.", &
976 if (len_trim(tmpstr) > 0)
then
984 call mom_error(fatal,
"opacity_init: #DEFINE OPACITY_SCHEME "//&
985 trim(tmpstr) //
"in input file is invalid.")
987 call mom_mesg(
'opacity_init: opacity scheme set to "'//trim(tmpstr)//
'".', 5)
990 call mom_error(warning,
"opacity_init: No scheme has successfully "//&
991 "been specified for the opacity. Using the default MANIZZA_05.")
995 call get_param(param_file, mdl,
"BLUE_FRAC_SW", cs%blue_frac, &
996 "The fraction of the penetrating shortwave radiation "//&
997 "that is in the blue band.", default=0.5, units=
"nondim")
999 call get_param(param_file, mdl,
"EXP_OPACITY_SCHEME", tmpstr, &
1000 "This character string specifies which exponential "//&
1001 "opacity scheme to utilize. Currently "//&
1002 "valid options include:\n"//&
1003 " \t\t SINGLE_EXP - Single Exponent decay. \n"//&
1004 " \t\t DOUBLE_EXP - Double Exponent decay.", &
1006 if (len_trim(tmpstr) > 0)
then
1008 select case (tmpstr)
1014 call mom_mesg(
'opacity_init: opacity scheme set to "'//trim(tmpstr)//
'".', 5)
1016 call get_param(param_file, mdl,
"PEN_SW_SCALE", cs%pen_sw_scale, &
1017 "The vertical absorption e-folding depth of the "//&
1018 "penetrating shortwave radiation.", units=
"m", default=0.0)
1021 call get_param(param_file, mdl,
"PEN_SW_SCALE_2ND", cs%pen_sw_scale_2nd, &
1022 "The (2nd) vertical absorption e-folding depth of the "//&
1023 "penetrating shortwave radiation "//&
1024 "(use if SW_EXP_MODE==double.)",&
1025 units=
"m", default=0.0)
1026 call get_param(param_file, mdl,
"SW_1ST_EXP_RATIO", cs%sw_1st_exp_ratio, &
1027 "The fraction of 1st vertical absorption e-folding depth "//&
1028 "penetrating shortwave radiation if SW_EXP_MODE==double.",&
1029 units=
"m", default=0.0)
1030 elseif (cs%OPACITY_SCHEME ==
single_exp)
then
1032 cs%pen_sw_scale_2nd = 0.0
1033 cs%sw_1st_exp_ratio = 1.0
1035 call get_param(param_file, mdl,
"PEN_SW_FRAC", cs%pen_sw_frac, &
1036 "The fraction of the shortwave radiation that penetrates "//&
1037 "below the surface.", units=
"nondim", default=0.0)
1040 call get_param(param_file, mdl,
"PEN_SW_NBANDS", optics%nbands, &
1041 "The number of bands of penetrating shortwave radiation.", &
1044 if (optics%nbands /= 2)
call mom_error(fatal, &
1045 "set_opacity: \Cannot use a double_exp opacity scheme with nbands!=2.")
1046 elseif (cs%Opacity_scheme ==
single_exp )
then
1047 if (optics%nbands /= 1)
call mom_error(fatal, &
1048 "set_opacity: \Cannot use a single_exp opacity scheme with nbands!=1.")
1051 call get_param(param_file, mdl,
"DEFAULT_2018_ANSWERS", default_2018_answers, &
1052 "This sets the default value for the various _2018_ANSWERS parameters.", &
1054 call get_param(param_file, mdl,
"OPTICS_2018_ANSWERS", optics%answers_2018, &
1055 "If true, use the order of arithmetic and expressions that recover the "//&
1056 "answers from the end of 2018. Otherwise, use updated expressions for "//&
1057 "handling the absorption of small remaining shortwave fluxes.", &
1058 default=default_2018_answers)
1060 call get_param(param_file, mdl,
"PEN_SW_FLUX_ABSORB", optics%PenSW_flux_absorb, &
1061 "A minimum remaining shortwave heating rate that will be simply absorbed in "//&
1062 "the next sufficiently thick layers for computational efficiency, instead of "//&
1063 "continuing to penetrate. The default, 2.5e-11 degC m s-1, is about 1e-4 W m-2 "//&
1064 "or 0.08 degC m century-1, but 0 is also a valid value.", &
1065 default=2.5e-11, units=
"degC m s-1", scale=gv%m_to_H*us%T_to_s)
1067 if (optics%answers_2018)
then ; pensw_minthick_dflt = 0.001 ;
else ; pensw_minthick_dflt = 1.0 ;
endif
1068 call get_param(param_file, mdl,
"PEN_SW_ABSORB_MINTHICK", pensw_absorb_minthick, &
1069 "A thickness that is used to absorb the remaining penetrating shortwave heat "//&
1070 "flux when it drops below PEN_SW_FLUX_ABSORB.", &
1071 default=pensw_minthick_dflt, units=
"m", scale=gv%m_to_H)
1072 optics%PenSW_absorb_Invlen = 1.0 / (pensw_absorb_minthick + gv%H_subroundoff)
1074 if (.not.
associated(optics%min_wavelength_band)) &
1075 allocate(optics%min_wavelength_band(optics%nbands))
1076 if (.not.
associated(optics%max_wavelength_band)) &
1077 allocate(optics%max_wavelength_band(optics%nbands))
1080 optics%min_wavelength_band(1) =0
1081 optics%max_wavelength_band(1) =550
1082 if (optics%nbands >= 2)
then
1083 optics%min_wavelength_band(2)=550
1084 optics%max_wavelength_band(2)=700
1086 if (optics%nbands > 2)
then
1087 do n=3,optics%nbands
1088 optics%min_wavelength_band(n) =700
1089 optics%max_wavelength_band(n) =2800
1094 call get_param(param_file, mdl,
"OPACITY_LAND_VALUE", cs%opacity_land_value, &
1095 "The value to use for opacity over land. The default is "//&
1096 "10 m-1 - a value for muddy water.", units=
"m-1", default=10.0)
1098 if (.not.
associated(optics%opacity_band)) &
1099 allocate(optics%opacity_band(optics%nbands,isd:ied,jsd:jed,nz))
1100 if (.not.
associated(optics%sw_pen_band)) &
1101 allocate(optics%sw_pen_band(optics%nbands,isd:ied,jsd:jed))
1102 allocate(cs%id_opacity(optics%nbands)) ; cs%id_opacity(:) = -1
1105 'Penetrating shortwave radiation flux into ocean',
'W m-2')
1107 'Visible penetrating shortwave radiation flux into ocean',
'W m-2')
1108 do n=1,optics%nbands
1109 write(bandnum,
'(i3)') n
1110 shortname =
'opac_'//trim(adjustl(bandnum))
1111 longname =
'Opacity for shortwave radiation in band '//trim(adjustl(bandnum)) &
1112 //
', saved as L^-1 tanh(Opacity * L) for L = 10^-10 m'
1124 if (
associated(cs%id_opacity))
deallocate(cs%id_opacity)
1125 if (
associated(cs))
deallocate(cs)
1127 if (
present(optics))
then ;
if (
associated(optics))
then
1128 if (
associated(optics%opacity_band))
deallocate(optics%opacity_band)
1129 if (
associated(optics%sw_pen_band))
deallocate(optics%sw_pen_band)