From b51c5b314d3af2d6ee01aa2cd7545c25dd8eccb0 Mon Sep 17 00:00:00 2001 From: Ben Green Date: Mon, 26 Apr 2021 18:45:22 +0000 Subject: [PATCH 1/2] Removing GSD_SURFACE_FLUXES_BUGFIX and replacing with flag thsfc_loc --- physics/GFS_surface_composites.F90 | 24 ++++++---- physics/GFS_surface_composites.meta | 8 ++++ physics/sfc_diag.f | 14 +++--- physics/sfc_diag.meta | 8 ++++ physics/sfc_diff.f | 71 +++++++++++++++++++++-------- physics/sfc_diff.meta | 8 ++++ physics/sfc_nst.f | 41 ++++++++++------- physics/sfc_nst.meta | 8 ++++ physics/sfc_sice.f | 46 +++++++++++-------- physics/sfc_sice.meta | 8 ++++ 10 files changed, 164 insertions(+), 72 deletions(-) diff --git a/physics/GFS_surface_composites.F90 b/physics/GFS_surface_composites.F90 index d5dc67f54..81d9ebf60 100644 --- a/physics/GFS_surface_composites.F90 +++ b/physics/GFS_surface_composites.F90 @@ -373,7 +373,7 @@ end subroutine GFS_surface_composites_post_finalize !! \htmlinclude GFS_surface_composites_post_run.html !! subroutine GFS_surface_composites_post_run ( & - im, kice, km, cplflx, cplwav2atm, frac_grid, flag_cice, islmsk, dry, wet, icy, wind, t1, q1, prsl1, & + im, kice, km, cplflx, cplwav2atm, frac_grid, flag_cice, thsfc_loc, islmsk, dry, wet, icy, wind, t1, q1, prsl1, & rd, rvrdm1, landfrac, lakefrac, oceanfrac, zorl, zorlo, zorll, zorli, & cd, cd_wat, cd_lnd, cd_ice, cdq, cdq_wat, cdq_lnd, cdq_ice, rb, rb_wat, rb_lnd, rb_ice, stress, stress_wat, stress_lnd, & stress_ice, ffmm, ffmm_wat, ffmm_lnd, ffmm_ice, ffhh, ffhh_wat, ffhh_lnd, ffhh_ice, uustar, uustar_wat, uustar_lnd, & @@ -410,6 +410,7 @@ subroutine GFS_surface_composites_post_run ( real(kind=kind_phys), dimension(im, km), intent(inout) :: stc ! Additional data needed for calling "stability" + logical, intent(in ) :: thsfc_loc real(kind=kind_phys), intent(in ) :: grav real(kind=kind_phys), dimension(:), intent(in ) :: prslki, z1, ztmax_wat, ztmax_lnd, ztmax_ice @@ -420,7 +421,7 @@ subroutine GFS_surface_composites_post_run ( integer :: i, k real(kind=kind_phys) :: txl, txi, txo, wfrac, q0, rho ! For calling "stability" - real(kind=kind_phys) :: tsurf, virtfac, thv1, tvs, z0max, ztmax + real(kind=kind_phys) :: tsurf, virtfac, tv1, thv1, tvs, z0max, ztmax ! Initialize CCPP error handling variables errmsg = '' @@ -475,24 +476,27 @@ subroutine GFS_surface_composites_post_run ( q0 = max( q1(i), qmin ) virtfac = one + rvrdm1 * q0 -#ifdef GSD_SURFACE_FLUXES_BUGFIX - thv1 = t1(i) / prslk1(i) * virtfac ! Theta-v at lowest level - tvs = half * (tsfc(i)+tsurf)/prsik1(i) * virtfac -#else - thv1 = t1(i) * prslki(i) * virtfac ! Theta-v at lowest level - tvs = half * (tsfc(i)+tsurf) * virtfac -#endif + tv1 = t1(i) * virtfac + + if(thsfc_loc) then ! Use local potential temperature + thv1 = t1(i) * prslki(i) * virtfac ! Theta-v at lowest level + tvs = half * (tsfc(i)+tsurf) * virtfac + else ! Use potential temperature referenced to 1000 hPa + thv1 = t1(i) / prslk1(i) * virtfac ! Theta-v at lowest level + tvs = half * (tsfc(i)+tsurf)/prsik1(i) * virtfac + endif zorl(i) = exp(txl*log(zorll(i)) + txi*log(zorli(i)) + txo*log(zorlo(i))) z0max = 0.01_kind_phys * zorl(i) ztmax = exp(txl*log(ztmax_lnd(i)) + txi*log(ztmax_ice(i)) + txo*log(ztmax_wat(i))) call stability(z1(i), snowd(i), thv1, wind(i), z0max, ztmax, tvs, grav, & ! inputs + tv1, thsfc_loc, & ! inputs rb(i), ffmm(i), ffhh(i), fm10(i), fh2(i), cd(i), cdq(i), & ! outputs stress(i), uustar(i)) ! BWG, 2021/02/25: cmm=cd*wind, chh=cdq*wind, so use composite cd, cdq - rho = prsl1(i) / (rd*t1(i)*(one + rvrdm1*q0)) + rho = prsl1(i) / (rd*t1(i)*virtfac) cmm(i) = cd(i)*wind(i) !txl*cmm_lnd(i) + txi*cmm_ice(i) + txo*cmm_wat(i) chh(i) = rho*cdq(i)*wind(i) !txl*chh_lnd(i) + txi*chh_ice(i) + txo*chh_wat(i) diff --git a/physics/GFS_surface_composites.meta b/physics/GFS_surface_composites.meta index 34766e9cb..1ad173852 100644 --- a/physics/GFS_surface_composites.meta +++ b/physics/GFS_surface_composites.meta @@ -71,6 +71,14 @@ type = logical intent = in optional = F +[thsfc_loc] + standard_name = flag_for_reference_pressure_theta + long_name = flag for reference pressure in theta calculation + units = flag + dimensions = () + type = logical + intent = in + optional = F [cplflx] standard_name = flag_for_flux_coupling long_name = flag controlling cplflx collection (default off) diff --git a/physics/sfc_diag.f b/physics/sfc_diag.f index b78c9b2f7..ceeaf1be8 100644 --- a/physics/sfc_diag.f +++ b/physics/sfc_diag.f @@ -23,7 +23,7 @@ end subroutine sfc_diag_finalize !! @{ subroutine sfc_diag_run & & (im,grav,cp,eps,epsm1,ps,u1,v1,t1,q1,prslki, & - & evap,fm,fh,fm10,fh2,tskin,qsurf, & + & evap,fm,fh,fm10,fh2,tskin,qsurf,thsfc_loc, & & f10m,u10m,v10m,t2m,q2m,errmsg,errflg & & ) ! @@ -32,6 +32,7 @@ subroutine sfc_diag_run & implicit none ! integer, intent(in) :: im + logical, intent(in) :: thsfc_loc ! Flag for reference pot. temp. real(kind=kind_phys), intent(in) :: grav,cp,eps,epsm1 real(kind=kind_phys), dimension(im), intent(in) :: & & ps, u1, v1, t1, q1, tskin, & @@ -74,11 +75,12 @@ subroutine sfc_diag_run & ! t2m(i) = t2m(i) * sig2k wrk = 1.0 - fhi -#ifdef GSD_SURFACE_FLUXES_BUGFIX - t2m(i) = tskin(i)*wrk + t1(i)*fhi - (grav+grav)/cp -#else - t2m(i) = tskin(i)*wrk + t1(i)*prslki(i)*fhi - (grav+grav)/cp -#endif + + if(thsfc_loc) then ! Use local potential temperature + t2m(i) = tskin(i)*wrk + t1(i)*prslki(i)*fhi - (grav+grav)/cp + else ! Use potential temperature referenced to 1000 hPa + t2m(i) = tskin(i)*wrk + t1(i)*fhi - (grav+grav)/cp + endif if(evap(i) >= 0.) then ! for evaporation>0, use inferred qsurf to deduce q2m q2m(i) = qsurf(i)*wrk + max(qmin,q1(i))*fhi diff --git a/physics/sfc_diag.meta b/physics/sfc_diag.meta index deebf23df..9c1e72433 100644 --- a/physics/sfc_diag.meta +++ b/physics/sfc_diag.meta @@ -168,6 +168,14 @@ kind = kind_phys intent = in optional = F +[thsfc_loc] + standard_name = flag_for_reference_pressure_theta + long_name = flag for reference pressure in theta calculation + units = flag + dimensions = () + type = logical + intent = in + optional = F [f10m] standard_name = ratio_of_wind_at_lowest_model_layer_and_wind_at_10m long_name = ratio of fm10 and fm diff --git a/physics/sfc_diff.f b/physics/sfc_diff.f index 669262982..12441e23a 100644 --- a/physics/sfc_diff.f +++ b/physics/sfc_diff.f @@ -69,6 +69,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) & flag_iter,redrag, & !intent(in) & u10m,v10m,sfc_z0_type, & !hafs,z0 type !intent(in) & wet,dry,icy, & !intent(in) + & thsfc_loc, & !intent(in) & tskin_wat, tskin_lnd, tskin_ice, & !intent(in) & tsurf_wat, tsurf_lnd, tsurf_ice, & !intent(in) & snwdph_wat,snwdph_lnd,snwdph_ice, & !intent(in) @@ -97,6 +98,8 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) logical, intent(in) :: redrag ! reduced drag coeff. flag for high wind over sea (j.han) logical, dimension(im), intent(in) :: flag_iter, wet, dry, icy + logical, intent(in) :: thsfc_loc ! Flag for reference pressure in theta calculation + real(kind=kind_phys), dimension(im), intent(in) :: u10m,v10m real(kind=kind_phys), intent(in) :: rvrdm1, eps, epsm1, grav real(kind=kind_phys), dimension(im), intent(in) :: & @@ -133,6 +136,9 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) real(kind=kind_phys) :: rat, thv1, restar, wind10m, & czilc, tem1, tem2, virtfac ! + + real(kind=kind_phys) :: tv1 + real(kind=kind_phys) :: tvs, z0, z0max ! real(kind=kind_phys), parameter :: @@ -178,18 +184,26 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) ztmax_wat(i) = 1. ! log(1) = 0 virtfac = one + rvrdm1 * max(q1(i),qmin) - thv1 = t1(i) * prslki(i) * virtfac + + tv1 = t1(i) * virtfac ! Virtual temperature in middle of lowest layer + if(thsfc_loc) then ! Use local potential temperature + thv1 = t1(i) * prslki(i) * virtfac + else ! Use potential temperature reference to 1000 hPa + thv1 = t1(i) / prslk1(i) * virtfac + endif ! compute stability dependent exchange coefficients ! this portion of the code is presently suppressed ! if (dry(i)) then ! Some land -#ifdef GSD_SURFACE_FLUXES_BUGFIX - tvs = half * (tsurf_lnd(i)+tskin_lnd(i))/prsik1(i) - & * virtfac -#else - tvs = half * (tsurf_lnd(i)+tskin_lnd(i)) * virtfac -#endif + + if(thsfc_loc) then ! Use local potential temperature + tvs = half * (tsurf_lnd(i)+tskin_lnd(i)) * virtfac + else ! Use potential temperature referenced to 1000 hPa + tvs = half * (tsurf_lnd(i)+tskin_lnd(i))/prsik1(i) + & * virtfac + endif + z0max = max(zmin, min(0.01_kp * z0rl_lnd(i), z1(i))) !** xubin's new z0 over land tem1 = one - shdmax(i) @@ -253,14 +267,21 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) call stability ! --- inputs: & (z1(i), snwdph_lnd(i), thv1, wind(i), - & z0max, ztmax_lnd(i), tvs, grav, + & z0max, ztmax_lnd(i), tvs, grav, tv1, thsfc_loc, ! --- outputs: & rb_lnd(i), fm_lnd(i), fh_lnd(i), fm10_lnd(i), fh2_lnd(i), & cm_lnd(i), ch_lnd(i), stress_lnd(i), ustar_lnd(i)) endif ! Dry points if (icy(i)) then ! Some ice - tvs = half * (tsurf_ice(i)+tskin_ice(i)) * virtfac + + if(thsfc_loc) then ! Use local potential temperature + tvs = half * (tsurf_ice(i)+tskin_ice(i)) * virtfac + else ! Use potential temperature referenced to 1000 hPa + tvs = half * (tsurf_ice(i)+tskin_ice(i))/prsik1(i) + & * virtfac + endif + z0max = max(zmin, min(0.01_kp * z0rl_ice(i), z1(i))) !** xubin's new z0 over land and sea ice tem1 = one - shdmax(i) @@ -288,7 +309,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) call stability ! --- inputs: & (z1(i), snwdph_ice(i), thv1, wind(i), - & z0max, ztmax_ice(i), tvs, grav, + & z0max, ztmax_ice(i), tvs, grav, tv1, thsfc_loc, ! --- outputs: & rb_ice(i), fm_ice(i), fh_ice(i), fm10_ice(i), fh2_ice(i), & cm_ice(i), ch_ice(i), stress_ice(i), ustar_ice(i)) @@ -298,7 +319,14 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) ! the stuff now put into "stability" if (wet(i)) then ! Some open ocean - tvs = half * (tsurf_wat(i)+tskin_wat(i)) * virtfac + + if(thsfc_loc) then ! Use local potential temperature + tvs = half * (tsurf_wat(i)+tskin_wat(i)) * virtfac + else + tvs = half * (tsurf_wat(i)+tskin_wat(i))/prsik1(i) + & * virtfac + endif + z0 = 0.01_kp * z0rl_wat(i) z0max = max(zmin, min(z0,z1(i))) ustar_wat(i) = sqrt(grav * z0 / charnock) @@ -332,7 +360,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) call stability ! --- inputs: & (z1(i), snwdph_wat(i), thv1, wind(i), - & z0max, ztmax_wat(i), tvs, grav, + & z0max, ztmax_wat(i), tvs, grav, tv1, thsfc_loc, ! --- outputs: & rb_wat(i), fm_wat(i), fh_wat(i), fm10_wat(i), fh2_wat(i), & cm_wat(i), ch_wat(i), stress_wat(i), ustar_wat(i)) @@ -392,6 +420,7 @@ end subroutine sfc_diff_run subroutine stability & ! --- inputs: & ( z1, snwdph, thv1, wind, z0max, ztmax, tvs, grav, & + & tv1, thsfc_loc, & ! --- outputs: & rb, fm, fh, fm10, fh2, cm, ch, stress, ustar) !----- @@ -400,6 +429,8 @@ subroutine stability & ! --- inputs: real(kind=kind_phys), intent(in) :: & & z1, snwdph, thv1, wind, z0max, ztmax, tvs, grav + real(kind=kind_phys), intent(in) :: tv1 + logical, intent(in) :: thsfc_loc ! --- outputs: real(kind=kind_phys), intent(out) :: & @@ -435,13 +466,15 @@ subroutine stability & dtv = thv1 - tvs adtv = max(abs(dtv),0.001_kp) dtv = sign(1.,dtv) * adtv -#ifdef GSD_SURFACE_FLUXES_BUGFIX - rb = max(-5000.0_kp, grav * dtv * z1 - & / (thv1 * wind * wind)) -#else - rb = max(-5000.0_kp, (grav+grav) * dtv * z1 - & / ((thv1 + tvs) * wind * wind)) -#endif + + if(thsfc_loc) then ! Use local potential temperature + rb = max(-5000.0_kp, (grav+grav) * dtv * z1 + & / ((thv1 + tvs) * wind * wind)) + else ! Use potential temperature referenced to 1000 hPa + rb = max(-5000.0_kp, grav * dtv * z1 + & / (tv1 * wind * wind)) + endif + tem1 = one / z0max tem2 = one / ztmax fm = log((z0max+z1) * tem1) diff --git a/physics/sfc_diff.meta b/physics/sfc_diff.meta index 63935ac11..17a30f28c 100644 --- a/physics/sfc_diff.meta +++ b/physics/sfc_diff.meta @@ -250,6 +250,14 @@ type = logical intent = in optional = F +[thsfc_loc] + standard_name = flag_for_reference_pressure_theta + long_name = flag for reference pressure in theta calculation + units = flag + dimensions = () + type = logical + intent = in + optional = F [tskin_wat] standard_name = surface_skin_temperature_over_water_interstitial long_name = surface skin temperature over water (temporary use as interstitial) diff --git a/physics/sfc_nst.f b/physics/sfc_nst.f index 517aa7ff0..99aab7dd0 100644 --- a/physics/sfc_nst.f +++ b/physics/sfc_nst.f @@ -32,7 +32,7 @@ subroutine sfc_nst_run & & sinlat, stress, & & sfcemis, dlwflx, sfcnsw, rain, timestep, kdt, solhr,xcosz, & & wind, flag_iter, flag_guess, nstf_name1, nstf_name4, & - & nstf_name5, lprnt, ipr, & + & nstf_name5, lprnt, ipr, thsfc_loc, & & tskin, tsurf, xt, xs, xu, xv, xz, zm, xtts, xzts, dt_cool, & ! --- input/output: & z_c, c_0, c_d, w_0, w_d, d_conv, ifd, qrain, & & qsurf, gflux, cmm, chh, evap, hflx, ep, errmsg, errflg & ! --- outputs: @@ -50,7 +50,7 @@ subroutine sfc_nst_run & ! prsl1, prslki, wet, use_flake, xlon, sinlat, stress, ! ! sfcemis, dlwflx, sfcnsw, rain, timestep, kdt,solhr,xcosz, ! ! wind, flag_iter, flag_guess, nstf_name1, nstf_name4, ! -! nstf_name5, lprnt, ipr, ! +! nstf_name5, lprnt, ipr, thsfc_loc, ! ! input/outputs: ! ! tskin, tsurf, xt, xs, xu, xv, xz, zm, xtts, xzts, dt_cool, ! ! z_c, c_0, c_d, w_0, w_d, d_conv, ifd, qrain, ! @@ -123,6 +123,7 @@ subroutine sfc_nst_run & ! nstf_name5 : zsea2 in mm 1 ! ! lprnt - logical, control flag for check print out 1 ! ! ipr - integer, grid index for check print out 1 ! +! thsfc_loc- logical, flag for reference pressure in theta 1 ! ! ! ! input/outputs: ! li added for oceanic components @@ -199,6 +200,7 @@ subroutine sfc_nst_run & & use_flake ! &, icy logical, intent(in) :: lprnt + logical, intent(in) :: thsfc_loc ! --- input/outputs: ! control variables of dtl system (5+2) and sl (2) and coefficients for d(tz)/d(ts) calculation @@ -297,11 +299,13 @@ subroutine sfc_nst_run & wndmag(i) = sqrt(u1(i)*u1(i) + v1(i)*v1(i)) q0(i) = max(q1(i), 1.0e-8_kp) -#ifdef GSD_SURFACE_FLUXES_BUGFIX - theta1(i) = t1(i) / prslk1(i) ! potential temperature at the middle of lowest model layer -#else - theta1(i) = t1(i) * prslki(i) -#endif + + if(thsfc_loc) then ! Use local potential temperature + theta1(i) = t1(i) * prslki(i) + else ! Use potential temperature referenced to 1000 hPa + theta1(i) = t1(i) / prslk1(i) ! potential temperature at the middle of lowest model layer + endif + tv1(i) = t1(i) * (one + rvrdm1*q0(i)) rho_a(i) = prsl1(i) / (rd*tv1(i)) qss(i) = fpvs(tsurf(i)) ! pa @@ -322,11 +326,12 @@ subroutine sfc_nst_run & ! at previous time step evap(i) = elocp * rch(i) * (qss(i) - q0(i)) qsurf(i) = qss(i) -#ifdef GSD_SURFACE_FLUXES_BUGFIX - hflx(i) = rch(i) * (tsurf(i)/prsik1(i) - theta1(i)) -#else - hflx(i) = rch(i) * (tsurf(i) - theta1(i)) -#endif + + if(thsfc_loc) then ! Use local potential temperature + hflx(i) = rch(i) * (tsurf(i) - theta1(i)) + else ! Use potential temperature referenced to 1000 hPa + hflx(i) = rch(i) * (tsurf(i)/prsik1(i) - theta1(i)) + endif ! if (lprnt .and. i == ipr) print *,' tskin=',tskin(i),' theta1=', ! & theta1(i),' hflx=',hflx(i),' t1=',t1(i),'prslki=',prslki(i) @@ -621,11 +626,13 @@ subroutine sfc_nst_run & qss(i) = eps*qss(i) / (ps(i) + epsm1*qss(i)) qsurf(i) = qss(i) evap(i) = elocp*rch(i) * (qss(i) - q0(i)) -#ifdef GSD_SURFACE_FLUXES_BUGFIX - hflx(i) = rch(i) * (tskin(i)/prsik1(i) - theta1(i)) -#else - hflx(i) = rch(i) * (tskin(i) - theta1(i)) -#endif + + if(thsfc_loc) then ! Use local potential temperature + hflx(i) = rch(i) * (tskin(i) - theta1(i)) + else ! Use potential temperature referenced to 1000 hPa + hflx(i) = rch(i) * (tskin(i)/prsik1(i) - theta1(i)) + endif + endif enddo endif ! if ( nstf_name1 > 1 ) then diff --git a/physics/sfc_nst.meta b/physics/sfc_nst.meta index a29f10f90..dc0056aeb 100644 --- a/physics/sfc_nst.meta +++ b/physics/sfc_nst.meta @@ -410,6 +410,14 @@ type = integer intent = in optional = F +[thsfc_loc] + standard_name = flag_for_reference_pressure_theta + long_name = flag for reference pressure in theta calculation + units = flag + dimensions = () + type = logical + intent = in + optional = F [tskin] standard_name = surface_skin_temperature_for_nsst long_name = ocean surface skin temperature diff --git a/physics/sfc_sice.f b/physics/sfc_sice.f index 081bbf48e..7b40c9d25 100644 --- a/physics/sfc_sice.f +++ b/physics/sfc_sice.f @@ -44,7 +44,7 @@ subroutine sfc_sice_run & & t0c, rd, ps, t1, q1, delt, & & sfcemis, dlwflx, sfcnsw, sfcdsw, srflag, & & cm, ch, prsl1, prslki, prsik1, prslk1, wind, & - & flag_iter, lprnt, ipr, & + & flag_iter, lprnt, ipr, thsfc_loc, & & hice, fice, tice, weasd, tskin, tprcp, tiice, ep, & ! --- input/outputs: & snwdph, qsurf, snowmt, gflux, cmm, chh, evap, hflx, & ! & frac_grid, icy, islmsk_cice, & @@ -110,6 +110,7 @@ subroutine sfc_sice_run & ! islimsk - integer, sea/land/ice mask (=0/1/2) im ! ! wind - real, im ! ! flag_iter- logical, im ! +! thsfc_loc- logical, reference pressure for potential temp im ! ! ! ! input/outputs: ! ! hice - real, sea-ice thickness im ! @@ -152,6 +153,7 @@ subroutine sfc_sice_run & ! --- inputs: integer, intent(in) :: im, kice, ipr logical, intent(in) :: lprnt + logical, intent(in) :: thsfc_loc logical, intent(in) :: frac_grid real (kind=kind_phys), intent(in) :: sbc, hvap, tgice, cp, eps, & @@ -276,11 +278,13 @@ subroutine sfc_sice_run & q0 = max(q1(i), qmin) ! tsurf(i) = tskin(i) -#ifdef GSD_SURFACE_FLUXES_BUGFIX - theta1(i) = t1(i) / prslk1(i) ! potential temperature in middle of lowest atm. layer -#else - theta1(i) = t1(i) * prslki(i) -#endif + + if(thsfc_loc) then ! Use local potential temperature + theta1(i) = t1(i) * prslki(i) + else ! Use potential temperature referenced to 1000 hPa + theta1(i) = t1(i) / prslk1(i) ! potential temperature in middle of lowest atm. layer + endif + rho(i) = prsl1(i) / (rd*t1(i)*(one+rvrdm1*q0)) qs1 = fpvs(t1(i)) qs1 = max(eps*qs1 / (prsl1(i) + epsm1*qs1), qmin) @@ -333,13 +337,14 @@ subroutine sfc_sice_run & !> - Calculate net non-solar and upir heat flux @ ice surface \a hfi. -#ifdef GSD_SURFACE_FLUXES_BUGFIX - hfi(i) = -dlwflx(i) + sfcemis(i)*sbc*t14 + evapi(i) & - & + rch(i)*(tice(i)/prsik1(i) - theta1(i)) -#else - hfi(i) = -dlwflx(i) + sfcemis(i)*sbc*t14 + evapi(i) & - & + rch(i)*(tice(i) - theta1(i)) -#endif + if(thsfc_loc) then ! Use local potential temperature + hfi(i) = -dlwflx(i) + sfcemis(i)*sbc*t14 + evapi(i) & + & + rch(i)*(tice(i) - theta1(i)) + else ! Use potential temperature referenced to 1000 hPa + hfi(i) = -dlwflx(i) + sfcemis(i)*sbc*t14 + evapi(i) & + & + rch(i)*(tice(i)/prsik1(i) - theta1(i)) + endif + !> - Calculate heat flux derivative at surface \a hfd. hfd(i) = 4.0_kind_phys*sfcemis(i)*sbc*tice(i)*t12 & & + (one + elocp*eps*hvap*qs1/(rd*t12)) * rch(i) @@ -415,13 +420,14 @@ subroutine sfc_sice_run & if (flag(i)) then ! --- ... calculate sensible heat flux (& evap over sea ice) -#ifdef GSD_SURFACE_FLUXES_BUGFIX - hflxi = rch(i) * (tice(i)/prsik1(i) - theta1(i)) - hflxw = rch(i) * (tgice / prsik1(i) - theta1(i)) -#else - hflxi = rch(i) * (tice(i) - theta1(i)) - hflxw = rch(i) * (tgice - theta1(i)) -#endif + if(thsfc_loc) then ! Use local potential temperature + hflxi = rch(i) * (tice(i) - theta1(i)) + hflxw = rch(i) * (tgice - theta1(i)) + else ! Use potential temperature referenced to 1000 hPa + hflxi = rch(i) * (tice(i)/prsik1(i) - theta1(i)) + hflxw = rch(i) * (tgice / prsik1(i) - theta1(i)) + endif + hflx(i) = fice(i)*hflxi + ffw(i)*hflxw evap(i) = fice(i)*evapi(i) + ffw(i)*evapw(i) ! diff --git a/physics/sfc_sice.meta b/physics/sfc_sice.meta index 4ce931bac..b256d54ff 100644 --- a/physics/sfc_sice.meta +++ b/physics/sfc_sice.meta @@ -281,6 +281,14 @@ type = integer intent = in optional = F +[thsfc_loc] + standard_name = flag_for_reference_pressure_theta + long_name = flag for reference pressure in theta calculation + units = flag + dimensions = () + type = logical + intent = in + optional = F [hice] standard_name = sea_ice_thickness long_name = sea-ice thickness From b6d0ede38a4d0c231999b0cc3d9241452e5b6fa3 Mon Sep 17 00:00:00 2001 From: Ben Green Date: Mon, 26 Apr 2021 18:52:56 +0000 Subject: [PATCH 2/2] Bugfix --- physics/GFS_surface_composites.F90 | 2 +- physics/GFS_surface_composites.meta | 16 ++++++++-------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/physics/GFS_surface_composites.F90 b/physics/GFS_surface_composites.F90 index 81d9ebf60..70515cf9b 100644 --- a/physics/GFS_surface_composites.F90 +++ b/physics/GFS_surface_composites.F90 @@ -476,8 +476,8 @@ subroutine GFS_surface_composites_post_run ( q0 = max( q1(i), qmin ) virtfac = one + rvrdm1 * q0 - tv1 = t1(i) * virtfac + tv1 = t1(i) * virtfac ! Virtual temperature in middle of lowest layer if(thsfc_loc) then ! Use local potential temperature thv1 = t1(i) * prslki(i) * virtfac ! Theta-v at lowest level tvs = half * (tsfc(i)+tsurf) * virtfac diff --git a/physics/GFS_surface_composites.meta b/physics/GFS_surface_composites.meta index 1ad173852..10e19ec4c 100644 --- a/physics/GFS_surface_composites.meta +++ b/physics/GFS_surface_composites.meta @@ -71,14 +71,6 @@ type = logical intent = in optional = F -[thsfc_loc] - standard_name = flag_for_reference_pressure_theta - long_name = flag for reference pressure in theta calculation - units = flag - dimensions = () - type = logical - intent = in - optional = F [cplflx] standard_name = flag_for_flux_coupling long_name = flag controlling cplflx collection (default off) @@ -916,6 +908,14 @@ type = logical intent = in optional = F +[thsfc_loc] + standard_name = flag_for_reference_pressure_theta + long_name = flag for reference pressure in theta calculation + units = flag + dimensions = () + type = logical + intent = in + optional = F [islmsk] standard_name = sea_land_ice_mask long_name = sea/land/ice mask (=0/1/2)