From 63e5196236317600052e1ec48dc59dc0362cef22 Mon Sep 17 00:00:00 2001 From: ShanSunNOAA Date: Tue, 30 Jul 2024 03:24:28 +0000 Subject: [PATCH] Creating skinsst_20240726: have SkinSST operate on ocean points, while NSST on lake points --- physics/SFC_Layer/UFS/sfc_nst_lake.f90 | 684 +++++++++++++++++++ physics/SFC_Layer/UFS/sfc_nst_lake.meta | 668 ++++++++++++++++++ physics/SFC_Layer/UFS/sfc_nst_lake_post.f90 | 89 +++ physics/SFC_Layer/UFS/sfc_nst_lake_post.meta | 205 ++++++ physics/SFC_Layer/UFS/skinsst.f90 | 578 ++++++++++++++++ physics/SFC_Layer/UFS/skinsst.meta | 368 ++++++++++ 6 files changed, 2592 insertions(+) create mode 100644 physics/SFC_Layer/UFS/sfc_nst_lake.f90 create mode 100644 physics/SFC_Layer/UFS/sfc_nst_lake.meta create mode 100644 physics/SFC_Layer/UFS/sfc_nst_lake_post.f90 create mode 100644 physics/SFC_Layer/UFS/sfc_nst_lake_post.meta create mode 100644 physics/SFC_Layer/UFS/skinsst.f90 create mode 100644 physics/SFC_Layer/UFS/skinsst.meta diff --git a/physics/SFC_Layer/UFS/sfc_nst_lake.f90 b/physics/SFC_Layer/UFS/sfc_nst_lake.f90 new file mode 100644 index 000000000..945f7c064 --- /dev/null +++ b/physics/SFC_Layer/UFS/sfc_nst_lake.f90 @@ -0,0 +1,684 @@ +!>\file sfc_nst.f90 +!! This file contains the GFS NSST model. + +!> This module contains the CCPP-compliant GFS near-surface sea temperature scheme. +module sfc_nst + + use machine , only : kind_phys, kp => kind_phys + use funcphys , only : fpvs + use module_nst_parameters , only : one, zero, half + use module_nst_parameters , only : t0k, cp_w, omg_m, omg_sh, sigma_r, solar_time_6am, sst_max + use module_nst_parameters , only : ri_c, z_w_max, delz, wd_max, rad2deg, const_rot, tau_min, tw_max + use module_nst_water_prop , only : get_dtzm_point, density, rhocoef, grv, sw_ps_9b + use nst_module , only : cool_skin, dtm_1p, cal_w, cal_ttop, convdepth, dtm_1p_fca + use nst_module , only : dtm_1p_tla, dtm_1p_mwa, dtm_1p_mda, dtm_1p_mta, dtl_reset + ! + implicit none +contains + + !>\defgroup gfs_nst_main_mod GFS Near-Surface Sea Temperature Module + !! This module contains the CCPP-compliant GFS near-surface sea temperature scheme. + !> @{ + !! This subroutine calls the Thermal Skin-layer and Diurnal Thermocline models to update the NSST profile. + !! \section arg_table_sfc_nst_run Argument Table + !! \htmlinclude sfc_nst_run.html + !! + !> \section NSST_general_algorithm GFS Near-Surface Sea Temperature Scheme General Algorithm + subroutine sfc_nst_run & + ( im, hvap, cp, hfus, jcal, eps, epsm1, rvrdm1, rd, rhw0, & ! --- inputs: + pi, tgice, sbc, ps, u1, v1, usfco, vsfco, icplocn2atm, t1, & + q1, tref, cm, ch, lseaspray, fm, fm10, & + prsl1, prslki, prsik1, prslk1, wet, use_lake_model, xlon, & + sinlat, stress, oceanfrac, & + sfcemis, dlwflx, sfcnsw, rain, timestep, kdt, solhr,xcosz, & + wind, flag_iter, flag_guess, nstf_name1, nstf_name4, & + 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: + ) + ! + ! ===================================================================== ! + ! description: ! + ! ! + ! ! + ! usage: ! + ! ! + ! call sfc_nst ! + ! inputs: ! + ! ( im, ps, u1, v1, t1, q1, tref, cm, ch, ! + ! lseaspray, fm, fm10, ! + ! prsl1, prslki, wet, use_lake_model, xlon, sinlat, stress, ! + ! sfcemis, dlwflx, sfcnsw, rain, timestep, kdt,solhr,xcosz, ! + ! wind, flag_iter, flag_guess, nstf_name1, nstf_name4, ! + ! 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, ! + ! -- outputs: + ! qsurf, gflux, cmm, chh, evap, hflx, ep ! + ! ) + ! ! + ! ! + ! subprogram/functions called: fpvs, density, rhocoef, cool_skin ! + ! ! + ! program history log: ! + ! 2007 -- xu li createad original code ! + ! 2008 -- s. moorthi adapted to the parallel version ! + ! may 2009 -- y.-t. hou modified to include input lw surface ! + ! emissivity from radiation. also replaced the ! + ! often comfusing combined sw and lw suface ! + ! flux with separate sfc net sw flux (defined ! + ! as dn-up) and lw flux. added a program doc block. ! + ! sep 2009 -- s. moorthi removed rcl and additional reformatting ! + ! and optimization + made pa as input pressure unit.! + ! 2009 -- xu li recreatead the code ! + ! feb 2010 -- s. moorthi added some changes made to the previous ! + ! version ! + ! Jul 2016 -- X. Li, modify the diurnal warming event reset ! + ! ! + ! ! + ! ==================== definition of variables ==================== ! + ! ! + ! inputs: size ! + ! im - integer, horiz dimension 1 ! + ! ps - real, surface pressure (pa) im ! + ! u1, v1 - real, u/v component of surface layer wind (m/s) im ! + ! usfco, vsfco - real, u/v component of surface current (m/s) im ! + ! icplocn2atm - integer, option to include ocean surface 1 ! + ! current in the computation of flux ! + ! t1 - real, surface layer mean temperature ( k ) im ! + ! q1 - real, surface layer mean specific humidity im ! + ! tref - real, reference/foundation temperature ( k ) im ! + ! cm - real, surface exchange coeff for momentum (m/s) im ! + ! ch - real, surface exchange coeff heat & moisture(m/s) im ! + ! lseaspray- logical, .t. for parameterization for sea spray 1 ! + ! fm - real, a stability profile function for momentum im ! + ! fm10 - real, a stability profile function for momentum im ! + ! at 10m ! + ! prsl1 - real, surface layer mean pressure (pa) im ! + ! prslki - real, im ! + ! prsik1 - real, im ! + ! prslk1 - real, im ! + ! wet - logical, =T if any ocn/lake water (F otherwise) im ! + ! use_lake_model- logical, =T if flake model is used for lake im ! + ! icy - logical, =T if any ice im ! + ! xlon - real, longitude (radians) im ! + ! sinlat - real, sin of latitude im ! + ! stress - real, wind stress (n/m**2) im ! + ! sfcemis - real, sfc lw emissivity (fraction) im ! + ! dlwflx - real, total sky sfc downward lw flux (w/m**2) im ! + ! sfcnsw - real, total sky sfc netsw flx into ocean (w/m**2) im ! + ! rain - real, rainfall rate (kg/m**2/s) im ! + ! timestep - real, timestep interval (second) 1 ! + ! kdt - integer, time step counter 1 ! + ! solhr - real, fcst hour at the end of prev time step 1 ! + ! xcosz - real, consine of solar zenith angle 1 ! + ! wind - real, wind speed (m/s) im ! + ! flag_iter- logical, execution or not im ! + ! when iter = 1, flag_iter = .true. for all grids im ! + ! when iter = 2, flag_iter = .true. when wind < 2 im ! + ! for both land and ocean (when nstf_name1 > 0) im ! + ! flag_guess-logical, .true.= guess step to get CD et al im ! + ! when iter = 1, flag_guess = .true. when wind < 2 im ! + ! when iter = 2, flag_guess = .false. for all grids im ! + ! nstf_name - integers , NSST related flag parameters 1 ! + ! nstf_name1 : 0 = NSSTM off 1 ! + ! 1 = NSSTM on but uncoupled 1 ! + ! 2 = NSSTM on and coupled 1 ! + ! nstf_name4 : zsea1 in mm 1 ! + ! 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 + ! tskin - real, ocean surface skin temperature ( k ) im ! + ! tsurf - real, the same as tskin ( k ) but for guess run im ! + ! xt - real, heat content in dtl im ! + ! xs - real, salinity content in dtl im ! + ! xu - real, u-current content in dtl im ! + ! xv - real, v-current content in dtl im ! + ! xz - real, dtl thickness im ! + ! zm - real, mxl thickness im ! + ! xtts - real, d(xt)/d(ts) im ! + ! xzts - real, d(xz)/d(ts) im ! + ! dt_cool - real, sub-layer cooling amount im ! + ! d_conv - real, thickness of free convection layer (fcl) im ! + ! z_c - sub-layer cooling thickness im ! + ! c_0 - coefficient1 to calculate d(tz)/d(ts) im ! + ! c_d - coefficient2 to calculate d(tz)/d(ts) im ! + ! w_0 - coefficient3 to calculate d(tz)/d(ts) im ! + ! w_d - coefficient4 to calculate d(tz)/d(ts) im ! + ! ifd - real, index to start dtlm run or not im ! + ! qrain - real, sensible heat flux due to rainfall (watts) im ! + + ! outputs: ! + + ! qsurf - real, surface air saturation specific humidity im ! + ! gflux - real, soil heat flux (w/m**2) im ! + ! cmm - real, im ! + ! chh - real, im ! + ! evap - real, evaperation from latent heat flux im ! + ! hflx - real, sensible heat flux im ! + ! ep - real, potential evaporation im ! + ! ! + ! ===================================================================== ! + + + + ! --- inputs: + integer, intent(in) :: im, kdt, ipr, nstf_name1, nstf_name4, nstf_name5 + integer, intent(in) :: icplocn2atm + + real (kind=kind_phys), intent(in) :: hvap, cp, hfus, jcal, eps, & + epsm1, rvrdm1, rd, rhw0, sbc, pi, tgice + real (kind=kind_phys), dimension(:), intent(in) :: ps, u1, v1, & + usfco, vsfco, t1, q1, cm, ch, fm, fm10, oceanfrac, & + prsl1, prslki, prsik1, prslk1, xlon, xcosz, & + sinlat, stress, sfcemis, dlwflx, sfcnsw, rain, wind + real (kind=kind_phys), dimension(:), intent(in), optional :: & + tref + real (kind=kind_phys), intent(in) :: timestep + real (kind=kind_phys), intent(in) :: solhr + + ! For sea spray effect + logical, intent(in) :: lseaspray + ! + logical, dimension(:), intent(in) :: flag_iter, flag_guess, wet + integer, dimension(:), intent(in) :: use_lake_model + 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 + real (kind=kind_phys), dimension(:), intent(inout) :: tskin, & + tsurf + real (kind=kind_phys), dimension(:), intent(inout), optional :: & + xt, xs, xu, xv, xz, zm, xtts, xzts, dt_cool, & + z_c, c_0, c_d, w_0, w_d, d_conv, ifd, qrain + + ! --- outputs: + real (kind=kind_phys), dimension(:), intent(inout) :: qsurf, gflux, cmm, chh, evap, hflx, ep + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! + ! locals + ! + integer :: k,i + ! + real (kind=kind_phys), dimension(im) :: q0, qss, rch, rho_a, theta1, tv1, wndmag + + real(kind=kind_phys) :: elocp,tem,cpinv,hvapi + ! + ! nstm related prognostic fields + ! + logical :: flag(im) + real (kind=kind_phys), dimension(im) :: xt_old, xs_old, xu_old, xv_old, xz_old, & + zm_old,xtts_old, xzts_old, ifd_old, tref_old, tskin_old, dt_cool_old,z_c_old + + real(kind=kind_phys) :: ulwflx(im), nswsfc(im) + ! real(kind=kind_phys) rig(im), + ! & ulwflx(im),dlwflx(im), + ! & slrad(im),nswsfc(im) + real(kind=kind_phys) :: alpha,beta,rho_w,f_nsol,sss,sep, cosa,sina,taux,tauy, & + grav,dz,t0,ttop0,ttop + + real(kind=kind_phys) :: le,fc,dwat,dtmp,wetc,alfac,ustar_a,rich + real(kind=kind_phys) :: rnl_ts,hs_ts,hl_ts,rf_ts,q_ts + real(kind=kind_phys) :: fw,q_warm + real(kind=kind_phys) :: t12,alon,tsea,sstc,dta,dtz + real(kind=kind_phys) :: zsea1,zsea2,soltim + logical :: do_nst + ! + ! parameters for sea spray effect + ! + real (kind=kind_phys) :: f10m, u10m, v10m, ws10, ru10, qss1, & + bb1, hflxs, evaps, ptem + ! + ! real (kind=kind_phys), parameter :: alps=0.5, bets=0.5, gams=0.1, + ! real (kind=kind_phys), parameter :: alps=0.5, bets=0.5, gams=0.0, + ! real (kind=kind_phys), parameter :: alps=1.0, bets=1.0, gams=0.2, + real (kind=kind_phys), parameter :: alps=0.75,bets=0.75,gams=0.15, & + ws10cr=30., conlf=7.2e-9, consf=6.4e-8 + real (kind=kind_phys) :: windrel + ! + !====================================================================================================== + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + if (nstf_name1 == 0) return ! No NSST model used + + cpinv = one/cp + hvapi = one/hvap + elocp = hvap/cp + + sss = 34.0_kp ! temporarily, when sea surface salinity data is not ready + ! + ! flag for open water and where the iteration is on + ! + do_nst = .false. + do i = 1, im + ! flag(i) = wet(i) .and. .not.icy(i) .and. flag_iter(i) + flag(i) = wet(i) .and. flag_iter(i) .and. use_lake_model(i)/=1 & + .and. oceanfrac(i)==0. + do_nst = do_nst .or. flag(i) + enddo + if (.not. do_nst) return + ! + ! save nst-related prognostic fields for guess run + ! + do i=1, im + ! if(wet(i) .and. .not.icy(i) .and. flag_guess(i)) then + if(wet(i) .and. flag_guess(i) .and. use_lake_model(i)/=1 .and. & + oceanfrac(i)==0.) then + xt_old(i) = xt(i) + xs_old(i) = xs(i) + xu_old(i) = xu(i) + xv_old(i) = xv(i) + xz_old(i) = xz(i) + zm_old(i) = zm(i) + xtts_old(i) = xtts(i) + xzts_old(i) = xzts(i) + ifd_old(i) = ifd(i) + tskin_old(i) = tskin(i) + dt_cool_old(i) = dt_cool(i) + z_c_old(i) = z_c(i) + endif + enddo + + + ! --- ... initialize variables. all units are m.k.s. unless specified. + ! ps is in pascals, wind is wind speed, theta1 is surface air + ! estimated from level 1 temperature, rho_a is air density and + ! qss is saturation specific humidity at the water surface + !! + do i = 1, im + if ( flag(i) ) then + + nswsfc(i) = sfcnsw(i) ! net solar radiation at the air-sea surface (positive downward) + wndmag(i) = sqrt(u1(i)*u1(i) + v1(i)*v1(i)) + + q0(i) = max(q1(i), 1.0e-8_kp) + + 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 + qss(i) = eps*qss(i) / (ps(i) + epsm1*qss(i)) ! pa + ! + evap(i) = zero + hflx(i) = zero + gflux(i) = zero + ep(i) = zero + + ! --- ... rcp = rho cp ch v + + if (icplocn2atm ==0) then + rch(i) = rho_a(i) * cp * ch(i) * wind(i) + cmm(i) = cm (i) * wind(i) + chh(i) = rho_a(i) * ch(i) * wind(i) + else if (icplocn2atm ==1) then + windrel= sqrt( (u1(i)-usfco(i))**2 + (v1(i)-vsfco(i))**2 ) + rch(i) = rho_a(i) * cp * ch(i) * windrel + cmm(i) = cm (i) * windrel + chh(i) = rho_a(i) * ch(i) * windrel + endif + + !> - Calculate latent and sensible heat flux over open water with tskin. + ! at previous time step + evap(i) = elocp * rch(i) * (qss(i) - q0(i)) + qsurf(i) = qss(i) + + 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) + ! &,' tsurf=',tsurf(i) + endif + enddo + + ! run nst model: dtm + slm + ! + zsea1 = 0.001_kp*real(nstf_name4) + zsea2 = 0.001_kp*real(nstf_name5) + + !> - Call module_nst_water_prop::density() to compute sea water density. + !> - Call module_nst_water_prop::rhocoef() to compute thermal expansion + !! coefficient (\a alpha) and saline contraction coefficient (\a beta). + do i = 1, im + if ( flag(i) ) then + tsea = tsurf(i) + t12 = tsea*tsea + ulwflx(i) = sfcemis(i) * sbc * t12 * t12 + alon = xlon(i)*rad2deg + grav = grv(sinlat(i)) + soltim = mod(alon/15.0_kp + solhr, 24.0_kp)*3600.0_kp + call density(tsea,sss,rho_w) ! sea water density + call rhocoef(tsea,sss,rho_w,alpha,beta) ! alpha & beta + ! + !> - Calculate sensible heat flux (\a qrain) due to rainfall. + ! + le = (2.501_kp-0.00237_kp*tsea)*1.0e6_kp + dwat = 2.11e-5_kp*(t1(i)/t0k)**1.94_kp ! water vapor diffusivity + dtmp = (one+3.309e-3_kp*(t1(i)-t0k)-1.44e-6_kp*(t1(i)-t0k) & + * (t1(i)-t0k))*0.02411_kp/(rho_a(i)*cp) ! heat diffusivity + wetc = 622.0_kp*le*qss(i)/(rd*t1(i)*t1(i)) + alfac = one / (one + (wetc*le*dwat)/(cp*dtmp)) ! wet bulb factor + tem = (1.0e3_kp * rain(i) / rho_w) * alfac * cp_w + qrain(i) = tem * (tsea-t1(i)+1.0e3_kp*(qss(i)-q0(i))*le/cp) + + !> - Calculate input non solar heat flux as upward = positive to models here + + f_nsol = hflx(i) + evap(i) + ulwflx(i) - dlwflx(i) + omg_sh*qrain(i) + + ! if (lprnt .and. i == ipr) print *,' f_nsol=',f_nsol,' hflx=', + ! &hflx(i),' evap=',evap(i),' ulwflx=',ulwflx(i),' dlwflx=',dlwflx(i) + ! &,' omg_sh=',omg_sh,' qrain=',qrain(i) + + sep = sss*(evap(i)/le-rain(i))/rho_w + ustar_a = sqrt(stress(i)/rho_a(i)) ! air friction velocity + ! + ! sensitivities of heat flux components to ts + ! + rnl_ts = 4.0_kp*sfcemis(i)*sbc*tsea*tsea*tsea ! d(rnl)/d(ts) + hs_ts = rch(i) + hl_ts = rch(i)*elocp*eps*hvap*qss(i)/(rd*t12) + rf_ts = tem * (one+rch(i)*hl_ts) + q_ts = rnl_ts + hs_ts + hl_ts + omg_sh*rf_ts + ! + !> - Call cool_skin(), which is the sub-layer cooling parameterization + !! (Fairfall et al. (1996) \cite fairall_et_al_1996). + ! & calculate c_0, c_d + ! + call cool_skin(ustar_a,f_nsol,nswsfc(i),evap(i),sss,alpha,beta, & + rho_w,rho_a(i),tsea,q_ts,hl_ts,grav,le, & + dt_cool(i),z_c(i),c_0(i),c_d(i)) + + tem = one / wndmag(i) + cosa = u1(i)*tem + sina = v1(i)*tem + taux = max(stress(i),tau_min)*cosa + tauy = max(stress(i),tau_min)*sina + fc = const_rot*sinlat(i) + ! + ! Run DTM-1p system. + ! + if ( (soltim > solar_time_6am .and. ifd(i) == zero) ) then + else + ifd(i) = one + ! + ! calculate fcl thickness with current forcing and previous time's profile + ! + ! if (lprnt .and. i == ipr) print *,' beg xz=',xz(i) + + !> - Call convdepth() to calculate depth for convective adjustments. + if ( f_nsol > zero .and. xt(i) > zero ) then + call convdepth(kdt,timestep,nswsfc(i),f_nsol,sss,sep,rho_w, & + alpha,beta,xt(i),xs(i),xz(i),d_conv(i)) + else + d_conv(i) = zero + endif + + ! if (lprnt .and. i == ipr) print *,' beg xz1=',xz(i) + ! + ! determine rich: wind speed dependent (right now) + ! + ! if ( wind(i) < 1.0 ) then + ! rich = 0.25 + 0.03*wind(i) + ! elseif ( wind(i) >= 1.0 .and. wind(i) < 1.5 ) then + ! rich = 0.25 + 0.1*wind(i) + ! elseif ( wind(i) >= 1.5 .and. wind(i) < 6.0 ) then + ! rich = 0.25 + 0.6*wind(i) + ! elseif ( wind(i) >= 6.0 ) then + ! rich = 0.25 + min(0.8*wind(i),0.50) + ! endif + + rich = ri_c + + !> - Call the diurnal thermocline layer model dtm_1p(). + call dtm_1p(kdt,timestep,rich,taux,tauy,nswsfc(i), & + f_nsol,sss,sep,q_ts,hl_ts,rho_w,alpha,beta,alon, & + sinlat(i),soltim,grav,le,d_conv(i), & + xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i)) + + ! if (lprnt .and. i == ipr) print *,' beg xz2=',xz(i) + + ! apply mda + if ( xt(i) > zero ) then + !> - If \a dtl heat content \a xt > 0.0, call dtm_1p_mda() to apply + !! minimum depth adjustment (mda). + call dtm_1p_mda(xt(i),xtts(i),xz(i),xzts(i)) + if ( xz(i) >= z_w_max ) then + !> - If \a dtl thickness >= module_nst_parameters::z_w_max, call dtl_reset() + !! to reset xt/xs/x/xv to zero, and xz to module_nst_parameters::z_w_max. + call dtl_reset(xt(i),xs(i),xu(i),xv(i),xz(i),xtts(i), xzts(i)) + + ! if (lprnt .and. i == ipr) print *,' beg xz3=',xz(i),' z_w_max=' + ! &,z_w_max + endif + + ! apply fca + if ( d_conv(i) > zero ) then + !> - If thickness of free convection layer > 0.0, call dtm_1p_fca() + !! to apply free convection adjustment. + !> - If \a dtl thickness >= module_nst_parameters::z_w_max(), call dtl_reset() + !! to reset xt/xs/x/xv to zero, and xz to module_nst_parameters::z_w_max(). + call dtm_1p_fca(d_conv(i),xt(i),xtts(i),xz(i),xzts(i)) + if ( xz(i) >= z_w_max ) then + call dtl_reset (xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i)) + endif + endif + + ! if (lprnt .and. i == ipr) print *,' beg xz4=',xz(i) + + ! apply tla + dz = min(xz(i),max(d_conv(i),delz)) + ! + !> - Call sw_ps_9b() to compute the fraction of the solar radiation + !! absorbed by the depth \a delz (Paulson and Simpson (1981) \cite paulson_and_simpson_1981). + !! And calculate the total heat absorbed in warm layer. + call sw_ps_9b(delz,fw) + q_warm = fw*nswsfc(i)-f_nsol !total heat absorbed in warm layer + + !> - Call cal_ttop() to calculate the diurnal warming amount at the top layer with + !! thickness of \a dz. + if ( q_warm > zero ) then + call cal_ttop(kdt,timestep,q_warm,rho_w,dz, xt(i),xz(i),ttop0) + + ! if (lprnt .and. i == ipr) print *,' d_conv=',d_conv(i),' delz=', + ! &delz,' kdt=',kdt,' timestep=',timestep,' nswsfc=',nswsfc(i), + ! &' f_nsol=',f_nsol,' rho_w=',rho_w,' dz=',dz,' xt=',xt(i), + ! &' xz=',xz(i),' qrain=',qrain(i) + + ttop = ((xt(i)+xt(i))/xz(i))*(one-dz/((xz(i)+xz(i)))) + + ! if (lprnt .and. i == ipr) print *,' beg xz4a=',xz(i) + ! &,' ttop=',ttop,' ttop0=',ttop0,' xt=',xt(i),' dz=',dz + ! &,' xznew=',(xt(i)+sqrt(xt(i)*(xt(i)-dz*ttop0)))/ttop0 + + !> - Call dtm_1p_tla() to apply top layer adjustment. + if ( ttop > ttop0 ) then + call dtm_1p_tla(dz,ttop0,xt(i),xtts(i),xz(i),xzts(i)) + + ! if (lprnt .and. i == ipr) print *,' beg xz4b=',xz(i),'z_w_max=', + ! &z_w_max + if ( xz(i) >= z_w_max ) then + call dtl_reset (xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i)) + endif + endif + endif ! if ( q_warm > 0.0 ) then + + ! if (lprnt .and. i == ipr) print *,' beg xz5=',xz(i) + + ! apply mwa + !> - Call dt_1p_mwa() to apply maximum warming adjustment. + t0 = (xt(i)+xt(i))/xz(i) + if ( t0 > tw_max ) then + call dtm_1p_mwa(xt(i),xtts(i),xz(i),xzts(i)) + if ( xz(i) >= z_w_max ) then + call dtl_reset (xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i)) + endif + endif + + ! if (lprnt .and. i == ipr) print *,' beg xz6=',xz(i) + + ! apply mta + !> - Call dtm_1p_mta() to apply maximum temperature adjustment. + sstc = tref(i) + (xt(i)+xt(i))/xz(i) - dt_cool(i) + + if ( sstc > sst_max ) then + dta = sstc - sst_max + call dtm_1p_mta(dta,xt(i),xtts(i),xz(i),xzts(i)) + ! write(*,'(a,f3.0,7f8.3)') 'mta, sstc,dta :',islimsk(i), + ! & sstc,dta,tref(i),xt(i),xz(i),2.0*xt(i)/xz(i),dt_cool(i) + if ( xz(i) >= z_w_max ) then + call dtl_reset (xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i)) + endif + endif + ! + endif ! if ( xt(i) > 0.0 ) then + ! reset dtl at midnight and when solar zenith angle > 89.994 degree + if ( abs(soltim) < 2.0_kp*timestep ) then + call dtl_reset (xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i)) + endif + + endif ! if (solar_time > solar_time_6am .and. ifd(i) == 0.0 ) then: too late to start the first day + + ! if (lprnt .and. i == ipr) print *,' beg xz7=',xz(i) + + ! update tsurf (when flag(i) .eqv. .true. ) + !> - Call get_dtzm_point() to computes \a dtz and \a tsurf. + call get_dtzm_point(xt(i),xz(i),dt_cool(i),z_c(i), zsea1,zsea2,dtz) + tsurf(i) = max(tgice, tref(i) + dtz ) + + ! if (lprnt .and. i == ipr) print *,' tsurf=',tsurf(i),' tref=', + ! &tref(i),' xz=',xz(i),' dt_cool=',dt_cool(i) + + !> - Call cal_w() to calculate \a w_0 and \a w_d. + if ( xt(i) > zero ) then + call cal_w(kdt,xz(i),xt(i),xzts(i),xtts(i),w_0(i),w_d(i)) + else + w_0(i) = zero + w_d(i) = zero + endif + + ! if ( xt(i) > 0.0 ) then + ! rig(i) = grav*xz(i)*xz(i)*(alpha*xt(i)-beta*xs(i)) + ! & /(2.0*(xu(i)*xu(i)+xv(i)*xv(i))) + ! else + ! rig(i) = 0.25 + ! endif + + ! qrain(i) = rig(i) + zm(i) = wind(i) + + endif + enddo + + ! restore nst-related prognostic fields for guess run + do i=1, im + ! if (wet(i) .and. .not.icy(i)) then + if (wet(i) .and. use_lake_model(i)/=1 .and. oceanfrac(i)==0.) & + then + if (flag_guess(i)) then ! when it is guess of + xt(i) = xt_old(i) + xs(i) = xs_old(i) + xu(i) = xu_old(i) + xv(i) = xv_old(i) + xz(i) = xz_old(i) + zm(i) = zm_old(i) + xtts(i) = xtts_old(i) + xzts(i) = xzts_old(i) + ifd(i) = ifd_old(i) + tskin(i) = tskin_old(i) + dt_cool(i) = dt_cool_old(i) + z_c(i) = z_c_old(i) + else + ! + ! update tskin when coupled and not guess run + ! (all other NSST variables have been updated in this case) + ! + if ( nstf_name1 > 1 ) then + tskin(i) = tsurf(i) + endif ! if nstf_name1 > 1 then + endif ! if flag_guess(i) then + endif ! if wet(i) .and. .not.icy(i) then + enddo + + ! if (lprnt .and. i == ipr) print *,' beg xz8=',xz(i) + + if ( nstf_name1 > 1 ) then + !> - Calculate latent and sensible heat flux over open water with updated tskin + !! for the grids of open water and the iteration is on. + do i = 1, im + if ( flag(i) ) then + qss(i) = fpvs( tskin(i) ) + qss(i) = eps*qss(i) / (ps(i) + epsm1*qss(i)) + qsurf(i) = qss(i) + evap(i) = elocp*rch(i) * (qss(i) - q0(i)) + + 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 + ! + !> - Include sea spray effects + ! + do i=1,im + if(lseaspray .and. flag(i)) then + f10m = fm10(i) / fm(i) + u10m = f10m * u1(i) + v10m = f10m * v1(i) + ws10 = sqrt(u10m*u10m + v10m*v10m) + ws10 = max(ws10,1.) + ws10 = min(ws10,ws10cr) + tem = .015 * ws10 * ws10 + ru10 = 1. - .087 * log(10./tem) + qss1 = fpvs(t1(i)) + qss1 = eps * qss1 / (prsl1(i) + epsm1 * qss1) + tem = rd * cp * t1(i) * t1(i) + tem = 1. + eps * hvap * hvap * qss1 / tem + bb1 = 1. / tem + evaps = conlf * (ws10**5.4) * ru10 * bb1 + evaps = evaps * rho_a(i) * hvap * (qss1 - q0(i)) + evap(i) = evap(i) + alps * evaps + hflxs = consf * (ws10**3.4) * ru10 + hflxs = hflxs * rho_a(i) * cp * (tskin(i) - t1(i)) + ptem = alps - gams + hflx(i) = hflx(i) + bets * hflxs - ptem * evaps + endif + enddo + ! + do i=1,im + if ( flag(i) ) then + tem = one / rho_a(i) + hflx(i) = hflx(i) * tem * cpinv + evap(i) = evap(i) * tem * hvapi + endif + enddo + ! + ! if (lprnt) print *,' tskin=',tskin(ipr) + + return + end subroutine sfc_nst_run + !> @} +end module sfc_nst diff --git a/physics/SFC_Layer/UFS/sfc_nst_lake.meta b/physics/SFC_Layer/UFS/sfc_nst_lake.meta new file mode 100644 index 000000000..55d351e65 --- /dev/null +++ b/physics/SFC_Layer/UFS/sfc_nst_lake.meta @@ -0,0 +1,668 @@ +[ccpp-table-properties] + name = sfc_nst + type = scheme + dependencies = date_def.f,../../tools/funcphys.f90,../../hooks/machine.F,module_nst_model.f90,module_nst_parameters.f90,module_nst_water_prop.f90 + +######################################################################## +[ccpp-arg-table] + name = sfc_nst_run + type = scheme +[im] + standard_name = horizontal_loop_extent + long_name = horizontal loop extent + units = count + dimensions = () + type = integer + intent = in +[hvap] + standard_name = latent_heat_of_vaporization_of_water_at_0C + long_name = latent heat of evaporation/sublimation + units = J kg-1 + dimensions = () + type = real + kind = kind_phys + intent = in +[cp] + standard_name = specific_heat_of_dry_air_at_constant_pressure + long_name = specific heat of dry air at constant pressure + units = J kg-1 K-1 + dimensions = () + type = real + kind = kind_phys + intent = in +[hfus] + standard_name = latent_heat_of_fusion_of_water_at_0C + long_name = latent heat of fusion + units = J kg-1 + dimensions = () + type = real + kind = kind_phys + intent = in +[jcal] + standard_name = joules_per_calorie_constant + long_name = joules per calorie constant + units = J cal-1 + dimensions = () + type = real + kind = kind_phys + intent = in +[eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys + intent = in +[epsm1] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants_minus_one + long_name = (rd/rv) - 1 + units = none + dimensions = () + type = real + kind = kind_phys + intent = in +[rvrdm1] + standard_name = ratio_of_vapor_to_dry_air_gas_constants_minus_one + long_name = (rv/rd) - 1 (rv = ideal gas constant for water vapor) + units = none + dimensions = () + type = real + kind = kind_phys + intent = in +[rd] + standard_name = gas_constant_of_dry_air + long_name = ideal gas constant for dry air + units = J kg-1 K-1 + dimensions = () + type = real + kind = kind_phys + intent = in +[rhw0] + standard_name = sea_water_reference_density + long_name = sea water reference density + units = kg m-3 + dimensions = () + type = real + kind = kind_phys + intent = in +[pi] + standard_name = pi + long_name = ratio of a circle's circumference to its diameter + units = none + dimensions = () + type = real + kind = kind_phys + intent = in +[tgice] + standard_name = freezing_point_temperature_of_seawater + long_name = freezing point temperature of seawater + units = K + dimensions = () + type = real + kind = kind_phys + intent = in +[sbc] + standard_name = stefan_boltzmann_constant + long_name = Stefan-Boltzmann constant + units = W m-2 K-4 + dimensions = () + type = real + kind = kind_phys + intent = in +[ps] + standard_name = surface_air_pressure + long_name = surface pressure + units = Pa + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[u1] + standard_name = x_wind_at_surface_adjacent_layer + long_name = x component of surface layer wind + units = m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[v1] + standard_name = y_wind_at_surface_adjacent_layer + long_name = y component of surface layer wind + units = m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[usfco] + standard_name = x_ocean_current + long_name = zonal current at ocean surface + units = m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[vsfco] + standard_name = y_ocean_current + long_name = meridional current at ocean surface + units = m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[icplocn2atm] + standard_name = control_for_air_sea_flux_computation_over_water + long_name = air-sea flux option + units = 1 + dimensions = () + type = integer + intent = in +[t1] + standard_name = air_temperature_at_surface_adjacent_layer + long_name = surface layer mean temperature + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[q1] + standard_name = specific_humidity_at_surface_adjacent_layer + long_name = surface layer mean specific humidity + units = kg kg-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[tref] + standard_name = reference_sea_surface_temperature + long_name = reference/foundation temperature + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = True +[cm] + standard_name = surface_drag_coefficient_for_momentum_in_air_over_water + long_name = surface exchange coeff for momentum over water + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[ch] + standard_name = surface_drag_coefficient_for_heat_and_moisture_in_air_over_water + long_name = surface exchange coeff heat surface exchange coeff heat & moisture over ocean moisture over water + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[lseaspray] + standard_name = flag_for_sea_spray + long_name = flag for sea spray parameterization + units = flag + dimensions = () + type = logical + intent = in +[fm] + standard_name = Monin_Obukhov_similarity_function_for_momentum_over_water + long_name = Monin-Obukhov similarity function for momentum over water + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[fm10] + standard_name = Monin_Obukhov_similarity_function_for_momentum_at_10m_over_water + long_name = Monin-Obukhov similarity parameter for momentum at 10m over water + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[prsl1] + standard_name = air_pressure_at_surface_adjacent_layer + long_name = surface layer mean pressure + units = Pa + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[prslki] + standard_name = ratio_of_exner_function_between_midlayer_and_interface_at_lowest_model_layer + long_name = Exner function ratio bt midlayer and interface at 1st layer + units = ratio + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[prsik1] + standard_name = surface_dimensionless_exner_function + long_name = dimensionless Exner function at the ground surface + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[prslk1] + standard_name = dimensionless_exner_function_at_surface_adjacent_layer + long_name = dimensionless Exner function at the lowest model layer + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[wet] + standard_name = flag_nonzero_wet_surface_fraction + long_name = flag indicating presence of some ocean or lake surface area fraction + units = flag + dimensions = (horizontal_loop_extent) + type = logical + intent = in +[oceanfrac] + standard_name = sea_area_fraction + long_name = fraction of horizontal grid area occupied by ocean + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[use_lake_model] + standard_name = flag_for_using_lake_model + long_name = flag indicating lake points using a lake model + units = flag + dimensions = (horizontal_loop_extent) + type = integer + intent = in +[xlon] + standard_name = longitude + long_name = longitude + units = radian + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[sinlat] + standard_name = sine_of_latitude + long_name = sine of latitude + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[stress] + standard_name = surface_wind_stress_over_water + long_name = surface wind stress over water + units = m2 s-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[sfcemis] + standard_name = surface_longwave_emissivity_over_water + long_name = surface lw emissivity in fraction over water + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[dlwflx] + standard_name = surface_downwelling_longwave_flux_absorbed_by_ground_over_water + long_name = total sky surface downward longwave flux absorbed by the ground over water + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[sfcnsw] + standard_name = surface_net_downwelling_shortwave_flux + long_name = total sky sfc net sw flx into ocean + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[rain] + standard_name = nonnegative_lwe_thickness_of_precipitation_amount_on_dynamics_timestep_over_water + long_name = total precipitation amount in each time step over water + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[timestep] + standard_name = timestep_for_dynamics + long_name = timestep interval + units = s + dimensions = () + type = real + kind = kind_phys + intent = in +[kdt] + standard_name = index_of_timestep + long_name = current time step index + units = index + dimensions = () + type = integer + intent = in +[solhr] + standard_name = forecast_utc_hour + long_name = time in hours after 00z at the current timestep + units = h + dimensions = () + type = real + kind = kind_phys + intent = in +[xcosz] + standard_name = instantaneous_cosine_of_zenith_angle + long_name = cosine of solar zenith angle + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[wind] + standard_name = wind_speed_at_lowest_model_layer + long_name = wind speed at lowest model level + units = m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[flag_iter] + standard_name = flag_for_iteration + long_name = flag for iteration + units = flag + dimensions = (horizontal_loop_extent) + type = logical + intent = in +[flag_guess] + standard_name = flag_for_guess_run + long_name = flag for guess run + units = flag + dimensions = (horizontal_loop_extent) + type = logical + intent = in +[nstf_name1] + standard_name = control_for_nsstm + long_name = NSSTM flag: off/uncoupled/coupled=0/1/2 + units = flag + dimensions = () + type = integer + intent = in +[nstf_name4] + standard_name = lower_bound_for_depth_of_sea_temperature_for_nsstm + long_name = zsea1 + units = mm + dimensions = () + type = integer + intent = in +[nstf_name5] + standard_name = upper_bound_for_depth_of_sea_temperature_for_nsstm + long_name = zsea2 + units = mm + dimensions = () + type = integer + intent = in +[lprnt] + standard_name = flag_print + long_name = flag for printing diagnostics to output + units = flag + dimensions = () + type = logical + intent = in +[ipr] + standard_name = horizontal_index_of_printed_column + long_name = horizontal index of printed column + units = index + dimensions = () + type = integer + intent = in +[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 +[tskin] + standard_name = surface_skin_temperature_for_nsst + long_name = ocean surface skin temperature + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[tsurf] + standard_name = surface_skin_temperature_after_iteration_over_water + long_name = surface skin temperature after iteration over water + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[xt] + standard_name = heat_content_in_diurnal_thermocline + long_name = heat content in diurnal thermocline layer + units = K m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[xs] + standard_name = sea_water_salinity_in_diurnal_thermocline + long_name = salinity content in diurnal thermocline layer + units = ppt m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[xu] + standard_name = x_current_in_diurnal_thermocline + long_name = u-current content in diurnal thermocline layer + units = m2 s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[xv] + standard_name = y_current_in_diurnal_thermocline + long_name = v-current content in diurnal thermocline layer + units = m2 s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[xz] + standard_name = diurnal_thermocline_layer_thickness + long_name = diurnal thermocline layer thickness + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[zm] + standard_name = ocean_mixed_layer_thickness + long_name = mixed layer thickness + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[xtts] + standard_name = derivative_of_heat_content_in_diurnal_thermocline_wrt_surface_skin_temperature + long_name = d(xt)/d(ts) + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[xzts] + standard_name = derivative_of_diurnal_thermocline_layer_thickness_wrt_surface_skin_temperature + long_name = d(xz)/d(ts) + units = m K-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[dt_cool] + standard_name = molecular_sublayer_temperature_correction_in_sea_water + long_name = sub-layer cooling amount + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[z_c] + standard_name = molecular_sublayer_thickness_in_sea_water + long_name = sub-layer cooling thickness + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[c_0] + standard_name = coefficient_c_0 + long_name = coefficient1 to calculate d(tz)/d(ts) + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[c_d] + standard_name = coefficient_c_d + long_name = coefficient2 to calculate d(tz)/d(ts) + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[w_0] + standard_name = coefficient_w_0 + long_name = coefficient3 to calculate d(tz)/d(ts) + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[w_d] + standard_name = coefficient_w_d + long_name = coefficient4 to calculate d(tz)/d(ts) + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[d_conv] + standard_name = free_convection_layer_thickness_in_sea_water + long_name = thickness of free convection layer + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[ifd] + standard_name = control_for_diurnal_thermocline_calculation + long_name = index to start dtlm run or not + units = index + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[qrain] + standard_name = surface_sensible_heat_due_to_rainfall + long_name = sensible heat flux due to rainfall + units = W + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[qsurf] + standard_name = surface_specific_humidity_over_water + long_name = surface air saturation specific humidity over water + units = kg kg-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[gflux] + standard_name = upward_heat_flux_in_soil_over_water + long_name = soil heat flux over water + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[cmm] + standard_name = surface_drag_wind_speed_for_momentum_in_air_over_water + long_name = momentum exchange coefficient over water + units = m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[chh] + standard_name = surface_drag_mass_flux_for_heat_and_moisture_in_air_over_water + long_name = thermal exchange coefficient over water + units = kg m-2 s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[evap] + standard_name = kinematic_surface_upward_latent_heat_flux_over_water + long_name = kinematic surface upward latent heat flux over water + units = kg kg-1 m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[hflx] + standard_name = kinematic_surface_upward_sensible_heat_flux_over_water + long_name = kinematic surface upward sensible heat flux over water + units = K m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[ep] + standard_name = surface_upward_potential_latent_heat_flux_over_water + long_name = surface upward potential latent heat flux over water + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out diff --git a/physics/SFC_Layer/UFS/sfc_nst_lake_post.f90 b/physics/SFC_Layer/UFS/sfc_nst_lake_post.f90 new file mode 100644 index 000000000..74c553440 --- /dev/null +++ b/physics/SFC_Layer/UFS/sfc_nst_lake_post.f90 @@ -0,0 +1,89 @@ +!> \file sfc_nst_post.f90 +!! This file contains code to be executed after the GFS NSST model. + +module sfc_nst_post + + use machine , only : kind_phys, kp => kind_phys + use module_nst_water_prop , only : get_dtzm_2d + + implicit none + +contains + + ! \defgroup GFS_NSST_POST GFS Near-Surface Sea Temperature Post + + !> \section arg_table_sfc_nst_post_run Argument Table + !! \htmlinclude sfc_nst_post_run.html + !! + ! \section NSST_general_post_algorithm General Algorithm + ! + ! \section NSST_detailed_post_algorithm Detailed Algorithm + ! @{ + subroutine sfc_nst_post_run & + ( im, kdt, rlapse, tgice, wet, use_lake_model, icy, oro, & + oro_uf, nstf_name1, oceanfrac, & + nstf_name4, nstf_name5, xt, xz, dt_cool, z_c, tref, xlon, & + tsurf_wat, tsfc_wat, nthreads, dtzm, errmsg, errflg & + ) + ! --- inputs: + integer, intent(in) :: im, kdt, nthreads + logical, dimension(:), intent(in) :: wet, icy + integer, dimension(:), intent(in) :: use_lake_model + real (kind=kind_phys), intent(in) :: rlapse, tgice + real (kind=kind_phys), dimension(:), intent(in) :: oro, oro_uf + integer, intent(in) :: nstf_name1, nstf_name4, nstf_name5 + real (kind=kind_phys), dimension(:), intent(in) :: xlon, oceanfrac + real (kind=kind_phys), dimension(:), intent(in), optional :: xt, xz, dt_cool, z_c, tref + + ! --- input/outputs: + real (kind=kind_phys), dimension(:), intent(inout) :: tsurf_wat, tsfc_wat + + ! --- outputs: + real (kind=kind_phys), dimension(:), intent(out) :: dtzm + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! --- locals + integer :: i + real(kind=kind_phys) :: zsea1, zsea2 + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + ! if (lprnt) print *,' tseaz2=',tseal(ipr),' tref=',tref(ipr), + ! & ' dt_cool=',dt_cool(ipr),' dt_warm=',2.0*xt(ipr)/xz(ipr), + ! & ' kdt=',kdt + + ! do i = 1, im + ! if (wet(i) .and. .not. icy(i)) then + ! tsurf_wat(i) = tsurf_wat(i) - (oro(i)-oro_uf(i)) * rlapse + ! endif + ! enddo + + ! --- ... run nsst model ... --- + + if (nstf_name1 > 1) then + zsea1 = 0.001_kp*real(nstf_name4) + zsea2 = 0.001_kp*real(nstf_name5) + call get_dtzm_2d (xt, xz, dt_cool, z_c, wet, zsea1, zsea2, im, 1, nthreads, dtzm) + do i = 1, im + ! if (wet(i) .and. .not.icy(i)) then + ! if (wet(i) .and. (frac_grid .or. .not. icy(i))) then + if (wet(i) .and. use_lake_model(i) /=1 .and. oceanfrac(i)==0.)& + then + tsfc_wat(i) = max(tgice, tref(i) + dtzm(i)) + ! tsfc_wat(i) = max(271.2, tref(i) + dtzm(i)) - & + ! (oro(i)-oro_uf(i))*rlapse + endif + enddo + endif + + ! if (lprnt) print *,' tseaz2=',tsea(ipr),' tref=',tref(ipr), & + ! & ' dt_cool=',dt_cool(ipr),' dt_warm=',dt_warm(ipr),' kdt=',kdt + + return + end subroutine sfc_nst_post_run + +end module sfc_nst_post diff --git a/physics/SFC_Layer/UFS/sfc_nst_lake_post.meta b/physics/SFC_Layer/UFS/sfc_nst_lake_post.meta new file mode 100644 index 000000000..c4734162f --- /dev/null +++ b/physics/SFC_Layer/UFS/sfc_nst_lake_post.meta @@ -0,0 +1,205 @@ +######################################################################## +[ccpp-table-properties] + name = sfc_nst_post + type = scheme + dependencies = ../../hooks/machine.F,module_nst_parameters.f90,module_nst_water_prop.f90 + +######################################################################## +[ccpp-arg-table] + name = sfc_nst_post_run + type = scheme +[im] + standard_name = horizontal_loop_extent + long_name = horizontal loop extent + units = count + dimensions = () + type = integer + intent = in +[kdt] + standard_name = index_of_timestep + long_name = current time step index + units = index + dimensions = () + type = integer + intent = in +[rlapse] + standard_name = air_temperature_lapse_rate_constant + long_name = environmental air temperature lapse rate constant + units = K m-1 + dimensions = () + type = real + kind = kind_phys + intent = in +[tgice] + standard_name = freezing_point_temperature_of_seawater + long_name = freezing point temperature of seawater + units = K + dimensions = () + type = real + kind = kind_phys + intent = in +[wet] + standard_name = flag_nonzero_wet_surface_fraction + long_name = flag indicating presence of some ocean or lake surface area fraction + units = flag + dimensions = (horizontal_loop_extent) + type = logical + intent = in +[oceanfrac] + standard_name = sea_area_fraction + long_name = fraction of horizontal grid area occupied by ocean + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[use_lake_model] + standard_name = flag_for_using_lake_model + long_name = flag indicating lake points using a lake model + units = flag + dimensions = (horizontal_loop_extent) + type = integer + intent = in +[icy] + standard_name = flag_nonzero_sea_ice_surface_fraction + long_name = flag indicating presence of some sea ice surface area fraction + units = flag + dimensions = (horizontal_loop_extent) + type = logical + intent = in +[oro] + standard_name = height_above_mean_sea_level + long_name = height_above_mean_sea_level + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[oro_uf] + standard_name = unfiltered_height_above_mean_sea_level + long_name = unfiltered height_above_mean_sea_level + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[nstf_name1] + standard_name = control_for_nsstm + long_name = NSSTM flag: off/uncoupled/coupled=0/1/2 + units = flag + dimensions = () + type = integer + intent = in +[nstf_name4] + standard_name = lower_bound_for_depth_of_sea_temperature_for_nsstm + long_name = zsea1 + units = mm + dimensions = () + type = integer + intent = in +[nstf_name5] + standard_name = upper_bound_for_depth_of_sea_temperature_for_nsstm + long_name = zsea2 + units = mm + dimensions = () + type = integer + intent = in +[xt] + standard_name = heat_content_in_diurnal_thermocline + long_name = heat content in diurnal thermocline layer + units = K m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = True +[xz] + standard_name = diurnal_thermocline_layer_thickness + long_name = diurnal thermocline layer thickness + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = True +[dt_cool] + standard_name = molecular_sublayer_temperature_correction_in_sea_water + long_name = sub-layer cooling amount + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = True +[z_c] + standard_name = molecular_sublayer_thickness_in_sea_water + long_name = sub-layer cooling thickness + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = True +[tref] + standard_name = reference_sea_surface_temperature + long_name = reference/foundation temperature + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = True +[xlon] + standard_name = longitude + long_name = longitude + units = radian + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[tsurf_wat] + standard_name = surface_skin_temperature_after_iteration_over_water + long_name = surface skin temperature after iteration over water + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[tsfc_wat] + standard_name = surface_skin_temperature_over_water + long_name = surface skin temperature over water + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[nthreads] + standard_name = number_of_openmp_threads + long_name = number of OpenMP threads available for physics schemes + units = count + dimensions = () + type = integer + intent = in +[dtzm] + standard_name = mean_change_over_depth_in_sea_water_temperature + long_name = mean of dT(z) (zsea1 to zsea2) + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out diff --git a/physics/SFC_Layer/UFS/skinsst.f90 b/physics/SFC_Layer/UFS/skinsst.f90 new file mode 100644 index 000000000..d517def2e --- /dev/null +++ b/physics/SFC_Layer/UFS/skinsst.f90 @@ -0,0 +1,578 @@ +!>\file skinsst.f90 +!! This file contains Rainer's skin temperature scheme. + +module state_eqn + implicit none +! --- coefficients for sigma-0 (based on Brydon & Sun fit, JGR 1999) + real, parameter, dimension(7) :: coef = (/ & + -1.36471E-01, 4.68181E-02, 8.07004E-01,-7.45353E-03,-2.94418E-03, & + 3.43570E-05, 3.48658E-05 /) +contains + + real function sig(t,s) +! --- sea water density (sigma = density - 1000) at p=0 + real, intent(in) :: t, s +! sig = coef(1)+s*coef(3)+ & +! t*(coef(2)+s*coef(5)+ & +! t*(coef(4)+s*coef(7)+t*coef(6))) + sig = t*(t*(s*coef(7)+t*coef(6)+coef(4)) & ! alt.grouping + +s*coef(5) +coef(2)) & + +s*coef(3) +coef(1) + return + end function sig + + real function dsigdt(t,s) +! --- thermal expansion coefficient + real, intent(in) :: t, s +! dsigdt = coef(2)+s*coef(5)+2.*t*(coef(4)+s*coef(7)+1.5*t*coef(6)) + dsigdt = 2.*t*(1.5*t*coef(6)+s*coef(7)+coef(4))+s*coef(5)+coef(2) ! alt.grouping + return + end + + real function dsigds(t,s) +! --- saline contraction coefficient + real, intent(in) :: t, s +! dsigds = coef(3)+t*(coef(5)+t*coef(7)) + dsigds = t*(t*coef(7)+coef(5))+coef(3) ! alt.grouping + return + end +end module state_eqn + + +module skinsst + implicit none + private + public :: skinsst_init, skinsst_run, skinsst_finalize + + contains + + subroutine skinsst_init() + end subroutine skinsst_init + + subroutine skinsst_finalize() + end subroutine skinsst_finalize + +!>\defgroup gfs_ocean_main GFS Simple Ocean Scheme Module +!! This subroutine calculates thermodynamical properties over +!! open water. +!! \section arg_table_skinsst_run Argument Table +!! \htmlinclude skinsst_run.html +!! + + subroutine skinsst_run( & + im, & ! horiz. loop extent in + iter, & ! ccpp loop counter in + wet, & ! .true. at ocean & lake points in + oceanfrac, & ! cell portion covered by ocean in + timestep, & ! model timestep in + xlon, xlat, & ! longitude, latitude in + sfcemis, & ! sea surface emissivity in + dlwflx, & ! absorbed downwelling LW flux in + sfcnsw, & ! net SW flux in + tsfco, & ! ocean/lake top layer temperature in + psfc, & ! surface pressure in + stress, & ! wind stress (N/m^2) in + wind, & ! atm. mid-layer 1 wind in + plyr1, & ! atm. mid-layer 1 pressure in + tlyr1, & ! atm. mid-layer 1 temperature in + qlyr1, & ! atm. mid-layer 1 humidity in + ulyr1, & ! atm. mid-layer 1 zonal wind in + vlyr1, & ! atm. mid-layer 1 meridional wind in + compres, & ! adiabat.compression factor, layer 1 => sfc in + cm, & ! drag coeff for momentum in + ch, & ! drag coeff for heat and moisture in + hvap, & ! latent heat of evaporation in + cp, & ! specif.heat of air in + rd, & ! gas constant of dry air in + eps, & ! ratio of gas constants, rd/rv in + sbc, & ! stefan-boltzmann constant in + ulwflx, & ! upwelling LW flux inout + tskin, & ! skin temp inout + xzts, & ! holding place for tskin inout + qsat, & ! saturation specif. humidity out + z_c, & ! sub-layer cooling thickness out + dt_warm, & ! warm-layer surface warming amount out + dt_cool, & ! skin layer cooling amount inout + evap, & ! kinematic latent heat flux, pos.up out + hflx, & ! kinematic sensible heat flux, pos.up out + ep, & ! potential latent heat flux, pos.up out + gflux, & ! heat flux over ground out + cmm, & ! momentum exchange coeff out + chh, & ! thermal exchange coeff out + lseaspray, & ! sea spray flag in + fm, & ! Monin-Obukhov function at surface in + fm10, & ! Monin-Obukhov function at 10m in + errmsg, errflg) + + use machine , only : kind_phys + use state_eqn + use funcphys, only : fpvs ! vapor pressure + + implicit none + +! --- input: + integer, intent(in) :: im,iter + logical, dimension(:), intent(in) :: wet + logical, intent(in) :: lseaspray + real (kind=kind_phys), dimension(:), intent(in) :: xlon,xlat, & + sfcemis, dlwflx, sfcnsw, tsfco, wind, psfc, plyr1, tlyr1, qlyr1, & + ulyr1, vlyr1, cm, ch, compres, stress, fm, fm10, oceanfrac + real (kind=kind_phys), intent(in) :: timestep, hvap, cp, rd, eps, & + sbc + +! --- inout: + real (kind=kind_phys), dimension(:), intent(inout) :: ulwflx, & + tskin, dt_cool, xzts + +! --- output: + real (kind=kind_phys), dimension(:), intent(out) :: evap, hflx, & + ep, qsat, gflux, cmm, chh, dt_warm, z_c + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + +! --- locals: + integer :: i, n, loop + real :: alon, alat, virt, rho_air, rho_wat, pvap, tsq, piston, vel, & + nonsol, & ! sum of nonsolar air-sea fluxes + spcifh = 3990., & ! seawater specific heat + grav = 9.806, & ! gravity + sss = 34.7, & ! sea surface salinity + penetr = 6., & ! shortwave penetration depth (m) + grnblu + integer,parameter :: itmax = 5 ! regula falsi iterations + real :: rnl_ts, hs_ts, rf_ts, alpha, beta, rch, ustar, & + hist(0:itmax) = 0., x1, x2, x3, y1, y2, y3, dif1, dif2, dif3 + +! variables for sea spray effect + real (kind=kind_phys) :: f10m, u10m, v10m, ws10, ru10, qss1, & + bb1, hflxs, evaps, ptem, tem + real (kind=kind_phys), parameter :: alps=0.75,bets=0.75,gams=0.15, & + ws10cr=30., conlf=7.2e-9, consf=6.4e-8 + + logical :: doprint, details + real(kind=kind_phys) :: frz=273.15, small=.05, testlon, testlat + real,parameter :: rad2deg = 57.2957795 + doprint(alon,alat)=abs(testlon-alon).lt.small .and. & + abs(testlat-alat).lt.small + +! --- latitude-dependent green/blue water transition factor affecting pendep +! grnblu(alat)= 1. - (alat/90.)**2 * (1. - .5*(alat/90.)**2) +! grnblu(alat)= 1. - .25*(alat/90.)**2 * (3. - (alat/90.)**2 * (alat/90.)**2) + grnblu(alat)= 1. - (alat/90.)**2 * (alat/90.)**2 * (1.5 - (alat/90.)**2) + +! --- piston velocity: at 0 m/s, molecular diffusion only. +! --- at 8 m/s, destroy warm layer during 600 sec time step +! --- based on 1 m thickness scale. +! piston(vel)= 1.4e-7 + 1.66667e-3*(vel/8.)**2 ! quadratic interpolation + piston(vel)= 1.4e-7 + 1.66667e-3*(vel/8.) ! linear interpolation + + if (iter.gt.1) return + + call get_testpt(testlon,testlat) +! --- temporary +! print '(a,2f8.2)','entering skinsst_run, testpt =',testlon,testlat + + do i = 1,im + if (wet(i) .and. oceanfrac(i) .gt. 0.) then + + alon=xlon(i)*rad2deg + alat=xlat(i)*rad2deg + if (doprint(alon,alat)) then + print 99,'entering skinsst_run lon,lat=',alon,alat, & + 'ocnfrac',oceanfrac(i), & ! ocean fraction + 'stress',stress(i), & ! wind stress (N/m^2) +! 'sfcemis',sfcemis(i), & ! sfc emissivity + 'wind',wind(i), & ! surface wind + 'pstonE3',piston(wind(i))*1.e3, & ! piston velocity + 'sfcnsw',sfcnsw(i), & ! total sky net SW flx into ocean + 'dlwflx',dlwflx(i), & ! absorbed downwelling LW flux + 'ulwflx',ulwflx(i), & ! upwelling LW flux + 'psfc',psfc(i)*.01, & ! surface pressure (mb) + 'plyr1',plyr1(i)*.01, & ! atm.layer 1 presure + 'tlyr1',tlyr1(i)-frz, & ! atm.layer 1 air temp + 'qlyr1',qlyr1(i)*1.e3, & ! atm.layer 1 humidity (g/kg) +! 'sigma_t',sig(tlyr1(i)-frz,sss), & ! sea water density - 1000 +! 'compres',compres(i), & ! midlyr-to-sfc adiab.compression + 'xzts',xzts(i)-frz, & ! previous tskin + 'dcoolE2',dt_cool(i)*100., & ! previous dtcool + 'tsfco',tsfco(i)-frz ! ocean top layer temperature + print '(5(a13,"=",l2))','lseaspray',lseaspray + end if + 99 format (/a,2f7.2/(5(a8,"=",f7.2))) + 98 format (/a,2f7.2/(4(a8,"=",es11.4))) + + virt = tlyr1(i) * (1. + (eps-1.)*qlyr1(i)) + rho_air = plyr1(i) / (rd*virt) + rch = rho_air * cp * ch(i) * wind(i) ! W/m^2/deg + cmm(i) = cm(i) * wind(i) + chh(i) = rho_air * ch(i) * wind(i) + ep(i) = 0. + rho_wat = 1000. + sig(tsfco(i)-frz,sss) + alpha = -dsigdt(tsfco(i)-frz,sss)/rho_wat + beta = dsigds(tsfco(i)-frz,sss)/rho_wat + ustar = sqrt(stress(i)/rho_air) ! air friction velocity + +! --- apply warm layer correction +! - - - - - - - - - - - - - - - - - - - - - - - - - - - +! --- tskin from last time step has been saved in xzts +! - - - - - - - - - - - - - - - - - - - - - - - - - - - + if (xzts(i).eq.0.) then ! initial xzts assumed to be 0 + tskin(i) = tsfco(i) + else + tskin(i) = xzts(i) ! tskin from previous time step + end if + dt_warm(i)=0. + + if (tskin(i)-frz.lt.-5.) print '(a,2f8.2,es13.3)', & + 'excessively cold tskin at lon,lat',alon,alat,tskin(i)-frz + + if (sfcnsw(i).lt.2.) then ! 2 W/m^2 cutoff +!! if (1.gt.0) then ! uncomment to disable warm lyr + tskin(i)=tsfco(i) +!! end if ! 1 > 0 + else ! SW flux > 0 + +! --- evaluate warmr-layer increment during current time step +! --- (curretly disabled for lake points) + + dt_warm(i) = sfcnsw(i) * timestep & + * 2./(penetr * grnblu(alat) * rho_wat * spcifh) +! --- note: dt_warm is cumulative. + tskin(i) = tskin(i) + dt_warm(i) & +! --- subtract old dt_cool (since dt_cool is not cumulative) + + dt_cool(i) ! sign convention: dt_cool > 0 +! --- cooling by heat diffusion + tskin(i) = tskin(i) + (tsfco(i)-tskin(i)) & + * min(1.,timestep*piston(wind(i))) + end if ! SW flux > 0 and ocnfrac > 0 + +! --- save (tskin - top layer T) for diagnostic purposes + dt_warm(i) = tskin(i) - tsfco(i) + +! --- start cool-skin iteration, using regula falsi (aiming for x_n = y_n) +! --- x1,x2,x3,y1,y2,y3 are consecutive dt_cool approximations. +! --- (curretly disabled for lake points) + + details = doprint(alon,alat) + + x1 = -.5 + x2 = +.5 + call surflx(nonsol, tskin(i)+x1, tlyr1(i)*compres(i), qlyr1(i), & + psfc(i), hflx(i), qsat(i), evap(i), hvap/cp, eps, rch, sbc, & + sfcemis(i), dlwflx(i), ulwflx(i), alon, alat, doprint(alon,alat)) + + call coolskin(ustar,nonsol,sfcnsw(i),evap(i),sss,alpha,beta, & + rho_wat,rho_air,tskin(i)+x1,grav,hvap, & + y1,z_c(i),alon,alat,doprint(alon,alat)) + dif1 = y1 - x1 + + if (y1.ne.0.) then + + do loop = 1,itmax + + call surflx(nonsol, tskin(i)+x2, tlyr1(i)*compres(i), qlyr1(i), & + psfc(i), hflx(i), qsat(i), evap(i), hvap/cp, eps, rch, sbc, & + sfcemis(i), dlwflx(i), ulwflx(i), alon, alat, doprint(alon,alat)) + + call coolskin(ustar,nonsol,sfcnsw(i),evap(i),sss,alpha,beta, & + rho_wat,rho_air,tskin(i)+x2,grav,hvap, & + y2,z_c(i),alon,alat,doprint(alon,alat)) + dif2 = y2 - x2 + + if (details) print '(a,3es11.3,i7)','(skinsst) x1,y1,y1-x1 =', & + x1,y1,dif1,loop + if (details) print '(a,3es11.3,i7)','(skinsst) x2,y2,y2-x2 =', & + x2,y2,dif2,loop + + x3 = (x1*dif2-x2*dif1)/(dif2-dif1) ! regula falsi + + if (abs(dif2).gt.1.e-5) then + + if (abs(dif1).gt.abs(dif2)) then + x1 = x2 + y1 = y2 + dif1 = dif2 + end if + x2 = x3 + + else ! we have convergence + dt_cool(i) = x3 + exit ! all done + end if ! convergence + hist(loop) = y2 + + end do ! iteration loop + + if (loop.eq.itmax) then + if (abs(hist(loop)).gt..5) then + print '(a,3f8.2/(11f7.2))','tskin not converging at lon,lat', & + alon,alat,hist(loop),(hist(n),n=1,loop) + end if + end if + tskin(i) = tskin(i)-dt_cool(i) + end if ! y1 nonzero + +! --- without lake model, use foundation teperature as sfe temp in lakes +! if (oceanfrac(i).eq.0.) tskin(i) = max(-2.,tref(i)) + +! --- save tskin and dt_cool for next call to skinsst + + xzts(i) = tskin(i) ! for use at next time step + +! --- according to GFS_surface_composites_inter.F90, +! --- dlwflx is the absorbed portion of downwelling LW flux. +! --- hence, the total downwelling flux is dlwflx/sfcemis +! --- and the reflected part is (1-sfcemis)*dlwflx/sfcemis + + ulwflx(i) = ulwflx(i) + dlwflx(i)*(1.-sfcemis(i))/sfcemis(i) + + if (lseaspray) then + f10m = fm10(i) / fm(i) + u10m = f10m * ulyr1(i) + v10m = f10m * vlyr1(i) + ws10 = sqrt(u10m*u10m + v10m*v10m) + ws10 = max(ws10,1.) + ws10 = min(ws10,ws10cr) + tem = .015 * ws10 * ws10 + ru10 = 1. - .087 * log(10./tem) + qss1 = fpvs(tlyr1(i)) + qss1 = eps * qss1 / (plyr1(i) + (eps-1.) * qss1) + tem = rd * cp * tlyr1(i) * tlyr1(i) + tem = 1. + eps * hvap * hvap * qss1 / tem + bb1 = 1. / tem + evaps = conlf * (ws10**5.4) * ru10 * bb1 + evaps = evaps * rho_air * hvap * (qss1 - qlyr1(i)) + evap(i) = evap(i) + alps * evaps + hflxs = consf * (ws10**3.4) * ru10 + hflxs = hflxs * rho_air * cp * (tskin(i) - tlyr1(i)) + ptem = alps - gams + hflx(i) = hflx(i) + bets * hflxs - ptem * evaps + endif + + if (doprint(alon,alat)) then + print 99,'exiting skinsst_run lon,lat=',alon,alat, & + 'virt',virt-frz, & ! virtual temp + 'rho_air',rho_air, & ! air density + 'pvap',pvap, & ! satur. vapor pressure (mb) + 'qsat',qsat(i), & ! satur. specif.humidity + 'hflx',hflx(i), & ! sensible heat flux + 'evap',evap(i), & ! latent heat flux + 'nonsol',nonsol, & ! net non-solar surface flux + 'sfcnsw',sfcnsw(i), & ! net solar surface flux + 'ulwflx',ulwflx(i), & ! upwelling LW flux + 'dwarmE2',dt_warm(i)*100., & ! temperature increment due to SW + 'dcoolE2',dt_cool(i)*100., & ! cool-skin temperature correction + 'tskin',tskin(i)-frz, & ! skin temperature + 'cmm',cmm(i), & + 'chh',chh(i), & + 'spray',bets*hflxs-ptem*evaps ! spray contrib. to heat flux + end if + +! --- convert fluxes from W/m^2 to "kinematic" i.e. velocity x fluxed variable + hflx(i) = hflx(i)/(rho_air * cp) ! deg m/sec + evap(i) = evap(i)/(rho_air * hvap) ! m/sec + + end if ! wet + end do ! im loop + + return + end subroutine skinsst_run + + + subroutine coolskin(ustar_a,f_nsol,f_sol_0,evap,sss,alpha,beta, & + rho_w,rho_a,ts,grav,latnt,deltat_c,z_c,alon,alat,doprint) + +! upper ocean cool-skin parameterizaion, Fairall et al, 1996. +! code extracted from NSST package + +! input: +! ustar_a : atmosphreic friction velocity at the air-sea interface (m/s) +! f_nsol : the "nonsolar" part of the surface heat flux (w/m^s) +! f_sol_0 : solar radiation at the ocean surface (w/m^2) +! evap : latent heat flux (w/m^2) +! sss : ocean upper mixed layer salinity (ppu) +! alpha : thermal expansion coefficient +! beta : saline contraction coefficient +! rho_w : oceanic density +! rho_a : atmospheric density +! ts : oceanic surface temperature +! grav : gravity +! latnt : latent heat of evaporation + +! output: +! deltat_c: cool-skin temperature correction (>0 means cooling) +! z_c : molecular sublayer (cool-skin) thickness (m) + + use machine , only : kind_phys + use module_nst_parameters, only: kw => tc_w,visw,cp_w, z_c_max, & + z_c_ini,ustar_a_min + implicit none + logical,intent(in) :: doprint + real(kind=kind_phys), intent(in) :: ustar_a,f_nsol,f_sol_0,evap, & + sss,alpha,beta,rho_w,rho_a,ts,grav,latnt,alon,alat + real(kind=kind_phys), intent(out):: deltat_c,z_c +! local variables: + real(kind=kind_phys), parameter :: frz=273.15 + real(kind=kind_phys) :: xi,hb,ustar1_a,bigc,deltaf,fxp + + if (doprint) print 99,'entering coolskin lon,lat=',alon,alat, & + 'ustar_a',ustar_a, & ! atmospheric friction velocity (m/s) + 'f_nsol',f_nsol, & ! "nonsolar" part of the surface heat flux + 'f_sol_0',f_sol_0, & ! solar radiation at the ocean surface + 'evap',evap, & ! latent heat flux (w/m^2) + 'rho_w',rho_w, & ! sea water density + 'rho_a',rho_a, & ! air density + 'ts',ts-frz, & ! ocean surface temperature + 'zc_ini',z_c_ini*1.e3 ! sublayer thickness (mm) + + 99 format (/a,2f7.2/(5(a8,"=",f7.2))) + 98 format (/a,2f7.2/(4(a8,"=",es11.4))) + + z_c = z_c_ini ! initial guess + + ustar1_a = max(ustar_a,ustar_a_min) + + call sw_rad_skin(z_c,fxp) + deltaf = f_sol_0 * fxp + + hb = alpha * (f_nsol-deltaf) + beta * sss * cp_w * evap/latnt + bigc = 16. * grav * cp_w * (rho_w * visw)**3 / (rho_a * kw)**2 + + if ( hb > 0 ) then + xi = 6./(1+(bigc * hb/ustar1_a**4)**0.75)**0.3333333 + else + xi = 6.0 + endif + z_c = min(z_c_max,xi * visw/(sqrt(rho_a/rho_w) * ustar1_a )) + + call sw_rad_skin(z_c,fxp) + + deltaf = f_sol_0 * fxp + deltaf = f_nsol - deltaf + if ( deltaf > 0 ) then + deltat_c = deltaf * z_c / kw + else + deltat_c = 0. + z_c = 0. + endif + + if (doprint) print 99,'exiting coolskin lon,lat=',alon,alat, & + 'fxp',fxp, & + 'deltaf',deltaf, & + 'hb',hb, & + 'bigc',bigc, & + 'xi',xi, & + 'delt_c',deltat_c, & ! skin layer temperature correction + 'z_c',z_c*1.e3 ! skin layer thickness (mm) + + return + end subroutine coolskin + + + subroutine sw_rad_skin(z,fxp) + ! original name: elemental subroutine sw_ohlmann_v1(z,fxp) + ! fraction of the solar radiation absorbed by the ocean at the depth z + + ! input: + ! z: depth (m) + + ! output: + ! fxp: fraction of solar radiation absorbed by the ocean at depth z (w/m^2) + + use machine , only : kind_phys + implicit none + real(kind=kind_phys),intent(in):: z + real(kind=kind_phys),intent(out):: fxp + + if(z>0) then + fxp=.065+11.*z-6.6e-5/z*(1.-exp(-z/8.0e-4)) + else + fxp=0. + endif + +! end subroutine sw_ohlmann_v1 + end subroutine sw_rad_skin + + + subroutine surflx( & + nonsol, & ! sum of nonsolar heat fluxes, pos.up + tsfc, & ! skin temperature + tlyr1, & ! temperature in lowest atmo layer + qlyr1, & ! spfc.humidity in lowest atmo layer + psfc, & ! surface pressure + hflx, & ! sensible heat flux, pos.up (out) + qsat, & ! satur.specf. humidity (out) + evap, & ! latent heat flux, pos.up (out) + elocp, & ! heat of evaporation over specif.heat, hvap/cp + eps, & ! ratio of air/vapor gas constants + rch, & ! rho * cp * ch * wind [W/deg] + sbc, & ! stefan-boltzmann constant + sfcemis, & ! sea surface emissivity + dlwflx, & ! absorbed downwelling LW flux, pos.down + ulwflx, & ! surface-emitted LW flux, pos.up (out) + alon,alat,doprint) + +! --- compute sum of nonsolar air-sea fluxes +! --- watch out for nonstandard sign convention: upwd fluxes are treated as > 0 + + use funcphys, only : fpvs ! vapor pressure + implicit none + real, intent(in) :: tsfc, tlyr1, qlyr1, psfc, elocp, eps, rch, sbc, & + sfcemis, dlwflx, alon, alat + logical, intent(in) :: doprint + real, intent(out) :: nonsol, qsat, evap, hflx, ulwflx + real :: pvap, frz=273.15 + + pvap = fpvs(tsfc) ! saturation vapor pressure (pa) + qsat = eps*pvap / (psfc + (eps-1.)*pvap) + evap = elocp * rch * (qsat - qlyr1) + hflx = rch * (tsfc - tlyr1) + ulwflx = sfcemis * sbc * tsfc**4 + nonsol = hflx + evap + ulwflx - dlwflx + + if (doprint) print 99,'exiting surflx lon,lat=',alon,alat, & + 'tsfc',tsfc-frz, & ! skin temperature + 'psfc',psfc*.01, & ! surface pressure (mb) + 'pvap',pvap, & ! saturation vapor pressure + 'qsat',qsat*1.e3, & !saturation specif. humidity (g/kg) + 'evap',evap, & ! latent heat flux + 'hflx',hflx, & ! sensible heat flux + 'ulwflx',ulwflx, & ! upwelling long-wave flux + 'dlwflx',dlwflx, & ! downwelling long-wave flux + 'nonsol',nonsol ! sum of nonsolar heat fluxes + 99 format (/a,2f7.2/(5(a8,"=",f7.2))) + 98 format (/a,2f7.2/(4(a8,"=",es11.4))) + + return + end subroutine surflx + +end module skinsst + + + subroutine get_testpt(testlon,testlat) + real,intent(OUT) :: testlon,testlat + +! testlon=145.56 ; testlat= -4.02 +! testlon=201.824 ; testlat= 21.538 +! testlat=169.215 ; testlon=-45.193 + +! testlon= 180.84 ; testlat= -0.51 +! testlon= 180.84 ; testlat= 0.51 +! testlon= 179.82 ; testlat= 0.51 +! testlon= 179.82 ; testlat= -0.51 + +! testlon= 215.00 ; testlat=-76.24 +! testlon= 273.39 ; testlat= 47.79 +! testlon= 273.76 ; testlat= 46.14 +! testlon= 274.14 ; testlat= 46.89 +! testlon= 274.53 ; testlat= 47.65 +! testlon= 274.84 ; testlat= 46.00 + +! testlon= 272.60 ; testlat= 48.70 + testlon= 273.39 ; testlat= 47.79 +! testlon= 292.34 ; testlat= 66.86 +! testlon= 245.31 ; testlat= 61.72 + +! print '(a,2f8.2)','(get_testpt) set test point location',testlon,testlat + return + end diff --git a/physics/SFC_Layer/UFS/skinsst.meta b/physics/SFC_Layer/UFS/skinsst.meta new file mode 100644 index 000000000..d9b6ee137 --- /dev/null +++ b/physics/SFC_Layer/UFS/skinsst.meta @@ -0,0 +1,368 @@ +[ccpp-table-properties] + name = skinsst + type = scheme + dependencies = ../../hooks/machine.F,../../tools/funcphys.f90 + +######################################################################## +[ccpp-arg-table] + name = skinsst_run + type = scheme +[im] + standard_name = horizontal_loop_extent + long_name = horizontal loop extent + units = count + dimensions = () + type = integer + intent = in +[iter] + standard_name = ccpp_loop_counter + long_name = loop counter for subcycling loops in CCPP + units = index + dimensions = () + type = integer + intent = in +[timestep] + standard_name = timestep_for_dynamics + long_name = timestep interval + units = s + dimensions = () + type = real + kind = kind_phys + intent = in +[hvap] + standard_name = latent_heat_of_vaporization_of_water_at_0C + long_name = latent heat of evaporation/sublimation + units = J kg-1 + dimensions = () + type = real + kind = kind_phys + intent = in +[cp] + standard_name = specific_heat_of_dry_air_at_constant_pressure + long_name = specific heat of dry air at constant pressure + units = J kg-1 K-1 + dimensions = () + type = real + kind = kind_phys + intent = in +[rd] + standard_name = gas_constant_of_dry_air + long_name = ideal gas constant for dry air + units = J kg-1 K-1 + dimensions = () + type = real + kind = kind_phys + intent = in +[eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys + intent = in +[psfc] + standard_name = surface_air_pressure + long_name = surface pressure + units = Pa + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[plyr1] + standard_name = air_pressure_at_surface_adjacent_layer + long_name = surface layer mean pressure + units = Pa + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[tlyr1] + standard_name = air_temperature_at_surface_adjacent_layer + long_name = surface layer mean temperature + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[qlyr1] + standard_name = specific_humidity_at_surface_adjacent_layer + long_name = surface layer mean specific humidity + units = kg kg-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[ulyr1] + standard_name = x_wind_at_surface_adjacent_layer + long_name = x component of surface layer wind + units = m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[vlyr1] + standard_name = y_wind_at_surface_adjacent_layer + long_name = y component of surface layer wind + units = m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[qsat] + standard_name = surface_specific_humidity_over_water + long_name = surface air saturation specific humidity over water + units = kg kg-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[tsfco] + standard_name = sea_surface_temperature + long_name = bulk sea surface temperature + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[xzts] + standard_name = derivative_of_diurnal_thermocline_layer_thickness_wrt_surface_skin_temperature + long_name = d(xz)/d(ts) + units = m K-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[tskin] + standard_name = surface_skin_temperature_over_water + long_name = surface skin temperature over water + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[dt_warm] + standard_name = derivative_of_heat_content_in_diurnal_thermocline_wrt_surface_skin_temperature + long_name = warm-skin surface temperature increment + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out + optional = True +[dt_cool] + standard_name = molecular_sublayer_temperature_correction_in_sea_water + long_name = sub-layer cooling amount (pos. means cooling) + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[z_c] + standard_name = molecular_sublayer_thickness_in_sea_water + long_name = sub-layer cooling thickness + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = True +[cm] + standard_name = surface_drag_coefficient_for_momentum_in_air_over_water + long_name = surface exchange coeff for momentum over water + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[ch] + standard_name = surface_drag_coefficient_for_heat_and_moisture_in_air_over_water + long_name = surface exchange coeff heat surface exchange coeff heat & moisture over ocean moisture over water + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[compres] + standard_name = ratio_of_exner_function_between_midlayer_and_interface_at_lowest_model_layer + long_name = Exner function ratio btw middle of lowest layer and ground + units = ratio + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[wet] + standard_name = flag_nonzero_wet_surface_fraction + long_name = flag indicating presence of some ocean or lake surface area fraction + units = flag + dimensions = (horizontal_loop_extent) + type = logical + intent = in +[oceanfrac] + standard_name = sea_area_fraction + long_name = fraction of horizontal grid area occupied by ocean + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[xlon] + standard_name = longitude + long_name = longitude of grid box + units = radian + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[xlat] + standard_name = latitude + long_name = latitude + units = radian + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[wind] + standard_name = wind_speed_at_lowest_model_layer + long_name = wind speed at lowest model level + units = m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[stress] + standard_name = surface_wind_stress_over_water + long_name = surface wind stress over water + units = m2 s-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[cmm] + standard_name = surface_drag_wind_speed_for_momentum_in_air_over_water + long_name = momentum exchange coefficient over water + units = m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[chh] + standard_name = surface_drag_mass_flux_for_heat_and_moisture_in_air_over_water + long_name = thermal exchange coefficient over water + units = kg m-2 s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[gflux] + standard_name = upward_heat_flux_in_soil_over_water + long_name = soil heat flux over water + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[evap] + standard_name = kinematic_surface_upward_latent_heat_flux_over_water + long_name = kinematic surface upward latent heat flux over water + units = kg kg-1 m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[hflx] + standard_name = kinematic_surface_upward_sensible_heat_flux_over_water + long_name = kinematic surface upward sensible heat flux over water + units = K m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[ep] + standard_name = surface_upward_potential_latent_heat_flux_over_water + long_name = surface upward potential latent heat flux over water + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[sbc] + standard_name = stefan_boltzmann_constant + long_name = Stefan-Boltzmann constant + units = W m-2 K-4 + dimensions = () + type = real + kind = kind_phys + intent = in +[sfcemis] + standard_name = surface_longwave_emissivity_over_water + long_name = surface lw emissivity in fraction over water + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[sfcnsw] + standard_name = surface_net_downwelling_shortwave_flux + long_name = total sky sfc net sw flx into ocean + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[ulwflx] + standard_name = surface_upwelling_longwave_flux_over_water + long_name = surface upwelling longwave flux at current time over water + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[dlwflx] + standard_name = surface_downwelling_longwave_flux + long_name = surface downwelling longwave flux at current time + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[lseaspray] + standard_name = flag_for_sea_spray + long_name = flag for sea spray parameterization + units = flag + dimensions = () + type = logical + intent = in +[fm] + standard_name = Monin_Obukhov_similarity_function_for_momentum_over_water + long_name = Monin-Obukhov similarity function for momentum over water + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[fm10] + standard_name = Monin_Obukhov_similarity_function_for_momentum_at_10m_over_water + long_name = Monin-Obukhov similarity parameter for momentum at 10m over water + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out