From a4779cc71125b40a7db3a4da8512247cbf2b0955 Mon Sep 17 00:00:00 2001 From: "David A. Bailey" Date: Tue, 14 Mar 2023 19:27:54 -0600 Subject: [PATCH] FSD bug fixes (#424) * FSD bug fixes for lateral melt and weld * put back declare sicen * another small welding fix * Change print statements to use icepack_warnings * Fix some comments * Fix some comments * Fix some comments * Check for wlat present when tr_fsd = .true. * Fix some syntax errors --------- Co-authored-by: cmbitz --- columnphysics/icepack_fsd.F90 | 14 ++- columnphysics/icepack_therm_itd.F90 | 125 +++++++++++++---------- columnphysics/icepack_therm_vertical.F90 | 29 +++--- configuration/driver/icedrv_flux.F90 | 1 + configuration/driver/icedrv_step.F90 | 6 +- 5 files changed, 96 insertions(+), 79 deletions(-) diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90 index bee46cd35..1454fdefc 100644 --- a/columnphysics/icepack_fsd.F90 +++ b/columnphysics/icepack_fsd.F90 @@ -699,7 +699,10 @@ subroutine fsd_add_new_ice (ncat, n, nfsd, & DO WHILE (elapsed_t.lt.dt) nsubt = nsubt + 1 - if (nsubt.gt.100) print *, 'latg not converging' + if (nsubt.gt.100) then + write(warnstr,*) subname,'latg not converging' + call icepack_warnings_add(warnstr) + endif ! finite differences df_flx(:) = c0 ! NB could stay zero if all in largest FS cat @@ -711,8 +714,6 @@ subroutine fsd_add_new_ice (ncat, n, nfsd, & df_flx(k) = f_flx(k+1) - f_flx(k) end do -! if (abs(sum(df_flx)) > puny) print*,'fsd_add_new ERROR df_flx /= 0' - dafsd_tmp(:) = c0 do k = 1, nfsd dafsd_tmp(k) = (-df_flx(k) + c2 * G_radial * afsdn_latg(k,n) & @@ -931,7 +932,6 @@ subroutine fsd_weld_thermo (ncat, nfsd, & gain, loss ! welding tendencies real(kind=dbl_kind) :: & - prefac , & ! multiplies kernel kern , & ! kernel subdt , & ! subcycling time step for stability (s) elapsed_t ! elapsed subcycling time @@ -942,7 +942,6 @@ subroutine fsd_weld_thermo (ncat, nfsd, & afsdn (:,:) = c0 afsd_init(:) = c0 stability = c0 - prefac = p5 do n = 1, ncat @@ -986,8 +985,7 @@ subroutine fsd_weld_thermo (ncat, nfsd, & if (k > i) then kern = c_weld * floe_area_c(i) * aicen(n) loss(i) = loss(i) + kern*afsd_tmp(i)*afsd_tmp(j) - if (i.eq.j) prefac = c1 ! otherwise 0.5 - gain(k) = gain(k) + prefac*kern*afsd_tmp(i)*afsd_tmp(j) + gain(k) = gain(k) + kern*afsd_tmp(i)*afsd_tmp(j) end if end do end do @@ -1011,11 +1009,11 @@ subroutine fsd_weld_thermo (ncat, nfsd, & END DO ! time + afsdn(:,n) = afsd_tmp(:) call icepack_cleanup_fsdn (nfsd, afsdn(:,n)) if (icepack_warnings_aborted(subname)) return do k = 1, nfsd - afsdn(k,n) = afsd_tmp(k) trcrn(nt_fsd+k-1,n) = afsdn(k,n) ! history/diagnostics d_afsdn_weld(k,n) = afsdn(k,n) - afsd_init(k) diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index 93ed8bebd..3aac1ac70 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -884,7 +884,7 @@ subroutine lateral_melt (dt, ncat, & fhocn, faero_ocn, & fiso_ocn, & rside, meltl, & - fside, sss, & + fside, wlat, & aicen, vicen, & vsnon, trcrn, & fzsal, flux_bio, & @@ -914,6 +914,9 @@ subroutine lateral_melt (dt, ncat, & real (kind=dbl_kind), intent(in) :: & rside , & ! fraction of ice that melts laterally + wlat ! lateral melt rate (m/s) + + real (kind=dbl_kind), intent(inout) :: & fside ! lateral heat flux (W/m^2) real (kind=dbl_kind), intent(inout) :: & @@ -955,7 +958,7 @@ subroutine lateral_melt (dt, ncat, & dfsalt , & ! change in fsalt dvssl , & ! snow surface layer volume dvint , & ! snow interior layer - cat1_arealoss, tmp ! + bin1_arealoss, tmp ! logical (kind=log_kind) :: & flag ! .true. if there could be lateral melting @@ -965,7 +968,6 @@ subroutine lateral_melt (dt, ncat, & vicen_init, & ! volume per unit area of ice (m) G_radialn , & ! rate of lateral melt (m/s) delta_an , & ! change in the ITD - qin , & ! enthalpy rsiden ! delta_an/aicen real (kind=dbl_kind), dimension (nfsd,ncat) :: & @@ -979,11 +981,9 @@ subroutine lateral_melt (dt, ncat, & real (kind=dbl_kind), dimension(nfsd+1) :: & f_flx ! -!echmod - for average qin - real (kind=dbl_kind), intent(in) :: & - sss real (kind=dbl_kind) :: & - Ti, Si0, qi0, sicen, & + sicen, & + etot, & ! column energy per itd cat, for FSD code elapsed_t, & ! FSD subcycling subdt ! FSD timestep (s) @@ -996,12 +996,11 @@ subroutine lateral_melt (dt, ncat, & dfsalt = c0 dvssl = c0 dvint = c0 - cat1_arealoss = c0 + bin1_arealoss = c0 tmp = c0 vicen_init = c0 G_radialn = c0 delta_an = c0 - qin = c0 rsiden = c0 afsdn = c1 afsdn_init = c0 @@ -1018,59 +1017,57 @@ subroutine lateral_melt (dt, ncat, & d_afsd_latm(:) = c0 end if - if (tr_fsd .and. fside < c0) then + if (tr_fsd .and. wlat > puny) then flag = .true. - - -! echmod - using category values would be preferable to the average value - ! Compute average enthalpy of ice (taken from add_new_ice) - if (sss > c2 * dSin0_frazil) then - Si0 = sss - dSin0_frazil - else - Si0 = sss**2 / (c4*dSin0_frazil) - endif - Ti = min(liquidus_temperature_mush(Si0/phi_init), -p1) - qi0 = enthalpy_mush(Ti, Si0) - + ! for FSD rside and fside not yet computed correctly, redo here + fside = c0 do n = 1, ncat - if (ktherm == 2) then ! mushy - do k = 1, nilyr - !qin(n) = qin(n) & - ! + trcrn(nt_qice+k-1,n)*vicen(n)/real(nilyr,kind=dbl_kind) - qin(n) = qi0 - enddo - else - qin(n) = -rhoi*Lfresh - endif - if (qin(n) < -puny) G_radialn(n) = -fside/qin(n) ! negative + G_radialn(n) = -wlat ! negative - if (G_radialn(n) < -puny) then + if (any(afsdn(:,n) < c0)) then + write(warnstr,*) subname, 'lateral_melt B afsd < 0 ',n + call icepack_warnings_add(warnstr) + endif + bin1_arealoss = -trcrn(nt_fsd+1-1,n) * aicen(n) * dt & + * G_radialn(n) / floe_binwidth(1) - if (any(afsdn(:,n) < c0)) print*,& - 'lateral_melt B afsd < 0',n + delta_an(n) = c0 + do k = 1, nfsd + delta_an(n) = delta_an(n) + ((c2/floe_rad_c(k))*aicen(n) & + * trcrn(nt_fsd+k-1,n)*G_radialn(n)*dt) ! delta_an < 0 + end do - cat1_arealoss = -trcrn(nt_fsd+1-1,n) * aicen(n) * dt & - * G_radialn(n) / floe_binwidth(1) + ! add negative area loss from fsd + delta_an(n) = delta_an(n) - bin1_arealoss - delta_an(n) = c0 - do k = 1, nfsd - delta_an(n) = delta_an(n) + ((c2/floe_rad_c(k))*aicen(n) & - * trcrn(nt_fsd+k-1,n)*G_radialn(n)*dt) ! delta_an < 0 - end do + if (delta_an(n) > c0) then + write(warnstr,*) subname, 'ERROR delta_an > 0 ',delta_an(n) + call icepack_warnings_add(warnstr) + endif - ! add negative area loss from fsd - delta_an(n) = delta_an(n) - cat1_arealoss + ! following original code, not necessary for fsd + if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n)/aicen(n),c1) - if (delta_an(n) > c0) print*,'ERROR delta_an > 0', delta_an(n) + if (rsiden(n) < c0) then + write(warnstr,*) subname, 'ERROR rsiden < 0 ',rsiden(n) + call icepack_warnings_add(warnstr) + endif - ! following original code, not necessary for fsd - if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n)/aicen(n),c1) + ! melting energy/unit area in each column, etot < 0 + etot = c0 + do k = 1, nslyr + etot = etot + trcrn(nt_qsno+k-1,n) * vsnon(n)/real(nslyr,kind=dbl_kind) + enddo - if (rsiden(n) < c0) print*,'ERROR rsiden < 0', rsiden(n) + do k = 1, nilyr + etot = etot + trcrn(nt_qice+k-1,n) * vicen(n)/real(nilyr,kind=dbl_kind) + enddo ! nilyr + + ! lateral heat flux, fside < 0 + fside = fside + rsiden(n)*etot/dt - end if ! G_radialn enddo ! ncat else if (rside > c0) then ! original, non-fsd implementation @@ -1132,8 +1129,10 @@ subroutine lateral_melt (dt, ncat, & DO WHILE (elapsed_t.lt.dt) nsubt = nsubt + 1 - if (nsubt.gt.100) & - print *, 'latm not converging' + if (nsubt.gt.100) then + write(warnstr,*) subname, 'latm not converging' + call icepack_warnings_add(warnstr) + endif ! finite differences df_flx(:) = c0 @@ -1146,8 +1145,10 @@ subroutine lateral_melt (dt, ncat, & df_flx(k) = f_flx(k+1) - f_flx(k) end do - if (abs(sum(df_flx(:))) > puny) & - print*,'sum(df_flx)/=0' + if (abs(sum(df_flx(:))) > puny) then + write(warnstr,*) subname, 'sum(df_flx) /= 0' + call icepack_warnings_add(warnstr) + endif ! this term ensures area conservation tmp = SUM(afsd_tmp(:)/floe_rad_c(:)) @@ -1243,6 +1244,8 @@ subroutine lateral_melt (dt, ncat, & if (tr_fsd) then + trcrn(nt_fsd:nt_fsd+nfsd-1,:) = afsdn + call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) ) if (icepack_warnings_aborted(subname)) return @@ -1970,7 +1973,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & Tf, sss, & salinz, & rside, meltl, & - fside, & + fside, wlat, & frzmlt, frazil, & frain, fpond, & fresh, fsalt, & @@ -2010,10 +2013,12 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & Tf , & ! freezing temperature (C) sss , & ! sea surface salinity (ppt) rside , & ! fraction of ice that melts laterally - fside , & ! lateral heat flux (W/m^2) frzmlt , & ! freezing/melting potential (W/m^2) wave_sig_ht ! significant height of waves in ice (m) + real (kind=dbl_kind), intent(in), optional :: & + wlat ! lateral melt rate (m/s) + real (kind=dbl_kind), dimension(:), intent(in) :: & wave_spectrum ! ocean surface wave spectrum E(f) (m^2 s) @@ -2052,6 +2057,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & real (kind=dbl_kind), intent(inout) :: & aice , & ! sea ice concentration aice0 , & ! concentration of open water + fside , & ! lateral heat flux (W/m^2) frain , & ! rainfall rate (kg/m^2 s) fpond , & ! fresh water flux to ponds (kg/m^2/s) fresh , & ! fresh water flux to ocean (kg/m^2/s) @@ -2122,6 +2128,13 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & return endif endif + if (tr_fsd) then + if (.not.(present(wlat))) then + call icepack_warnings_add(subname//' error in FSD arguments, tr_fsd=T') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif + endif endif !----------------------------------------------------------------- @@ -2221,7 +2234,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & fhocn, faero_ocn, & fiso_ocn, & rside, meltl, & - fside, sss, & + fside, wlat, & aicen, vicen, & vsnon, trcrn, & fzsal, flux_bio, & diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index ac7cf232a..32a9e8f4f 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -493,7 +493,7 @@ subroutine frzmlt_bottom_lateral (dt, ncat, & strocnxT, strocnyT, & Tbot, fbot, & rside, Cdn_ocn, & - fside) + fside, wlat) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories @@ -527,6 +527,7 @@ subroutine frzmlt_bottom_lateral (dt, ncat, & real (kind=dbl_kind), intent(out) :: & Tbot , & ! ice bottom surface temperature (deg C) fbot , & ! heat flux to ice bottom (W/m^2) + wlat , & ! lateral melt rate (m/s) rside , & ! fraction of ice that melts laterally fside ! lateral heat flux (W/m^2) @@ -543,7 +544,6 @@ subroutine frzmlt_bottom_lateral (dt, ncat, & real (kind=dbl_kind) :: & deltaT , & ! SST - Tbot >= 0 ustar , & ! skin friction velocity for fbot (m/s) - wlat , & ! lateral melt rate (m/s) xtmp ! temporary variable ! Parameters for bottom melting @@ -616,31 +616,24 @@ subroutine frzmlt_bottom_lateral (dt, ncat, & do n = 1, ncat etot = c0 - qavg = c0 - ! melting energy/unit area in each column, etot < 0 do k = 1, nslyr etot = etot + qsnon(k,n) * vsnon(n)/real(nslyr,kind=dbl_kind) - qavg = qavg + qsnon(k,n) enddo do k = 1, nilyr etot = etot + qicen(k,n) * vicen(n)/real(nilyr,kind=dbl_kind) - qavg = qavg + qicen(k,n) enddo ! nilyr ! lateral heat flux, fside < 0 - if (tr_fsd) then ! floe size distribution - fside = fside + wlat*qavg - else ! default floe size - fside = fside + rside*etot/dt - endif + fside = fside + rside*etot/dt enddo ! n !----------------------------------------------------------------- ! Limit bottom and lateral heat fluxes if necessary. + ! FYI: fside is not yet correct for fsd, may need to adjust fbot further !----------------------------------------------------------------- xtmp = frzmlt/(fbot + fside + puny) @@ -2119,7 +2112,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fbot , & Tbot , Tsnice , & frzmlt , rside , & - fside , & + fside , wlat , & fsnow , frain , & fpond , fsloss , & fsurf , fsurfn , & @@ -2260,6 +2253,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & mlt_onset , & ! day of year that sfc melting begins frz_onset ! day of year that freezing begins (congel or frazil) + real (kind=dbl_kind), intent(out), optional :: & + wlat ! lateral melt rate (m/s) + real (kind=dbl_kind), intent(inout), optional :: & fswthru_vdr , & ! vis dir shortwave penetrating to ocean (W/m^2) fswthru_vdf , & ! vis dif shortwave penetrating to ocean (W/m^2) @@ -2449,6 +2445,13 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & call icepack_warnings_setabort(.true.,__FILE__,__LINE__) return endif + if (tr_fsd) then + if (.not.(present(wlat))) then + call icepack_warnings_add(subname//' error in FSD arguments, tr_fsd=T') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif + endif endif !----------------------------------------------------------------- @@ -2543,7 +2546,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & strocnxT, strocnyT, & Tbot, fbot, & rside, Cdn_ocn, & - fside) + fside, wlat) if (icepack_warnings_aborted(subname)) return diff --git a/configuration/driver/icedrv_flux.F90 b/configuration/driver/icedrv_flux.F90 index 7e0e992ee..749e52396 100644 --- a/configuration/driver/icedrv_flux.F90 +++ b/configuration/driver/icedrv_flux.F90 @@ -276,6 +276,7 @@ module icedrv_flux real (kind=dbl_kind), dimension (nx), public :: & rside , & ! fraction of ice that melts laterally fside , & ! lateral heat flux (W/m^2) + wlat , & ! lateral melt rate (m/s) fsw , & ! incoming shortwave radiation (W/m^2) coszen , & ! cosine solar zenith angle, < 0 for sun below horizon rdg_conv, & ! convergence term for ridging (1/s) diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90 index bc6279d22..0b5e74dda 100644 --- a/configuration/driver/icedrv_step.F90 +++ b/configuration/driver/icedrv_step.F90 @@ -111,7 +111,7 @@ subroutine step_therm1 (dt) use icedrv_arrays_column, only: meltsliqn, meltsliq use icedrv_calendar, only: yday use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, n_iso, nx - use icedrv_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rside, fside, & + use icedrv_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rside, fside, wlat, & fbot, Tbot, Tsnice use icedrv_flux, only: meltsn, melttn, meltbn, congeln, snoicen, uatm, vatm use icedrv_flux, only: wind, rhoa, potT, Qa, Qa_iso, zlvl, strax, stray, flatn @@ -324,6 +324,7 @@ subroutine step_therm1 (dt) fbot = fbot(i), frzmlt = frzmlt(i), & Tbot = Tbot(i), Tsnice = Tsnice(i), & rside = rside(i), fside = fside(i), & + wlat = wlat(i), & fsnow = fsnow(i), frain = frain(i), & fpond = fpond(i), fsloss = fsloss(i), & fsurf = fsurf(i), fsurfn = fsurfn(i,:), & @@ -433,7 +434,7 @@ subroutine step_therm2 (dt) use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, nblyr, & nltrcr, nx, nfsd use icedrv_flux, only: fresh, frain, fpond, frzmlt, frazil, frz_onset - use icedrv_flux, only: update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside, fside + use icedrv_flux, only: update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside, fside, wlat use icedrv_flux, only: meltl, frazil_diag, flux_bio, faero_ocn, fiso_ocn use icedrv_flux, only: HDO_ocn, H2_16O_ocn, H2_18O_ocn use icedrv_init, only: tmask @@ -496,6 +497,7 @@ subroutine step_therm2 (dt) nt_strata=nt_strata(1:ntrcr,:), & Tf=Tf(i), sss=sss(i), & salinz=salinz(i,:), fside=fside(i), & + wlat=wlat(i), & rside=rside(i), meltl=meltl(i), & frzmlt=frzmlt(i), frazil=frazil(i), & frain=frain(i), fpond=fpond(i), &