From 840a09be9316720221c2e458591988dda4779b06 Mon Sep 17 00:00:00 2001 From: bzhao Date: Tue, 25 Oct 2022 17:33:38 -0400 Subject: [PATCH 01/30] added ifdef to protect additional arguments --- columnphysics/icepack_therm_bl99.F90 | 18 ++++++++++++++++++ columnphysics/icepack_therm_vertical.F90 | 17 +++++++++++++++++ 2 files changed, 35 insertions(+) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 1c3a3aabb..5f8d4e61f 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -69,6 +69,9 @@ subroutine temperature_changes (dt, & Tsf, Tbot, & fsensn, flatn, & flwoutn, fsurfn, & +#ifdef GEOSCOUPLED + dfsurfdt_in, & +#endif fcondtopn,fcondbot, & einit ) @@ -114,6 +117,11 @@ subroutine temperature_changes (dt, & flatn , & ! surface downward latent heat (W m-2) flwoutn ! upward LW at surface (W m-2) +#ifdef GEOSCOUPLED + real (kind=dbl_kind), intent(in):: & + dfsurfdt_in +#endif + real (kind=dbl_kind), intent(out):: & fcondbot ! downward cond flux at bottom surface (W m-2) @@ -226,6 +234,10 @@ subroutine temperature_changes (dt, & dflat_dT = c0 dflwout_dT = c0 einex = c0 +#ifdef GEOSCOUPLED + ! derivative information is passed by GEOS + dfsurf_dT = dfsurfdt_in +#endif dt_rhoi_hlyr = dt / (rhoi*hilyr) ! hilyr > 0 if (hslyr > hs_min/real(nslyr,kind=dbl_kind)) & l_snow = .true. @@ -340,7 +352,10 @@ subroutine temperature_changes (dt, & !----------------------------------------------------------------- converged = .true. +#ifndef GEOSCOUPLED + ! derivative information is passed by GEOS dfsurf_dT = c0 +#endif avg_Tsi = c0 enew = c0 einex = c0 @@ -366,6 +381,7 @@ subroutine temperature_changes (dt, & if (calc_Tsfc) then +#ifndef GEOSCOUPLED !----------------------------------------------------------------- ! Update radiative and turbulent fluxes and their derivatives ! with respect to Tsf. @@ -387,6 +403,8 @@ subroutine temperature_changes (dt, & dfsens_dT, dflat_dT ) if (icepack_warnings_aborted(subname)) return +#endif + !----------------------------------------------------------------- ! Compute conductive flux at top surface, fcondtopn. ! If fsurfn < fcondtopn and Tsf = 0, then reset Tsf to slightly less diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 67e4b4f16..4f3f8fa76 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -93,6 +93,9 @@ subroutine thermo_vertical (nilyr, nslyr, & fswsfc, fswint, & Sswabs, Iswabs, & fsurfn, fcondtopn, & +#ifdef GEOSCOUPLED + dfsurfdt_in, & +#endif fcondbotn, & fsensn, flatn, & flwoutn, evapn, & @@ -184,6 +187,11 @@ subroutine thermo_vertical (nilyr, nslyr, & fcondtopn, & ! downward cond flux at top surface (W m-2) fcondbotn ! downward cond flux at bottom surface (W m-2) +#ifdef GEOSCOUPLED + real (kind=dbl_kind), intent(inout):: & + dfsurfdt_in +#endif + ! coupler fluxes to ocean real (kind=dbl_kind), intent(out):: & freshn , & ! fresh water flux to ocean (kg/m^2/s) @@ -346,6 +354,9 @@ subroutine thermo_vertical (nilyr, nslyr, & Tsf, Tbot, & fsensn, flatn, & flwoutn, fsurfn, & +#ifdef GEOSCOUPLED + dfsurfdt_in, & +#endif fcondtopn, fcondbotn, & einit ) if (icepack_warnings_aborted(subname)) return @@ -2163,6 +2174,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fsnow , frain , & fpond , fsloss , & fsurf , fsurfn , & +#ifdef GEOSCOUPLED + dfsurfdt_in , & +#endif fcondtop , fcondtopn , & fcondbot , fcondbotn , & fswsfcn , fswintn , & @@ -2786,6 +2800,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fswsfc=fswsfcn (n), fswint=fswintn (n), & Sswabs=Sswabsn(:,n), Iswabs=Iswabsn(:,n), & fsurfn=fsurfn (n), fcondtopn=fcondtopn(n), & +#ifdef GEOSCOUPLED + dfsurfdt_in = dfsurfdt(n), & +#endif fcondbotn=fcondbotn(n), & fsensn=fsensn (n), flatn=flatn (n), & flwoutn=flwoutn, evapn=evapn, & From 0fea3b783e7ce869fd0c0ea1b5a8d6b2fa2ce246 Mon Sep 17 00:00:00 2001 From: bzhao Date: Mon, 7 Nov 2022 15:17:11 -0500 Subject: [PATCH 02/30] add extra arguments to various thermo routines; skip surface flux step and use passed in values to update Ts --- columnphysics/icepack_therm_bl99.F90 | 12 +++++++++--- columnphysics/icepack_therm_vertical.F90 | 18 +++++++++++++++--- 2 files changed, 24 insertions(+), 6 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 5f8d4e61f..3c82f5346 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -71,6 +71,8 @@ subroutine temperature_changes (dt, & flwoutn, fsurfn, & #ifdef GEOSCOUPLED dfsurfdt_in, & + flatn_f, & + dflatdt_in, & #endif fcondtopn,fcondbot, & einit ) @@ -119,7 +121,9 @@ subroutine temperature_changes (dt, & #ifdef GEOSCOUPLED real (kind=dbl_kind), intent(in):: & - dfsurfdt_in + dfsurfdt_in , & + flatn_f , & + dflatdt_in #endif real (kind=dbl_kind), intent(out):: & @@ -237,6 +241,8 @@ subroutine temperature_changes (dt, & #ifdef GEOSCOUPLED ! derivative information is passed by GEOS dfsurf_dT = dfsurfdt_in + dflat_dT = dflatdt_in + flatn = flatn_f #endif dt_rhoi_hlyr = dt / (rhoi*hilyr) ! hilyr > 0 if (hslyr > hs_min/real(nslyr,kind=dbl_kind)) & @@ -402,7 +408,6 @@ subroutine temperature_changes (dt, & dfsurf_dT, dflwout_dT, & dfsens_dT, dflat_dT ) if (icepack_warnings_aborted(subname)) return - #endif !----------------------------------------------------------------- @@ -811,10 +816,11 @@ subroutine temperature_changes (dt, & endif if (calc_Tsfc) then - +#ifndef GEOSCOUPLED ! update fluxes that depend on Tsf flwoutn = flwoutn + dTsf_prev * dflwout_dT fsensn = fsensn + dTsf_prev * dfsens_dT +#endif flatn = flatn + dTsf_prev * dflat_dT endif ! calc_Tsfc diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 4f3f8fa76..9bd05ef14 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -95,6 +95,8 @@ subroutine thermo_vertical (nilyr, nslyr, & fsurfn, fcondtopn, & #ifdef GEOSCOUPLED dfsurfdt_in, & + flatn_f, & + dflatdt_in, & #endif fcondbotn, & fsensn, flatn, & @@ -188,8 +190,10 @@ subroutine thermo_vertical (nilyr, nslyr, & fcondbotn ! downward cond flux at bottom surface (W m-2) #ifdef GEOSCOUPLED - real (kind=dbl_kind), intent(inout):: & - dfsurfdt_in + real (kind=dbl_kind), intent(in):: & + dfsurfdt_in, & + flatn_f, & + dflatdt_in #endif ! coupler fluxes to ocean @@ -356,6 +360,8 @@ subroutine thermo_vertical (nilyr, nslyr, & flwoutn, fsurfn, & #ifdef GEOSCOUPLED dfsurfdt_in, & + flatn_f, & + dflatdt_in, & #endif fcondtopn, fcondbotn, & einit ) @@ -2175,7 +2181,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fpond , fsloss , & fsurf , fsurfn , & #ifdef GEOSCOUPLED - dfsurfdt_in , & + dfsurfdt , dflatdt , & #endif fcondtop , fcondtopn , & fcondbot , fcondbotn , & @@ -2360,6 +2366,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fsurfn , & ! net flux to top surface, excluding fcondtop fcondtopn , & ! downward cond flux at top surface (W m-2) fcondbotn , & ! downward cond flux at bottom surface (W m-2) +#ifdef GEOSCOUPLED + dfsurfdt , & ! + dflatdt , & ! +#endif flatn , & ! latent heat flux (W m-2) fsensn , & ! sensible heat flux (W m-2) fsurfn_f , & ! net flux to top surface, excluding fcondtop @@ -2802,6 +2812,8 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fsurfn=fsurfn (n), fcondtopn=fcondtopn(n), & #ifdef GEOSCOUPLED dfsurfdt_in = dfsurfdt(n), & + flatn_f = flatn_f(n), & + dflatdt_in = dflatdt(n), & #endif fcondbotn=fcondbotn(n), & fsensn=fsensn (n), flatn=flatn (n), & From 98647f4a1f78dcc976f0b830eba5f7e4e63fc542 Mon Sep 17 00:00:00 2001 From: bzhao Date: Mon, 7 Nov 2022 19:23:02 -0500 Subject: [PATCH 03/30] mods for recomputing sublimation/condensation --- columnphysics/icepack_therm_vertical.F90 | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 9bd05ef14..8c30550b2 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -1281,7 +1281,11 @@ subroutine thickness_changes (nilyr, nslyr, & evapin = c0 ! initialize if (hsn > puny) then ! add snow with enthalpy zqsn(1) +#ifdef GEOSCOUPLED + dhs = econ / (-rhos*Lfresh - rhos*Lvap) ! econ < 0, dhs > 0 +#else dhs = econ / (zqsn(1) - rhos*Lvap) ! econ < 0, dhs > 0 +#endif mass = massice(1) + massliq(1) massi = c0 @@ -1293,7 +1297,11 @@ subroutine thickness_changes (nilyr, nslyr, & evapn = evapn + dhs*rhos evapsn = evapsn + dhs*rhos else ! add ice with enthalpy zqin(1) +#ifdef GEOSCOUPLED + dhi = econ / (-rhoi*Lfresh - rhoi*Lvap) ! econ < 0, dhi > 0 +#else dhi = econ / (qm(1) - rhoi*Lvap) ! econ < 0, dhi > 0 +#endif dzi(1) = dzi(1) + dhi evapn = evapn + dhi*rhoi evapin = evapin + dhi*rhoi @@ -1393,8 +1401,11 @@ subroutine thickness_changes (nilyr, nslyr, & !-------------------------------------------------------------- ! Sublimation of snow (evapn < 0) !-------------------------------------------------------------- - +#ifdef GEOSCOUPLED + qsub = -rhos*Lfresh - rhos*Lvap ! qsub < 0 +#else qsub = zqsn(k) - rhos*Lvap ! qsub < 0 +#endif dhs = max (-dzs(k), esub/qsub) ! esub > 0, dhs < 0 mass = massice(1) + massliq(1) @@ -1438,7 +1449,11 @@ subroutine thickness_changes (nilyr, nslyr, & ! Sublimation of ice (evapn < 0) !-------------------------------------------------------------- +#ifdef GEOSCOUPLED + qsub = -rhoi*Lfresh - rhoi*Lvap ! qsub < 0 +#else qsub = qm(k) - rhoi*Lvap ! qsub < 0 +#endif dhi = max (-dzi(k), esub/qsub) ! esub < 0, dhi < 0 dzi(k) = dzi(k) + dhi esub = esub - dhi*qsub From efc5483facb1fb69ff0f0ed38d194be705faf9af Mon Sep 17 00:00:00 2001 From: bzhao Date: Tue, 20 Dec 2022 15:12:35 -0500 Subject: [PATCH 04/30] properly account for extra energy due to inconsistent cond/subl treatment; add solar to fsurf --- columnphysics/icepack_therm_bl99.F90 | 12 +++++++ columnphysics/icepack_therm_vertical.F90 | 46 +++++++++++++++++++++++- 2 files changed, 57 insertions(+), 1 deletion(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 3c82f5346..243698c71 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -339,6 +339,10 @@ subroutine temperature_changes (dt, & endif +#ifdef GEOSCOUPLED + fsurfn = fsurfn + fswsfc ! this is the total heat flux +#endif + !----------------------------------------------------------------- ! Solve for new temperatures. ! Iterate until temperatures converge with minimal energy error. @@ -751,6 +755,14 @@ subroutine temperature_changes (dt, & call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'fsurf:', fsurfn call icepack_warnings_add(warnstr) + write(warnstr,*) subname, 'dfsurf_dT:', dfsurf_dT + call icepack_warnings_add(warnstr) + write(warnstr,*) subname, 'enew:', enew + call icepack_warnings_add(warnstr) + write(warnstr,*) subname, 'einit:', einit + call icepack_warnings_add(warnstr) + write(warnstr,*) subname, 'dt:', dt + call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'fcondtop, fcondbot, fswint', & fcondtopn, fcondbot, fswint call icepack_warnings_add(warnstr) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 8c30550b2..4796374fa 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -252,6 +252,7 @@ subroutine thermo_vertical (nilyr, nslyr, & efinal , & ! final energy of melting (J m-2) einter ! intermediate energy + real (kind=dbl_kind) :: & fadvocn ! advective heat flux to ocean @@ -1180,6 +1181,11 @@ subroutine thickness_changes (nilyr, nslyr, & qm , & ! energy of melting (J m-3) = zqin in BL99 formulation qmlt ! enthalpy of melted ice (J m-3) = zero in BL99 formulation +#ifdef GEOSCOUPLED + real (kind=dbl_kind) :: & + sblx ! flux error due to cond/sub inconsistency +#endif + real (kind=dbl_kind) :: & qbotm , & qbotp , & @@ -1198,6 +1204,10 @@ subroutine thickness_changes (nilyr, nslyr, & dhs = c0 hsn_new = c0 +#ifdef GEOSCOUPLED + sblx = c0 +#endif + do k = 1, nilyr dzi(k) = hilyr enddo @@ -1292,7 +1302,17 @@ subroutine thickness_changes (nilyr, nslyr, & if (dzs(1) > puny) massi = c1 + dhs/dzs(1) massice(1) = massice(1) * massi massliq(1) = max(c0, mass + rhos*dhs - massice(1)) ! conserve new total mass - +#ifdef GEOSCOUPLED + hstot = dzs(1) + dhs + ! adjust top layer snow enthalpy b.c. we added them at 0C + zqsnew = -rhos*Lfresh + if (hstot > puny) then + zqsn(1) = (dzs(1) * zqsn(1) & + + dhs * zqsnew) / hstot + ! avoid roundoff errors + zqsn(1) = min(zqsn(1), -rhos*Lfresh) + endif +#endif dzs(1) = dzs(1) + dhs evapn = evapn + dhs*rhos evapsn = evapsn + dhs*rhos @@ -1302,7 +1322,18 @@ subroutine thickness_changes (nilyr, nslyr, & #else dhi = econ / (qm(1) - rhoi*Lvap) ! econ < 0, dhi > 0 #endif + +#ifdef GEOSCOUPLED + ! adjust top layer ice enthalpy b.c. we added them at 0C + zqsnew = -rhoi*Lfresh + hqtot = dzi(1)*qm(1) + dhi*zqsnew +#endif + dzi(1) = dzi(1) + dhi +#ifdef GEOSCOUPLED + if (dzi(1) > puny) & + zqin(1) = hqtot / dzi(1) ! need to revisit for kterm = 2 +#endif evapn = evapn + dhi*rhoi evapin = evapin + dhi*rhoi ! enthalpy of melt water @@ -1408,6 +1439,10 @@ subroutine thickness_changes (nilyr, nslyr, & #endif dhs = max (-dzs(k), esub/qsub) ! esub > 0, dhs < 0 +#ifdef GEOSCOUPLED + sblx = sblx + (-dhs)*min(zqsn(k) - (-rhos*Lfresh), c0) ! sblx < 0 (J m-2) +#endif + mass = massice(1) + massliq(1) massi = c0 if (dzs(k) > puny) massi = c1 + dhs/dzs(k) @@ -1455,6 +1490,11 @@ subroutine thickness_changes (nilyr, nslyr, & qsub = qm(k) - rhoi*Lvap ! qsub < 0 #endif dhi = max (-dzi(k), esub/qsub) ! esub < 0, dhi < 0 + +#ifdef GEOSCOUPLED + sblx = sblx + (-dhi)*(qm(k) - (-rhoi*Lfresh)) ! sblx can be v+- (J m-2) +#endif + dzi(k) = dzi(k) + dhi esub = esub - dhi*qsub esub = max(esub, c0) @@ -1802,7 +1842,11 @@ subroutine thickness_changes (nilyr, nslyr, & ! sublimated/condensed ice. !----------------------------------------------------------------- +#ifdef GEOSCOUPLED + efinal = -evapn*Lvap + sblx +#else efinal = -evapn*Lvap +#endif evapn = evapn/dt evapsn = evapsn/dt evapin = evapin/dt From ace2808f37f1de8eeff3c8f537fa2e176193106f Mon Sep 17 00:00:00 2001 From: bzhao Date: Fri, 13 Jan 2023 15:06:12 -0500 Subject: [PATCH 05/30] move flatn update into convergence loop --- columnphysics/icepack_therm_bl99.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 243698c71..ac15d4a99 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -689,6 +689,7 @@ subroutine temperature_changes (dt, & !----------------------------------------------------------------- fsurfn = fsurfn + dTsf*dfsurf_dT + flatn = flatn + dTsf*dflat_dT if (l_snow) then fcondtopn = kh(1) * (Tsf-zTsn(1)) else @@ -832,8 +833,8 @@ subroutine temperature_changes (dt, & ! update fluxes that depend on Tsf flwoutn = flwoutn + dTsf_prev * dflwout_dT fsensn = fsensn + dTsf_prev * dfsens_dT -#endif flatn = flatn + dTsf_prev * dflat_dT +#endif endif ! calc_Tsfc From 72e6f7ac6d10c9dac2c6c5b1b532016037527503 Mon Sep 17 00:00:00 2001 From: bzhao Date: Fri, 13 Jan 2023 15:07:51 -0500 Subject: [PATCH 06/30] change enthalpy computation for condensed ice; use sblx to track differenc from fresh ice --- columnphysics/icepack_therm_vertical.F90 | 48 ++++++++++++++++++------ 1 file changed, 36 insertions(+), 12 deletions(-) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 4796374fa..9673ce825 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -111,7 +111,8 @@ subroutine thermo_vertical (nilyr, nslyr, & congel, snoice, & mlt_onset, frz_onset, & yday, dsnow, & - prescribed_ice) + prescribed_ice, & + my_tsk, my_i, my_j, my_blk) integer (kind=int_kind), intent(in) :: & nilyr , & ! number of ice layers @@ -137,6 +138,9 @@ subroutine thermo_vertical (nilyr, nslyr, & logical (kind=log_kind), intent(in), optional :: & prescribed_ice ! if .true., use prescribed ice instead of computed + integer (kind=int_kind), intent(in), optional :: & + my_tsk, my_i, my_j, my_blk + real (kind=dbl_kind), dimension (:), intent(inout) :: & zqsn , & ! snow layer enthalpy, zqsn < 0 (J m-3) zqin , & ! ice layer enthalpy, zqin < 0 (J m-3) @@ -460,7 +464,11 @@ subroutine thermo_vertical (nilyr, nslyr, & fsnow, einit, & einter, efinal, & fcondtopn, fcondbotn, & - fadvocn, fbot ) + fadvocn, fbot, & + my_tsk = my_tsk, & + my_i = my_i, & + my_j = my_j, & + my_blk = my_blk) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- @@ -1325,15 +1333,15 @@ subroutine thickness_changes (nilyr, nslyr, & #ifdef GEOSCOUPLED ! adjust top layer ice enthalpy b.c. we added them at 0C - zqsnew = -rhoi*Lfresh - hqtot = dzi(1)*qm(1) + dhi*zqsnew + !zqsnew = -rhoi*Lfresh + !hqtot = dzi(1)*qm(1) + dhi*zqsnew + sblx = sblx + (-dhi)*(qm(1) - (-rhoi*Lfresh)) ! sblx can be v+- (J m-2) #endif - dzi(1) = dzi(1) + dhi -#ifdef GEOSCOUPLED - if (dzi(1) > puny) & - zqin(1) = hqtot / dzi(1) ! need to revisit for kterm = 2 -#endif +!#ifdef GEOSCOUPLED + ! if (dzi(1) > puny) & + ! zqin(1) = hqtot / dzi(1) ! need to revisit for kterm = 2 +!#endif evapn = evapn + dhi*rhoi evapin = evapin + dhi*rhoi ! enthalpy of melt water @@ -2010,7 +2018,9 @@ subroutine conservation_check_vthermo(dt, & einit, einter, & efinal, & fcondtopn,fcondbotn, & - fadvocn, fbot ) + fadvocn, fbot, & + my_tsk, my_i, & + my_j, my_blk) real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -2031,6 +2041,9 @@ subroutine conservation_check_vthermo(dt, & efinal , & ! final energy of melting (J m-2) fcondbotn + integer(kind=int_kind), intent(in), optional :: & + my_tsk, my_i, my_j, my_blk + ! local variables real (kind=dbl_kind) :: & @@ -2060,6 +2073,10 @@ subroutine conservation_check_vthermo(dt, & write(warnstr,*) subname, 'Thermo energy conservation error' call icepack_warnings_add(warnstr) + write(warnstr,*) subname, 'at task, i, j, iblk' + call icepack_warnings_add(warnstr) + write(warnstr,*) subname, my_tsk, my_i, my_j, my_blk + call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Flux error (W/m^2) =', ferr call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Energy error (J) =', ferr*dt @@ -2285,7 +2302,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & lmask_n , lmask_s , & mlt_onset , frz_onset , & yday , prescribed_ice, & - zlvs) + zlvs , my_tsk, & + my_i, my_j, & + my_blk ) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories @@ -2307,6 +2326,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & logical (kind=log_kind), intent(in), optional :: & prescribed_ice ! if .true., use prescribed ice instead of computed + integer (kind=int_kind), intent(in), optional :: & + my_tsk, my_i, my_j, my_blk + real (kind=dbl_kind), intent(inout) :: & aice , & ! sea ice concentration vice , & ! volume per unit area of ice (m) @@ -2887,7 +2909,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & congel=congeln (n), snoice=snoicen (n), & mlt_onset=mlt_onset, frz_onset=frz_onset, & yday=yday, dsnow=dsnown (n), & - prescribed_ice=prescribed_ice) + prescribed_ice=prescribed_ice, & + my_tsk = my_tsk, my_i = my_i, & + my_j = my_j, my_blk = my_blk ) if (icepack_warnings_aborted(subname)) then write(warnstr,*) subname, ' ice: Vertical thermo error, cat ', n From 7b3ddb12abc9dfece477ca7f92b1fba058f76078 Mon Sep 17 00:00:00 2001 From: bzhao Date: Tue, 17 Jan 2023 10:55:33 -0500 Subject: [PATCH 07/30] refactor addtional arguments passing and protect against potential errors --- columnphysics/icepack_therm_bl99.F90 | 27 +++++--- columnphysics/icepack_therm_vertical.F90 | 80 ++++++++++++++++-------- 2 files changed, 74 insertions(+), 33 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index ac15d4a99..6be4c8027 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -69,13 +69,11 @@ subroutine temperature_changes (dt, & Tsf, Tbot, & fsensn, flatn, & flwoutn, fsurfn, & -#ifdef GEOSCOUPLED + fcondtopn,fcondbot, & + einit, & dfsurfdt_in, & flatn_f, & - dflatdt_in, & -#endif - fcondtopn,fcondbot, & - einit ) + dflatdt_in) integer (kind=int_kind), intent(in) :: & nilyr , & ! number of ice layers @@ -119,12 +117,10 @@ subroutine temperature_changes (dt, & flatn , & ! surface downward latent heat (W m-2) flwoutn ! upward LW at surface (W m-2) -#ifdef GEOSCOUPLED - real (kind=dbl_kind), intent(in):: & + real (kind=dbl_kind), intent(in), optional:: & dfsurfdt_in , & flatn_f , & dflatdt_in -#endif real (kind=dbl_kind), intent(out):: & fcondbot ! downward cond flux at bottom surface (W m-2) @@ -239,6 +235,21 @@ subroutine temperature_changes (dt, & dflwout_dT = c0 einex = c0 #ifdef GEOSCOUPLED + if (.not. present(dfsurfdt_in)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing dfsurfdt_in" ) + return + endif + if (.not. present(dflatdt_in)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing dflatdt_in" ) + return + endif + if (.not. present(flatn_f)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing flatn_f" ) + return + endif ! derivative information is passed by GEOS dfsurf_dT = dfsurfdt_in dflat_dT = dflatdt_in diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 9673ce825..cda15db00 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -93,11 +93,6 @@ subroutine thermo_vertical (nilyr, nslyr, & fswsfc, fswint, & Sswabs, Iswabs, & fsurfn, fcondtopn, & -#ifdef GEOSCOUPLED - dfsurfdt_in, & - flatn_f, & - dflatdt_in, & -#endif fcondbotn, & fsensn, flatn, & flwoutn, evapn, & @@ -110,6 +105,9 @@ subroutine thermo_vertical (nilyr, nslyr, & smliq, massliq, & congel, snoice, & mlt_onset, frz_onset, & + dfsurfdt_in, & + flatn_f, & + dflatdt_in, & yday, dsnow, & prescribed_ice, & my_tsk, my_i, my_j, my_blk) @@ -193,12 +191,10 @@ subroutine thermo_vertical (nilyr, nslyr, & fcondtopn, & ! downward cond flux at top surface (W m-2) fcondbotn ! downward cond flux at bottom surface (W m-2) -#ifdef GEOSCOUPLED - real (kind=dbl_kind), intent(in):: & + real (kind=dbl_kind), intent(in), optional:: & dfsurfdt_in, & flatn_f, & dflatdt_in -#endif ! coupler fluxes to ocean real (kind=dbl_kind), intent(out):: & @@ -285,6 +281,24 @@ subroutine thermo_vertical (nilyr, nslyr, & zTin(:) = c0 meltsliq= c0 +#ifdef GEOSCOUPLED + if (.not. present(dfsurfdt_in)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing dfsurfdt_in" ) + return + endif + if (.not. present(dflatdt_in)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing dflatdt_in" ) + return + endif + if (.not. present(flatn_f)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing flatn_f" ) + return + endif +#endif + if (calc_Tsfc) then fsensn = c0 flatn = c0 @@ -363,13 +377,15 @@ subroutine thermo_vertical (nilyr, nslyr, & Tsf, Tbot, & fsensn, flatn, & flwoutn, fsurfn, & + fcondtopn, fcondbotn, & #ifdef GEOSCOUPLED - dfsurfdt_in, & - flatn_f, & - dflatdt_in, & -#endif - fcondtopn, fcondbotn, & + einit, & + dfsurfdt_in = dfsurfdt_in,& + flatn_f = flatn_f, & + dflatdt_in = dflatdt_in) +#else einit ) +#endif if (icepack_warnings_aborted(subname)) return endif ! ktherm @@ -2256,9 +2272,6 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fsnow , frain , & fpond , fsloss , & fsurf , fsurfn , & -#ifdef GEOSCOUPLED - dfsurfdt , dflatdt , & -#endif fcondtop , fcondtopn , & fcondbot , fcondbotn , & fswsfcn , fswintn , & @@ -2301,6 +2314,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & smicen , smliqn , & lmask_n , lmask_s , & mlt_onset , frz_onset , & + dfsurfdt , dflatdt , & yday , prescribed_ice, & zlvs , my_tsk, & my_i, my_j, & @@ -2447,10 +2461,6 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fsurfn , & ! net flux to top surface, excluding fcondtop fcondtopn , & ! downward cond flux at top surface (W m-2) fcondbotn , & ! downward cond flux at bottom surface (W m-2) -#ifdef GEOSCOUPLED - dfsurfdt , & ! - dflatdt , & ! -#endif flatn , & ! latent heat flux (W m-2) fsensn , & ! sensible heat flux (W m-2) fsurfn_f , & ! net flux to top surface, excluding fcondtop @@ -2471,6 +2481,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & snoicen , & ! snow-ice growth (m) dsnown ! change in snow thickness (m/step-->cm/day) + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & + dfsurfdt , & ! + dflatdt ! + real (kind=dbl_kind), optional, dimension(:), intent(inout) :: & fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2) fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2) @@ -2573,6 +2587,22 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & character(len=*),parameter :: subname='(icepack_step_therm1)' + + + +#ifdef GEOSCOUPLED + if (.not. present(dfsurfdt)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing dfsurfdt" ) + return + endif + if (.not. present(dflatdt)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing dflatdt" ) + return + endif +#endif + !----------------------------------------------------------------- ! allocate local optional arguments !----------------------------------------------------------------- @@ -2891,11 +2921,6 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fswsfc=fswsfcn (n), fswint=fswintn (n), & Sswabs=Sswabsn(:,n), Iswabs=Iswabsn(:,n), & fsurfn=fsurfn (n), fcondtopn=fcondtopn(n), & -#ifdef GEOSCOUPLED - dfsurfdt_in = dfsurfdt(n), & - flatn_f = flatn_f(n), & - dflatdt_in = dflatdt(n), & -#endif fcondbotn=fcondbotn(n), & fsensn=fsensn (n), flatn=flatn (n), & flwoutn=flwoutn, evapn=evapn, & @@ -2908,6 +2933,11 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & smliq=l_smliq (:,n), massliq=massliqn(:,n), & congel=congeln (n), snoice=snoicen (n), & mlt_onset=mlt_onset, frz_onset=frz_onset, & +#ifdef GEOSCOUPLED + dfsurfdt_in = dfsurfdt(n), & + flatn_f = flatn_f(n), & + dflatdt_in = dflatdt(n), & +#endif yday=yday, dsnow=dsnown (n), & prescribed_ice=prescribed_ice, & my_tsk = my_tsk, my_i = my_i, & From b9b1dade6acb2e4d69f0b82b2046543373c07c21 Mon Sep 17 00:00:00 2001 From: bzhao Date: Tue, 24 Jan 2023 20:42:36 -0500 Subject: [PATCH 08/30] add sblx to fhocn to conserve energy --- columnphysics/icepack_therm_vertical.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index cda15db00..d3f03daac 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -1868,6 +1868,7 @@ subroutine thickness_changes (nilyr, nslyr, & #ifdef GEOSCOUPLED efinal = -evapn*Lvap + sblx + fhocnn = fhocnn + sblx/dt #else efinal = -evapn*Lvap #endif From 05d53f69f41c5b0a0a50b81996b23ca893e2e22b Mon Sep 17 00:00:00 2001 From: bzhao Date: Thu, 26 Jan 2023 17:01:04 -0500 Subject: [PATCH 09/30] fixed a bug (sblx is added to fhocn instead of efinal) --- columnphysics/icepack_therm_vertical.F90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index d3f03daac..183e321d7 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -1867,11 +1867,9 @@ subroutine thickness_changes (nilyr, nslyr, & !----------------------------------------------------------------- #ifdef GEOSCOUPLED - efinal = -evapn*Lvap + sblx fhocnn = fhocnn + sblx/dt -#else - efinal = -evapn*Lvap #endif + efinal = -evapn*Lvap evapn = evapn/dt evapsn = evapsn/dt evapin = evapin/dt From 8b853582fe780eb5752599b0f1a077f3dd9c4f23 Mon Sep 17 00:00:00 2001 From: bzhao Date: Fri, 27 Jan 2023 16:55:10 -0500 Subject: [PATCH 10/30] allow addtional band of penetrative sw flux to be computed --- columnphysics/icepack_shortwave.F90 | 128 ++++++++++++++++++++++++++++ 1 file changed, 128 insertions(+) diff --git a/columnphysics/icepack_shortwave.F90 b/columnphysics/icepack_shortwave.F90 index d1d86621e..14829a158 100644 --- a/columnphysics/icepack_shortwave.F90 +++ b/columnphysics/icepack_shortwave.F90 @@ -104,6 +104,8 @@ subroutine shortwave_ccsm3 (aicen, vicen, & fswthru_vdf, & fswthru_idr, & fswthru_idf, & + druvr, dfuvr, & + drpar, dfpar, & fswpenl, & Iswabs, SSwabs, & albin, albsn, & @@ -126,6 +128,12 @@ subroutine shortwave_ccsm3 (aicen, vicen, & swidr , & ! sw down, near IR, direct (W/m^2) swidf ! sw down, near IR, diffuse (W/m^2) + real (kind=dbl_kind), intent(in), optional :: & + druvr , & ! + dfuvr , & ! + drpar , & ! + dfpar ! + ! baseline albedos for ccsm3 shortwave, set in namelist real (kind=dbl_kind), intent(in) :: & albicev , & ! visible ice albedo for h > ahmax @@ -190,6 +198,29 @@ subroutine shortwave_ccsm3 (aicen, vicen, & character(len=*),parameter :: subname='(shortwave_ccsm3)' +#ifdef GEOSCOUPLED + if (.not. present(druvr)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing druvr" ) + return + endif + if (.not. present(dfuvr)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing dfuvr" ) + return + endif + if (.not. present(drpar)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing drpar" ) + return + endif + if (.not. present(dfpar)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing dfpar" ) + return + endif +#endif + !----------------------------------------------------------------- ! Solar radiation: albedo and absorbed shortwave !----------------------------------------------------------------- @@ -304,6 +335,12 @@ subroutine shortwave_ccsm3 (aicen, vicen, & fswthru_vdf=l_fswthru_vdf(n),& fswthru_idr=l_fswthru_idr(n),& fswthru_idf=l_fswthru_idf(n),& +#ifdef GEOSCOUPLED + druvr = druvr, & + dfuvr = dfuvr, & + drpar = drpar, & + dfpar = dfpar, & +#endif fswpenl=fswpenl(:,n), & Iswabs=Iswabs(:,n)) @@ -580,6 +617,8 @@ subroutine absorbed_solar (heat_capacity, & fswthru_vdf, & fswthru_idr, & fswthru_idf, & + druvr, dfuvr, & + drpar, dfpar, & fswpenl, & Iswabs) @@ -606,6 +645,12 @@ subroutine absorbed_solar (heat_capacity, & alvdfns , & ! visible, diffuse albedo, snow alidfns ! near-ir, diffuse albedo, snow + real (kind=dbl_kind), intent(in), optional :: & + druvr , & ! + dfuvr , & ! + drpar , & ! + dfpar ! + real (kind=dbl_kind), intent(out):: & fswsfc , & ! SW absorbed at ice/snow surface (W m-2) fswint , & ! SW absorbed in ice interior, below surface (W m-2) @@ -645,11 +690,43 @@ subroutine absorbed_solar (heat_capacity, & hilyr , & ! ice layer thickness asnow ! fractional area of snow cover +#ifdef GEOSCOUPLED + real (kind=dbl_kind) :: & + druvrpen , & ! + dfuvrpen , & ! + drparpen , & ! + dfparpen ! +#endif + + + character(len=*),parameter :: subname='(absorbed_solar)' !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- +#ifdef GEOSCOUPLED + if (.not. present(druvr)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing druvr" ) + return + endif + if (.not. present(dfuvr)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing dfuvr" ) + return + endif + if (.not. present(drpar)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing drpar" ) + return + endif + if (.not. present(dfpar)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing dfpar" ) + return + endif +#endif trantop = c0 tranbot = c0 @@ -688,6 +765,13 @@ subroutine absorbed_solar (heat_capacity, & fswpenvdr = swvdr * (c1-alvdrni) * (c1-asnow) * i0vis fswpenvdf = swvdf * (c1-alvdfni) * (c1-asnow) * i0vis +#ifdef GEOSCOUPLED + druvrpen = druvr * (c1-alvdrni) * (c1-asnow) * i0vis + dfuvrpen = dfuvr * (c1-alvdfni) * (c1-asnow) * i0vis + drparpen = drpar * (c1-alvdrni) * (c1-asnow) * i0vis + dfparpen = dfpar * (c1-alvdfni) * (c1-asnow) * i0vis +#endif + ! no penetrating radiation in near IR ! fswpenidr = swidr * (c1-alidrni) * (c1-asnow) * i0nir ! fswpenidf = swidf * (c1-alidfni) * (c1-asnow) * i0nir @@ -724,10 +808,17 @@ subroutine absorbed_solar (heat_capacity, & ! SW penetrating thru ice into ocean fswthru = fswpen * tranbot +#ifdef GEOSCOUPLED + fswthru_vdr = druvrpen * tranbot + fswthru_vdf = dfuvrpen * tranbot + fswthru_idr = drparpen * tranbot + fswthru_idf = dfparpen * tranbot +#else fswthru_vdr = fswpenvdr * tranbot fswthru_vdf = fswpenvdf * tranbot fswthru_idr = c0 fswthru_idf = c0 +#endif ! SW absorbed in ice interior fswint = fswpen - fswthru @@ -4046,6 +4137,8 @@ subroutine icepack_step_radiation (dt, ncat, & modal_aero, & swvdr, swvdf, & swidr, swidf, & + druvr, dfuvr, & + drpar, dfpar, & coszen, fsnow, & alvdrn, alvdfn, & alidrn, alidfn, & @@ -4155,6 +4248,12 @@ subroutine icepack_step_radiation (dt, ncat, & fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2) fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) + real (kind=dbl_kind), intent(in), optional :: & + druvr , & ! + dfuvr , & ! + drpar , & ! + dfpar ! + real (kind=dbl_kind), dimension(:,:), intent(inout) :: & fswpenln , & ! visible SW entering ice layers (W m-2) Iswabsn , & ! SW radiation absorbed in ice layers (W m-2) @@ -4196,6 +4295,29 @@ subroutine icepack_step_radiation (dt, ncat, & character(len=*),parameter :: subname='(icepack_step_radiation)' +#ifdef GEOSCOUPLED + if (.not. present(druvr)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing druvr" ) + return + endif + if (.not. present(dfuvr)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing dfuvr" ) + return + endif + if (.not. present(drpar)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing drpar" ) + return + endif + if (.not. present(dfpar)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing dfpar" ) + return + endif +#endif + allocate(l_fswthrun_vdr(ncat)) allocate(l_fswthrun_vdf(ncat)) allocate(l_fswthrun_idr(ncat)) @@ -4321,6 +4443,12 @@ subroutine icepack_step_radiation (dt, ncat, & fswthru_vdf=l_fswthrun_vdf,& fswthru_idr=l_fswthrun_idr,& fswthru_idf=l_fswthrun_idf,& +#ifdef GEOSCOUPLED + druvr = druvr, & + dfuvr = dfuvr, & + drpar = drpar, & + dfpar = dfpar, & +#endif fswpenl=fswpenln, & Iswabs=Iswabsn, & Sswabs=Sswabsn, & From a588900b47122e576da429edcf9e2f69c5f7f28c Mon Sep 17 00:00:00 2001 From: bzhao Date: Mon, 30 Jan 2023 18:28:33 -0500 Subject: [PATCH 11/30] refactored code to get addtional penertrative sw fluxes through ice; remove ifdefs and use optional arguments --- columnphysics/icepack_shortwave.F90 | 254 +++++++++++++++++----------- 1 file changed, 158 insertions(+), 96 deletions(-) diff --git a/columnphysics/icepack_shortwave.F90 b/columnphysics/icepack_shortwave.F90 index 14829a158..744fc3d4a 100644 --- a/columnphysics/icepack_shortwave.F90 +++ b/columnphysics/icepack_shortwave.F90 @@ -104,6 +104,10 @@ subroutine shortwave_ccsm3 (aicen, vicen, & fswthru_vdf, & fswthru_idr, & fswthru_idf, & + fswthru_uvrdr, & + fswthru_uvrdf, & + fswthru_pardr, & + fswthru_pardf, & druvr, dfuvr, & drpar, dfpar, & fswpenl, & @@ -165,6 +169,12 @@ subroutine shortwave_ccsm3 (aicen, vicen, & fswthru_idr , & ! nir dir SW through ice to ocean (W m-2) fswthru_idf ! nir dif SW through ice to ocean (W m-2) + real (kind=dbl_kind), dimension(:), intent(out), optional :: & + fswthru_uvrdr , & ! vis dir SW through ice to ocean (W/m^2) + fswthru_uvrdf , & ! vis dif SW through ice to ocean (W/m^2) + fswthru_pardr , & ! nir dir SW through ice to ocean (W/m^2) + fswthru_pardf ! nir dif SW through ice to ocean (W/m^2) + real (kind=dbl_kind), intent(inout) :: & coszen ! cosine(zenith angle) @@ -196,30 +206,40 @@ subroutine shortwave_ccsm3 (aicen, vicen, & l_fswthru_idr , & ! nir dir SW through ice to ocean (W m-2) l_fswthru_idf ! nir dif SW through ice to ocean (W m-2) + real (kind=dbl_kind), dimension(:), allocatable :: & + l_fswthru_uvrdr , & ! vis uvr dir SW through ice to ocean (W/m^2) + l_fswthru_uvrdf , & ! vis uvr dif SW through ice to ocean (W/m^2) + l_fswthru_pardr , & ! vis par dir SW through ice to ocean (W/m^2) + l_fswthru_pardf ! vis par dif SW through ice to ocean (W/m^2) + + real(kind=dbl_kind) :: & + l_druvr, & ! + l_dfuvr, & ! + l_drpar, & ! + l_dfpar ! + character(len=*),parameter :: subname='(shortwave_ccsm3)' -#ifdef GEOSCOUPLED - if (.not. present(druvr)) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//": missing druvr" ) - return - endif - if (.not. present(dfuvr)) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//": missing dfuvr" ) - return - endif - if (.not. present(drpar)) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//": missing drpar" ) - return - endif - if (.not. present(dfpar)) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//": missing dfpar" ) - return - endif -#endif + if (present(druvr)) then + l_druvr = druvr + else + l_druvr = c0 + endif + if (present(dfuvr)) then + l_dfuvr = dfuvr + else + l_dfuvr = c0 + endif + if (present(drpar)) then + l_drpar = drpar + else + l_drpar = c0 + endif + if (present(dfpar)) then + l_dfpar = dfpar + else + l_dfpar = c0 + endif !----------------------------------------------------------------- ! Solar radiation: albedo and absorbed shortwave @@ -230,6 +250,16 @@ subroutine shortwave_ccsm3 (aicen, vicen, & allocate(l_fswthru_idr(ncat)) allocate(l_fswthru_idf(ncat)) + allocate(l_fswthru_uvrdr(ncat)) + allocate(l_fswthru_uvrdf(ncat)) + allocate(l_fswthru_pardr(ncat)) + allocate(l_fswthru_pardf(ncat)) + + l_fswthru_uvrdr = c0 + l_fswthru_uvrdf = c0 + l_fswthru_pardr = c0 + l_fswthru_pardf = c0 + ! For basic shortwave, set coszen to a constant between 0 and 1. coszen = p5 ! sun above the horizon @@ -335,12 +365,14 @@ subroutine shortwave_ccsm3 (aicen, vicen, & fswthru_vdf=l_fswthru_vdf(n),& fswthru_idr=l_fswthru_idr(n),& fswthru_idf=l_fswthru_idf(n),& -#ifdef GEOSCOUPLED - druvr = druvr, & - dfuvr = dfuvr, & - drpar = drpar, & - dfpar = dfpar, & -#endif + fswthru_uvrdr=l_fswthru_uvrdr(n),& + fswthru_uvrdf=l_fswthru_uvrdf(n),& + fswthru_pardr=l_fswthru_pardr(n),& + fswthru_pardf=l_fswthru_pardf(n),& + druvr = l_druvr, & + dfuvr = l_dfuvr, & + drpar = l_drpar, & + dfpar = l_dfpar, & fswpenl=fswpenl(:,n), & Iswabs=Iswabs(:,n)) @@ -360,6 +392,17 @@ subroutine shortwave_ccsm3 (aicen, vicen, & deallocate(l_fswthru_idr) deallocate(l_fswthru_idf) + if (present(fswthru_uvrdr)) fswthru_uvrdr = l_fswthru_uvrdr + if (present(fswthru_uvrdf)) fswthru_uvrdf = l_fswthru_uvrdf + if (present(fswthru_pardr)) fswthru_pardr = l_fswthru_pardr + if (present(fswthru_pardf)) fswthru_pardf = l_fswthru_pardf + + deallocate(l_fswthru_uvrdr) + deallocate(l_fswthru_uvrdf) + deallocate(l_fswthru_pardr) + deallocate(l_fswthru_pardf) + + end subroutine shortwave_ccsm3 !======================================================================= @@ -617,6 +660,10 @@ subroutine absorbed_solar (heat_capacity, & fswthru_vdf, & fswthru_idr, & fswthru_idf, & + fswthru_uvrdr, & + fswthru_uvrdf, & + fswthru_pardr, & + fswthru_pardf, & druvr, dfuvr, & drpar, dfpar, & fswpenl, & @@ -645,7 +692,7 @@ subroutine absorbed_solar (heat_capacity, & alvdfns , & ! visible, diffuse albedo, snow alidfns ! near-ir, diffuse albedo, snow - real (kind=dbl_kind), intent(in), optional :: & + real (kind=dbl_kind), intent(in) :: & druvr , & ! dfuvr , & ! drpar , & ! @@ -662,6 +709,12 @@ subroutine absorbed_solar (heat_capacity, & fswthru_idr , & ! nir dir SW through ice to ocean (W m-2) fswthru_idf ! nir dif SW through ice to ocean (W m-2) + real (kind=dbl_kind), intent(out) :: & + fswthru_uvrdr , & ! vis dir SW through ice to ocean (W/m^2) + fswthru_uvrdf , & ! vis dif SW through ice to ocean (W/m^2) + fswthru_pardr , & ! nir dir SW through ice to ocean (W/m^2) + fswthru_pardf ! nir dif SW through ice to ocean (W/m^2) + real (kind=dbl_kind), dimension (:), intent(out) :: & Iswabs , & ! SW absorbed in particular layer (W m-2) fswpenl ! visible SW entering ice layers (W m-2) @@ -690,43 +743,17 @@ subroutine absorbed_solar (heat_capacity, & hilyr , & ! ice layer thickness asnow ! fractional area of snow cover -#ifdef GEOSCOUPLED real (kind=dbl_kind) :: & druvrpen , & ! dfuvrpen , & ! drparpen , & ! dfparpen ! -#endif - - character(len=*),parameter :: subname='(absorbed_solar)' !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- -#ifdef GEOSCOUPLED - if (.not. present(druvr)) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//": missing druvr" ) - return - endif - if (.not. present(dfuvr)) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//": missing dfuvr" ) - return - endif - if (.not. present(drpar)) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//": missing drpar" ) - return - endif - if (.not. present(dfpar)) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//": missing dfpar" ) - return - endif -#endif trantop = c0 tranbot = c0 @@ -765,16 +792,14 @@ subroutine absorbed_solar (heat_capacity, & fswpenvdr = swvdr * (c1-alvdrni) * (c1-asnow) * i0vis fswpenvdf = swvdf * (c1-alvdfni) * (c1-asnow) * i0vis -#ifdef GEOSCOUPLED + ! no penetrating radiation in near IR +! fswpenidr = swidr * (c1-alidrni) * (c1-asnow) * i0nir +! fswpenidf = swidf * (c1-alidfni) * (c1-asnow) * i0nir + druvrpen = druvr * (c1-alvdrni) * (c1-asnow) * i0vis dfuvrpen = dfuvr * (c1-alvdfni) * (c1-asnow) * i0vis drparpen = drpar * (c1-alvdrni) * (c1-asnow) * i0vis dfparpen = dfpar * (c1-alvdfni) * (c1-asnow) * i0vis -#endif - - ! no penetrating radiation in near IR -! fswpenidr = swidr * (c1-alidrni) * (c1-asnow) * i0nir -! fswpenidf = swidf * (c1-alidfni) * (c1-asnow) * i0nir fswpen = fswpenvdr + fswpenvdf @@ -808,17 +833,14 @@ subroutine absorbed_solar (heat_capacity, & ! SW penetrating thru ice into ocean fswthru = fswpen * tranbot -#ifdef GEOSCOUPLED - fswthru_vdr = druvrpen * tranbot - fswthru_vdf = dfuvrpen * tranbot - fswthru_idr = drparpen * tranbot - fswthru_idf = dfparpen * tranbot -#else fswthru_vdr = fswpenvdr * tranbot fswthru_vdf = fswpenvdf * tranbot fswthru_idr = c0 fswthru_idf = c0 -#endif + fswthru_uvrdr = druvrpen * tranbot + fswthru_uvrdf = dfuvrpen * tranbot + fswthru_pardr = drparpen * tranbot + fswthru_pardf = dfparpen * tranbot ! SW absorbed in ice interior fswint = fswpen - fswthru @@ -4148,6 +4170,10 @@ subroutine icepack_step_radiation (dt, ncat, & fswthrun_vdf, & fswthrun_idr, & fswthrun_idf, & + fswthrun_uvrdr, & + fswthrun_uvrdf, & + fswthrun_pardr, & + fswthrun_pardf, & fswpenln, & Sswabsn, Iswabsn, & albicen, albsnon, & @@ -4248,6 +4274,12 @@ subroutine icepack_step_radiation (dt, ncat, & fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2) fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & + fswthrun_uvrdr , & ! vis dir SW through ice to ocean (W/m^2) + fswthrun_uvrdf , & ! vis dif SW through ice to ocean (W/m^2) + fswthrun_pardr , & ! nir dir SW through ice to ocean (W/m^2) + fswthrun_pardf ! nir dif SW through ice to ocean (W/m^2) + real (kind=dbl_kind), intent(in), optional :: & druvr , & ! dfuvr , & ! @@ -4290,38 +4322,57 @@ subroutine icepack_step_radiation (dt, ncat, & l_fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2) l_fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) + real (kind=dbl_kind), dimension(:), allocatable :: & + l_fswthrun_uvrdr , & ! vis uvr dir SW through ice to ocean (W/m^2) + l_fswthrun_uvrdf , & ! vis uvr dif SW through ice to ocean (W/m^2) + l_fswthrun_pardr , & ! vis par dir SW through ice to ocean (W/m^2) + l_fswthrun_pardf ! vis par dif SW through ice to ocean (W/m^2) + + real(kind=dbl_kind) :: & + l_druvr, & ! + l_dfuvr, & ! + l_drpar, & ! + l_dfpar ! + real (kind=dbl_kind), dimension(:,:), allocatable :: & l_rsnow ! snow grain radius tracer (10^-6 m) character(len=*),parameter :: subname='(icepack_step_radiation)' -#ifdef GEOSCOUPLED - if (.not. present(druvr)) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//": missing druvr" ) - return - endif - if (.not. present(dfuvr)) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//": missing dfuvr" ) - return - endif - if (.not. present(drpar)) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//": missing drpar" ) - return - endif - if (.not. present(dfpar)) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//": missing dfpar" ) - return - endif -#endif + if (present(druvr)) then + l_druvr = druvr + else + l_druvr = c0 + endif + if (present(dfuvr)) then + l_dfuvr = dfuvr + else + l_dfuvr = c0 + endif + if (present(drpar)) then + l_drpar = drpar + else + l_drpar = c0 + endif + if (present(dfpar)) then + l_dfpar = dfpar + else + l_dfpar = c0 + endif allocate(l_fswthrun_vdr(ncat)) allocate(l_fswthrun_vdf(ncat)) allocate(l_fswthrun_idr(ncat)) allocate(l_fswthrun_idf(ncat)) + allocate(l_fswthrun_uvrdr(ncat)) + allocate(l_fswthrun_uvrdf(ncat)) + allocate(l_fswthrun_pardr(ncat)) + allocate(l_fswthrun_pardf(ncat)) + + l_fswthrun_uvrdr = c0 + l_fswthrun_uvrdf = c0 + l_fswthrun_pardr = c0 + l_fswthrun_pardf = c0 hin = c0 hbri = c0 @@ -4443,12 +4494,14 @@ subroutine icepack_step_radiation (dt, ncat, & fswthru_vdf=l_fswthrun_vdf,& fswthru_idr=l_fswthrun_idr,& fswthru_idf=l_fswthrun_idf,& -#ifdef GEOSCOUPLED - druvr = druvr, & - dfuvr = dfuvr, & - drpar = drpar, & - dfpar = dfpar, & -#endif + fswthru_uvrdr=l_fswthrun_uvrdr,& + fswthru_uvrdf=l_fswthrun_uvrdf,& + fswthru_pardr=l_fswthrun_pardr,& + fswthru_pardf=l_fswthrun_pardf,& + druvr = l_druvr, & + dfuvr = l_dfuvr, & + drpar = l_drpar, & + dfpar = l_dfpar, & fswpenl=fswpenln, & Iswabs=Iswabsn, & Sswabs=Sswabsn, & @@ -4507,10 +4560,19 @@ subroutine icepack_step_radiation (dt, ncat, & if (present(fswthrun_idr)) fswthrun_idr = l_fswthrun_idr if (present(fswthrun_idf)) fswthrun_idf = l_fswthrun_idf + if (present(fswthrun_uvrdr)) fswthrun_uvrdr = l_fswthrun_uvrdr + if (present(fswthrun_uvrdf)) fswthrun_uvrdf = l_fswthrun_uvrdf + if (present(fswthrun_pardr)) fswthrun_pardr = l_fswthrun_pardr + if (present(fswthrun_pardf)) fswthrun_pardf = l_fswthrun_pardf + deallocate(l_fswthrun_vdr) deallocate(l_fswthrun_vdf) deallocate(l_fswthrun_idr) deallocate(l_fswthrun_idf) + deallocate(l_fswthrun_uvrdr) + deallocate(l_fswthrun_uvrdf) + deallocate(l_fswthrun_pardr) + deallocate(l_fswthrun_pardf) deallocate(l_rsnow) end subroutine icepack_step_radiation From 711bcafb0f9ae8595309fe1fcb5c05e8f172f014 Mon Sep 17 00:00:00 2001 From: bzhao Date: Thu, 2 Feb 2023 18:19:43 -0500 Subject: [PATCH 12/30] add additional penetrating sw flux to merge_fluxes call --- columnphysics/icepack_flux.F90 | 23 ++++++++ columnphysics/icepack_therm_vertical.F90 | 72 ++++++++++++++++++++++++ 2 files changed, 95 insertions(+) diff --git a/columnphysics/icepack_flux.F90 b/columnphysics/icepack_flux.F90 index 4f951fff2..95482bbc7 100644 --- a/columnphysics/icepack_flux.F90 +++ b/columnphysics/icepack_flux.F90 @@ -43,6 +43,8 @@ subroutine merge_fluxes (aicen, & fhocnn, fswthrun, & fswthrun_vdr, fswthrun_vdf,& fswthrun_idr, fswthrun_idf,& + fswthrun_uvrdr, fswthrun_uvrdf,& + fswthrun_pardr, fswthrun_pardf,& strairxT, strairyT, & Cdn_atm_ratio, & fsurf, fcondtop, & @@ -56,6 +58,8 @@ subroutine merge_fluxes (aicen, & fhocn, fswthru, & fswthru_vdr, fswthru_vdf,& fswthru_idr, fswthru_idf,& + fswthru_uvrdr, fswthru_uvrdf,& + fswthru_pardr, fswthru_pardf,& melttn, meltsn, meltbn, congeln, snoicen, & meltt, melts, & meltb, dsnow, dsnown,& @@ -93,6 +97,10 @@ subroutine merge_fluxes (aicen, & fswthrun_vdf, & ! vis dif sw radiation through ice bot (W/m**2) fswthrun_idr, & ! nir dir sw radiation through ice bot (W/m**2) fswthrun_idf, & ! nir dif sw radiation through ice bot (W/m**2) + fswthrun_uvrdr, & ! vis dir sw radiation through ice bot (W/m**2) + fswthrun_uvrdf, & ! vis dir sw radiation through ice bot (W/m**2) + fswthrun_pardr, & ! vis dir sw radiation through ice bot (W/m**2) + fswthrun_pardf, & ! vis dir sw radiation through ice bot (W/m**2) melttn , & ! top ice melt (m) meltbn , & ! bottom ice melt (m) meltsn , & ! snow melt (m) @@ -139,6 +147,12 @@ subroutine merge_fluxes (aicen, & fswthru_idr , & ! nir dir sw radiation through ice bot (W/m**2) fswthru_idf ! nir dif sw radiation through ice bot (W/m**2) + real (kind=dbl_kind), intent(inout), optional :: & + fswthru_uvrdr , & ! vis dir sw radiation through ice bot (W/m**2) + fswthru_uvrdf , & ! vis dif sw radiation through ice bot (W/m**2) + fswthru_pardr , & ! nir dir sw radiation through ice bot (W/m**2) + fswthru_pardf ! nir dif sw radiation through ice bot (W/m**2) + real (kind=dbl_kind), optional, intent(inout):: & Uref ! air speed reference level (m/s) @@ -212,6 +226,15 @@ subroutine merge_fluxes (aicen, & if (present(fswthru_idf)) & fswthru_idf = fswthru_idf + fswthrun_idf * aicen + if (present(fswthru_uvrdr)) & + fswthru_uvrdr = fswthru_uvrdr + fswthrun_uvrdr * aicen + if (present(fswthru_uvrdf)) & + fswthru_uvrdf = fswthru_uvrdf + fswthrun_uvrdf * aicen + if (present(fswthru_pardr)) & + fswthru_pardr = fswthru_pardr + fswthrun_pardr * aicen + if (present(fswthru_pardf)) & + fswthru_pardf = fswthru_pardf + fswthrun_pardf * aicen + ! ice/snow thickness meltt = meltt + melttn * aicen diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 183e321d7..a5cf99497 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -2279,6 +2279,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fswthrun_vdf, & fswthrun_idr, & fswthrun_idf, & + fswthrun_uvrdr, & + fswthrun_uvrdf, & + fswthrun_pardr, & + fswthrun_pardf, & fswabs , & flwout , & Sswabsn , Iswabsn , & @@ -2294,6 +2298,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fswthru_vdf , & fswthru_idr , & fswthru_idf , & + fswthru_uvrdr , & + fswthru_uvrdf , & + fswthru_pardr , & + fswthru_pardf , & flatn_f , fsensn_f , & fsurfn_f , fcondtopn_f , & faero_atm , faero_ocn , & @@ -2423,6 +2431,12 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & meltsliq , & ! mass of snow melt (kg/m^2) fsloss ! rate of snow loss to leads (kg/m^2/s) + real (kind=dbl_kind), intent(inout), optional :: & + fswthru_uvrdr , & ! vis dir shortwave penetrating to ocean (W/m^2) + fswthru_uvrdf , & ! vis dif shortwave penetrating to ocean (W/m^2) + fswthru_pardr , & ! nir dir shortwave penetrating to ocean (W/m^2) + fswthru_pardf ! nir dif shortwave penetrating to ocean (W/m^2) + real (kind=dbl_kind), dimension(:), optional, intent(inout) :: & Qa_iso , & ! isotope specific humidity (kg/kg) Qref_iso , & ! isotope 2m atm reference spec humidity (kg/kg) @@ -2490,6 +2504,12 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2) fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) + real (kind=dbl_kind), optional, dimension(:), intent(inout) :: & + fswthrun_uvrdr , & ! vis dir SW through ice to ocean (W/m^2) + fswthrun_uvrdf , & ! vis dif SW through ice to ocean (W/m^2) + fswthrun_pardr , & ! nir dir SW through ice to ocean (W/m^2) + fswthrun_pardf ! nir dif SW through ice to ocean (W/m^2) + real (kind=dbl_kind), dimension(:,:), intent(inout) :: & zqsn , & ! snow layer enthalpy (J m-3) zqin , & ! ice layer enthalpy (J m-3) @@ -2575,12 +2595,24 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & l_fswthru_idr , & ! nir dir SW through ice to ocean (W/m^2) l_fswthru_idf ! nir dif SW through ice to ocean (W/m^2) + real (kind=dbl_kind) :: & + l_fswthru_uvrdr , & ! vis dir SW through ice to ocean (W/m^2) + l_fswthru_uvrdf , & ! vis dir SW through ice to ocean (W/m^2) + l_fswthru_pardr , & ! vis dir SW through ice to ocean (W/m^2) + l_fswthru_pardf ! vis dir SW through ice to ocean (W/m^2) + real (kind=dbl_kind), dimension(:), allocatable :: & l_fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2) l_fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2) l_fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2) l_fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) + real (kind=dbl_kind), dimension(:), allocatable :: & + l_fswthrun_uvrdr , & ! vis dir SW through ice to ocean (W/m^2) + l_fswthrun_uvrdf , & ! vis dir SW through ice to ocean (W/m^2) + l_fswthrun_pardr , & ! vis dir SW through ice to ocean (W/m^2) + l_fswthrun_pardf ! vis dir SW through ice to ocean (W/m^2) + real (kind=dbl_kind) :: & pond ! water retained in ponds (m) @@ -2702,6 +2734,26 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & l_fswthrun_idf = c0 if (present(fswthrun_idf)) l_fswthrun_idf = fswthrun_idf + allocate(l_fswthrun_vdf(ncat)) + l_fswthrun_vdf = c0 + if (present(fswthrun_vdf)) l_fswthrun_vdf = fswthrun_vdf + + allocate(l_fswthrun_uvrdr(ncat)) + l_fswthrun_uvrdr = c0 + if (present(fswthrun_uvrdr)) l_fswthrun_uvrdr = fswthrun_uvrdr + + allocate(l_fswthrun_uvrdf(ncat)) + l_fswthrun_uvrdf = c0 + if (present(fswthrun_uvrdf)) l_fswthrun_uvrdf = fswthrun_uvrdf + + allocate(l_fswthrun_pardr(ncat)) + l_fswthrun_pardr = c0 + if (present(fswthrun_pardr)) l_fswthrun_pardr = fswthrun_pardr + + allocate(l_fswthrun_pardf(ncat)) + l_fswthrun_pardf = c0 + if (present(fswthrun_pardf)) l_fswthrun_pardf = fswthrun_pardf + allocate(l_meltsliqn(ncat)) l_meltsliqn = c0 if (present(meltsliqn)) l_meltsliqn = meltsliqn @@ -3113,6 +3165,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fswthrun_vdf=l_fswthrun_vdf(n), & fswthrun_idr=l_fswthrun_idr(n), & fswthrun_idf=l_fswthrun_idf(n), & + fswthrun_uvrdr=l_fswthrun_uvrdr(n), & + fswthrun_uvrdf=l_fswthrun_uvrdf(n), & + fswthrun_pardr=l_fswthrun_pardr(n), & + fswthrun_pardf=l_fswthrun_pardf(n), & strairxT=strairxT, strairyT=strairyT,& Cdn_atm_ratio=Cdn_atm_ratio, & fsurf=fsurf, fcondtop=fcondtop,& @@ -3129,6 +3185,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fswthru_vdf=l_fswthru_vdf, & fswthru_idr=l_fswthru_idr, & fswthru_idf=l_fswthru_idf, & + fswthru_uvrdr=l_fswthru_uvrdr, & + fswthru_uvrdf=l_fswthru_uvrdf, & + fswthru_pardr=l_fswthru_pardr, & + fswthru_pardf=l_fswthru_pardf, & melttn=melttn (n), meltsn=meltsn(n), & meltbn=meltbn (n), congeln=congeln(n),& meltt=meltt, melts=melts, & @@ -3185,10 +3245,18 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & if (present(fswthrun_vdf)) fswthrun_vdf = l_fswthrun_vdf if (present(fswthrun_idr)) fswthrun_idr = l_fswthrun_idr if (present(fswthrun_idf)) fswthrun_idf = l_fswthrun_idf + if (present(fswthrun_uvrdr)) fswthrun_uvrdr = l_fswthrun_uvrdr + if (present(fswthrun_uvrdf)) fswthrun_uvrdf = l_fswthrun_uvrdf + if (present(fswthrun_pardr)) fswthrun_pardr = l_fswthrun_pardr + if (present(fswthrun_pardf)) fswthrun_pardf = l_fswthrun_pardf if (present(fswthru_vdr )) fswthru_vdr = l_fswthru_vdr if (present(fswthru_vdf )) fswthru_vdf = l_fswthru_vdf if (present(fswthru_idr )) fswthru_idr = l_fswthru_idr if (present(fswthru_idf )) fswthru_idf = l_fswthru_idf + if (present(fswthru_uvrdr )) fswthru_uvrdr = l_fswthru_uvrdr + if (present(fswthru_uvrdf )) fswthru_uvrdf = l_fswthru_uvrdf + if (present(fswthru_pardr )) fswthru_pardr = l_fswthru_pardr + if (present(fswthru_pardf )) fswthru_pardf = l_fswthru_pardf if (present(fsloss )) fsloss = l_fsloss if (present(meltsliqn )) meltsliqn = l_meltsliqn if (present(meltsliq )) meltsliq = l_meltsliq @@ -3206,6 +3274,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & deallocate(l_fswthrun_vdf) deallocate(l_fswthrun_idr) deallocate(l_fswthrun_idf) + deallocate(l_fswthrun_uvrdr) + deallocate(l_fswthrun_uvrdf) + deallocate(l_fswthrun_pardr) + deallocate(l_fswthrun_pardf) deallocate(l_meltsliqn) deallocate(l_rsnw) deallocate(l_smice) From 200b1b4035313d86b874308661ef267c89d9b046 Mon Sep 17 00:00:00 2001 From: bzhao Date: Fri, 3 Feb 2023 19:34:56 -0500 Subject: [PATCH 13/30] fixed a bug --- columnphysics/icepack_therm_vertical.F90 | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index a5cf99497..5ff24625b 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -2718,6 +2718,18 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & l_fswthru_idf = c0 if (present(fswthru_idf)) l_fswthru_idf = fswthru_idf + l_fswthru_uvrdr = c0 + if (present(fswthru_uvrdr)) l_fswthru_uvrdr = fswthru_uvrdr + + l_fswthru_uvrdf = c0 + if (present(fswthru_uvrdf)) l_fswthru_uvrdf = fswthru_uvrdf + + l_fswthru_pardr = c0 + if (present(fswthru_pardr)) l_fswthru_pardr = fswthru_pardr + + l_fswthru_pardf = c0 + if (present(fswthru_pardf)) l_fswthru_pardf = fswthru_pardf + allocate(l_fswthrun_vdr(ncat)) l_fswthrun_vdr = c0 if (present(fswthrun_vdr)) l_fswthrun_vdr = fswthrun_vdr @@ -2734,10 +2746,6 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & l_fswthrun_idf = c0 if (present(fswthrun_idf)) l_fswthrun_idf = fswthrun_idf - allocate(l_fswthrun_vdf(ncat)) - l_fswthrun_vdf = c0 - if (present(fswthrun_vdf)) l_fswthrun_vdf = fswthrun_vdf - allocate(l_fswthrun_uvrdr(ncat)) l_fswthrun_uvrdr = c0 if (present(fswthrun_uvrdr)) l_fswthrun_uvrdr = fswthrun_uvrdr From 800d9f5aea980bafa70a7249b9b52447a9bdaf1f Mon Sep 17 00:00:00 2001 From: bzhao Date: Mon, 6 Feb 2023 18:43:28 -0500 Subject: [PATCH 14/30] fixed a bug causing fswthru_*dr/df to accumulate values over time steps --- columnphysics/icepack_therm_vertical.F90 | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 5ff24625b..9f3814288 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -2431,7 +2431,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & meltsliq , & ! mass of snow melt (kg/m^2) fsloss ! rate of snow loss to leads (kg/m^2/s) - real (kind=dbl_kind), intent(inout), optional :: & + real (kind=dbl_kind), intent(out), optional :: & fswthru_uvrdr , & ! vis dir shortwave penetrating to ocean (W/m^2) fswthru_uvrdf , & ! vis dif shortwave penetrating to ocean (W/m^2) fswthru_pardr , & ! nir dir shortwave penetrating to ocean (W/m^2) @@ -2719,16 +2719,12 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & if (present(fswthru_idf)) l_fswthru_idf = fswthru_idf l_fswthru_uvrdr = c0 - if (present(fswthru_uvrdr)) l_fswthru_uvrdr = fswthru_uvrdr l_fswthru_uvrdf = c0 - if (present(fswthru_uvrdf)) l_fswthru_uvrdf = fswthru_uvrdf l_fswthru_pardr = c0 - if (present(fswthru_pardr)) l_fswthru_pardr = fswthru_pardr l_fswthru_pardf = c0 - if (present(fswthru_pardf)) l_fswthru_pardf = fswthru_pardf allocate(l_fswthrun_vdr(ncat)) l_fswthrun_vdr = c0 From 49bb2cf91fd9038db47bf066267a3cc8d129a5f2 Mon Sep 17 00:00:00 2001 From: bzhao Date: Tue, 21 Feb 2023 21:10:23 -0500 Subject: [PATCH 15/30] do not add frain to fresh when coupled to GEOS --- columnphysics/icepack_therm_itd.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index 4af3848dd..28d204b7b 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -2127,8 +2127,9 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & !----------------------------------------------------------------- ! Let rain drain through to the ocean. !----------------------------------------------------------------- - +#ifndef GEOSCOUPLED fresh = fresh + frain * aice +#endif !----------------------------------------------------------------- ! Given thermodynamic growth rates, transport ice between From b35edc03e4f9ffb8c7dd55d947c4e8db0e7c7158 Mon Sep 17 00:00:00 2001 From: bzhao Date: Wed, 22 Mar 2023 16:50:58 -0400 Subject: [PATCH 16/30] not allowing fsurfn to be reset to zero when coupled to GEOS --- columnphysics/icepack_therm_vertical.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 9f3814288..526edbc25 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -302,7 +302,9 @@ subroutine thermo_vertical (nilyr, nslyr, & if (calc_Tsfc) then fsensn = c0 flatn = c0 +#ifndef GEOSCOUPLED fsurfn = c0 +#endif fcondtopn = c0 endif From 53a333238983dfbe491955634f9b26923e951022 Mon Sep 17 00:00:00 2001 From: bzhao Date: Fri, 5 May 2023 10:57:32 -0400 Subject: [PATCH 17/30] refactor thermo routines; use module vars to pass cpl fluxes and minimize routine interface changes --- columnphysics/icepack_therm_bl99.F90 | 42 ++++++++---------------- columnphysics/icepack_therm_shared.F90 | 8 +++++ columnphysics/icepack_therm_vertical.F90 | 38 +++++++++++++-------- 3 files changed, 45 insertions(+), 43 deletions(-) diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 92c6cf88b..3338b221a 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -24,6 +24,12 @@ module icepack_therm_bl99 use icepack_therm_shared, only: ferrmax, l_brine use icepack_therm_shared, only: surface_heat_flux, dsurface_heat_flux_dTsf +#ifdef GEOSCOUPLED + use icepack_therm_shared, only: dfsurfdts_cpl, & ! + dflatdts_cpl, & ! + fsurf_cpl, & ! + flat_cpl ! +#endif implicit none @@ -70,10 +76,7 @@ subroutine temperature_changes (dt, & fsensn, flatn, & flwoutn, fsurfn, & fcondtopn,fcondbot, & - einit, & - dfsurfdt_in, & - flatn_f, & - dflatdt_in) + einit ) integer (kind=int_kind), intent(in) :: & nilyr , & ! number of ice layers @@ -117,11 +120,6 @@ subroutine temperature_changes (dt, & flatn , & ! surface downward latent heat (W m-2) flwoutn ! upward LW at surface (W m-2) - real (kind=dbl_kind), intent(in), optional:: & - dfsurfdt_in , & - flatn_f , & - dflatdt_in - real (kind=dbl_kind), intent(out):: & fcondbot ! downward cond flux at bottom surface (W m-2) @@ -235,25 +233,10 @@ subroutine temperature_changes (dt, & dflwout_dT = c0 einex = c0 #ifdef GEOSCOUPLED - if (.not. present(dfsurfdt_in)) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//": missing dfsurfdt_in" ) - return - endif - if (.not. present(dflatdt_in)) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//": missing dflatdt_in" ) - return - endif - if (.not. present(flatn_f)) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//": missing flatn_f" ) - return - endif - ! derivative information is passed by GEOS - dfsurf_dT = dfsurfdt_in - dflat_dT = dflatdt_in - flatn = flatn_f + dfsurf_dT = dfsurfdts_cpl + dflat_dT = dflatdts_cpl + fsurfn = fsurf_cpl + flatn = flat_cpl #endif dt_rhoi_hlyr = dt / (rhoi*hilyr) ! hilyr > 0 if (hslyr > hs_min/real(nslyr,kind=dbl_kind)) & @@ -374,7 +357,6 @@ subroutine temperature_changes (dt, & converged = .true. #ifndef GEOSCOUPLED - ! derivative information is passed by GEOS dfsurf_dT = c0 #endif avg_Tsi = c0 @@ -700,7 +682,9 @@ subroutine temperature_changes (dt, & !----------------------------------------------------------------- fsurfn = fsurfn + dTsf*dfsurf_dT +#ifdef GEOSCOUPLED flatn = flatn + dTsf*dflat_dT +#endif if (l_snow) then fcondtopn = kh(1) * (Tsf-zTsn(1)) else diff --git a/columnphysics/icepack_therm_shared.F90 b/columnphysics/icepack_therm_shared.F90 index ad00db2ac..99027173d 100644 --- a/columnphysics/icepack_therm_shared.F90 +++ b/columnphysics/icepack_therm_shared.F90 @@ -51,6 +51,14 @@ module icepack_therm_shared real (kind=dbl_kind), public :: & hi_min ! minimum ice thickness allowed (m) +#ifdef GEOSCOUPLED + real (kind=dbl_kind), public :: & + dfsurfdts_cpl, & ! + dflatdts_cpl, & ! + fsurf_cpl, & ! + flat_cpl ! +#endif + !======================================================================= contains diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 67610def6..2750f40d9 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -37,6 +37,12 @@ module icepack_therm_vertical use icepack_therm_shared, only: calculate_tin_from_qin, Tmin use icepack_therm_shared, only: hi_min use icepack_therm_shared, only: adjust_enthalpy +#ifdef GEOSCOUPLED + use icepack_therm_shared, only: dfsurfdts_cpl, & ! + dflatdts_cpl, & ! + fsurf_cpl, & ! + flat_cpl ! +#endif use icepack_therm_bl99, only: temperature_changes use icepack_therm_mushy, only: temperature_changes_salinity @@ -104,8 +110,9 @@ subroutine thermo_vertical (nilyr, nslyr, & smliq, massliq, & congel, snoice, & mlt_onset, frz_onset, & + fsurfn_in, & + flatn_in, & dfsurfdt_in, & - flatn_f, & dflatdt_in, & yday, dsnow, & prescribed_ice, & @@ -191,8 +198,9 @@ subroutine thermo_vertical (nilyr, nslyr, & fcondbotn ! downward cond flux at bottom surface (W m-2) real (kind=dbl_kind), intent(in), optional:: & + fsurfn_in, & + flatn_in, & dfsurfdt_in, & - flatn_f, & dflatdt_in ! coupler fluxes to ocean @@ -291,19 +299,27 @@ subroutine thermo_vertical (nilyr, nslyr, & call icepack_warnings_add(subname//": missing dflatdt_in" ) return endif - if (.not. present(flatn_f)) then + if (.not. present(flatn_in)) then + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//": missing flatn_in" ) + return + endif + if (.not. present(fsurfn_in)) then call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//": missing flatn_f" ) + call icepack_warnings_add(subname//": missing fsurfn_in" ) return endif + + dfsurfdts_cpl = dfsurfdt_in + dflatdts_cpl = dflatdt_in + fsurf_cpl = fsurfn_in + flat_cpl = flatn_in #endif if (calc_Tsfc) then fsensn = c0 flatn = c0 -#ifndef GEOSCOUPLED fsurfn = c0 -#endif fcondtopn = c0 endif @@ -385,14 +401,7 @@ subroutine thermo_vertical (nilyr, nslyr, & fsensn, flatn, & flwoutn, fsurfn, & fcondtopn, fcondbotn, & -#ifdef GEOSCOUPLED - einit, & - dfsurfdt_in = dfsurfdt_in,& - flatn_f = flatn_f, & - dflatdt_in = dflatdt_in) -#else einit ) -#endif if (icepack_warnings_aborted(subname)) return endif ! ktherm @@ -2952,8 +2961,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & congel=congeln (n), snoice=snoicen (n), & mlt_onset=mlt_onset, frz_onset=frz_onset, & #ifdef GEOSCOUPLED + fsurfn_in = fsurfn_f(n), & + flatn_in = flatn_f(n), & dfsurfdt_in = dfsurfdt(n), & - flatn_f = flatn_f(n), & dflatdt_in = dflatdt(n), & #endif yday=yday, dsnow=dsnown (n), & From 7eeba49785c1f376322afefaf4a1d52180561399 Mon Sep 17 00:00:00 2001 From: bzhao Date: Mon, 8 May 2023 10:34:30 -0400 Subject: [PATCH 18/30] change account code to be able to submit jobs --- configuration/scripts/machines/env.discover_intel | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/configuration/scripts/machines/env.discover_intel b/configuration/scripts/machines/env.discover_intel index 4f13d1a5b..b8bdec18f 100755 --- a/configuration/scripts/machines/env.discover_intel +++ b/configuration/scripts/machines/env.discover_intel @@ -47,7 +47,7 @@ setenv ICE_MACHINE_WKDIR /discover/nobackup/$user/ICEPACK_RUNS setenv ICE_MACHINE_INPUTDATA /discover/nobackup/sakella/icepack_data setenv ICE_MACHINE_BASELINE /discover/nobackup/$user/ICEPACK_BASELINE setenv ICE_MACHINE_SUBMIT "sbatch" -setenv ICE_MACHINE_ACCT g0613 +setenv ICE_MACHINE_ACCT g0609 setenv ICE_MACHINE_QUEUE "share" setenv ICE_MACHINE_TPNODE 36 setenv ICE_MACHINE_BLDTHRDS 1 From f240d679e9d24f1a3ca15ecd24e081ceb546b6d3 Mon Sep 17 00:00:00 2001 From: bzhao Date: Mon, 8 May 2023 12:49:57 -0400 Subject: [PATCH 19/30] change batching test jobs to running them on interactive nodes because debug queue has limit --- configuration/scripts/icepack.test.setup.csh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/configuration/scripts/icepack.test.setup.csh b/configuration/scripts/icepack.test.setup.csh index b5d4ee433..ae79769fa 100755 --- a/configuration/scripts/icepack.test.setup.csh +++ b/configuration/scripts/icepack.test.setup.csh @@ -63,7 +63,8 @@ chmod +x ${jobfile} cat >! ${subfile} << EOFS #!/bin/csh -f -${ICE_MACHINE_SUBMIT} ./${jobfile} +#${ICE_MACHINE_SUBMIT} ./${jobfile} +./${jobfile} echo "\`date\` \${0}: ${ICE_CASENAME} job submitted" >> ${ICE_CASEDIR}/README.case EOFS From 81a7628405d49b7b91b2a704a74afbbd05be1252 Mon Sep 17 00:00:00 2001 From: bzhao Date: Thu, 11 May 2023 13:11:23 -0400 Subject: [PATCH 20/30] on discover, use interactive node to run tests; fix tsdir when not using suite; add a small_suite and a 15 day run set --- configuration/scripts/icepack.test.setup.csh | 15 ++++++++++++++- configuration/scripts/options/set_nml.run15day | 4 ++++ configuration/scripts/tests/small_suite.ts | 4 ++++ icepack.setup | 2 +- 4 files changed, 23 insertions(+), 2 deletions(-) create mode 100644 configuration/scripts/options/set_nml.run15day create mode 100644 configuration/scripts/tests/small_suite.ts diff --git a/configuration/scripts/icepack.test.setup.csh b/configuration/scripts/icepack.test.setup.csh index ae79769fa..6055f5357 100755 --- a/configuration/scripts/icepack.test.setup.csh +++ b/configuration/scripts/icepack.test.setup.csh @@ -60,14 +60,27 @@ cat >> ${jobfile} < ${ICE_SCRIPTS}/tests/baseline.script chmod +x ${jobfile} +if ($ICE_MACHINE == 'discover') then + cat >! ${subfile} << EOFS #!/bin/csh -f -#${ICE_MACHINE_SUBMIT} ./${jobfile} ./${jobfile} echo "\`date\` \${0}: ${ICE_CASENAME} job submitted" >> ${ICE_CASEDIR}/README.case EOFS +else + +cat >! ${subfile} << EOFS +#!/bin/csh -f + +${ICE_MACHINE_SUBMIT} ./${jobfile} +echo "\`date\` \${0}: ${ICE_CASENAME} job submitted" >> ${ICE_CASEDIR}/README.case + +EOFS + +endif + chmod +x ${subfile} exit 0 diff --git a/configuration/scripts/options/set_nml.run15day b/configuration/scripts/options/set_nml.run15day new file mode 100644 index 000000000..2ed171930 --- /dev/null +++ b/configuration/scripts/options/set_nml.run15day @@ -0,0 +1,4 @@ +npt = 360 +dumpfreq = 'd' + + diff --git a/configuration/scripts/tests/small_suite.ts b/configuration/scripts/tests/small_suite.ts new file mode 100644 index 000000000..ef3f2aa8b --- /dev/null +++ b/configuration/scripts/tests/small_suite.ts @@ -0,0 +1,4 @@ +# Test Grid PEs Sets BFB-compare +smoke col 1x1 diag1,run1year +smoke col 1x1 debug,run1year + diff --git a/icepack.setup b/icepack.setup index 4fee72006..b2ca4c97e 100755 --- a/icepack.setup +++ b/icepack.setup @@ -567,7 +567,7 @@ foreach envname ( $nenvnames ) if (${dosuite} == 1) then set testname = "${tsdir}/${testname_base}" else - set testname = "${testname_base}" + set testname = "${tsdir}/${testname_base}" endif set case = ${testname} endif From e7b4551ec39450ab2289e7b2024bcebbfae8fd62 Mon Sep 17 00:00:00 2001 From: bzhao Date: Tue, 22 Aug 2023 16:39:56 -0400 Subject: [PATCH 21/30] uncomment code recomputing zqin for BL99 --- columnphysics/icepack_therm_vertical.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 2750f40d9..9731802bc 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -1030,12 +1030,12 @@ subroutine init_vertical_profile(nilyr, nslyr, & endif ! echmod: is this necessary? -! if (ktherm == 1) then -! if (zTin(k)>= -zSin(k)*depressT) then -! zTin(k) = -zSin(k)*depressT - puny -! zqin(k) = -rhoi*cp_ocn*zSin(k)*depressT -! endif -! endif + if (ktherm == 1) then + if (zTin(k)>= -zSin(k)*depressT) then + zTin(k) = -zSin(k)*depressT - puny + zqin(k) = -rhoi*cp_ocn*zSin(k)*depressT + endif + endif !----------------------------------------------------------------- ! initial energy per unit area of ice/snow, relative to 0 C From abbbd626c5c48121ad02ceff125d341ae9f37f62 Mon Sep 17 00:00:00 2001 From: Bin Zhao Date: Thu, 31 Aug 2023 11:44:58 -0400 Subject: [PATCH 22/30] enable mushy layer thermodynamics scheme (#2) * first attempt at enabling mushy layer thero scheme in GEOS coupling; built ok; needs tests * fixed a tice_high bug in init_vertical; add updating fsurf/fsurf_cpl when Tsf gets reset to Tmlt * protect ismyturn with GEOSCOUPLED macro * protect ismyturn with another GEOSCOUPLED macro * remove printout i,j,tsk --- columnphysics/icepack_therm_mushy.F90 | 38 ++++++++++++++++++++++++ columnphysics/icepack_therm_shared.F90 | 20 +++++++++++++ columnphysics/icepack_therm_vertical.F90 | 32 ++++++++++++++++++++ 3 files changed, 90 insertions(+) diff --git a/columnphysics/icepack_therm_mushy.F90 b/columnphysics/icepack_therm_mushy.F90 index 80abd03f5..7c1526012 100644 --- a/columnphysics/icepack_therm_mushy.F90 +++ b/columnphysics/icepack_therm_mushy.F90 @@ -21,6 +21,15 @@ module icepack_therm_mushy use icepack_therm_shared, only: ferrmax use icepack_warnings, only: warnstr, icepack_warnings_add use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted +#ifdef GEOSCOUPLED + use icepack_therm_shared, only: dfsurfdts_cpl, & ! + dflatdts_cpl, & ! + flat_cpl0, & ! + fsurf_cpl0, & ! + fsurf_cpl, & ! + flat_cpl ! + use icepack_therm_shared, only: ismyturn +#endif implicit none @@ -866,6 +875,11 @@ subroutine two_stage_solver_nosnow(nilyr, nslyr, & else ! initially melting +#ifdef GEOSCOUPLED + fsurf_cpl = fsurf_cpl + dfsurfdts_cpl * (Tmlt - Tsf) + flat_cpl = flat_cpl + dflatdts_cpl * (Tmlt - Tsf) +#endif + ! solve the system for melt and no snow Tsf = Tmlt @@ -911,6 +925,11 @@ subroutine two_stage_solver_nosnow(nilyr, nslyr, & fcondtop1 = fcondtop fsurfn1 = fsurfn +#ifdef GEOSCOUPLED + fsurf_cpl = fsurf_cpl0 + flat_cpl = flat_cpl0 +#endif + ! reset the solution to initial values Tsf = Tsf0 zqin = zqin0 @@ -1266,9 +1285,20 @@ subroutine picard_solver(nilyr, nslyr, & zTsn_prev = zTsn zTin_prev = zTin +#ifdef GEOSCOUPLED + dfsurfn_dTsf = dfsurfdts_cpl + dflatn_dTsf = dflatdts_cpl + fsurfn = fsurf_cpl + flatn = flat_cpl + fsurfn = fsurfn + fswsfc + flwoutn = c0 !prevent compiler warning + fsensn = c0 !prevent compiler warning +#endif + ! picard iteration picard: do nit = 1, nit_max +#ifndef GEOSCOUPLED ! surface heat flux call surface_heat_flux(Tsf, fswsfc, & rhoa, flw, & @@ -1284,6 +1314,7 @@ subroutine picard_solver(nilyr, nslyr, & dfsurfn_dTsf, dflwoutn_dTsf, & dfsensn_dTsf, dflatn_dTsf) if (icepack_warnings_aborted(subname)) return +#endif ! tridiagonal solve of new temperatures call solve_heat_conduction(lsnow, lcold, & @@ -1334,6 +1365,11 @@ subroutine picard_solver(nilyr, nslyr, & fadvheat_nit) if (icepack_warnings_aborted(subname)) return +#ifdef GEOSCOUPLED + fsurfn = fsurfn + (Tsf - Tsf_prev)*dfsurfn_dTsf + flatn = flatn + (Tsf - Tsf_prev)*dflatn_dTsf +#endif + if (lconverged) exit Tsf_prev = Tsf @@ -1357,6 +1393,7 @@ subroutine picard_solver(nilyr, nslyr, & dt, nilyr) if (icepack_warnings_aborted(subname)) return +#ifndef GEOSCOUPLED ! final surface heat flux call surface_heat_flux(Tsf, fswsfc, & rhoa, flw, & @@ -1365,6 +1402,7 @@ subroutine picard_solver(nilyr, nslyr, & flwoutn, fsensn, & flatn, fsurfn) if (icepack_warnings_aborted(subname)) return +#endif ! if not converged if (.not. lconverged) then diff --git a/columnphysics/icepack_therm_shared.F90 b/columnphysics/icepack_therm_shared.F90 index 99027173d..29a26c015 100644 --- a/columnphysics/icepack_therm_shared.F90 +++ b/columnphysics/icepack_therm_shared.F90 @@ -36,6 +36,9 @@ module icepack_therm_shared icepack_liquidus_temperature, & icepack_sea_freezing_temperature, & icepack_enthalpy_snow, & +#ifdef GEOSCOUPLED + ismyturn, & +#endif adjust_enthalpy real (kind=dbl_kind), parameter, public :: & @@ -55,8 +58,13 @@ module icepack_therm_shared real (kind=dbl_kind), public :: & dfsurfdts_cpl, & ! dflatdts_cpl, & ! + fsurf_cpl0, & ! + flat_cpl0, & ! fsurf_cpl, & ! flat_cpl ! + + integer(kind=int_kind), public :: & + local_tsk, local_i, local_j, local_blk #endif !======================================================================= @@ -550,6 +558,18 @@ subroutine adjust_enthalpy (nlyr, & end subroutine adjust_enthalpy +#ifdef GEOSCOUPLED + function ismyturn() result(ret) + + logical(kind=log_kind) :: ret + + ret = (local_tsk == 38 .and. local_i == 5 .and. & + local_j == 53 .and. local_blk == 1) + + end function ismyturn +#endif + + !======================================================================= end module icepack_therm_shared diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 9731802bc..63c81e2ba 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -40,8 +40,15 @@ module icepack_therm_vertical #ifdef GEOSCOUPLED use icepack_therm_shared, only: dfsurfdts_cpl, & ! dflatdts_cpl, & ! + flat_cpl0, & ! + fsurf_cpl0, & ! fsurf_cpl, & ! flat_cpl ! + use icepack_therm_shared, only: local_tsk, & ! + local_i, & ! + local_j, & ! + ismyturn, & + local_blk #endif use icepack_therm_bl99, only: temperature_changes use icepack_therm_mushy, only: temperature_changes_salinity @@ -314,6 +321,22 @@ subroutine thermo_vertical (nilyr, nslyr, & dflatdts_cpl = dflatdt_in fsurf_cpl = fsurfn_in flat_cpl = flatn_in + fsurf_cpl0 = fsurfn_in + flat_cpl0 = flatn_in + + if(present(my_tsk)) then + local_tsk = my_tsk + endif + if(present(my_i)) then + local_i = my_i + endif + if(present(my_j)) then + local_j = my_j + endif + if(present(my_blk)) then + local_blk = my_blk + endif + #endif if (calc_Tsfc) then @@ -995,6 +1018,7 @@ subroutine init_vertical_profile(nilyr, nslyr, & call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'zTin=',zTin(k) call icepack_warnings_add(warnstr) + tice_high = .false. else call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//" init_vertical_profile: Starting thermo, zTin > Tmax, layer" ) @@ -2974,6 +2998,14 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & if (icepack_warnings_aborted(subname)) then write(warnstr,*) subname, ' ice: Vertical thermo error, cat ', n call icepack_warnings_add(warnstr) + write(warnstr,*) subname, ' ice: Vertical thermo error, my_tsk ', my_tsk + call icepack_warnings_add(warnstr) + write(warnstr,*) subname, ' ice: Vertical thermo error, my_i ', my_i + call icepack_warnings_add(warnstr) + write(warnstr,*) subname, ' ice: Vertical thermo error, my_j ', my_j + call icepack_warnings_add(warnstr) + write(warnstr,*) subname, ' ice: Vertical thermo error, my_blk ', my_blk + call icepack_warnings_add(warnstr) return endif From 73467f62940bfaa7ce3eaeb53cf93c0b0fbc8416 Mon Sep 17 00:00:00 2001 From: bzhao Date: Wed, 4 Oct 2023 16:16:12 -0400 Subject: [PATCH 23/30] fixed a bug of correcting zqsn(1) prematually when ktherm=2 --- columnphysics/icepack_therm_vertical.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 63c81e2ba..cdd793330 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -1344,8 +1344,10 @@ subroutine thickness_changes (nilyr, nslyr, & if (hstot > puny) then zqsn(1) = (dzs(1) * zqsn(1) & + dhs * zqsnew) / hstot + if (ktherm < 2) then ! avoid roundoff errors - zqsn(1) = min(zqsn(1), -rhos*Lfresh) + zqsn(1) = min(zqsn(1), -rhos*Lfresh) + endif endif #endif dzs(1) = dzs(1) + dhs From 3122eef78a469da4fca663e29afd4883bcb44994 Mon Sep 17 00:00:00 2001 From: bzhao Date: Wed, 13 Dec 2023 19:02:25 -0500 Subject: [PATCH 24/30] first attempt at enabling delta Eddington shortwave scheme in GEOS --- columnphysics/icepack_shortwave.F90 | 76 +++++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) diff --git a/columnphysics/icepack_shortwave.F90 b/columnphysics/icepack_shortwave.F90 index 84378d464..07f67bbae 100644 --- a/columnphysics/icepack_shortwave.F90 +++ b/columnphysics/icepack_shortwave.F90 @@ -79,6 +79,18 @@ module icepack_shortwave real (kind=dbl_kind), parameter :: & exp_argmax = c10 ! maximum argument of exponential +#ifdef GEOSCOUPLED + real (kind=dbl_kind), private :: & + fswthru_uvrdr_out, & + fswthru_uvrdf_out, & + fswthru_pardr_out, & + fswthru_pardf_out, & + druvr_in, & + dfuvr_in, & + drpar_in, & + dfpar_in +#endif + !======================================================================= contains @@ -885,6 +897,10 @@ subroutine run_dEdd(dt, ncat, & fswthrun_vdf, & fswthrun_idr, & fswthrun_idf, & + fswthrun_uvrdr, & + fswthrun_uvrdf, & + fswthrun_pardr, & + fswthrun_pardf, & fswpenln, & Sswabsn, Iswabsn, & albicen, albsnon, & @@ -992,6 +1008,12 @@ subroutine run_dEdd(dt, ncat, & fswthrun_idr, & ! nir dir SW through ice to ocean (W/m^2) fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) + real(kind=dbl_kind), dimension(:), intent(out), optional :: & + fswthrun_uvrdr, & ! uvr dir SW through ice to ocean (W/m^2) + fswthrun_uvrdf, & ! uvr dif SW through ice to ocean (W/m^2) + fswthrun_pardr, & ! par dir SW through ice to ocean (W/m^2) + fswthrun_pardf ! par dif SW through ice to ocean (W/m^2) + real(kind=dbl_kind), dimension(:,:), intent(inout) :: & rsnow , & ! snow grain radius tracer (10^-6 m) Sswabsn , & ! SW radiation absorbed in snow layers (W m-2) @@ -1061,6 +1083,10 @@ subroutine run_dEdd(dt, ncat, & endif ! cosine of the zenith angle + +#ifdef GEOSCOUPLED + ! coszen already computed by GEOS +#else #ifdef CESMCOUPLED call compute_coszen (tlat, tlon, & yday, sec, coszen, & @@ -1068,6 +1094,7 @@ subroutine run_dEdd(dt, ncat, & #else call compute_coszen (tlat, tlon, & yday, sec, coszen) +#endif #endif if (icepack_warnings_aborted(subname)) return @@ -1256,6 +1283,21 @@ subroutine run_dEdd(dt, ncat, & enddo endif +#ifdef GEOSCOUPLED + if(present(fswthrun_uvrdr)) then + fswthrun_uvrdr(n) = fswthru_uvrdr_out + endif + if(present(fswthrun_uvrdf)) then + fswthrun_uvrdf(n) = fswthru_uvrdf_out + endif + if(present(fswthrun_pardr)) then + fswthrun_pardr(n) = fswthru_pardr_out + endif + if(present(fswthrun_pardf)) then + fswthrun_pardf(n) = fswthru_pardf_out + endif +#endif + endif ! aicen > puny enddo ! ncat @@ -1469,6 +1511,12 @@ subroutine shortwave_dEdd (dEdd_algae, & fswthru_vdf = c0 fswthru_idr = c0 fswthru_idf = c0 +#ifdef GEOSCOUPLED + fswthru_uvrdr_out = c0 + fswthru_uvrdf_out = c0 + fswthru_pardr_out = c0 + fswthru_pardf_out = c0 +#endif ! compute fraction of nir down direct to total over all points: fnidr = c0 if( swidr + swidf > puny ) then @@ -1953,6 +2001,14 @@ subroutine compute_dEdd (nilyr, nslyr, klev, klevp, & fthruidr, & ! nir dir shortwave through snow/bare ice/ponded ice to ocean (W/m^2) fthruidf ! nir dif shortwave through snow/bare ice/ponded ice to ocean (W/m^2) +#ifdef GEOSCOUPLED + real (kind=dbl_kind) :: & + fthru_uvrdr, & + fthru_uvrdf, & + fthru_pardr, & + fthru_pardf +#endif + real (kind=dbl_kind), dimension(nslyr) :: & Sabs ! shortwave absorbed in snow layer (W m-2) @@ -3043,6 +3099,13 @@ subroutine compute_dEdd (nilyr, nslyr, klev, klevp, & fthruvdr = fthruvdr + dfdir(klevp)*swdr fthruvdf = fthruvdf + dfdif(klevp)*swdf +#ifdef GEOSCOUPLED + fthru_uvrdr = dfdir(klevp)*druvr_in + fthru_uvrdf = dfdif(klevp)*dfuvr_in + fthru_pardr = dfdir(klevp)*drpar_in + fthru_pardf = dfdif(klevp)*dfpar_in +#endif + ! if snow covered ice, set snow internal absorption; else, Sabs=0 if( srftyp == 1 ) then ki = 0 @@ -3160,6 +3223,12 @@ subroutine compute_dEdd (nilyr, nslyr, klev, klevp, & fswthru_vdf = fswthru_vdf + fthruvdf*fi fswthru_idr = fswthru_idr + fthruidr*fi fswthru_idf = fswthru_idf + fthruidf*fi +#ifdef GEOSCOUPLED + fswthru_uvrdr_out = fswthru_uvrdr_out + fthru_uvrdr*fi + fswthru_uvrdf_out = fswthru_uvrdf_out + fthru_uvrdf*fi + fswthru_pardr_out = fswthru_pardr_out + fthru_pardr*fi + fswthru_pardf_out = fswthru_pardf_out + fthru_pardf*fi +#endif do k = 1, nslyr Sswabs(k) = Sswabs(k) + Sabs(k)*fi @@ -4302,6 +4371,13 @@ subroutine icepack_step_radiation (dt, ncat, & l_fswthrun_pardr = c0 l_fswthrun_pardf = c0 +#ifdef GEOSCOUPLED + druvr_in = l_druvr + dfuvr_in = l_dfuvr + drpar_in = l_drpar + dfpar_in = l_dfpar +#endif + hin = c0 hbri = c0 linitonly = .false. From 9c0b042ee8fa974ff264023065ed7025478a8d11 Mon Sep 17 00:00:00 2001 From: bzhao Date: Tue, 23 Apr 2024 21:09:47 -0400 Subject: [PATCH 25/30] fix errors due to merging --- columnphysics/icepack_flux.F90 | 4 ---- columnphysics/icepack_shortwave.F90 | 1 - columnphysics/icepack_therm_vertical.F90 | 1 + 3 files changed, 1 insertion(+), 5 deletions(-) diff --git a/columnphysics/icepack_flux.F90 b/columnphysics/icepack_flux.F90 index c10381d6b..e57f17eed 100644 --- a/columnphysics/icepack_flux.F90 +++ b/columnphysics/icepack_flux.F90 @@ -93,10 +93,6 @@ subroutine merge_fluxes (aicen, & fsaltn , & ! salt flux to ocean (kg/m2/s) fhocnn , & ! actual ocn/ice heat flx (W/m**2) fswthrun, & ! sw radiation through ice bot (W/m**2) - fswthrun_vdr, & ! vis dir sw radiation through ice bot (W/m**2) - fswthrun_vdf, & ! vis dif sw radiation through ice bot (W/m**2) - fswthrun_idr, & ! nir dir sw radiation through ice bot (W/m**2) - fswthrun_idf, & ! nir dif sw radiation through ice bot (W/m**2) fswthrun_uvrdr, & ! vis dir sw radiation through ice bot (W/m**2) fswthrun_uvrdf, & ! vis dir sw radiation through ice bot (W/m**2) fswthrun_pardr, & ! vis dir sw radiation through ice bot (W/m**2) diff --git a/columnphysics/icepack_shortwave.F90 b/columnphysics/icepack_shortwave.F90 index 597a4c5b5..63d615fd8 100644 --- a/columnphysics/icepack_shortwave.F90 +++ b/columnphysics/icepack_shortwave.F90 @@ -4095,7 +4095,6 @@ subroutine icepack_step_radiation (dt, & endif #endif endif ->>>>>>> main hin = c0 hbri = c0 diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 54570d87f..704648796 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -1337,6 +1337,7 @@ subroutine thickness_changes (nilyr, nslyr, & zqsn(1) = min(zqsn(1), -rhos*Lfresh) endif endif +#endif dzs(1) = dzs(1) + dhs evapn = evapn + dhs*rhos From a65d6204f62808e78f83aea40f19c39f390d4c7a Mon Sep 17 00:00:00 2001 From: bzhao Date: Wed, 5 Jun 2024 15:48:13 -0400 Subject: [PATCH 26/30] remove debugging facility code --- columnphysics/icepack_therm_mushy.F90 | 1 - columnphysics/icepack_therm_shared.F90 | 18 ------- columnphysics/icepack_therm_vertical.F90 | 61 ++---------------------- 3 files changed, 5 insertions(+), 75 deletions(-) diff --git a/columnphysics/icepack_therm_mushy.F90 b/columnphysics/icepack_therm_mushy.F90 index dc294ed38..b24fb151c 100644 --- a/columnphysics/icepack_therm_mushy.F90 +++ b/columnphysics/icepack_therm_mushy.F90 @@ -28,7 +28,6 @@ module icepack_therm_mushy fsurf_cpl0, & ! fsurf_cpl, & ! flat_cpl ! - use icepack_therm_shared, only: ismyturn #endif implicit none diff --git a/columnphysics/icepack_therm_shared.F90 b/columnphysics/icepack_therm_shared.F90 index c23b08ae4..f6e57bf75 100644 --- a/columnphysics/icepack_therm_shared.F90 +++ b/columnphysics/icepack_therm_shared.F90 @@ -35,9 +35,6 @@ module icepack_therm_shared icepack_snow_temperature, & icepack_liquidus_temperature, & icepack_sea_freezing_temperature, & -#ifdef GEOSCOUPLED - ismyturn, & -#endif adjust_enthalpy real (kind=dbl_kind), parameter, public :: & @@ -58,9 +55,6 @@ module icepack_therm_shared flat_cpl0, & ! fsurf_cpl, & ! flat_cpl ! - - integer(kind=int_kind), public :: & - local_tsk, local_i, local_j, local_blk #endif !======================================================================= @@ -590,18 +584,6 @@ subroutine adjust_enthalpy (nlyr, & end subroutine adjust_enthalpy -#ifdef GEOSCOUPLED - function ismyturn() result(ret) - - logical(kind=log_kind) :: ret - - ret = (local_tsk == 38 .and. local_i == 5 .and. & - local_j == 53 .and. local_blk == 1) - - end function ismyturn -#endif - - !======================================================================= end module icepack_therm_shared diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 704648796..910df6aec 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -44,11 +44,6 @@ module icepack_therm_vertical fsurf_cpl0, & ! fsurf_cpl, & ! flat_cpl ! - use icepack_therm_shared, only: local_tsk, & ! - local_i, & ! - local_j, & ! - ismyturn, & - local_blk #endif use icepack_therm_bl99, only: temperature_changes use icepack_therm_mushy, only: temperature_changes_salinity @@ -120,8 +115,7 @@ subroutine thermo_vertical (nilyr, nslyr, & dfsurfdt_in, & dflatdt_in, & yday, dsnow, & - prescribed_ice, & - my_tsk, my_i, my_j, my_blk) + prescribed_ice ) integer (kind=int_kind), intent(in) :: & nilyr , & ! number of ice layers @@ -147,9 +141,6 @@ subroutine thermo_vertical (nilyr, nslyr, & logical (kind=log_kind), intent(in), optional :: & prescribed_ice ! if .true., use prescribed ice instead of computed - integer (kind=int_kind), intent(in), optional :: & - my_tsk, my_i, my_j, my_blk - real (kind=dbl_kind), dimension (:), intent(inout) :: & zqsn , & ! snow layer enthalpy, zqsn < 0 (J m-3) zqin , & ! ice layer enthalpy, zqin < 0 (J m-3) @@ -322,20 +313,6 @@ subroutine thermo_vertical (nilyr, nslyr, & flat_cpl = flatn_in fsurf_cpl0 = fsurfn_in flat_cpl0 = flatn_in - - if(present(my_tsk)) then - local_tsk = my_tsk - endif - if(present(my_i)) then - local_i = my_i - endif - if(present(my_j)) then - local_j = my_j - endif - if(present(my_blk)) then - local_blk = my_blk - endif - #endif if (calc_Tsfc) then @@ -489,11 +466,7 @@ subroutine thermo_vertical (nilyr, nslyr, & fsnow, einit, & einter, efinal, & fcondtopn, fcondbotn, & - fadvocn, fbot, & - my_tsk = my_tsk, & - my_i = my_i, & - my_j = my_j, & - my_blk = my_blk) + fadvocn, fbot ) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- @@ -2022,9 +1995,7 @@ subroutine conservation_check_vthermo(dt, & einit, einter, & efinal, & fcondtopn,fcondbotn, & - fadvocn, fbot, & - my_tsk, my_i, & - my_j, my_blk) + fadvocn, fbot ) real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -2045,9 +2016,6 @@ subroutine conservation_check_vthermo(dt, & efinal , & ! final energy of melting (J m-2) fcondbotn - integer(kind=int_kind), intent(in), optional :: & - my_tsk, my_i, my_j, my_blk - ! local variables real (kind=dbl_kind) :: & @@ -2077,10 +2045,6 @@ subroutine conservation_check_vthermo(dt, & write(warnstr,*) subname, 'Thermo energy conservation error' call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'at task, i, j, iblk' - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, my_tsk, my_i, my_j, my_blk - call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Flux error (W/m^2) =', ferr call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Energy error (J) =', ferr*dt @@ -2312,9 +2276,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & mlt_onset , frz_onset , & dfsurfdt , dflatdt , & yday , prescribed_ice, & - zlvs , my_tsk, & - my_i, my_j, & - my_blk ) + zlvs ) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories @@ -2336,9 +2298,6 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & logical (kind=log_kind), intent(in), optional :: & prescribed_ice ! if .true., use prescribed ice instead of computed - integer (kind=int_kind), intent(in), optional :: & - my_tsk, my_i, my_j, my_blk - real (kind=dbl_kind), intent(inout) :: & aice , & ! sea ice concentration vice , & ! volume per unit area of ice (m) @@ -2904,21 +2863,11 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & dflatdt_in = dflatdt(n), & #endif yday=yday, dsnow=dsnown (n), & - prescribed_ice=prescribed_ice, & - my_tsk = my_tsk, my_i = my_i, & - my_j = my_j, my_blk = my_blk ) + prescribed_ice=prescribed_ice ) if (icepack_warnings_aborted(subname)) then write(warnstr,*) subname, ' ice: Vertical thermo error, cat ', n call icepack_warnings_add(warnstr) - write(warnstr,*) subname, ' ice: Vertical thermo error, my_tsk ', my_tsk - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, ' ice: Vertical thermo error, my_i ', my_i - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, ' ice: Vertical thermo error, my_j ', my_j - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, ' ice: Vertical thermo error, my_blk ', my_blk - call icepack_warnings_add(warnstr) return endif From 86cf3b39b109f58ec5792e6d2bc7e5323e0cf94a Mon Sep 17 00:00:00 2001 From: bzhao Date: Thu, 6 Jun 2024 09:55:07 -0400 Subject: [PATCH 27/30] remove some files --- configuration/scripts/options/set_nml.run15day | 4 ---- configuration/scripts/tests/small_suite.ts | 4 ---- 2 files changed, 8 deletions(-) delete mode 100644 configuration/scripts/options/set_nml.run15day delete mode 100644 configuration/scripts/tests/small_suite.ts diff --git a/configuration/scripts/options/set_nml.run15day b/configuration/scripts/options/set_nml.run15day deleted file mode 100644 index 2ed171930..000000000 --- a/configuration/scripts/options/set_nml.run15day +++ /dev/null @@ -1,4 +0,0 @@ -npt = 360 -dumpfreq = 'd' - - diff --git a/configuration/scripts/tests/small_suite.ts b/configuration/scripts/tests/small_suite.ts deleted file mode 100644 index ef3f2aa8b..000000000 --- a/configuration/scripts/tests/small_suite.ts +++ /dev/null @@ -1,4 +0,0 @@ -# Test Grid PEs Sets BFB-compare -smoke col 1x1 diag1,run1year -smoke col 1x1 debug,run1year - From 7edcd0509c7bf39e996ebe63cf2ef7c5426fb4f0 Mon Sep 17 00:00:00 2001 From: bzhao Date: Thu, 6 Jun 2024 12:38:07 -0400 Subject: [PATCH 28/30] restore to consortium version; not needed to be speific to geos --- configuration/scripts/icepack.test.setup.csh | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/configuration/scripts/icepack.test.setup.csh b/configuration/scripts/icepack.test.setup.csh index 6055f5357..b5d4ee433 100755 --- a/configuration/scripts/icepack.test.setup.csh +++ b/configuration/scripts/icepack.test.setup.csh @@ -60,18 +60,6 @@ cat >> ${jobfile} < ${ICE_SCRIPTS}/tests/baseline.script chmod +x ${jobfile} -if ($ICE_MACHINE == 'discover') then - -cat >! ${subfile} << EOFS -#!/bin/csh -f - -./${jobfile} -echo "\`date\` \${0}: ${ICE_CASENAME} job submitted" >> ${ICE_CASEDIR}/README.case - -EOFS - -else - cat >! ${subfile} << EOFS #!/bin/csh -f @@ -80,7 +68,5 @@ echo "\`date\` \${0}: ${ICE_CASENAME} job submitted" >> ${ICE_CASEDIR}/README.c EOFS -endif - chmod +x ${subfile} exit 0 From f31b5a9feb2eca9e542b15f8829fce3dc124e8e9 Mon Sep 17 00:00:00 2001 From: bzhao Date: Thu, 6 Jun 2024 17:12:05 -0400 Subject: [PATCH 29/30] revert back to consortium version --- icepack.setup | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/icepack.setup b/icepack.setup index ec40e1c70..788a273c8 100755 --- a/icepack.setup +++ b/icepack.setup @@ -570,7 +570,7 @@ foreach envname ( $nenvnames ) if (${dosuite} == 1) then set testname = "${tsdir}/${testname_base}" else - set testname = "${tsdir}/${testname_base}" + set testname = "${testname_base}" endif set case = ${testname} endif From 572ea3d30c4d2c9c6b5d4f6d83168ebce06d4007 Mon Sep 17 00:00:00 2001 From: bzhao Date: Tue, 11 Jun 2024 10:14:47 -0400 Subject: [PATCH 30/30] add/correct some comments to geos specific SW fluxes --- columnphysics/icepack_shortwave.F90 | 64 ++++++++++++++--------------- 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/columnphysics/icepack_shortwave.F90 b/columnphysics/icepack_shortwave.F90 index 63d615fd8..71c9eb72e 100644 --- a/columnphysics/icepack_shortwave.F90 +++ b/columnphysics/icepack_shortwave.F90 @@ -123,14 +123,14 @@ module icepack_shortwave #ifdef GEOSCOUPLED real (kind=dbl_kind), private :: & - fswthru_uvrdr_out, & - fswthru_uvrdf_out, & - fswthru_pardr_out, & - fswthru_pardf_out, & - druvr_in, & - dfuvr_in, & - drpar_in, & - dfpar_in + fswthru_uvrdr_out, & ! penetrative uvr flux, direct + fswthru_uvrdf_out, & ! penetrative uvr flux, diffuse + fswthru_pardr_out, & ! penetrative par flux, direct + fswthru_pardf_out, & ! penetrative par flux, diffuse + druvr_in, & ! uvr flux, direct + dfuvr_in, & ! uvr flux, diffuse + drpar_in, & ! par flux, direct + dfpar_in ! par flux, diffuse #endif ! dEdd tuning parameters, set in namelist ! R_ice ! sea ice tuning parameter; +1 > 1sig increase in albedo @@ -297,10 +297,10 @@ subroutine shortwave_ccsm3 (aicen, vicen, & fswthrun_idf ! nir dif SW through ice to ocean (W m-2) real (kind=dbl_kind), dimension(:), intent(out), optional :: & - fswthru_uvrdr , & ! vis dir SW through ice to ocean (W/m^2) - fswthru_uvrdf , & ! vis dif SW through ice to ocean (W/m^2) - fswthru_pardr , & ! nir dir SW through ice to ocean (W/m^2) - fswthru_pardf ! nir dif SW through ice to ocean (W/m^2) + fswthru_uvrdr , & ! vis dir uvr SW through ice to ocean (W m-2) + fswthru_uvrdf , & ! vis dif uvr SW through ice to ocean (W m-2) + fswthru_pardr , & ! vis dir par SW through ice to ocean (W m-2) + fswthru_pardf ! vis dif par SW through ice to ocean (W m-2) real (kind=dbl_kind), intent(inout) :: & coszen ! cosine(zenith angle) @@ -335,16 +335,16 @@ subroutine shortwave_ccsm3 (aicen, vicen, & l_fswthru_idf ! nir dif SW through ice to ocean (W m-2) real (kind=dbl_kind) :: & - l_fswthru_uvrdr , & ! vis uvr dir SW through ice to ocean (W/m^2) - l_fswthru_uvrdf , & ! vis uvr dif SW through ice to ocean (W/m^2) - l_fswthru_pardr , & ! vis par dir SW through ice to ocean (W/m^2) - l_fswthru_pardf ! vis par dif SW through ice to ocean (W/m^2) + l_fswthru_uvrdr , & ! vis uvr dir SW through ice to ocean (W m-2) + l_fswthru_uvrdf , & ! vis uvr dif SW through ice to ocean (W m-2) + l_fswthru_pardr , & ! vis par dir SW through ice to ocean (W m-2) + l_fswthru_pardf ! vis par dif SW through ice to ocean (W m-2) real(kind=dbl_kind) :: & - l_druvr, & ! - l_dfuvr, & ! - l_drpar, & ! - l_dfpar ! + l_druvr, & ! uvr flux, direct (W m-2) + l_dfuvr, & ! uvr flux, diffuse (W m-2) + l_drpar, & ! par flux, direct (W m-2) + l_dfpar ! par flux, diffuse (W m-2) character(len=*),parameter :: subname='(shortwave_ccsm3)' @@ -785,10 +785,10 @@ subroutine absorbed_solar (aicen, & alidfns ! near-ir, diffuse albedo, snow real (kind=dbl_kind), intent(in) :: & - druvr , & ! - dfuvr , & ! - drpar , & ! - dfpar ! + druvr , & ! uvr flux, direct (W m-2) + dfuvr , & ! uvr flux, diffuse (W m-2) + drpar , & ! par flux, direct (W m-2) + dfpar ! par flux, diffuse (W m-2) real (kind=dbl_kind), intent(out):: & fswsfc , & ! SW absorbed at ice/snow surface (W m-2) @@ -802,10 +802,10 @@ subroutine absorbed_solar (aicen, & fswthru_idf ! nir dif SW through ice to ocean (W m-2) real (kind=dbl_kind), intent(out) :: & - fswthru_uvrdr , & ! vis dir SW through ice to ocean (W/m^2) - fswthru_uvrdf , & ! vis dif SW through ice to ocean (W/m^2) - fswthru_pardr , & ! nir dir SW through ice to ocean (W/m^2) - fswthru_pardf ! nir dif SW through ice to ocean (W/m^2) + fswthru_uvrdr , & ! vis dir uvr SW through ice to ocean (W m-2) + fswthru_uvrdf , & ! vis dif uvr SW through ice to ocean (W m-2) + fswthru_pardr , & ! vis dir par SW through ice to ocean (W m-2) + fswthru_pardf ! vis dif par SW through ice to ocean (W m-2) real (kind=dbl_kind), dimension (:), intent(out) :: & Iswabs , & ! SW absorbed in particular layer (W m-2) @@ -3986,10 +3986,10 @@ subroutine icepack_step_radiation (dt, & fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) real (kind=dbl_kind), dimension(:), intent(inout), optional :: & - fswthrun_uvrdr , & ! vis dir SW through ice to ocean (W/m^2) - fswthrun_uvrdf , & ! vis dif SW through ice to ocean (W/m^2) - fswthrun_pardr , & ! nir dir SW through ice to ocean (W/m^2) - fswthrun_pardf ! nir dif SW through ice to ocean (W/m^2) + fswthrun_uvrdr , & ! vis dir uvr SW through ice to ocean (W/m^2) + fswthrun_uvrdf , & ! vis dif uvr SW through ice to ocean (W/m^2) + fswthrun_pardr , & ! vis dir par SW through ice to ocean (W/m^2) + fswthrun_pardf ! vis dif par SW through ice to ocean (W/m^2) real (kind=dbl_kind), intent(in), optional :: & druvr , & !